如何在 Python 上以数值方式求解这个微分方程系统？答案

【问题标题】：How can I solve this system of differential equations numerically on Python?如何在 Python 上以数值方式求解这个微分方程系统？
【发布时间】：2022-08-24 06:41:04
【问题描述】：

我想用 Python 解决这个system of differential equations，我想知道该怎么做。自从我对这些系统不熟悉以来，我什么都没尝试过，尽管我已经解决了一些个人差异。 eqs.. 现在，我对solve_ivp 很熟悉，实际上就是这样。

我建议从阅读 scipy 教程开始
你的问题让我想起了The Simpsons。请先do some research，这将允许您制定一个具体的这里的问题是on-topic。请使用tour，并阅读How to Ask 和question checklist。完成一些研究并尝试后，请记住在问题中包含您的minimal reproducible example。欢迎来到堆栈溢出！
谢谢你们俩：我终于解决了。我发誓我在 Youtube 和这里搜索过如何解决这些系统之一，但我找不到使用 solve_ivp 的任何系统。我终于找到了一个。我会发布解决方案。

标签： python numerical-methods differential-equations

【解决方案1】：

这是代码：

import numpy as np
from numpy import sin,cos,sign,sqrt,pi
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
from matplotlib.pyplot import figure

g = 9.807
h = 1.8
l = 0.56 * h
R_sc = 14
mu = 0.1
k = 0.3
m = 80
v_0 = 12
phi_0 = np.pi/4

def velphi(t,y):
    dy=np.zeros([3])
    dy[1] = y[2]
    dy[0] = -mu*np.sqrt(g**2 + ((y[0])**2/(R_sc*np.cos(y[1])))**2)-(k/m)*(y[0])**2
    dy[2] = (g/l)*np.sin(y[1])-((y[0])**2/(l*R_sc))*np.sign(y[1])

    return dy

y0 = np.array([v_0, phi_0, 0])
time = np.linspace(0, 10, 10000)
sol = solve_ivp(velphi, (0,10), y0, method='RK45', t_eval=time, dense_output=True, rtol=1e-8, atol=1e-10)
t, pos, vel, phi, omega = sol.t, sol.y[0], sol.y[1], sol.y[2], sol.y[3]

plt.plot(t.T, phi.T,"b",linewidth = 0.65, label = "${\Phi}$")
plt.plot(t.T, omega.T,"g",linewidth = 0.65, label = "${\Omega}$")
s = 'Condicions inicials: ${\Phi}_{o}$ ='+str(np.degrees(phi_0))+'º,  ${\Omega}_{o}$ = 0 rad/s'
plt.title('Pèndol centrífug invertit | '+ s)
plt.xlabel('Temps (s)')
plt.ylabel(u'${\Phi}$ (rad) | ${\Omega}$ (rad / s)')
plt.grid(False)
plt.legend(loc = "upper right")
plt.show()
plt.plot(t.T, vel.T,"r",linewidth = 0.65, label = "v")
s = 'Condicions inicials: ${\Phi}_{o}$ ='+str(np.degrees(phi_0))+'º,  ${\Omega}_{o}$ = 0 rad/s'
plt.title('Pèndol centrífug invertit | '+ s)
plt.xlabel('Temps (s)')
plt.ylabel('v (m/s)')
plt.grid(False)
plt.legend(loc = "upper right")
plt.show()

这里的解决方案：

Velocity

Phi and the angular velocity

【讨论】：