数值分析——数值微分（欧拉法、改进欧拉格式、双校正预测模型、四阶Runge-Kutta、Adams等）

小刀来啦

已于 2022-05-23 15:35:10 修改

阅读量4.4k

点赞数 6

分类专栏：数值分析笔记文章标签： python 算法数据分析

于 2022-04-12 11:34:56 首次发布

本文链接：https://blog.csdn.net/weixin_40519529/article/details/124116413

版权

数值分析-数值微分

前言
敲黑板
梯形平均方法

前言

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

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

敲黑板

一阶方程的初值问题：

$\begin{cases} y^{'}=f(x,y)\\ y(x_0)=y_0 \end{cases}$

数值解法就是求出 $y (x)$ 在一系列等距节点 $x_1,x_2,...,x_n$ 处的值。对于初值问题数值解法一般采取步进式求解过程，即通过前面求解的节点结果推出后续节点的值。

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

设h为等距节点分隔间距，使用向前差分代替该点的导数值，即：
$y^{'}(x_0) \approx \frac{y(x_1)-y(x_0)}{h}$

即在已知初值点 $y(x_0)$ 的情况下，后续节点的计算方式为：
$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)$

代码实现

# 前向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]$ 范围内还比较拟合，之后显著偏离原来的曲线，即前向差分欧拉方法的误差还是很大的。

局部截断误差

记 $y (x)$ 是给定初值问题的准确解，假设在第n步计算准确的情况下，即 $y_n=y(x_n)$ ，考虑截断误差 $T_{n+1}=y(x_{n+1})-y_{n+1}$ ，称之为局部截断误差。
若某算法的局部截断误差为 $O(h^{p+1})$ ，则称该算法具有p阶精度。
一般在节点 $x_{n+1}$ 处采用Taylor展开分析局部截断误差，来看下欧拉前向差分的局部截断误差：

$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)$
即前向欧拉方法具有1阶精度。

后向差分欧拉

对于前向差分我们有， $y_{n+1}=y_n+hf(x_n,y_n)$ ，自然地，后向差分的Euler公式即为 $y_{n+1}=y_n+hf(x_n,y_{n+1})$ ，但是这类计算属于隐式计算，需要不断迭代，不易求解。若迭代收敛，则 $y_{n+1}=\lim_{k \rightarrow \infty} y_{n+1}^k$ 。
对于后向欧拉公式，可以证得（教材P108）其局部截断误差为 $T_{n+1}=-\frac{h^2}{2} y^{''}(x_n)+O(h^3)$

代码实现

# 前向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))