python_常微分方程的求解

11月10日计算物理实验——求解中心力场记

调用

scipy.integrate.odeint(func, y0, t, args=(),· · · )

代码说明：

• func：微分方程的描写函数。要求微分方程化为标准形式，即$\frac{dy}{dt}=f(y,t)$，接受数组形式；

• y0：初值；

• t：需要求解的时间点；

• args=()：函数func中其他参数

步骤

• 列出方程；

• 将方程化为一阶微分方程，比如:
  $$\{\begin{array}{c}{y}_1=x\\ y_2=\dot{x}\\ y_3=y\\ y_4=\dot{y}\end{array}$$
  由上式得到：
  $$\{\begin{array}{c}\dot{y_1}\\ \dot{y_2}\\ \dot{y_3}\\ \dot{y_4}\end{array}$$
  关于$y_1,y_2,y_3,y_4$的表达式，放入一个矩阵就是输入的func；

• 设置初值与时间；

• 代入求解。

例子：

解：

列出方程
$$F=-\frac{GMm}{r^2}$$
如上假设，写成将方程写成：
$$\{\begin{array}{l}\dot{y_1}=y_2\\ \dot{y_2}=-GMm(y_1^2+y_3^2{)}^{-\frac{3}{2}}\cdot y_1\\ \dot{y_3}=y_4\\ \dot{y_4}=-GMm(y_1^2+y_3^2{)}^{-\frac{3}{2}}\cdot y_3\end{array}$$

代码

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

a = 1#G*M
m = 1

#get func
def orb(d_list, t):
    y1, y2, y3, y4 = d_list
    return np.array([y2, -a*m*(y1**2 + y3**2)**(-3/2)*y1, y4, -a*m*(y1**2 + y3**2)**(-3/2)*y3])

#change E and p
t1 = np.linspace(0, 10, 2000)
t1_1 = np.linspace(0, 100, 2000)
t2 = np.linspace(0, 10, 2000)
t3 = np.linspace(0, 2, 2000)
result1 = odeint(orb, [1, 0, 0, 1], t1)
result1_1 = odeint(orb, [1, -0.5, 0, 1], t1_1)
result2 = odeint(orb, [1, 0, 0, np.sqrt(2)], t2)
result3 = odeint(orb, [1, 0, 0, 1.8], t3)

#plot
fig = plt.figure()
axe = fig.add_subplot(221)
axe.plot(result1[:, 0], result1[:, 2])
axe.set_title('E<0')
axe = fig.add_subplot(222)
axe.plot(result1_1[:, 0], result1_1[:, 2])
axe.set_title('E<0,p2')
axe = fig.add_subplot(223)
axe.plot(result2[:, 0], result2[:, 2])
axe.set_title('E=0')
axe = fig.add_subplot(224)
axe.plot(result3[:, 0], result3[:, 2], 'b-', result3[:, 0], -result3[:, 2], 'b-')
axe.set_title('E>0')
plt.show()

运行结果：