数值分析(11):常微分方程的数值解法之Euler法

本文介绍了常微分方程初值问题的存在唯一性定理及其充分条件,并详细阐述了Euler方法,包括显式Euler、隐式Euler和梯形方法。此外,讨论了隐式方程的迭代求解策略和预估-校正方法。同时,对显式和隐式单步方法进行了误差分析,探讨了整体截断误差和局部截断误差的概念及阶的定义。
摘要由CSDN通过智能技术生成

1. 引言

我们通常遇到的一阶常微分方程初值问题,如下所示:

在这里插入图片描述

如何证明上面的微分方程存在解,且解唯一呢?则需要以下的存在唯一性定理:

在这里插入图片描述

上面定理中的Lipschitz条件,其定义如下:
在这里插入图片描述
这个定义很难验证,通常使用以下的充分条件来判断函数是否满足Lipschitz条件,即:
在这里插入图片描述

那么对于上面存在唯一性定理中,同样可以用充分条件来判断 f f f是否满足对 y y y的Lipschitz条件,即:
在这里插入图片描述

讲述了存在唯一性定理后,接下来来看如何解决无法写出解析表达式的微分方程。一种最朴素的方式就是用直线去逼近,把积分区间分成很多段,每一段都使用直线逼近原曲线,用图像表示即为:

在这里插入图片描述

为了方便起见,把 y ( x ) y(x) y(x) x n x_n xn处的精确值记为 y ( x n ) y(x_n) y(xn), 其近似值 y n y_n yn表示。

在求微分方程的数值解时,做以下的归类:

在这里插入图片描述

由于是近似,自然要分析误差,误差分为两种,整体截断误差局部截断误差,这两个定义如下:

在这里插入图片描述

两者的区别是:整体截断误差是存在误差累积的,累积前n步的误差;局部截断误差是只计算当前第n步的误差。

2. Euler方法

Euler方法就是引言中用直线近似来求解微分方程的方法。

2.1 显式Euler方法

所谓“显式”,就是说 y n + 1 y_{n+1} yn+1只在等式左边出现。只要能够求出 y n + 1 y_{n+1} yn+1,根据区间分段结果,直接用直线相连即可。
在这里插入图片描述
当然上面这个计算公式可以用三种方法得到,分别是:数值积分方法、Taylor展开方法、数值微分方法,即:
在这里插入图片描述
在这里插入图片描述

上面只用到了 y n y_n yn即可求解 y n + 1 y_{n+1} yn+1,因此是单步法。

2.2 隐式Euler方法

