如何使用 scipy 解决数组的非线性系统

Question

我写了一门课，目的是求解微分方程组（以 numpy.array 形式给出），为了求解非线性系统，我使用scipy.optimize.fsolve 使用此处的示例在一篇文章中找到，如果我尝试用于微分方程组，该方法适用于单个方程而失败！我写了一个最小、完整且可验证的示例 通过这种方式，您可以验证并深入了解该类的工作原理！

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

编辑 不，我很抱歉，但这不是我想要的……基本上，当我有一个方程组时，我必须得到相同维度的 K1 和 K2 self.u[i]

Answer 1

对于您的两个方程，您的 self.f() 函数返回一个长度为 2 的数组。因此，f1 和 f2 都是长度为 2 的数组。 结果是 equations() 返回一个长度为 4 的扁平数组。 但是，fsolve 期望 equations（它是 func 参数）返回一个长度为 2 的数组，因为您最初的猜测 x0=(2,2) 是一个长度为 2 的数组。

我不确定以下是否是您在您的上下文中尝试执行的操作，但您可以通过替换来修复此错误：

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)

使用以下代码（在行尾添加引用[0]，[1]）：

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]