Euler方法及其改进方法在线视频

大家好

我们在上节课学习了如何离散

微分方程的初值问题

今天我们来学习Euler方法及其改进方法

首先 我们来看一下Euler方法的求解公式

所谓Euler方法

是以这样一个差分方程的初值问题

yn+1等于yn加上hf(xn yn)

y₀=y(a）

它的解作为微分方程初值问题的数值解

即y(xn)近似的等于yn

我们在前面讲了单步法

它的一般形式为

yn+1等于yn加上h乘上Φ(xn yn

xn+1 yn+1)

对于Euler方法来说

我们这里的Φ取的即为

f(xn yn)

即xn yxn这一点

它的切线斜率的一个近似

我们用差分方程的解

作为微分方程解的一个近似

从图形上来看

就是对于微方程解的这条曲线

用一条分段的折线来逼近它

我们来看一个例子

用Euler方法求解下列方程的数值解

即dy dx等于x根下y

当x大于等于0 小于等于1的时候

初值条件是y₀=1

这里我们的步长分别取0.5和0.25

我们知道Euler公式

它的迭代格式为yn+1等于yn加上hf(xn yn)

我们把f的表达式代进来

即为yn加上hxn根下yn

当步长h取0.5的时候

我们

可以得到y0.5

近似的就等于y₁等于y₀加上

hx₀根下y₀

把y₀=1代进去

我们就可以得到y0.5近似的等于1.0

求出y₁之后

我们再把它代入迭代格式的右端

可以计算y₂等于y₁加上hx₁根下y₁

把y₁的值1.0代进来之后 我们可以求出

y₂是1.25

因此 这时候y在1.0的

值 我们就近似的等于1.25

这是步长h取0.5的时候

我们可以求出在两个节点0.5和1.0它的近似值

当步长h取0.25的时候

我们可以分别代入这个公式求出

y(0.25)近似的等于

y₁=1

y(0.5) 近似的等于y₂

等于1.0625

y的0.75

近似的等于y₃

计算出来 是1.191347

而y的1.0近似的是等于y₄

这个值是1.39600

如果我们把它的解析解

和它的数值解画在图形上 我们可以看出

这条黄色的曲线 就是它的精确解

即y=(1+x²/4)²

而当步长h等于0.5的时候

近似解的曲线 就为这条绿色的曲线

当步长h等于0.25的时候

我们的数值解的曲线就为这条黑色的曲线

从这个图形上我们可以看出

当步长变小即由0.5变为0.25的时候

我们这个数值解对真解的逼近更精确

接下来 我们来看向后Euler方法 它的

计算公式为 yn+1等于yn加上hf(xn+1 yn+1)

初值 y₀=ya

我们知道 这个向后Euler方法是一个隐式的公式

当f是一个非线性函数的时候

我们通常需要来求解一个非线性方程

所以在实际计算当中

一般是要用迭代法来进行求解

这时候我们选的初值

一般用显格式 Euler格式给出初值

然后代入向后Euler格式的右端

不断进行迭代

即可以得到yn+1

它的k+1次迭代值等于 yn加上hf(xn+1 yn+1 k

第k次的迭代值

对于这样的一个迭代求解格式

我们这里有两个问题

向后Euler迭代序列是否是收敛的

向后Euler迭代序列在什么条件下收敛

我们来看一下

我们要讨论这个

迭代序列的收敛性

我们就来比较一下yn+1的(k+1)次减掉yn+1的(k) 次

它的差的绝对值

我们把这个迭代格式代入

就可以得到

它等于h乘上f(xn+1 yn+1(k)) 减掉

f(xn+1 yn+1(k-1))

我们知道右端函数f关于y满足Lipschitz条件

因此它小于等于h乘上这个Lipschitz常数大L

再乘以yn+1的(k)减yn+1(k-1)的绝对值

由这样的一个递推的不等式

我们继续往前递推即可以得到

yn+1(k+1)减yn+1(k)

它的绝对值小于等于(hL)的k次方

乘上yn+1(1)减掉yn+1(0)的绝对值

因此 如果(hL)是严格小于1的

当k趋向于无穷大的时候

我们这个迭代序列 就是收敛

所以向后Euler迭代法收敛的条件即为hL小于1

接下来我们再来看

梯形公式

梯形公式为yn+1等于yn

加上2分之h

[f(xn yn)+f(xn+1 yn+1)]

y₀=a

梯形公式也是一个隐式格式

所以我们在求解的时候 也需要迭代求解

它的求解公式初值也是由Euler格式给出

然后

进行迭代计算

yn+1(k+1)等于 yn加2分之h

[f(xn yn)+f(xn+1 yn+1(k))]

对于这个迭代序列

它的收敛性是什么样子的

我们来分析一下

yn+1(k+1)减yn+1(k)它的绝对值

把迭代公式代入

即得到它等于2分之h

f(xn+1 yn+1(k))减f(xn+1 yn+1(k-1)

还是由Lipschitz条件

可以得到 它小于等于2分之hL

乘上yn+1(k)减yn+1(k-1)

我们继续往前递推

即可以得到小于等于2分之hL的k次方

yn+1(1)减掉yn+1(0)

当2分之hL是大于0小于1的时候

我们可以证明这个迭代格式就是收敛的

因此梯形公式迭代法收敛的条件即为

2分之hL严格小于1

接下来 我们来学习改进Euler法

改进Euler法

它是用显格式Euler格式 先来做一个预测

即yn+1拔等于yn加上hf(xn yn)

我们把yn+1拔

去替代梯形格式右端的yn+1

即可以得到yn+1等于yn加上2分之h f(xn yn)

加f(xn+1 yn+1)拔

这样得到的格式

称为Euler公式与梯形公式组成的预测校正系统

又称改进Euler法

我们可以看出

改进Euler法是一个显格式

为了便于编制程序计算

对改进Euler法 我们通常改写成以下形式

即先用Euler格式做一个预测

yp等于yn加上hf(xn yn)

把它代入向后Euler的右端

得到yq等于yn加上hf(xn+1 yp)

然后我们把yp和yq做一个平均

作为我们的数值解

对于改进Euler法我们可以这样来理解

分为预测和校正两步

在预测的时候 我们是用Euler法来预测

这个时候我们选的这种步进的斜率

是(xn yxn)这一点它的切线的斜率

在校正的时候

我们选的斜率是

(xn yxn)和(xn+1 yxn+1)这两个点

它的切线斜率的一个平均

下面我们来看一个例子

用改进Euler法计算初值问题

y'=-y-y²sinx

初始条件是当x=1的时候y=1

它的解在1.2和1.4处的近似值

取步长h是等于0.2

要求计算过程保留小数点后5位

对于这个问题我们的右端项f(x y)=-y-y²sinx

我们把它代入改进Euler格式的计算公式得

yp=yn+h(-yn-yn²sinxn)

yq就等于

yn+h(-yp-yp²sinxn+1)

然后 我们把yp和yq做一下平均记为yn+1

初始值y₀是等于1

由这个公式即可以得到

当n等于0的时候

我们的yp是0.631706

yq是0.799272

因此我们算出来1.2这一点

它的真解的近似值y₁就是

0.71549

当n取1的时候我们算出yp是

0.476964

yq是0.575259

由此我们可以得出解在1.4的值近似的等于y₂

是0.52611

接下来 我们给出单步法的误差估计

我们在前面给出单步法的一般形式为

yn+1等于yn加上hΦ(xn yn xn+1 yn+1)

这里 我们首先给出局部截断误差的定义

局部截断误差它等于真解y(xn+1)减掉yn+1拔

其中这个yn+1拔

等于y(xn)加上h

我们把Φ里面的yn和yn+1都换成真解

即 Φ(xn y(xn) xn+1 y(xn+1))

我们这里得到局部截断误差的技巧 就是Taylor公式

我们再给出整体截断误差的定义

在节点xn+1处

整体截断误差en+1就等于y(xn+1)

减掉yn+1

其中y(xn+1)是这个节点的精确值

而yn+1是这个节点的数值解

我们这两张图上

左边给出来的这张图即为

局部截断误差

而右边我们给出来的

这个真解和近似解的误差就是整体截断误差

如果某种数值方法

它的局部截断误差是

大Oh的p+1次方

我们称这个方法是一个p阶方法

如果我们把它的局部截断误差做一个更精确的展开

把它的hp+1次方这一项精确的表述出来

即写为ψ(xn y(xn))乘上h的p+1次方

加上大Oh的p+2次方

那这时候我们就称ψ(xn y(xn))h的p+1次方

为局部截断误差的主项

下面我们给出Euler方法的局部截断误差估计

我们知道Euler方法的计算公式为

yn+1等于yn加上hf(xn yn)

那所以yn+1拔就等于

y(xn)加上hf(xn y(xn))

而f(xn y(xn))就是y'(xn)

所以yn+1-拔就等于y(xn)加上hy'(xn)

然后我们把y(xn+1)这个精确解

在xn点做Taylor展开

即得到

它等于y(xn)加上hy'(xn)加2分之h的平方y''(ξ)

这个ξ大于xn小于xn+1

我们把两个式子相减即可以得到

局部截断误差Rn+1等于2分之h平方y''(ξ)

ξ是大于xn小于xn+1

Euler方法的局部截断误差为

2分之h平方y''(ξ)

根据定义 我们可以知道

Euler方法是一个一阶方法

下面我们来看一下Euler方法它的整体截断误差

我们这里记 yn+1- 等于y(xn)加上hf(xn y(xn))

en+1就等于真解减掉数值解

我们这里把这个yn+1-给它插进去利用

这个绝对值的三角不等式可以得到

en+1的绝对值

小于等于y(xn+1)减yn+1- 再加上

yn+1- 减掉yn+1的绝对值

这两部分我们再用绝对值给它分开

即可以得到en+1的绝对值小于等于

Rn+1的绝对值加上y(xn)减yn的绝对值

加上hf(xn y(xn))减掉f(xn yn)

我们这里的第二项实际上就是

xn这个节点的整体截断误差

而第三部分 我们可以利用f

它关于y的 条件

把它放大到hL乘上y(xn)减掉yn

那因此啊

在x(n+1)这个节点的整体截断误差

en+1它的绝对值

就小于等于它的局部截断误差Rn+1的绝对值

再加上(1+hL)乘上xn这一点的

整体截断误差的绝对值

那这样的话我们就可以得到

整体截断误差的一个递推关系

如果我们假设右端函数充分光滑

使得这个解是二阶导数连续的

这样话就可以保证它的二阶导数

在这个区间上是有界的

因此局部截断误差Rn+1它的绝对值

就小于等于2分之h的平方乘上大M

那所以整体截断误差en+1它的绝对值

就小于等于2分之h平方大M

加上1+hL的en的绝对值

然后把en它就应该小于等于2分之h平方大M

加上1+hL的en-1的绝对值

然后我们一直递推下去

注意到这个e0是等于零的

那这样话我们即可以得到

en+1的绝对值小于等于2分之h的平方大M

乘上1加上1+hL

一直加到1+hL的n次方

括号里边这个我们可以利用等比数列的求和公式

可以得到en+1的绝对值就小于等于

2L分之hM乘上(1+hL)的n+1次方减1

我们注意到n+1乘上步长是小于等于b-a的

那因此这个(1+hL)的n+1次方

它小于等于(1+hL)的b-a除以h次幂

小于e的大L乘上b-a这样一个常数

那最后我们就可以得到这样一个估计

en+1的绝对值

小于等于2L分之hM乘上e的L(b-a)减1

那因此它是大O(h)的

所以Euler方法它的整体截断误差

是与h同阶的 是一阶方法

由这个误差估计我们可以看出

当步长h趋向于0的时候 在所有的节点

这个整体截断误差它的极限都是0

也即我们的数值解是收敛于真解的

对于向后Euler方法

我们可以完全类似的处理

可以得到它的局部截断误差

即y(xn+1)减掉yn+1-是

负的2分之h²y''(xn)加上大O(h)³

因此它的局部截断误差的主项就为

负的2分之h²y''(xn)

它是一个一阶的方法

对于梯形公式我们也可以做类似的推导

这个时候大家注意对于梯形公式yn+1- 它是

y(xn)加2分之h(f(xn y(xn))加f(xn+1 y(xn+1))

而f(xn) y(xn)即y'(xn)

f(xn+1) y(xn+1)

为y'(xn+1)

我们再把y'(xn+1)

在xn点继续做Taylor展开

就可以得到这个表达式

然后我们把y(xn+1)在xn点也做Taylor展开

我们把y(xn+1)的展开式和y(n+1)-的展开式

两式相减就可以得到

局部截断误差Rn+1等于负的12分之h³y'''(xn)加上大O(h⁴)

因此它的局部截断误差主项为

负的12分之h³y'''(xn)

因此梯形公式是一个二阶的方法

我们总结一下

在本节我们

主要介绍了Euler方法

向后Euler方法

梯形公式 改进Euler方法

以及它们的局部截断误差

和整体截断误差

谢谢大家