隐式法就是等式右边出现 y n + 1 y_{n+1} yn+1,可以看到显式和隐式的区别在于切线的斜率选取的问题,选取点 ( 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.3 梯形方法

如果把显式隐式的切线斜率公式相加取均值,就得到了梯形公式,即:

在这里插入图片描述

3. 隐式方程的迭代求解

梯形方法和隐式Euler方法都是隐式公式,这就涉及到一个问题,等式左右两边都有未知数 y n + 1 y_{n+1} yn+1,这个等式怎么解?

如果是线性的等式,自然可以通过移项转换变成显式的等式,如果是非线性的等式,可以通过前面《数值分析(3):线性代数方程组的迭代解法》中所述的迭代解法求解。

对于隐式Euler法,即用下式迭代:

在这里插入图片描述
上面的迭代公式中,第一步迭代还是需要用显式公式进行计算,后面的迭代采用隐式Euler法。很自然地,需要考虑上述迭代法的收敛性,即:

在这里插入图片描述

对于梯形公式,有迭代公式:

在这里插入图片描述

迭代收敛条件为:
在这里插入图片描述

4. 预估-校正方法

为了消除迭代,出现了预估-校正的方法,先给出粗糙估计然后再给出稍精确的求解,这是微分方程数值解常用方法。

比如改进的Euler方法,即:

在这里插入图片描述

这公式称为改进的Euler公式,其精度比 Euler 公式好,比梯形公式稍差些,但是,它是显式方法。

5. 误差分析

5.1 显式单步方法误差分析

显式单步方法的通式如下:
在这里插入图片描述
正如前述,其整体截断误差为通过 x 0 , y 0 x_0,y_0 x0,y0,计算 y 1 , y 2 , . . . , y n y_1,y_2,...,y_n y1,y2,...,yn,最后再计算整体截断误差 e n = y ( x n ) − y n e_n=y(x_n)-y_n en=y(xn)yn,注意带括号的是精确值,带下标的是估计值

当然整体截断误差与整个计算中每步情况有关,但一般求得较为困难,因此先考虑一步的误差。

显式单步法的局部截断误差为:

在这里插入图片描述

单步法的阶,用来衡量计算方法的精度:
在这里插入图片描述

主局部截断误差,它就相当于泰勒展开中的提取余项的操作,对于 p p p阶方法中的 p + 1 p+1 p+1阶无穷小,它进一步将其细化为 p + 1 p+1 p+1次的项,和 p + 2 p+2 p+2阶无穷小,而 p + 1 p+1 p+1次的项 就是主局部截断误差

在这里插入图片描述

5.2 隐式单步方法误差分析

在显式单步方法误差分析中,我们分析了局部截断误差表达式单步法的阶主局部截断误差。那么对于隐式法,完全可以类似的定义。

比如对于隐式Euler法:
在这里插入图片描述
局部截断误差表达式为:

T n + 1 = y ( x n + 1 ) − [ y ( x n ) + h f ( x n + 1 , y ( x n + 1 ) ) ] T_{n+1}=y\left( x_{n+1} \right) -\left[ y\left( x_n \right) +hf\left( x_{n+1},y\left( x_{n+1} \right) \right) \right] Tn+1=y(xn+1)[y(xn)+hf(xn+1,y(xn+1))]

对其进行泰勒展开化简后得:
T n + 1 = y ( x n + 1 ) − [ y ( x n ) + h f ( x n + 1 , y ( x n + 1 ) ) ] = y ( x n + h ) − [ y ( x n ) + h y ′ ( x n + h ) ] = y ( x n ) + h y ′ ( x n ) + h 2 2 ! y ′ ′ ( x n ) + o ( h 3 ) − y ( x n ) − h ( y ′ ( x n ) + h y ′ ′ ( x n ) + o ( h 2 ) ) = − h 2 2 y ′ ′ ( x n ) + o ( h 3 ) T_{n+1}=y\left( x_{n+1} \right) -\left[ y\left( x_n \right) +hf\left( x_{n+1},y\left( x_{n+1} \right) \right) \right] \\ =y\left( x_n+h \right) -\left[ y\left( x_n \right) +hy'\left( x_n+h \right) \right] \\ =y\left( x_n \right) +hy'\left( x_n \right) +\frac{h^2}{2!}y''\left( x_n \right) +o\left( h^3 \right) -y\left( x_n \right) -h\left( y'\left( x_n \right) +hy''\left( x_n \right) +o\left( h^2 \right) \right) \\ =-\frac{h^2}{2}y''\left( x_n \right) +o\left( h^3 \right) Tn+1=y(xn+1)[y(xn)+hf(xn+1,y(xn+1))]=y(xn+h)[y(xn)+hy(xn+h)]=y(xn)+hy(xn)+2!h2y(xn)+o(h3)y(xn)h(y(xn)+hy(xn)+o(h2))=2h2y(xn)+o(h3)

因此单步法的阶为:一阶方法
主局部截断误差为: − h 2 2 y ′ ′ ( x n ) -\frac{h^2}{2}y''\left( x_n \right) 2h2y(xn)

其他的隐式方法完全可以类比推导,这里就省略不写了。

求积公式的余项和数值解的局部截断误差主局部截断误差有如下关系,以下关系也是很容易理解的,因为微分方程的数值解就是可以通过求积公式导出,自然它们的余项具有非常相近的形式:
在这里插入图片描述

参考文献:

关治,陆金甫《数值方法》

  • 7
    点赞
  • 65
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
好的,您可以参考以下的步骤: Step 1: 准备工作 在使用Euler求解常微分方程的初值问题之前,您需要先了解一些预备知识,包括: - 常微分方程的基本概念; - 常微分方程初值问题的定义; - 数值解法的基本概念和误差分析方。 另外,您需要熟悉Matlab编程环境,并了解基本的Matlab语和函数库。 Step 2: 写出Euler数值格式 Euler常微分方程数值解法中最简单的一种,其数值格式为: y_{n+1}=y_n+f(y_n,t_n)\times h 其中,y_n和t_n分别表示自变量和因变量的取值,h为步长,f(y_n,t_n)为方程左侧y'的函数值,即: y'=f(y,t) Step 3: 编写Matlab代码 基于以上的Euler数值格式,可以编写出Matlab代码,如下所示: function [t,y]=euler(f,tspan,y0,h) % f为y'的函数句柄 % tspan为区间[t0,t1] % y0为初值 % h为步长 t=tspan(1):h:tspan(2); % 根据步长h生成时间节点t y=zeros(size(t)); % 初始化y向量 y(1)=y0; % 设定初值 for i=1:length(t)-1 y(i+1)=y(i)+f(y(i),t(i))*h; % 根据Euler更新y(i+1)的值 end Step 4: 编写示例程序并测试 下面是一个求解微分方程y'=-2y的示例程序,并使用Euler求解其初值问题: f=@(y,t)-2*y; % 定义函数句柄f(y,t)=-2*y tspan=[0,3]; % 区间[t0,t1]=[0,3] y0=1; % 初值y(0)=1 h=0.1; % 步长h=0.1 [t,y]=euler(f,tspan,y0,h); % 调用euler函数求解 % 绘制y关于t的图像 plot(t,y,'-o'); xlabel('t'); ylabel('y(t)'); title('Euler Method for y''=-2y'); 运行这个程序后,可以得到Euler数值解结果,并绘制出相应的y-t图像。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值