常微分方程初值问题: 一种改进的 Runge-Kutta-Fehlberg 方法
前言
Runge-Kutta-Fehlberg (RKF) 方法是一类非常经典且有趣的数值方法, 是由 Erwin Fehlberg 在十九世纪 60 年代为 NASA 工作时提出的一系列误差控制方法的统称. 我们介绍一种针对 RKF45 的改进方法, 是由 Bu 等人在 2014 年提出的 [1].
考虑初值问题\[ y’=f(t,y),\quad 0\leq t\leq L,\quad y(0)=y_0,\tag{1}\]RKF45 方法是一种显示地单步方法, 即利用已知的条件 \(y_i, t_i\) 通过计算得到 \(y(t_{i+1})\) 的近似值 \(y_{i+1}\). 基本的计算步骤可以简单地描述为\[ \begin{align*} y_{i+1} &= RK4(y_i),\\ e_{i+1} &= RK5(y_i)-y_{i+1} , \end{align*}\]其中 RK4, RK5 分别为四阶和五阶的 Runge-Kutta 方法, 具体可以参考 [2]. 我们知道, RKF45 方法的一个很重要的方面是可以动态调节时间步长的大小, 并且出于计算复杂度的考虑, 两种方法共享部分函数的计算值.
方法
在描述方法之前, 我们有必要对新方法的动机作一个介绍. 首先, 在经典的 RKF45 方法中, 由于 RK5 是高阶方法, 所以 \(RK5(y_i)\) 相比于 \(RK4(y_i)\) 更接近于真实解 \(y(t_{i+1})\). 如果我们定义绝对误差 \(E_i=y(t_i)-y_i\), 那么可以近似地成立\[ e_{i}\approx E_{i},\quad \text{i.e., }\quad RK5(y_{i-1})-y_i \approx y(t_{i})-y_{i}.\]由 RKF45 方法的定义可知, 局部误差满足 \(E_i=O(h^4)\) 和 \(E_i-e_i=y(t_{i})-RK5(y_{i-1})=O(h^5)\). 因此\[ y(t_i) = y_i+E_i = (y_i+e_i) + (E_i-e_i) = y_i+e_i +O(h^5), \tag{2}\]从而, 我们可以猜想 \(y_i+e_i\) 是一个比 \(y_i\) 更精确的近似值.
因此, 我们可以得到第 \(i\) 步我们需要求解的微分方程\[ y’=f(t,y),\quad t_i\leq t\leq t_{i+1},\quad y(t_i)=y_i+e_i,\tag{3}\]为了表述算法更加清楚, 我们把改进的算法 (见 [1]) 抄写如下:\[ \begin{cases} y_{i+1}=y_i+e_i+h \sum_{j=1}^6 b_j k_j,\\ e_{i+1} = h \sum_{j=1}^6(\hat{b}_j-b_j)k_j \\ \end{cases} \tag{4}\]其中, \(k_j=f \left( t_i+c_j h, y_i+e_i+h \sum_{l=1}^{i-1}\alpha_{j,l}k_l \right)\), \(b_j, \hat{b}_j\) 分别是 RK4 和 RK5 的系数, 具体可以参考 [1, 2].
等价形式
上一节中, 我们简述了改进的 RKF45 方法, 事实上它是 [2, 常微分方程初值问题: Runge-Kutta-Fehlberg 方法] 中介绍的改进方法的等价形式. 具体地, 由定义 \(e_i=RK5(y_{i-1})-y_i\) 可知 \(y_i+e_i=RK5(y_{i-1})=:w_i\), 算法 (4) 等价为\[ \begin{cases} y_{i+1}=w_i+h \sum_{j=1}^6 b_j k_j,\\ e_{i+1} = h \sum_{j=1}^6(\hat{b}_j-b_j)k_j, \\ \end{cases} \tag{4}\]其中, \(k_j=f \left( t_i+c_j h, w_i+h \sum_{l=1}^{i-1}\alpha_{j,l}k_l \right)\).
数值实验
例 1. 考虑描述刚体旋转运动的 Euler 方程 [3]:\[ \begin{align*} I_1 y’_1 &= (I_2-I_3)y_2y_3,\\ I_2 y’_2 &= (I_3-I_1)y_3y_1,\\ I_3 y’_3 &= (I_1-I_2)y_1y_2 + f(t), \end{align*}\]其中, \(t\in [0,T]\), 初值 \(y(0)=[1,0,0.9]\), 外力\[ f(t)= \begin{cases} 0.25 \sin^2 t, & t\in [3\pi,4\pi],\\ 0,& t\in [0,3\pi)\cup (4\pi, \infty ). \end{cases}\]在数值计算中, 我们取常数 \(I=[0.5, 2, 3]\) 和时间上界 \(L=20\). 首先, 我们利用 scipy.integrate 中的八阶显式 Runge-Kutta 方法 [3] DOP853 (精度设置为浮点数的机器精度 eps) 得到的值作为精确解:\[ y(20)\approx [0.9877945603406532, 0.1231409420171716, 1.2625251695845630],\]然后再利用经典的 RKF45 方法和改进的 RKF45 方法与精确解比较. 我们在下面的表格中列出了基于不同精度要求得到的数值解. 可以看出两种方法需要的时间迭代步数几乎完全相同. 另外, 在下图中, 我们也画出了误差-迭代步数间的相关图像, 其中每一个节点对应一种数值精度. 详细的代码实现见附录.
# results for rkf45:
10e-9: 1355 0.9877945589174362 0.1231409531491341 1.2625251693740960
10e-10: 2411 0.9877945602003257 0.1231409429644286 1.2625251695525301
10e-11: 4285 0.9877945603239571 0.1231409421165484 1.2625251695815041
10e-12: 7608 0.9877945603386764 0.1231409420282025 1.2625251695844755
10e-13: 13530 0.9877945603402347 0.1231409420193917 1.2625251695847712
# results for rkf45 improved:
10e-9: 1355 0.9877945602243190 0.1231409438915988 1.2625251696891766
10e-10: 2411 0.9877945603570603 0.1231409420073311 1.2625251695853295
10e-11: 4284 0.9877945603451845 0.1231409419985611 1.2625251695836004
10e-12: 7608 0.9877945603405285 0.1231409420180147 1.2625251695848110
10e-13: 13532 0.9877945603404181 0.1231409420184413 1.2625251695848039
参考
[1] Sunyoung Bu, Wonkyu Chung and Philsu Kim. A method for improving the embedded Runge Kutta Fehlberg 4(5). International Journal of Mathematical and Computational Science, 8(8), 2014.
[2] 常微分方程初值问题: Runge-Kutta-Fehlberg 方法.
[3] E. Hairer, S. P. Norsett G. Wanner, “Solving Ordinary Differential Equations I: Nonstiff Problems”, Sec. II.
附录
RKF45 方法的 Python 实现
# Runge-Kutta-Fehlberg 45 (improved) method
# for Euler’s equation of rotation of a rigid body.
#
# --- by numanal.com, 2021/21/17
#
import numpy as np
import matplotlib.pyplot as plt
class my_ode_ivp:
def __init__(self, fun, t_span, y0, yfun=None):
self.fun = fun
self.t_span = t_span
self.y0 = y0
self.tall = None
self.yall = None
def rkf45(self, TOL=1e-6, hmin=1e-7, hmax=0.2, improved_flag=False):
# improved_flag indicate to use improved rkf45 or classical rkf45
#
flag = True
h = hmax
yi = self.y0
ti = self.t_span[0]
N, Ny = 1000, 1 # Ny is the total number of computed y
yall = np.zeros((yi.shape[0], N))
tall = np.zeros(N)
yall[:,0] = yi
tall[0] = ti
while flag:
k1=h*self.fun(ti, yi)
k2=h*self.fun(ti+.25*h, yi+.25*k1)
k3=h*self.fun(ti+3*h/8, yi+3*k1/32+9*k2/32)
k4=h*self.fun(ti+12*h/13, yi+1932*k1/2197-7200*k2/2197+7296*k3/2197)
k5=h*self.fun(ti+h, yi+439*k1/216-8*k2+3680*k3/513-845*k4/4104)
k6=h*self.fun(ti+.5*h, yi-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40)
# compute |y_{i+1}-w_{i+1}|/h
D = k1/360-128*k3/4275-2197*k4/75240+k5/50+2*k6/55
R = np.linalg.norm(D)/h
if R<TOL: # accurate is acceptable
ti += h
yi += 25*k1/216+1408*k3/2565+2197*k4/4104-k5/5 # rkf45
if improved_flag: # imaproved rkf45
yi += D
# save new yi
if Ny>=yall.shape[1]:
yall = np.hstack((yall, np.zeros((yi.shape[0],1000))))
tall = np.concatenate((tall, np.zeros(1000)))
tall[Ny]=ti
yall[:,Ny]=yi
Ny += 1
# compute q
q = .84*(TOL/R)**(.25)
# new size for step i-th or (i+1)-th
if q<=.2:
h=.2*h
elif q>=4:
h=4*h
else:
h *= q
if h>hmax:
h=hmax
if ti >= self.t_span[1]:
flag=False
elif ti+h > self.t_span[1]:
h = self.t_span[1]-ti
if h<hmin:
flag=False
self.tall = tall[:Ny]
self.yall = yall[:,:Ny]
def plot(self):
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
# for j in range(3):
ax.plot(self.tall, self.yall[0,:], 'r-', label='$y_1$')
ax.plot(self.tall, self.yall[1,:], 'g-', label='$y_2$')
ax.plot(self.tall, self.yall[2,:], 'b-', label='$y_3$')
plt.xlabel('t')
plt.ylabel('y')
plt.legend()
plt.show()
# setting for ode: f(t,y), y0, t_span
def f(t, y):
I1, I2, I3 = 0.5, 2., 3.
ans = np.zeros(3)
ans[0] = (I2-I3)*y[1]*y[2]/I1
ans[1] = (I3-I1)*y[2]*y[0]/I2
ans[2] = (I1-I2)*y[0]*y[1]/I3
if t>= 3*np.pi and t<=4*np.pi:
ans[2] += 0.25*np.sin(t)**2
return(ans)
y0 = np.array([1., 0., 0.9])
t_span = np.array([0., 20.])
# solve and print sol
ode = my_ode_ivp(fun=f, t_span=t_span, y0=y0)
tol = 1.0e-11
ode.rkf45(TOL=tol, hmin=.00001, hmax=.25, improved_flag=False)
yend = ode.yall[:,-1]
print("{:16d} {:.16f} {:.16f} {:.16f} ".format(ode.tall.shape[0], yend[0], yend[1],yend[2]))
ode.plot()
计算真实解的 Python 代码
# Solve EULR Example in [1]:
# - Euler’s equation of rotation of a rigid body.
#
# [1] E. Hairer, S. P. Norsett G. Wanner,
# “Solving Ordinary Differential Equations I: Nonstiff Problems”, Sec. II.
#
# --- by numanal.com, 2021/21/17
#
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp, DOP853
# Define derivative function
def f(t, y, I1, I2, I3):
dydt = [(I2-I3)*y[1]*y[2]/I1, (I3-I1)*y[2]*y[0]/I2, (I1-I2)*y[0]*y[1]/I3]
if t>= 3*np.pi and t<=4*np.pi:
dydt[2] += 0.25*np.sin(t)**2
return dydt
# Define time spans, initial values, and constants
yinit = [1.,0.,0.9]
I1, I2, I3 = 0.5, 2., 3.
# define nodes to compute
t_eval = np.linspace(0, 20, 200, endpoint=True)
eps = np.finfo(float).eps
#
# Solve with tol=eps
#
sol = solve_ivp(f, t_span=[0, 20], y0=yinit, args=(I1, I2, I3), t_eval=t_eval, atol=eps, rtol=100*eps, method='DOP853')
# print sol on the endpoint
yend = sol.y[:,-1]
print("{:.7f} {:.16f} {:.16f} {:.16f} ".format(sol.t[-1], yend[0], yend[1],yend[2]))
#
# Plot
#
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(sol.t, sol.y[0],'-b')
ax.plot(sol.t, sol.y[1],'-r')
ax.plot(sol.t, sol.y[2],'-g')
ax.legend(['$y1$','y2','y3'])
plt.show()
误差图像的程序实现
#
# --- by numanal.com, 2021/21/17
#
import numpy as np
import matplotlib.pyplot as plt
y = np.array([[0.9877945603406532, 0.1231409420171716, 1.2625251695845630]])
yrkf = np.array([0.9877945589174362, 0.1231409531491341, 1.2625251693740960,
0.9877945602003257, 0.1231409429644286, 1.2625251695525301,
0.9877945603239571, 0.1231409421165484, 1.2625251695815041,
0.9877945603386764, 0.1231409420282025, 1.2625251695844755,
0.9877945603402347, 0.1231409420193917, 1.2625251695847712,]).reshape(5,3)
yrkf2 = np.array([0.9877945602243190, 0.1231409438915988, 1.2625251696891766,
0.9877945603570603, 0.1231409420073311, 1.2625251695853295 ,
0.9877945603451845, 0.1231409419985611, 1.2625251695836004,
0.9877945603405285, 0.1231409420180147, 1.2625251695848110 ,
0.9877945603404181, 0.1231409420184413, 1.2625251695848039 ,]).reshape(5,3)
e1 = np.linalg.norm(yrkf-y,axis=1,ord=2)
e2 = np.linalg.norm(yrkf2-y,axis=1)
n1 = np.array([1355, 2411 ,4285 ,7608 ,13530,])
n2 = np.array([1355, 2411 ,4284 ,7608 ,13532,])
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.plot(n1, e1, 'x--', label='RKF45')
ax.plot(n2, e2, '^--', label='RKF45 Improved')
ax.set_yscale('log')
plt.xlabel('Times of iteration')
plt.ylabel('Error')
plt.legend()
plt.show()
print()
