python - 如何用python求解非线性DE的9方程系统?
问题描述
我拼命地试图解决(并显示图表)一个由九个非线性微分方程组成的系统,这些方程模拟了回旋镖的路径。系统如下:
左边的所有字母都是变量,其他的要么是常数,要么是已知函数,取决于 v_G 和 w_z
我尝试过但scipy.odeint
没有确凿的结果(我遇到了这个问题,但解决方法不起作用。)
我开始认为这个问题与这些方程是非线性的或者分母中的函数可能导致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)
但是,我不断遇到与此处所示相同的问题,尤其是错误消息:
Excess work done on this call (perhaps wrong Dfun type)
我不太确定这意味着什么,但看起来求解器在解决系统时遇到了麻烦。无论如何,当我尝试通过 XYZ 坐标显示 3D 路径时,我只得到 3 或 4 个点,其中应该有 2000 之类的东西。
所以我的问题是: - 我在我的代码中做错了吗?- 如果没有,是否有其他可能更复杂的工具来解决这个系统?- 如果没有,是否有可能从这个 ODE 系统中得到我想要的?
提前致谢
解决方案
有几个问题:
- 如果我复制代码,它不会运行
- 您提到的解决方法不适用于 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})
推荐阅读
- css - 使用VScode,试图链接到一个资产,但链接中的两个字符是黄色而不是橙色......为什么?
- flutter - DataCell 溢出带有多个文本小部件的列 [dataTable]
- java - RecyclerView如何在每个项目上设置OnClickListener?
- java - 是否值得使用适配器设计模式而不是对现有接口或类进行更改?
- events - 事件处理程序中 QRCode 中的 Itext 7 getNumberOfPages 从第 2 页开始,带有 END_PAGE
- r - spinplot 如何绘制类似 ggplot2 的堆积条形图?
- javascript - 未经检查的 runtime.lastError:无法访问 url "" 的内容。扩展清单必须请求访问此主机的权限。在清单 3
- pine-script - 了解这行代码中的nz函数
- node.js - NodeJS writeFileSync “data”参数必须是字符串类型或 Buffer、TypedArray 或 DataView 的实例。收到一个对象实例
- php - 调用上下文中缺少 Codeigniter PHP SoapClient 端点名称(poolAlias)