Solve a system of non-linear equations in Python (scipy.optimize.fsolve)

[drudox]

I'm trying to solve this system of non linear equations using scipy.optimize.fsolve , I took this from an example in one other post [here][1]

my system of equation is the follow :

    for i in range(len(self.time)-1):
            
                def equations(variable):
                    k1,k2 = variable 
                    f1 = -k1 + self.f(self.time[i]+ (0.5+np.sqrt(3)/6)* self.dt , self.u[i]+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2) 
                    f2 = -k2 + self.f(self.time[i]+ (0.5-np.sqrt(3)/6)* self.dt , self.u[i]+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
                    return (f1,f2)
                            
                    
                k1,k2 = fsolve(equations,(5,5))

when I run the code I got :

    TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'equations'.Shape should be (2,) but it is (2, 1).

I don't know why this mismatch and how to fix it ..
I tried :

 for i in range(len(self.time)-1):
                ui = self.u[i]
                ti = self.time[i]
    
                def equations(variable):
                    k1,k2 = variable 
                    f1 = -k1 + self.f(ti+ (0.5+np.sqrt(3)/6)* self.dt , ui+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2) 
                    f2 = -k2 + self.f(ti+ (0.5-np.sqrt(3)/6)* self.dt , ui+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
                    return (f1,f2)
                            
                    
                k1,k2 = fsolve(equations,(1,1))

how can I fix this mismatch ?

[scidam]

It seems to be an implicit Runge-Kutta method implementation...
You are declaring equations in the loop, that is definitely inefficient.
equations should return ndarray, so, change return (f1,f2) to return np.array([f1, f2]).

[drudox]

yes it is an Implicit Runge Kutta method !! I've already done this ! and actually now my solve methods is this :

d

ef solve(self):
        '''
              perform Implicit RK4 scheme 
        '''
        self.time , self.u = self.dydt.createArray()
        print("Running Implicit RK 4th order") 
        
        for i in range(len(self.time)-1):

            def equations(variable):
                k1, k2 = variable
                f1 = -k1 + self.f(self.time[i]+ (0.5+np.sqrt(3)/6)* self.dt , self.u[i]+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2) 
                f2 = -k2 + self.f(self.time[i]+ (0.5-np.sqrt(3)/6)* self.dt , self.u[i]+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
                return np.array([f1,f2]).ravel() #.reshape(2,)  
                        
                
            k1 , k2 = fsolve(equations,(self.u[i],self.u[i]))
            
            self.u[i+1] = self.u[i] + self.dt/2* (k1 + k2)
 
        print('... Done!')

this works fine when I deal with a ODE single equations ... but doesn't works when I deal with system of ODE ... may you help me to find the right way to proceed ?

---

I declare the equation inside the for loop (where else ?) because k1 , k2 are function of

self.time[i], self.u[i] 

[drudox]

Here I report the whole class (I have cut the irrelevant part) in order to be testable for who want to try to give me help !

import numpy as np
from scipy.optimize import fsolve , newton_krylov
import matplotlib.pyplot as plt

class ImpRK4 :

    def __init__(self, fun , t0, tf, dt , y0):
        self.func = fun
        self.t0=t0
        self.tf=tf
        self.dt=dt
        self.u0=y0
        self.n = round((tf-t0)/dt)
        self.time  = np.linspace(self.t0, self.tf, self.n+1 )
        self.u     = np.array([self.u0  for i in range(self.n+1) ])

    def f(self,ti,ui):
         return  np.array([functions(ti,ui) for functions in self.func])     

    def solve(self): 


       for i in range(len(self.time)-1):

            def equations(variable):
                k1,k2 = variable
                f1 = -k1 + self.f(self.time[i]+ (0.5+np.sqrt(3)/6)* self.dt , self.u[i]+0.25*self.dt* k1+ (0.25+ np.sqrt(3)/6)*self.dt*k2) 
                f2 = -k2 + self.f(self.time[i]+ (0.5-np.sqrt(3)/6)* self.dt , self.u[i]+(0.25-np.sqrt(3)/6)*self.dt *k1 + 0.25*self.dt* k2)
                return np.array([f1,f2]).ravel() #.reshape(2,)  


            k1 , k2 = fsolve(equations,(2,2)) #(self.u[i],self.u[i]))
            self.u[i+1] = self.u[i] + self.dt/2* (k1 + k2)


       plt.plot(self.time,self.u)
       plt.show()    
def main():



func00 = lambda t,u : -10*(t-1)*u[0]

func01 = lambda t,u : u[1] 
func02 = lambda t,u : (1-u[0]**2)*u[1] - u[0]

func0x = np.array([func00])
func0 = np.array([func01,func02])



t0 = 0. 
tf = 2.      
u0 = y01   
dt = 0.008 

y01 = np.array([1.,1.])
diffeq = ImpRK4(func0,t0,tf,dt,y01)    


#y0  = np.array([np.exp(-5)])
#diffeq.solve()
#diffeq = ImpRK4(func0x,t0,tf,dt,y0) ## with single equations works
diffeq.solve()



if __name__ == '__main__': 
    main() 

[scidam]

In case of system of ordinary differential equations you will faced with necessity to solve
algebraic system of size m*s , where m -- the number of differential equations, s -- the number of stages in rk-method.

I slightly modified the code above to be able to handle systems of ODEs, but it still includes hardcoded
components (e.g. the number of stages, s=2), so improvements definitely needed:

import numpy as np
from scipy.optimize import fsolve , newton_krylov
import matplotlib.pyplot as plt



class ImpRK4 :
    def __init__(self, functions, t0, tf, dt, y0):
        self.functions = functions
        self.t0 = t0
        self.tf = tf
        self.dt = dt
        self.u0 = y0
        self.n = round((tf - t0) / dt)
        self.time  = np.linspace(self.t0, self.tf, self.n+1)
        self.u = np.array([self.u0  for i in range(self.n+1)])
        self.m = len(functions) # the number of equations (you can estimate it)
        
    def ode(self, t, u):
        return  np.array([func(t, u) for func in self.functions]) # u is an array of size `m`!
    
    def implicit_equations(self, variables, ti, ui):
        k1, k2 = variables[:self.m], variables[self.m:]
        f1 = -k1 + self.ode(ti + (0.5 + np.sqrt(3) / 6) * self.dt, ui + 0.25 * self.dt * k1 + (0.25 + np.sqrt(3) / 6) * self.dt * k2)
        f2 = -k2 + self.ode(ti + (0.5 - np.sqrt(3) / 6) * self.dt, ui + (0.25 - np.sqrt(3) / 6) * self.dt * k1 + 0.25 * self.dt * k2)
        return np.hstack([f1,f2]).ravel()
 
    def solve(self): 
        s = 2
        for i in range(len(self.time) - 1):
            res = fsolve(self.implicit_equations, np.random.rand(self.m * s), args=(self.time[i], self.u[i])) # s =2, the number of stadies of the num method
            k1, k2 = res[:self.m], res[self.m:]
            self.u[i + 1] = self.u[i] + self.dt / 2 * (k1 + k2)
        print('The system is solved.')



def main():
    func00 = lambda t,u : -10*(t-1)*u[0]
    func01 = lambda t,u : u[1] 
    func02 = lambda t,u : (1-u[0]**2)*u[1] - u[0]
    func0x = np.array([func00])
    func0 = np.array([func01,func02])
    t0 = 0. 
    tf = 2.      
    y01 = np.array([1., 1.])
    dt = 0.008 
    diffeq = ImpRK4(func0,t0,tf,dt,y01)    
    diffeq.solve()
    plt.plot(diffeq.u[:,1], diffeq.u[:,0])
    plt.show()
 
 
if __name__ == '__main__': 
    main()

[drudox]

Congratulations ! and thank you very much !! you are the best ;)

[drudox]

How can obtain the same result using another solver like newton_krylov ??

[scidam]

indeed, the newton_krylov method doesnt' support extra arguments, but you can wrap
the implicit-equations using lambda function, e.g.:

res = newton_krylov(lambda vars: self.implicit_equations(vars, self.time[i], self.u[i]), np.random.rand(self.m * s))