python初探常微分方程组数值解

预备知识

包括最常用的四阶Ronge-Kutta数值方法以及四阶Adams预估-校正格式

四阶R-K

之所以是四阶R-K，是因为三阶精度太低，在步长略大时无法满足正常求解精度要求，而五阶以上虽然精度很高，但算法耗时。四阶R-K较为折中，既可满足精度要求，对于复杂计算其计算速度也可接受。
四阶R-K公式的标准格式如下：
对于形如 $\frac{dy}{dx}=f(x,y)$ 的常微分方程：

$$\begin{gathered} y_{i+1}=y_i+\frac{h}{6}(K_1+2K_2+2K_3+K_4) \\ {K}_1=f(x_i,y_i) \\ {K}_2=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_1) \\ {K}_3=f(x_i+\frac{h}{2},y_i+\frac{h}{2}{K}_2) \\ {K}_4=f(x_i+h,y_i+h{K}_3)\text{\hspace{1em}(1)} \end{gathered}$$

对于常微分方程组来说，将$y$ 与$f$ 向量化即可，其中$h$ 为步长，应尽量取小。

四阶Adams预估-校正公式

四阶Adams预估-校正格式可以看做一种改进的Euler格式，其计算效率相较四阶R-K更高，但其精度不如四阶R-K，在各种工程应用中应做出各种尝试后选取适合实际的算法。
四阶Adams预估-校正公式标准格式如下：
对于形如 $\frac{dy}{dx}=f(x,y)$ 的常微分方程：

预估：$${\bar{y}}_{i+1}=y_i+\frac{h}{24}(55f_i-59f_{i-1}+37f_{i-2}-9f_{i-3})$$
校正：$$y_{i+1}=\frac{h}{24}(9f(x_{i+1},{\bar{y}}_{i+1})+19f_i-5f_{i-1}+f_{i-2})\text{\hspace{1em}(2)}$$
其中$i=3,4,5...,n-1$
由于在预估式中出现了$f_{i-3},f_{i-2},f_{i-1}$，因此在启动Adams预估校正算法之前需要用其他算法进行前三个时间步的运算，较为普遍的是采用四阶R-K公式进行启动计算。

实战演练

理论推导

本文采用难度中等的常微分方程组进行演练，其物理含义为Apollo卫星的运动轨迹$(x,y)$
其满足以下方程组：
$$\begin{gathered} \ddot{x}=2\dot{y}+x-\frac{{\mu }^{\ast }(x+\mu )}{{r_1}^3}-\frac{\mu (x+{\mu }^{\ast })}{{r_2}^3} \\ \ddot{y}=-2\dot{x}+y-\frac{{\mu }^{\ast }y}{{r_1}^3}-\frac{\mu y}{{r_2}^3}\text{\hspace{1em}(2)} \end{gathered}$$
其中$$\begin{gathered} \mu =1/82.45,{\mu }^{\ast }=1-\mu \\ {r}_1=\sqrt{(x+\mu {)}^2+y^2},r_2=\sqrt{(x-{\mu }^{\ast }{)}^2+y^2} \end{gathered}$$
之后对常微分方程组进行变量代换，取$x_1=x,x=\dot{x},x_3=y,x_4=\dot{y}$
则式$(2)$可化为下式$(3)$
$$\begin{gathered} \dot{x_1}=x_2 \\ \dot{x_2}=2x_4+x_1-{\mu }^{\ast }(x_1+\mu )/{r_1}^3-\mu (x_1-{\mu }^{\ast })/{r_2}^3 \\ \dot{x_3}=x_4 \\ \dot{x_4}=-2x_2+x_3-{\mu }^{\ast }{x}_3/{r_1}^3-\mu x_3/{r_2}^3\text{\hspace{1em}(3)} \end{gathered}$$
同时由于式$(2)$的初始条件为
$$x(0)=1.2,y(0)=0,\dot{x}(0)=0,\dot{y}(0)=-1.04935751$$
此时我们得到了常微分方程组，初始条件，利用四阶R-K公式，设定合适的步长即可对其进行求解。

python实现

下面对其python代码实现进行展示：

创建求解常微分方程组的简单类

