数值分析-数值微分
前言
上一篇在数值积分中我们介绍了插值型求积公式以及如何分析对应的代数精度、稳定性和收敛性。当然啦，有数值积分就有数值微分，在以往的学习中我们可能已经掌握了很多解析解法，但是解析解法只能用来求解一些具有特定类型的方程，大部分方程都是需要数值求解的。 按照教材内容，我们着重点是一阶常微分方程初值问题的数值微分方法。
具体内容详见 《数值分析》 李庆扬、王能超、易大义（华中科大出版社）第五章
敲黑板
-
一阶方程的初值问题：
{ y ′ = f ( x , y ) y ( x 0 ) = y 0 \begin{cases} y^{'}=f(x,y)\\ y(x_0)=y_0 \end{cases} {y′=f(x,y)y(x0)=y0
数值解法就是求出 y ( x ) y(x) y(x)在一系列等距节点 x 1 , x 2 , . . . , x n x_1,x_2,...,x_n x1,x2,...,xn处的值。对于初值问题数值解法一般采取步进式求解过程，即通过前面求解的节点结果推出后续节点的值。
欧拉前向差商法(欧拉折线法)
-
设h为等距节点分隔间距，使用向前差分代替该点的导数值，即：
y ′ ( x 0 ) ≈ y ( x 1 ) − y ( x 0 ) h y^{'}(x_0) \approx \frac{y(x_1)-y(x_0)}{h} y′(x0)≈hy(x1)−y(x0)即在已知初值点 y ( x 0 ) y(x_0) y(x0)的情况下，后续节点的计算方式为：
y 1 = y 0 + h f ( x 0 , y 0 ) , y 2 = y 1 + h f ( x 1 , y 1 ) , . . . , y n + 1 = y n + h f ( x n , y n ) y_1=y_0+hf(x_0,y_0), \,\, y_2=y_1+hf(x_1,y_1),..., \,\, y_{n+1}=y_n+hf(x_n,y_n) y1=y0+hf(x0,y0),y2=y1+hf(x1,y1),...,yn+1=yn+hf(xn,yn)
代码实现
# 前向Euler
# Author: Junno
# Date: 2022-04-14
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
# 1-order differential equation
F=lambda x,y: y-2*x/y
# parsing solution of y
Solve_Y=lambda x:np.sqrt(1+2*x)
# initial value
initial_x=0
initial_y=1
# sep width
h=0.1
x_min=0
x_max=5
x=np.arange(x_min,x_max,h)
Euler_polygonal_arc=[]
for i in range(len(x)):
if x[i]==initial_x:
Euler_polygonal_arc.append(initial_y)
else:
Euler_polygonal_arc.append(Euler_polygonal_arc[i-1]+h*F(x[i-1],Euler_polygonal_arc[i-1]))
plt.figure()
plt.scatter(x,Solve_Y(x),s=5, marker='*',c='black')
plt.scatter(x,Euler_polygonal_arc,s=5, marker='*',c='red')
plt.xlabel('x')
plt.ylabel('y')
Parsing_Solution=mpatches.Patch(color='black', label='Parsing_Solution')
Numerical_Solution=mpatches.Patch(color='red', label='Numerical_Solution')
plt.title('Euler_polygonal_arc method in range {}-{} with h={:.1f}'.format(x_min,x_max,h))
plt.legend(handles=[Parsing_Solution,Numerical_Solution],loc='best')
plt.show()
在0-5的范围内步长0.1查看结果：
可以看到欧拉前向差分仅在
[
0
,
2
]
[0,2]
[0,2]范围内还比较拟合，之后显著偏离原来的曲线，即前向差分欧拉方法的误差还是很大的。
局部截断误差
-
记 y ( x ) y(x) y(x)是给定初值问题的准确解，假设在第n步计算准确的情况下，即 y n = y ( x n ) y_n=y(x_n) yn=y(xn)，考虑截断误差 T n + 1 = y ( x n + 1 ) − y n + 1 T_{n+1}=y(x_{n+1})-y_{n+1} Tn+1=y(xn+1)−yn+1，称之为局部截断误差。
-
若某算法的局部截断误差为 O ( h p + 1 ) O(h^{p+1}) O(hp+1)，则称该算法具有p阶精度。
-
一般在节点 x n + 1 x_{n+1} xn+1处采用Taylor展开分析局部截断误差，来看下欧拉前向差分的局部截断误差：
T
n
+
1
=
y
(
x
n
+
1
)
−
y
n
+
1
=
[
y
(
x
n
)
+
h
y
′
(
x
n
)
+
h
2
2
y
′
′
(
x
n
)
+
O
(
h
3
)
]
−
[
y
n
+
h
f
(
x
n
,
y
n
)
]
=
h
2
2
y
′
′
(
x
n
)
+
O
(
h
3
)
≈
O
(
h
2
)
T_{n+1}=y(x_{n+1})-y_{n+1}=[y(x_n)+h y^{'}(x_n)+\frac{h^2}{2} y^{''}(x_n)+O(h^3)]-[y_n+h f(x_n,y_n)]=\frac{h^2}{2} y^{''}(x_n)+O(h^3) \approx O(h^2)
Tn+1=y(xn+1)−yn+1=[y(xn)+hy′(xn)+2h2y′′(xn)+O(h3)]−[yn+hf(xn,yn)]=2h2y′′(xn)+O(h3)≈O(h2)
即前向欧拉方法具有1阶精度。
后向差分欧拉
-
对于前向差分我们有， y n + 1 = y n + h f ( x n , y n ) y_{n+1}=y_n+hf(x_n,y_n) yn+1=yn+hf(xn,yn)，自然地，后向差分的Euler公式即为 y n + 1 = y n + h f ( x n , y n + 1 ) y_{n+1}=y_n+hf(x_n,y_{n+1}) yn+1=yn+hf(xn,yn+1)，但是这类计算属于隐式计算，需要不断迭代，不易求解。若迭代收敛，则 y n + 1 = lim k → ∞ y n + 1 k y_{n+1}=\lim_{k \rightarrow \infty} y_{n+1}^k yn+1=limk→∞yn+1k。
-
对于后向欧拉公式，可以证得（教材P108）其局部截断误差为 T n + 1 = − h 2 2 y ′ ′ ( x n ) + O ( h 3 ) T_{n+1}=-\frac{h^2}{2} y^{''}(x_n)+O(h^3) Tn+1=−2h2y′′(xn)+O(h3)
代码实现
# 前向Euler
# Author: Junno
# Date: 2022-04-24
def iter_back_Euler(x,y_,y,F,h,eps=1e-4,N=100):
v0=y-eps*10
v1=y
while N or abs(v1-v0)>eps:
v1=y_+h*F(x,v0)
v0=v1
N-=1
return y
Euler_polygonal_arc_back=[]
for i in range(len(x)):
if x[i]==initial_x:
Euler_polygonal_arc_back.append(initial_y)
else:
Euler_polygonal_arc_back.append(iter_back_Euler(x[i],Euler_polygonal_arc_back[i-1],Euler_polygonal_arc_back[i-1]+h*F(x[i-1],Euler_polygonal_arc_back[i-1]),\
F,h))
相应的求解结果：
梯形平均方法
-
由于前向和后向欧拉公式的误差主项互为相反数，自然想到两者结合来抵消误差，提高精度。这种方法称为梯形方法，即 y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] y_{n+1}=y_n+\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})] yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1)]
-
局部截断误差分析
T n + 1 = y ( x n + 1 ) − y n + 1 = y ( x n + 1 ) − y ( x n ) − h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ) ] = y ( x n + 1 ) − y ( x n ) − h 2 [ y ′ ( x n ) + y ′ ( x n + 1 ) ] = y ( x n ) + h y ′ ( x n ) + h 2 2 y ′ ′ ( x n ) + h 3 3 ! y ′ ′ ′ ( x n ) − y ( x n ) − h 2 [ y ′ ( x n ) + y ′ ( x n ) + h y ′ ′ ( x n ) + h 2 2 y ′ ′ ′ ( x n ) ] + O ( h 4 ) = − h 3 12 y ′ ′ ′ ( x n ) + O ( h 4 ) ≈ O ( h 3 ) \begin{aligned} T_{n+1}&=y(x_{n+1})-y_{n+1}\\ &=y(x_{n+1})-y(x_n)-\frac{h}{2}[f(x_n,y_n)+f(x_{n+1},y_{n+1})]\\ &=y(x_{n+1})-y(x_n)-\frac{h}{2}[y^{'}(x_n)+y^{'}(x_{n+1})]\\ &=y(x_n)+h y^{'}(x_n)+\frac{h^2}{2} y^{''}(x_n)+\frac{h^3}{3!} y^{'''}(x_n)-y(x_n)-\frac{h}{2}[y^{'}(x_n)+y^{'}(x_n)+hy^{''}(x_n)+\frac{h^2}{2} y^{'''}(x_n)]+O(h^4)\\ &=-\frac{h^3}{12}y^{'''}(x_n)+O(h^4)\\ & \approx O(h^3) \end{aligned} Tn+1=y(xn+1)−yn+1=y(xn+1)−y(xn)−2h[f(xn,yn)+f(xn+1,yn+1)]=y(xn+1)−y(xn)−2h[y′(xn)+y′(xn+1)]=y(xn)+hy′(xn)+2h2y′′(xn)+3!h3y′′′(xn)−y(xn)−2h[y′(xn)+y′(xn)+hy′′(xn)+2h2y′′′(xn)]+O(h4)=−12h3y′′′(xn)+O(h4)≈O(h3)
即梯形公式具有2阶精度。但梯形方法总体上来说仍属于隐式计算！
改进的Euler格式
-
结合前向欧拉和梯形公式可以得到显式计算下的Euler改进公式，即先用前向Euler计算一个较为粗糙的预测值，它的精度可能较差，然后再使用梯形公式进行校正，得到校正值。这样构成的预测—校正系统成为改进的Euler格式。过程如下：
预 测 ： y n + 1 ˉ = y n + h f ( x n , y n ) 校 正 ： y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ˉ ) ] \begin{aligned} &预测：\bar{y_{n+1}}=y_n+h f(x_n,y_n)\\ &校正：y_{n+1} =y_n+\frac{h}{2} [f(x_n,y_n)+f(x_{n+1},\bar{y_{n+1}})] \end{aligned} 预测：yn+1ˉ=yn+hf(xn,yn)校正：yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1ˉ)]
可以表示成以下平均化形式：
{ y p = y n + h f ( x n , y n ) y c = y n + h f ( x n + 1 , y p ) y n + 1 = 1 2 ( y p + y c ) \begin{cases} y_p=y_n+h f(x_n,y_n)\\ y_c=y_n+h f(x_{n+1},y_p)\\ y_{n+1}=\frac{1}{2} (y_p+y_c) \end{cases} ⎩⎪⎨⎪⎧yp=yn+hf(xn,yn)yc=yn+hf(xn+1,yp)yn+1=21(yp+yc)代码实现
# 改进的Euler格式
# Author: Junno
# Date: 2022-04-14
Euler_improved=[]
for i in range(len(x)):
if x[i]==initial_x:
Euler_improved.append(initial_y)
else:
yp=Euler_improved[i-1]+h*F(x[i-1],Euler_improved[i-1])
yc=Euler_improved[i-1]+h*F(x[i],yp)
Euler_improved.append(0.5*(yp+yc))
同样，在0-5的定义域内求解：
可以得到结论：改进的Euler格式精度比前向Euler高，但稳定性依旧不高。
Euler两步格式
-
无论是前向还是后向Euler，亦或是改进的Euler格式，都只用到了单步信息，即要得到 y n + 1 y_{n+1} yn+1的值，只需要求得前一步的结果即可。Euler两步格式即再向前考虑多一步，除 y n y_n yn外，还包括 y n − 1 y_{n-1} yn−1的信息。两步法的优点是调用了两个节点信息，能通过较少的计算量达到较高的精度。 我们可以得到预测—校正系统的Euler两步格式如下：
预 测 ： y n + 1 ˉ = y n − 1 + 2 h f ( x n , y n ) 校 正 ： y n + 1 = y n + h 2 [ f ( x n , y n ) + f ( x n + 1 , y n + 1 ˉ ) ] \begin{aligned} &预测：\bar{y_{n+1}}=y_{n-1}+2h f(x_n,y_n)\\ &校正：y_{n+1} =y_n+\frac{h}{2} [f(x_n,y_n)+f(x_{n+1},\bar{y_{n+1}})] \end{aligned} 预测：yn+1ˉ=yn−1+2hf(xn,yn)校正：yn+1=yn+2h[f(xn,yn)+f(xn+1,yn+1ˉ)]
此时的预测公式的局部截断误差是与校正公式匹配的，在改进的Euler格式中两者不匹配，分别是一阶精度和二阶精度。下面验证Euler两步格式下的预测公式的局部截断误差(将 y ( x n + 1 ) y(x_{n+1}) y(xn+1)与 y n − 1 y_{n-1} yn−1在 x n x_n xn处泰勒展开)：T n + 1 = y ( x n + 1 ) − y n − 1 − 2 h f ( x n , y n ) = y ( x n + 1 ) − y n − 1 − 2 h y n ′ = y n + h y n ′ + h 2 2 y n ′ ′ + h 3 6 y n ′ ′ ′ + O ( h 4 ) − ( y n − h y n ′ + 1 2 h 2 y n ′ ′ − 1 6 h 3 y n ′ ′ ′ + O ( h 4 ) ) − 2 h y n ′ ≈ 1 3 h 3 y ′ ′ ′ ( x n ) \begin{aligned} T_{n+1}&=y(x_{n+1})-y_{n-1}-2hf(x_n,y_n)\\ &= y(x_{n+1})-y_{n-1}-2h y^{'}_n \\ &= y_n+h y^{'}_n+\frac{h^2}{2} y^{''}_n+\frac{h^3}{6} y^{'''}_n+O(h^4)-(y_n-h y^{'}_n+\frac{1}{2}h^2 y^{''}_n-\frac{1}{6}h^3 y^{'''}_n+O(h^4))-2h y^{'}_n\\ &\approx \frac{1}{3}h^3y^{'''}(x_n) \end{aligned} Tn+1=y(xn+1)−yn−1−2hf(xn,yn)=y(xn+1)−yn−1−2hyn′=yn+hyn′+2h2yn′′+6h3yn′′′+O(h4)−(yn−hyn′+21h2yn′′−61h3yn′′′+O(h4))−2hyn′≈31h3y′′′(xn)
代码实现
# Euler 两步法
# Author: Junno
# Date: 2022-04-14
Euler_dual_step=[]
# use Euler_improved method to calculate y1
for i in range(len(x)):
if x[i]==initial_x:
Euler_dual_step.append(initial_y)
else:
if i==1:
yp=Euler_dual_step[i-1]+h*F(x[i-1],Euler_dual_step[i-1])
yc=Euler_dual_step[i-1]+h*F(x[i],yp)
Euler_dual_step.append(0.5*(yp+yc))
else:
yn_partial=F(x[i-1],Euler_dual_step[i-1])
yp=Euler_dual_step[i-2]+2*h*yn_partial
yc=Euler_dual_step[i-1]+h/2*(yn_partial+F(x[i],yp))
Euler_dual_step.append(yc)
同样，在0-5的定义域内求解：
预测-校正-双改进
- 对于预测公式和校正公式，其各自的局部截断误差为：
预 测 ： y ( x n + 1 ) − y n + 1 ˉ ≈ 1 3 h 3 y ′ ′ ′ ( x n ) 校 正 ： y ( x n + 1 ) − y n + 1 ≈ − 1 12 h 3 y ′ ′ ′ ( x n ) \begin{aligned} &预测：y(x_{n+1})-\bar{y_{n+1}}\approx \frac{1}{3}h^3y^{'''}(x_n)\\ &校正：y(x_{n+1})-y_{n+1} \approx -\frac{1}{12}h^3y^{'''}(x_n) \end{aligned} 预测：y(xn+1)−yn+1ˉ≈31h3y′′′(xn)校正：y(xn+1)−yn+1≈−121h3y′′′(xn)
即可以继续推得：
y
n
+
1
−
y
n
+
1
ˉ
≈
5
12
h
3
y
′
′
′
(
x
n
)
h
3
y
′
′
′
(
x
n
)
=
12
5
(
y
n
+
1
−
y
n
+
1
ˉ
)
\begin{aligned} y_{n+1}-\bar{y_{n+1}} & \approx \frac{5}{12}h^3y^{'''}(x_n)\\ h^3y^{'''}(x_n) &= \frac{12}{5}(y_{n+1}-\bar{y_{n+1}}) \end{aligned}
yn+1−yn+1ˉh3y′′′(xn)≈125h3y′′′(xn)=512(yn+1−yn+1ˉ)
导出事后估计式：
{
y
(
x
n
+
1
)
−
y
n
+
1
ˉ
≈
4
5
(
y
n
+
1
−
y
n
+
1
ˉ
)
y
(
x
n
+
1
)
−
y
n
+
1
≈
−
1
5
(
y
n
+
1
−
y
n
+
1
ˉ
)
\begin{cases} y(x_{n+1})-\bar{y_{n+1}}\approx \frac{4}{5}(y_{n+1}-\bar{y_{n+1}})\\ y(x_{n+1})-y_{n+1} \approx -\frac{1}{5}(y_{n+1}-\bar{y_{n+1}}) \end{cases}
{y(xn+1)−yn+1ˉ≈54(yn+1−yn+1ˉ)y(xn+1)−yn+1≈−51(yn+1−yn+1ˉ)
通过对得到的预测值和校正值进行误差补偿，可以进一步改善精度。
给出预测-校正-双改进的两步Euler格式：记
p
n
p_n
pn和
c
n
c_n
cn表示第n步的预测值和校正值
预
测
：
p
n
+
1
=
y
n
−
1
+
2
h
y
n
′
改
进
：
m
n
+
1
=
p
n
+
1
+
4
5
(
c
n
−
p
n
)
计
算
：
m
n
+
1
′
=
f
(
x
n
+
1
,
m
n
+
1
)
校
正
：
c
n
+
1
=
y
n
+
h
2
(
m
n
+
1
′
+
y
n
′
)
改
进
：
y
n
+
1
=
c
n
+
1
−
1
5
(
c
n
+
1
−
p
n
+
1
)
\begin{aligned} &预测：p_{n+1}=y_{n-1}+2h y^{'}_n\\ &改进：m_{n+1}=p_{n+1}+\frac{4}{5}(c_n-p_n)\\ &计算：m^{'}_{n+1}=f(x_{n+1},m_{n+1})\\ &校正：c_{n+1}=y_n+\frac{h}{2}(m^{'}_{n+1}+y^{'}_n)\\ &改进：y_{n+1}=c_{n+1}-\frac{1}{5}(c_{n+1}-p_{n+1})\\ \end{aligned}
预测：pn+1=yn−1+2hyn′改进：mn+1=pn+1+54(cn−pn)计算：mn+1′=f(xn+1,mn+1)校正：cn+1=yn+2h(mn+1′+yn′)改进：yn+1=cn+1−51(cn+1−pn+1)
代码实现
# 预测-校正-双改进的两步Euler格式
# Author: Junno
# Date: 2022-04-14
Euler_dual_step_dual_compensate=[]
# p1=c1=0
predictions=[0 for _ in range(2)]
corrections=[0 for _ in range(2)]
for i in range(len(x)):
if x[i]==initial_x:
Euler_dual_step_dual_compensate.append(initial_y)
else:
if i==1:
yp=Euler_dual_step_dual_compensate[i-1]+h*F(x[i-1],Euler_dual_step_dual_compensate[i-1])
yc=Euler_dual_step_dual_compensate[i-1]+h*F(x[i],yp)
Euler_dual_step_dual_compensate.append(0.5*(yp+yc))
else:
yn_partial=F(x[i-1],Euler_dual_step_dual_compensate[i-1])
# predict
p_=Euler_dual_step_dual_compensate[i-2]+2*h*yn_partial
predictions.append(p_)
# compensate
m_=p_+4/5*(corrections[i-1]-predictions[i-1])
# calculate
m_partial=F(x[i],m_)
# correct
c_=Euler_dual_step_dual_compensate[i-1]+h/2*(m_partial+yn_partial)
corrections.append(c_)
# compensate
y_=c_-1/5*(c_-p_)
Euler_dual_step_dual_compensate.append(y_)
同样，在0-5的定义域内求解：
果然，要想精度高，计算量必须要达到！
Runge-Kutta方法
- 单步法的思想是从点 ( x n , y n ) (x_n,y_n) (xn,yn)出发，以某一斜率达到点 ( x n + 1 , y n + 1 ) (x_{n+1},y_{n+1}) (xn+1,yn+1)，欧拉法及其各种变形改进最高精度为2阶。基于Taylor级数法，Runge和Kutta使用函数值的计算近似逼近高阶导数，具体过程是求解满足给定精度要求的各项系数。（详见P113）
四阶Runge-Kutta方法
{ y ( x n + 1 ) = y n + h 6 ( K 1 + 2 K 2 + 2 K 3 + k 4 ) K 1 = f ( x n , y n ) K 2 = f ( x n + h 2 , y n + h 2 K 1 ) K 3 = f ( x n + h 2 , y n + h 2 K 2 ) K 4 = f ( x n + h , y n + h K 3 ) \begin{cases} y(x_{n+1}) &=y_n+\frac{h}{6}(K_1+2K_2+2K_3+k_4)\\ K_1 &= f(x_n,y_n)\\ K_2 &= f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)\\ K_3 &= f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2)\\ K_4 &= f(x_n+h,y_n+h K_3) \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧y(xn+1)K1K2K3K4=yn+6h(K1+2K2+2K3+k4)=f(xn,yn)=f(xn+2h,yn+2hK1)=f(xn+2h,yn+2hK2)=f(xn+h,yn+hK3)
四阶Runge-Kutta方法在每一步中需要计算四次函数值f，其截断误差为 O ( h 5 ) O(h^5) O(h5)，属于一种高精度的单步数值微分法。
注意的是，Runge-Kutta方法是基于Taylor级数展开而来，需要所求解具有足够好的光滑性，否则效果会变差。
代码实现
# 四阶Runge-Kutta方法
# Author: Junno
# Date: 2022-04-14
Runge_Kutta_four_order=[]
for i in range(len(x)):
if x[i]==initial_x:
Runge_Kutta_four_order.append(initial_y)
else:
K_1=F(x[i-1],Runge_Kutta_four_order[i-1])
K_2=F(x[i-1]+h/2,Runge_Kutta_four_order[i-1]+h/2*K_1)
K_3=F(x[i-1]+h/2,Runge_Kutta_four_order[i-1]+h/2*K_2)
K_4=F(x[i-1]+h,Runge_Kutta_four_order[i-1]+h*K_3)
Runge_Kutta_four_order.append(Runge_Kutta_four_order[i-1]+h/6*(K_1+2*(K_2+K_3)+K_4))
同样比较0-5范围0.1步长下的结果：
可以看到四阶Runge-Kutta方法基本与解析解曲线重合，精度很高！
给出不同步长下四阶Runge-Kutta方法结果与解析解的无穷范数：
给出0.5步长下各方法的结果：
其他高精度方法
- 这里列出其他高精度的线性多步法，顾名思义，即用到了前面的多步信息，可以有效提高精度。不做证明，仅做代码实用性分析！按需取用即可
Adams 预测-校正模型
预 测 ： y ( x n + 1 ) ˉ = y n + h 24 ( 55 f n − 59 f n − 1 + 37 f n − 2 − 9 f n − 3 ) f n + 1 ˉ = f ( x n + 1 , y ( x n + 1 ) ˉ ) 校 正 ： y ( x n + 1 ) = y n + h 24 ( 9 f n + 1 ˉ + 19 f n − 5 f n − 1 + f n − 2 ) f n + 1 = f ( x n + 1 , y ( x n + 1 ) ) \begin{aligned} &预测：\bar{y(x_{n+1})}=y_n+\frac{h}{24}(55f_n-59f_{n-1}+37 f_{n-2}-9 f_{n-3})\\ &\bar{f_{n+1}}=f(x_{n+1},\bar{y(x_{n+1})})\\ &校正：y(x_{n+1})=y_n+\frac{h}{24}(9 \bar{f_{n+1}} +19f_n-5f_{n-1}+f_{n-2})\\ &f_{n+1}=f(x_{n+1},y(x_{n+1}))\\ \end{aligned} 预测：y(xn+1)ˉ=yn+24h(55fn−59fn−1+37fn−2−9fn−3)fn+1ˉ=f(xn+1,y(xn+1)ˉ)校正：y(xn+1)=yn+24h(9fn+1ˉ+19fn−5fn−1+fn−2)fn+1=f(xn+1,y(xn+1))
- 以上方法类似于改进的Euler格式，其精度为四阶，但由于需要前4步的f信息，所以可以用精度同为四阶的Runge-Kutto方法先计算出前四个f值。
代码实现
# Adams 预测-校正系统
# Author: Junno
# Date: 2022-04-14
def Runge_Kutta_4(x,y,h,F):
# x and y are results of the last step
K_1=F(x,y)
K_2=F(x+h/2,y+h/2*K_1)
K_3=F(x+h/2,y+h/2*K_2)
K_4=F(x+h,y+h*K_3)
return y+h/6*(K_1+2*(K_2+K_3)+K_4)
f_n=[]
Adams_pc=[]
for i in range(len(x)):
if x[i]==initial_x:
f_n=[F(initial_x,initial_y)]
Adams_pc.append(initial_y)
else:
if i<4:
yn=Runge_Kutta_4(x[i-1],Adams_pc[i-1],h,F)
Adams_pc.append(yn)
f_n.append(F(x[i],yn))
else:
p_=Adams_pc[i-1]+h/24*(55*f_n[i-1]-59*f_n[i-2]+37*f_n[i-3]-9*f_n[i-4])
f_n_=F(x[i],p_)
c_=Adams_pc[i-1]+h/24*(9*f_n_+19*f_n[i-1]-5*f_n[i-2]+f_n[i-3])
Adams_pc.append(c_)
f_n.append(F(x[i],c_))
在0-5范围步长0.1比较：
不同步长下四阶Runge-Kutto方法与Adams预测-校正系统的无穷范数误差对比：
Miline格式与Hamming格式
-
Miline格式： y n + 1 = y n − 3 + 4 h 3 ( 2 y n ′ − y n − 1 ′ + 2 y n − 2 ′ ) y_{n+1}=y_{n-3}+\frac{4h}{3}(2y^{'}_n-y^{'}_{n-1}+2y^{'}_{n-2}) yn+1=yn−3+34h(2yn′−yn−1′+2yn−2′)
其局部截断误差为： y ( x n + 1 ) − y n + 1 ≈ 14 45 h 5 y n ( 5 ) y(x_{n+1})-y_{n+1} \approx \frac{14}{45}h^5y^{(5)}_n y(xn+1)−yn+1≈4514h5yn(5)
-
Hamming格式： y n + 1 = 1 8 ( 9 y n − y n − 2 ) + 3 h 8 ( y n + 1 ′ + 2 y n ′ − y n − 1 ′ ) y_{n+1}=\frac{1}{8}(9 y_n-y_{n-2})+\frac{3h}{8}(y^{'}_{n+1}+2y^{'}_{n}-y^{'}_{n-1}) yn+1=81(9yn−yn−2)+83h(yn+1′+2yn′−yn−1′)
其局部截断误差为： y ( x n + 1 ) − y n + 1 ≈ − 1 40 h 5 y n ( 5 ) y(x_{n+1})-y_{n+1} \approx -\frac{1}{40}h^5y^{(5)}_n y(xn+1)−yn+1≈−401h5yn(5)
-
同样，可以将Miline与Hamming格式组成预测-校正系统：
预 测 ： y ( x n + 1 ) ˉ = y n − 3 + 4 h 3 ( 2 y n ′ − y n − 1 ′ + 2 y n − 2 ′ ) f n + 1 ˉ = f ( x n + 1 , y ( x n + 1 ) ˉ ) 校 正 ： y n + 1 = 1 8 ( 9 y n − y n − 2 ) + 3 h 8 ( f n + 1 ˉ + 2 y n ′ − y n − 1 ′ ) f n + 1 = f ( x n + 1 , y ( x n + 1 ) ) \begin{aligned} &预测：\bar{y(x_{n+1})}=y_{n-3}+\frac{4h}{3}(2y^{'}_n-y^{'}_{n-1}+2y^{'}_{n-2})\\ &\bar{f_{n+1}}=f(x_{n+1},\bar{y(x_{n+1})})\\ &校正：y_{n+1}=\frac{1}{8}(9 y_n-y_{n-2})+\frac{3h}{8}(\bar{f_{n+1}}+2y^{'}_{n}-y^{'}_{n-1})\\ &f_{n+1}=f(x_{n+1},y(x_{n+1}))\\ \end{aligned} 预测：y(xn+1)ˉ=yn−3+34h(2yn′−yn−1′+2yn−2′)fn+1ˉ=f(xn+1,y(xn+1)ˉ)校正：yn+1=81(9yn−yn−2)+83h(fn+1ˉ+2yn′−yn−1′)fn+1=f(xn+1,y(xn+1))
代码实现
# Miline-Hamming p-c
# Author: Junno
# Date: 2022-04-14
f_n=[]
Miline_Hamming_pc=[]
for i in range(len(x)):
if x[i]==initial_x:
f_n=[F(initial_x,initial_y)]
Miline_Hamming_pc.append(initial_y)
else:
if i<4:
yn=Runge_Kutta_4(x[i-1],Miline_Hamming_pc[i-1],h,F)
Miline_Hamming_pc.append(yn)
f_n.append(F(x[i],yn))
else:
p_=Miline_Hamming_pc[i-4]+4*h/3*(2*f_n[i-1]-f_n[i-2]+2*f_n[i-3])
f_n_=F(x[i],p_)
c_=1/8*(9*Miline_Hamming_pc[i-1]-Miline_Hamming_pc[i-3])+3*h/8*(f_n_+2*f_n[i-1]-f_n[i-2])
Miline_Hamming_pc.append(c_)
f_n.append(F(x[i],c_))
不同步长下三种四阶精度方法的无穷范数误差对比如下，可以看到随着步长增加，总体误差呈现抖动上升趋势。
Miline-Hamming 双改进预测-校正系统
- 类似于前面分析的，同样可以得到Miline-Hamming 双改进预测-校正系统的形式：
预 测 ： p n + 1 = y n − 3 + 4 h 3 ( 2 y n ′ − y n − 1 ′ + 2 y n − 2 ′ ) 改 进 ： m n + 1 = p n + 1 + 112 121 ( c n − p n ) 计 算 ： m n + 1 ′ = f ( x n + 1 , m n + 1 ) 校 正 ： c n + 1 = 1 8 ( 9 y n − y n − 2 ) + 3 h 8 ( m n + 1 ′ + 2 y n ′ − y n − 1 ′ ) 改 进 ： y n + 1 = c n + 1 − 9 121 ( c n + 1 − p n + 1 ) \begin{aligned} &预测：p_{n+1}=y_{n-3}+\frac{4h}{3}(2y^{'}_n-y^{'}_{n-1}+2y^{'}_{n-2})\\ &改进：m_{n+1}=p_{n+1}+\frac{112}{121}(c_n-p_n)\\ &计算：m^{'}_{n+1}=f(x_{n+1},m_{n+1})\\ &校正：c_{n+1}=\frac{1}{8}(9 y_n-y_{n-2})+\frac{3h}{8}(m^{'}_{n+1}+2y^{'}_{n}-y^{'}_{n-1})\\ &改进：y_{n+1}=c_{n+1}-\frac{9}{121}(c_{n+1}-p_{n+1})\\ \end{aligned} 预测：pn+1=yn−3+34h(2yn′−yn−1′+2yn−2′)改进：mn+1=pn+1+121112(cn−pn)计算：mn+1′=f(xn+1,mn+1)校正：cn+1=81(9yn−yn−2)+83h(mn+1′+2yn′−yn−1′)改进：yn+1=cn+1−1219(cn+1−pn+1)
代码实现
# Author: Junno
# Date: 2022-04-14
# Miline-Hamming 双改进预测-校正系统
Miline_Hamming_dual_compensate=[]
predictions=[0 for _ in range(4)]
corrections=[0 for _ in range(4)]
f_n=[]
for i in range(len(x)):
if x[i]==initial_x:
Miline_Hamming_dual_compensate.append(initial_y)
f_n.append(F(initial_x,initial_y))
else:
if i<4:
yn=Runge_Kutta_4(x[i-1],Miline_Hamming_dual_compensate[i-1],h,F)
Miline_Hamming_dual_compensate.append(yn)
f_n.append(F(x[i],yn))
else:
# predict
p_=Miline_Hamming_dual_compensate[i-4]+4*h/3*(2*f_n[i-1]-f_n[i-2]+2*f_n[i-3])
predictions.append(p_)
# compensate
m_=p_+112/121*(corrections[i-1]-predictions[i-1])
# calculate
m_partial=F(x[i],m_)
# correct
c_=1/8*(9*Miline_Hamming_dual_compensate[i-1]-Miline_Hamming_dual_compensate[i-3])+3*h/8*(m_partial+2*f_n[i-1]-f_n[i-2])
corrections.append(c_)
# compensate
y_=c_-9/121*(c_-p_)
Miline_Hamming_dual_compensate.append(y_)
f_n.append(F(x[i],y_))
不同步长下三种四阶精度方法的无穷范数误差对比：
几种方法也不是说越复杂精度越高，具体还是要尝试比较，如对于
y
′
=
y
y^{'}=y
y′=y，即
y
=
e
x
y=e^x
y=ex，上面几种四阶精度方法的误差如下：
2176
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
