基于python的常微分方程组数值解
- 预备知识
- 四阶R-K
- 四阶Adams预估-校正公式
- 实战演练
- 理论推导
- python实现
- 创建求解常微分方程组的简单类
- 之后将各种条件代入即可用指定算法进行运算:
- 附录
预备知识
包括最常用的四阶Ronge-Kutta数值方法以及四阶Adams预估-校正格式
四阶R-K
之所以是四阶R-K,是因为三阶精度太低,在步长略大时无法满足正常求解精度要求,而五阶以上虽然精度很高,但算法耗时。四阶R-K较为折中,既可满足精度要求,对于复杂计算其计算速度也可接受。
四阶R-K公式的标准格式如下:
对于形如
对于常微分方程组来说,将 与 向量化即可,其中
四阶Adams预估-校正公式
四阶Adams预估-校正格式可以看做一种改进的Euler格式,其计算效率相较四阶R-K更高,但其精度不如四阶R-K,在各种工程应用中应做出各种尝试后选取适合实际的算法。
四阶Adams预估-校正公式标准格式如下:
对于形如
预估:
校正:
其中
由于在预估式中出现了,因此在启动Adams预估校正算法之前需要用其他算法进行前三个时间步的运算,较为普遍的是采用四阶R-K公式进行启动计算。
实战演练
理论推导
本文采用难度中等的常微分方程组进行演练,其物理含义为Apollo卫星的运动轨迹
其满足以下方程组:
其中
之后对常微分方程组进行变量代换,取
则式可化为下式
同时由于式的初始条件为
此时我们得到了常微分方程组,初始条件,利用四阶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的结果了~
以上,请各位巨巨批评指正!