MATLAB第十三讲常微分方程初值问题数值解法ppt课件

xn1 xn
f (t, y(t))dt .
(2.4)
右端积分用左矩形公式hf(xn,y(xn))近似，再以yn代替 y(xn)，yn+1代替y(xn+1)也得到欧拉公式(2.1)，局部截
断误差也是(2.3).
如果右端积分用右矩形公式hf(xn+1,y(xn+1))近似， 则得到另一个公式
yn1 yn hf ( xn1, yn1 ),
基于上述几何解释，我们从初始点P0(x0, y0)出发, 先依方向场在该点的方向推进到x=x1上一点P1，然后 再从P1点依方向场在该点的方向推进到 x=x2 上一点 P2 , 循环前进做出一条折线P0 P1 P2.
上页 下页
一般地，设已做出该折线的顶点Pn,过Pn(xn, yn)依
方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个顶
y1 y0 hf ( x0, y0 ),
y2 y1 hf (x1, y1),
上页 下页
例1 用欧拉公式求解初值问题
y
y
2x y
(0 x 1),
y(0) 1.
(2.2)
解 取步长h=0.1，欧拉公式的具体形式为
yn1
yn
h( yn
2xn ) yn
其中xn=nh=0.1n (n=0,1,,10), 已知y0 =1, 由此式可得
)],
于是有
yn1
y(k1) n1
hL 2
yn1
y(k) n1
,
使得
hL 1,
2
则当k→∞时有yn(k)1 yn1 , 这说明迭代过程(2.8)是收敛的.
上页 下页
9.2.3 改进的欧拉公式
我们看到，梯形方法虽然提高了精度，但其算法 复杂，在应用迭代公式(2.8)进行实际计算时，每迭代 一次，都要重新计算函数f(x, y )的值，而迭代又要反 复进行若干次，计算量很大，而且往往难以预测. 为 了控制计算量，通常只迭代一两次就转入下一步的计 算，这就简化了算法.
具体地说，我们先用欧拉公式求得一个初步的近
似值 yn，1 称之为预测值，此预测值 的yn精1 度可能很差，
再用梯形公式(2.7)将它校正一次，即按(2.8)式迭代一 次得yn+1，这个结果称之为校正值.
上页 下页
这样建立的预测-校正系统通常称为改进的欧拉公式:
预测 yn1 yn hf (xn, yn ) , h
0.8 1.616476
1.0 1.737869
误差
0.002716 0.004024 0.05818
分方程
y y(t),
y(
t
0
)
y0 .
上页 下页
很多物理系统与时间有关，从卫星运行轨道到单
摆运动，从化学反应到物种竞争都是随时间的延续而
不断变化的. 解常微分方程是描述连续变化的数学语
言，微分方程的求解就是确定满足给定方程的可微函
数y(t)，研究它的数值方法是本章的主要目的. 考虑一
阶常微分方程的初值问题
右端积分用梯形求积公式近似, 并用yn代替y(xn), yn+1
代替y(xn+1)，则得
h
yn1 yn 2 f (xn, yn ) f (xn1, yn1) ,
(2.7)
称为矩形方法.
矩形方法是隐式单步法，用迭代法求解，同后退 的欧拉方法一样，仍用欧拉法提供迭代初值，则矩 形迭代公式为
上页 下页
所谓数值解法, 就是寻求解y(x)在一系列离散节点
x1 x2 xn xn1
上的近似值 y1,y2,,yn,yn+1,. 相邻两个节点的间距 hn=xn+1-xn称为步长. 今后如不特别说明，总是假定 hi=h(i=1,2,)为常数, 这时节点为xn=x0+nh(i=0,1,2,) (等距节点).
上页 下页
9.2 简单的数值方法
9.2.1 欧拉法与后退欧拉法 我们知道，在xy平面上，微分方程(1.1)式的解
y=f(x)称作它的积分曲线，积分曲线上一点(x, y)的切 线斜率等于函数f(x, y)的值. 如果按f(x, y)在xy平面上 建立一个方向场，那么，积分曲线上每一点的切线 方向均与方向场在该点的方向相一致.
替方程的解曲线，因而常称公式(2.1)为欧拉折线法.
还可以通过几何直观来考察欧拉方法的精度.假设 yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，
按欧拉方法做出的折线
PnPn+1便是y=y(x)过点Pn y
pn1
的切线.从图形上看,这样
y y(x)
定出的顶点Pn+1显著地 偏离了原来的积分曲线，
常微分方程初值问 题数值解法
上页 下页
9.1 引 言
科学技术中很多问题都可用微分方程的定解问题
来描述，主要有初值问题与边值问题两大类，本章只
考虑初值问题. 常微分方程初值问题中最简单的例子
是人口模型，设某特定区域在t0时刻人口为y(t0)=y0已
知的，该区域的人口自然增长率为，人口增长与人
口总数成正比，所以t时刻的人口总数y(t)满足以下微
(2.5)
称为(隐式)后退的欧拉公式.它也可以通过利用均差
近似导数y (xn+1)，即
y( xn1) y( xn ) 直接得到. xn1 xn
y( xn1 )
f ( xn , y( xn1 )),
上页 下页
后退的欧拉公式与欧拉公式有着本质的区别, 后 者是关于yn+1的一个直接计算公式，这类公式称作是 显式的；前者公式的右端含有未知的yn+1，它实际上 是关于yn+1的一个函数方程,这类方程称作是隐式的.
y1
y0
h( y0
2 x0 y0
)
1
0.1
1.1
y2
y1
h( y1
2 x1 y1
)
1.1
0.1(1.1
0.2 )
1.1
1.191818
上页 下页
依次计算下去，部分计算结果见下表.
xn 欧拉公式数值解yn 准确解y(xn)
0.2 1.191818
1.183216
0.4 1.358213
1.341641
点Pn,Pn+1的坐标有关系
yn1 yn xn1 xn
yn1 h
yn
f ( xn, yn )
y( xn )
f ( hf ( xn, yn )
(2 .1)
这就是著名的(显式)欧拉(Euler)公式. 若初值y0已知， 则依公式(2.1)可逐次逐步算出各点数值解.
y(0) n1
yn
hf ( xn ,
yn );
y
( k 1)
n1
yn
h 2
f
( xn ,
yn)
f (xn1,
y(k) n1
)
(2.8) (k 0,1,).
为了分析迭代过程的收敛性, 将(2.7)与(2.8)相减, 得
yn1
y(k1) n1
h[ 2
f
( xn1,
yn1)
f
( xn1,
y(k) n1
在yn=y(xn)的前提下，f(xn,yn )=f(xn,y(xn))=y(xn).于是
可得欧拉法(2.1)的公式误差为
y(xn1)
yn1
h2 2
y(n )
h2 2
y( xn ),
(2.3)
称为此方法的局部截断误差.
上页 下页
如果对方程(1.1)从xn到xn+1积分，得
y( xn1 ) y( xn )
pn
p xn1
可见欧拉方法是相当粗 糙的.
xn
xn 1
x
上页 下页
为了分析计算公式的精度，通常可用泰勒展开将
y(xn+1)在xn处展开，则有
y( xn1) y( xn h) y( xn )
hy( xn )
h2 2
y(n )
y(xn)
hf ( xn ,
yn)
h2 2
y(n )
n ( xn , xn1 ).
y f ( x, y),
y(
x0
)
y0 .
x [ x0 , b],
(1.1) (1.2)
如果存在实数L>0，使得
f ( x, y1 ) f ( x, y2 ) L y1 y2 , y1 , y2 R. (1.3)
则称f关于y满足利普希茨(Lipschitz)条件，L称为y的
利普希茨常数(简称Lips.常数).
函数f有关，当f的Lips.常数L比较小时，解对初值和
右端函数相对不敏感，可视为好条件. 若L较大则可
认为坏条件，即病态问题.
如果右端函数可导，由中值定理有
f ( x, y1 )
f (x, y2 )
f ( x, )
y
y1
y2 , 在y1 , y2之间.
若假定
f
( x, y
y)
在域D内有界,
设
f ( x, y) y
化为显式，直接计算得
y(1) n1
yn
hf
( xn1 ,
y(0) n1
),
然后再用 yn(代1)1 入(2.5)式，又有
y(2) n1
yn
hf
( xn1 ,
y(1) n1
).
如此反复进行，得
y(k1) n1
yn
hf (xn1,
y(k) n1
)
(k
0,1,).
(2.6)
上页 下页
由于f(x, y)对y满足Lipschitz条件(1.3). 由(2.6)减(2.5)得
(2.9)
校正 yn1 yn 2 [ f ( xn , yn ) f ( xn1, yn1 )] .
或表示为下列平均化形式
yp yc
yn yn
hf hf
( xn , ( xn1
yn ), , yp ),
yn1
1 2
(
yp
yc
).
上页 下页
例2 用改进的欧拉法解 例1中的初值问题(2.2).
定理2 设f在区域D (如定理1所定义) 上连续, 且 关于y满足利普希茨条件，设初值问题
y( x) f ( x, y), y( x0 ) s.
的解为y(x, s)，则
y( x, s1 ) y( x, s2 ) e L x x0 s1 s2 .
上页 下页
这个定理表明解对初值依赖的敏感性，它与右端
解 仍取步长h=0.1,改 进的欧拉公式为
部分计算结果见下表
y
p
yn
h( yn
2 xn yn
),
yc
yn
h( y p
2 xn1 ), yp
yn1
1 2
(
yp
yc ).
xn
yn
0 1.
0.2 1.184096
0.4 1.343360
误差
0 0.000088 0.001719
xn
yn
0.6 1.485956
上页 下页
定理1 设f在区域D={(x, y)|axb,yR}上连续, 关 于y满足利普希茨条件，则对任意x0[a, b], y0R，常 微分方程初值问题(1.1)式和(1.2)式当x[a, b]时存在唯 一的连续可微解y(x) .
解的存在唯一性定理是常微分方程理论的基本内 容，也是数值方法的出发点，此外还要考虑方程的解 对扰动的敏感性，它有以下结论.
y(k1) n1
yn1
h
f
(
xn1
,
y(k) n1
)
f (xn1, yn1)
hL
y(k) n1
yn1 .
由此可知，只要hL<1，迭代法(2.6)就收敛到解yn+1. 关于后退欧拉方法的公式误差，从积分公式看到它 与欧拉法是相似的.
上页 下页
9.2.2 梯形方法
为得到比欧拉法精度高的计算公式，在等式(2.4)
0.6 1.508966
1.483240
0.8 1.649783
1.612452
1.0 1.784770
1.732051
误差 0.008602 0.016572 0.025726 0.037331 0.052719
与准确解 y 1 2x 相比，可看出欧拉公式的计算结
果精度很差.
上页 下页
欧拉公式具有明显的几何意义, 就是用折线近似代
L, 则
f ( x, y1 ) f ( x, y2 ) L y1 y2 .
它表明f满足利普希茨条件，且L的大小反映了右端函 数f关于y变化的快慢，刻画了初值问(1.1)式和(1.2)式 是否为好条件. 这在数值求解中也是很重要的.
上页 下页
虽然求解常微分方程有各种各样的解析方法，但 解析方法只能用来求解一些特殊类型的方程，实际问 题中归结出来的微分方程主要靠数值解法.
显式与隐式两类方法各有特点，考虑到数值稳 定性等其他因素，人们有时需要选用隐式方法，但 使用显式算法远比隐式方便.
隐式方程(2.5)通常用迭代法求解，而迭代过程的 实质是逐步显式化.
上页 下页
设用欧拉公式
y(0) n1
yn
hf
( xn ,
yn )
给出迭代初值 yn(0)1，用它代入(2.5)式的右端，使之转
上页 下页
初值问题的数值解法有个基本特点，他们都采取 “步进式”，即求解过程顺着节点排列的次序一步一 步地向前推进. 描述这类算法，只要给出用已知信息 yn,yn-1,yn-2,计算yn+1的递推公式.
本章首先要对常微分方程(1.1)离散化，建立求解 数值解的递推公式. 一类是计算yn+1时只用到前一点的 值yn，称为单步法. 另一类是用到yn+1前面 k 点的值 yn,yn-1,, yn-k+1，称为k步法. 其次，要研究公式的局部 截断误差和阶，数值解yn与精确解y(xn)的误差估计及 收敛性，还有递推公式的计算稳定性等问题.