# -*- coding: utf-8 -*-
"""
Created on Mon Jan 25 11:12:42 2021

@author: xinghy
"""
# 1.2  0  0  -1.04935751'''initial'''
import numpy as np

class Ode_crew():
    '''equation crew'''
    def __init__(self, x1, x2, x3, x4, t_initial, t_last, t_iter):    
        self.x1 = [x1]
        self.x2 = [x2]
        self.x3 = [x3]
        self.x4 = [x4]
        self.time_pool = [t_initial]
        self.t_latest = t_initial
        self.t_iter = t_iter
        self.t_last = t_last
        self.x = [self.x1, self.x2, self.x3, self.x4]
        
    def f(self, x1, x2, x3, x4):
    '''方程组'''
        x1_ = x2
        x2_ = 2 * x4 + x1 - (1 - 1/82.45)*(x1 + 1/82.45)/np.sqrt((x1 + 1/82.45)**2+x3**2)**3 - (
            1/82.45)*(x1 -1 + 1/82.45)/np.sqrt((x1 - 1 + 1/82.45)**2 + x3**2)**3
        x3_ = x4
        x4_ = -2 * x2 + x3 - (1 - 1/82.45)*x3/np.sqrt((x1 + 1/82.45)**2+x3**2)**3 - (
            1/82.45)*x3/np.sqrt((x1- 1 + 1/82.45)**2 + x3**2)**3
        f = [x1_, x2_, x3_, x4_]
        return f
    
    def four_RK(self):
    '''四阶R-K'''
        K1 = self.f(self.x1[-1], self.x2[-1], self.x3[-1], self.x4[-1])
        K2 = self.f(self.x1[-1] + self.t_iter/2*K1[0], self.x2[-1] + self.t_iter/2*K1[1],
                    self.x3[-1] + self.t_iter/2*K1[2], self.x4[-1] + self.t_iter/2*K1[3])
        K3 = self.f(self.x1[-1] + self.t_iter/2*K2[0], self.x2[-1] + self.t_iter/2*K2[1],
                    self.x3[-1] + self.t_iter/2*K2[2], self.x4[-1] + self.t_iter/2*K2[3])
        K4 = self.f(self.x1[-1] + self.t_iter*K3[0], self.x2[-1] + self.t_iter*K3[1],
                    self.x3[-1] + self.t_iter*K3[2], self.x4[-1] + self.t_iter*K3[3])
        x_ori = np.array([self.x1[-1], self.x2[-1], self.x3[-1], self.x4[-1]])
        x_new = x_ori + np.array(self.t_iter/6 * (np.array(K1) + 2*np.array(K2) + 2*
                                                  np.array(K3) + np.array(K4)))
        self.x1.append(x_new[0])
        self.x2.append(x_new[1])
        self.x3.append(x_new[2])
        self.x4.append(x_new[3])
        self.t_latest += self.t_iter
        self.time_pool.append(self.t_latest)

之后将各种条件代入即可用指定算法进行运算：

import numpy as np
import matplotlib.pyplot as plt

from ode import Ode_crew

odecrew = Ode_crew(1.2, 0, 0, -1.04935751, 0, 20, 0.001)
while odecrew.t_latest < odecrew.t_last - 0.0001*0.1:
    odecrew.four_RK()
    # print(odecrew.time_pool[-1])
x_ = odecrew.x[0]
y_ = odecrew.x[2]
plt.scatter(x_, y_, c=y_, cmap=plt.cm.Blues, s=0.2)
plt.title('Apollo trajectory', fontsize=25)
plt.xlabel('x', fontsize=10)
plt.ylabel('y', fontsize=10)

运行后即可得到如下图所示的Apollo运行轨迹图：
由图可以看到卫星的运行轨迹结果，同时添加了一些颜色渐变的效果，比较简陋，之后有时间可以精进一下。

附录

前面也提到了四阶Adams预估-校正公式，是一个与四阶R-K算法总是同时出现的常用算法，但笔者运算之后发现四阶Adams所需步长要远小于四阶R-K所需步长，其精度总是与四阶R-K差几个量级，可能是方程组的固有性质所致，也可能是四阶Adams的固有弊端，在这里就不展开Adams的结果了~

以上，请各位巨巨批评指正！