如何用python求解非线性DE的9方程系统？

Question

我正在拼命地尝试求解（并显示图表）一个由九个非线性微分方程组成的系统，这些方程对回旋镖的路径进行建模。系统如下：

左边的所有字母都是变量，其他的要么是常数要么是已知函数，取决于 v_G 和 w_z

我尝试过使用scipy.odeint，但没有确定的结果（我有this issue，但解决方法无效。）

我开始认为问题与这些方程是非线性的或分母中的函数可能导致 scipy 求解器根本无法处理的奇点有关。但是，我对那种数学知识并不熟悉。 我必须解决这组方程组的哪些可能性？

编辑：对不起，如果我不够清楚。由于它模拟了回旋镖的路径，我的目标不是解析这个系统（即我不关心每个函数的数学表达式），而是获取每个函数在特定时间范围内的值（比如，从 t1 = 0s 到 t2 = 15s，每个值之间的间隔为 0.01s），以便显示每个函数的图形和回旋镖的质心图形（X，Y，Z 是它的坐标）。

这是我试过的代码：

import scipy.integrate as spi
import numpy as np

#Constants

I3 = 10**-3
lamb = 1
L = 5*10**-1
mu = I3
m = 0.1
Cz = 0.5
rho = 1.2
S = 0.03*0.4
Kz = 1/2*rho*S*Cz
g = 9.81

#Initial conditions

omega0 = 20*np.pi
V0 = 25
Psi0 = 0
theta0 = np.pi/2
phi0 = 0
psi0 = -np.pi/9
X0 = 0
Y0 = 0
Z0 = 1.8

INPUT = (omega0, V0, Psi0, theta0, phi0, psi0, X0, Y0, Z0) #initial conditions 


def diff_eqs(t, INP):  
    '''The main set of equations'''

    Y=np.zeros((9))

    Y[0] = (1/I3) * (Kz*L*(INP[1]**2+(L*INP[0])**2))
    Y[1] = -(lamb/m)*INP[1]
    Y[2] = -(1/(m * INP[1])) * ( Kz*L*(INP[1]**2+(L*INP[0])**2) + m*g) + (mu/I3)/INP[0]
    Y[3] = (1/(I3*INP[0]))*(-mu*INP[0]*np.sin(INP[6]))
    Y[4] = (1/(I3*INP[0]*np.sin(INP[3]))) * (mu*INP[0]*np.cos(INP[5]))
    Y[5] = -np.cos(INP[3])*Y[4]
    Y[6] = INP[1]*(-np.cos(INP[5])*np.cos(INP[4]) + np.sin(INP[5])*np.sin(INP[4])*np.cos(INP[3]))
    Y[7] = INP[1]*(-np.cos(INP[5])*np.sin(INP[4]) - np.sin(INP[5])*np.cos(INP[4])*np.cos(INP[3]))
    Y[8] = INP[1]*(-np.sin(INP[5])*np.sin(INP[3]))


    return Y   # For odeint



t_start = 0.0
t_end = 20
t_step = 0.01
t_range = np.arange(t_start, t_end, t_step)

RES = spi.odeint(diff_eqs, INPUT, t_range)

但是，我一直遇到与here 相同的问题，尤其是错误消息：

Excess work done on this call (perhaps wrong Dfun type)

我不太确定这意味着什么，但看起来求解器在求解系统时遇到了麻烦。无论如何，当我尝试通过 XYZ 坐标显示 3D 路径时，我只得到 3 或 4 个点，应该有 2000 之类的东西。

所以我的问题是： - 我在我的代码中做错了吗？ - 如果没有，是否有其他可能更复杂的工具来解决这个系统？ - 如果没有，是否有可能从这个 ODE 系统中得到我想要的？

提前致谢

Answer 1

有几个问题：

• 如果我复制代码，它不会运行
• 您提到的解决方法不适用于 odeint，给定 解决方案使用 ode
• odeint 的 scipy 参考说：“对于新代码，请使用 scipy.integrate.solve_ivp 求解微分方程。”
• 调用 RES = spi.odeint(diff_eqs, INPUT, t_range) 应该是 与函数头 def diff_eqs(t, INP) 一致。注意 顺序：RES = spi.odeint(diff_eqs,t_range, INPUT)

数学公式也存在一些问题：

• 看看你图片上的第三个公式。它没有趋势项，以零开头 - 这是什么意思？
• 很难检查您是否已将公式正确转换为代码，因为代码并未严格遵循公式。

下面我用 scipy solve_ivp 尝试了一个解决方案。如果 A 我能够运行钟摆，但如果 B 找不到回旋镖有意义的解决方案。所以检查一下数学，我猜数学表达式中有些错误。

对于图形，使用 pandas 将所有变量绘制在一起（参见下面的代码）。

import scipy.integrate as spi
import numpy as np
import pandas as pd

def diff_eqs_boomerang(t,Y):   
    INP = Y
    dY = np.zeros((9))
    dY[0] = (1/I3) * (Kz*L*(INP[1]**2+(L*INP[0])**2))
    dY[1] = -(lamb/m)*INP[1]
    dY[2] = -(1/(m * INP[1])) * ( Kz*L*(INP[1]**2+(L*INP[0])**2) + m*g) + (mu/I3)/INP[0]
    dY[3] = (1/(I3*INP[0]))*(-mu*INP[0]*np.sin(INP[6]))
    dY[4] = (1/(I3*INP[0]*np.sin(INP[3]))) * (mu*INP[0]*np.cos(INP[5]))
    dY[5] = -np.cos(INP[3])*INP[4]
    dY[6] = INP[1]*(-np.cos(INP[5])*np.cos(INP[4]) + np.sin(INP[5])*np.sin(INP[4])*np.cos(INP[3]))
    dY[7] = INP[1]*(-np.cos(INP[5])*np.sin(INP[4]) - np.sin(INP[5])*np.cos(INP[4])*np.cos(INP[3]))
    dY[8] = INP[1]*(-np.sin(INP[5])*np.sin(INP[3]))
    return dY   

def diff_eqs_pendulum(t,Y): 
    dY = np.zeros((3))
    dY[0] =  Y[1]
    dY[1] = -Y[0]
    dY[2] =  Y[0]*Y[1]
    return dY

t_start, t_end = 0.0, 12.0

case = 'A'

if case == 'A':         # pendulum
    Y = np.array([0.1, 1.0, 0.0]); 
    Yres = spi.solve_ivp(diff_eqs_pendulum, [t_start, t_end], Y, method='RK45', max_step=0.01)

if case == 'B':          # boomerang
    Y = np.array([omega0, V0, Psi0, theta0, phi0, psi0, X0, Y0, Z0])
    print('Y initial:'); print(Y); print()
    Yres = spi.solve_ivp(diff_eqs_boomerang, [t_start, t_end], Y, method='RK45', max_step=0.01)

#---- graphics ---------------------
yy = pd.DataFrame(Yres.y).T
tt = np.linspace(t_start,t_end,yy.shape[0])
with plt.style.context('fivethirtyeight'): 
    plt.figure(1, figsize=(20,5))
    plt.plot(tt,yy,lw=8, alpha=0.5);
    plt.grid(axis='y')
    for j in range(3):
        plt.fill_between(tt,yy[j],0, alpha=0.2, label='y['+str(j)+']')
    plt.legend(prop={'size':20})