使用 scipy 求解第一个 ODE 的运动方程

Posted 2023-02-18

【中文标题】使用 scipy 求解第一个 ODE 的运动方程

【英文标题】：Solve motion equations for first ODE using scipy

【发布时间】：2021-12-30 14:00:35

【问题描述】：

我想使用 scipy solve_ivp 函数求解运动一阶 ODE 方程。我可以看到我做错了什么，因为这应该是一个椭圆，但我只绘制了四个点。你能发现错误吗？

import math
import matplotlib.pyplot as plt 
import numpy as np 
import scipy.integrate


gim = 4*(math.pi**2)
x0 = 1 #x-position of the center or h
y0 = 0 #y-position of the center or k
vx0 = 0 #vx position
vy0 = 1.1* 2* math.pi #vy position
initial = [x0, y0, vx0, vy0] #initial state of the system
time = np.arange(0, 1000, 0.01) #period

def motion(t, Z): 
  dx = Z[2] # vx
  dy = Z[3] # vy
  dvx = -gim/(x**2+y**2)**(3/2) * x * Z[2]
  dvy = -gim/(x**2+y**2)**(3/2) * y * Z[3]
  return [dx, dy, dvx, dvy]

sol = scipy.integrate.solve_ivp(motion, t_span=time, y0= initial, method='RK45')
plt.plot(sol.y[0],sol.y[1],"x", label="Scipy RK45 solution")
plt.show()

【问题讨论】：

我认为您没有发布您正在使用的确切代码：x 和 y 未在 motion 中定义。 solve_ivp 参数 t_span 期望时间间隔的终点（即两个数字）。最后，这个模型是什么？让 dvx 依赖于 vx（对于 dvy 也是如此）引入了阻尼，我不认为会产生椭圆；不过方程是非线性的，所以很难说。

【参考方案1】：

代码不应从新的工作区运行。引力公式中使用的变量x,y 没有在任何地方声明。所以插入x,y = Z[:2]或类似的行。

引力公式通常不包含速度分量。删除Z[2]、Z[3]。

再次检查时间跨度和评估时间参数的期望值。时间跨度从数组中获取前两个值。所以更改为t_span=time[[0,-1]] 来构建第一个和最后一个时间值对。

第二个情节的评价点不足，使用的线段太大。使用您的 time 数组应该不是问题。