How to solve non linear system of array using scipy

Question

I wrote a class with the aims to solve the system of differential equations (given in numpy.array form) , in order to solve the non linear system I'm using the scipy.optimize.fsolve using an example find here in one post, the method works fine with single equation while fail if I'm try to using for a system of differential equations ! I wrote a Minimal, Complete, and Verifiable example in this way you can verify and deep understand how the class works!

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()

EDIT No I'm sorry but is not what I was looking for ... basically when I have a system of equations I have to get K1 and K2 of the same dimension of self.u[i]

Iddo Hanniel Iddo Hanniel · Accepted Answer · 2018-08-14T15:15:43

For your two equations, your self.f() function returns an array of length 2. Therefore, f1 and f2 are both arrays of length 2. The result is then that equations() returns a flattened array of length 4. However, fsolve expects equations (it's func argument) to return an array of length 2, since your initial guess x0=(2,2) is an array of length 2.

I'm not sure if the following is what you're trying to do in your context, but you can fix this error by replacing:

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)

with the following code (adding references [0], [1] at the line end):

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)[0]
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)[1]

How to solve non linear system of array using scipy

1 Answers