数值微分方法：欧拉法、改进欧拉、双校正预测、四阶Runge-Kutta、Adams等

数值分析-数值微分

• 前言
• 敲黑板
• 欧拉前向差商法(欧拉折线法) 代码实现 局部截断误差 后向差分欧拉 代码实现
• 梯形平均方法
• 改进的Euler格式 代码实现 Euler两步格式 代码实现 预测-校正-双改进 代码实现 Runge-Kutta方法 四阶Runge-Kutta方法 代码实现 其他高精度方法 Adams 预测-校正模型 代码实现 Miline格式与Hamming格式 代码实现 Miline-Hamming 双改进预测-校正系统 代码实现

前言

​   上一篇在数值积分中我们介绍了插值型求积公式以及如何分析对应的代数精度、稳定性和收敛性。当然啦，有数值积分就有数值微分，在以往的学习中我们可能已经掌握了很多解析解法，但是解析解法只能用来求解一些具有特定类型的方程，大部分方程都是需要数值求解的。 按照教材内容，我们着重点是一阶常微分方程初值问题的数值微分方法。

   具体内容详见 《数值分析》 李庆扬、王能超、易大义（华中科大出版社）第五章

敲黑板

• { y′=f(x,y)y(x0​)=y0​​

欧拉前向差商法(欧拉折线法)

• 设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​)+2h2​y′′(xn​)+O(h3)]−[yn​+hf(xn​,yn​)]=2h2​y′′(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​=−2h2​y′′(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​)]
• Tn+1​​=y(xn+1​)−yn+1​=y(xn+1​
  免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删