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

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 12x相比，可看出欧拉公式的计算结
果精度很差.
上页 下页
欧拉公式具有明显的几何意义, 就是用折线近似代
或表示为下列平均化形式
yp yc
yn yn
hf hf
( x n , y n ), ( x n1 , y p ),
f(x ,y 1 ) f(x ,y 2 ) L y 1 y 2 .
它表明f满足利普希茨条件，且L的大小反映了右端函 数f关于y变化的快慢，刻画了初值问(1.1)式和(1.2)式 是否为好条件. 这在数值求解中也是很重要的.
上页 下页
虽然求解常微分方程有各种各样的解析方法，但 解析方法只能用来求解一些特殊类型的方程，实际问 题中归结出来的微分方程主要靠数值解法.
利普希茨常数(简称Lips.常数).
上页 下页
定理1 设f在区域D={(x, y)|axb,yR}上连续, 关 于y满足利普希茨条件，则对任意x0[a, b], y0R，常 微分方程初值问题(1.1)式和(1.2)式当x[a, b]时存在 唯一的连续可微解y(x) .
解的存在唯一性定理是常微分方程理论的基本内 容，也是数值方法的出发点，此外还要考虑方程的解 对扰动的敏感性，它有以下结论.
可得欧拉法(2.1)的公式误差为
y (x n 1 ) y n 1 h 2 2y (n ) h 2 2y (x n ), (2 .3 )
称为此方法的局部截断误差.
上页 下页
如果对方程(1.1)从xn到xn+1积分，得
y (x n 1 ) y (x n )x x n n 1f(t,y (t)d ).t (2 .4 )
上页 下页
9.2.2 梯形方法 为得到比欧拉法精度高的计算公式，在等式(2.4)
右端积分用梯形求积公式近似, 并用yn代替y(xn), yn+1 代替y(xn+1)，则得
y n 1 y n h 2f(x n ,y n ) f(x n 1 ,y n 1 ), (2 .7 )
称为矩形方法. 矩形方法是隐式单步法，用迭代法求解，同后退
阶常微分方程的初值问题
y y(x 0)f (xy,0y.),x [x0,b],
(1.1) (1.2)
如果存在实数L>0，使得
f ( x ,y 1 ) f ( x ,y 2 ) L y 1 y 2 , y 1 ,y 2 R . ( 1 . 3 )
则称f关于y满足利普希茨(Lipschitz)条件，L称为y的
x y n n 1 1 x y n n y n 1 h y n f(x n ,y n ) y (x n ) f(x n ,y (x n )
即
y n 1 y n h (x n ,f y n )
( 2 .1 )
这就是著名的(显式)欧拉(Euler)公式. 若初值y0已知， 则依公式(2.1)可逐次逐步算出各点数值解.
人口总数成正比，所以t时刻的人口总数y(t)满足以下 微分方程
y y(t),
y(t0
)
y0 .
上页 下页
很多物理系统与时间有关，从卫星运行轨道到单
摆运动，从化学反应到物种竞争都是随时间的延续而
不断变化的. 解常微分方程是描述连续变化的数学语 言，微分方程的求解就是确定满足给定方程的可微函
数y(t)，研究它的数值方法是本章的主要目的. 考虑一
替方程的解曲线，因而常称公式(2.1)为欧拉折线法.
还可以通过几何直观来考察欧拉方法的精度.假设 yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，
按欧拉方法做出的折线
PnPn+1便是y=y(x)过点Pn y
p n1
的切线.从图形上看,这样
y y(x)
定出的顶点Pn+1显著地 偏离了原来的积分曲线，
所谓数值解法, 就是寻求解y(x)在一系列离散节点
x 1 x 2 x n x n 1
上的近似值 y1,y2,,yn,yn+1,. 相邻两个节点的间距 hn=xn+1-xn称为步长. 今后如不特别说明，总是假定 hi=h(i=1,2,)为常数, 这时节点为xn=x0+nh(i=0,1,2,) (等距节点).
的欧拉方法一样，仍用欧拉法提供迭代初值，则矩 形迭代公式为
上页 下页
y yn ((n 0 k 1 ) 11 ) yn yn hh 2(x ffn (,x yn n ,)yn ;)f(xn1,yn (k )1)(k0,1, ).(2.8)
为了分析迭代过程的收敛性, 将(2.7)与(2.8)相减, 得
矩形方法是隐式单步法用迭代法求解同后退的欧拉方法一样仍用欧拉法提供迭代初值则矩形迭代公式为为了分析迭代过程的收敛性将27与28相减hl则当k时有923改进的欧拉公式我们看到梯形方法虽然提高了精度但其算法复杂在应用迭代公式28进行实际计算时每迭代一次都要重新计算函数fx的值而迭代又要反复进行若干次计算量很大而且往往难以预测
显式与隐式两类方法各有特点，考虑到数值稳 定性等其他因素，人们有时需要选用隐式方法，但 使用显式算法远比隐式方便.
隐式方程(2.5)通常用迭代法求解，而迭代过程的 实质是逐步显式化.
上页 下页
设用欧拉公式
yn (0)1ynh(fxn,yn)
给出迭代初值
y
(0) n1
，用它代入(2.5)式的右端，使之转
y1y0h(fx0,y0),
y 2 y1h(fx1,y1),
上页 下页
例1 用欧拉公式求解初值问题
yy2yx (0x1), y(0)1.
(2.2)
解 取步长h=0.1，欧拉公式的具体形式为 yn1 ynh(yn2yxnn)
其中xn=nh=0.1n (n=0,1,,10), 已知y0 =1, 由此式可得
MATLAB第十三讲常微分方程初值问 题数值解法
常微分方程初值问 题数值解法
上页 下页
精品资料
• 你怎么称呼老师？ • 如果老师最后没有总结一节课的重点的难点，你
是否会认为老师的教学方法需要改进？ • 你所经历的课堂，是讲座式还是讨论式？ • 教师的教鞭 • “不怕太阳晒，也不怕那风雨狂，只怕先生骂我
pn
p xn1
可见欧拉方法是相当粗 糙的.
xn
x n1
x
上页 下页
为了分析计算公式的精度，通常可用泰勒展开将
y(yx(nx+n 1 )1在)xny处(x展n开h)，则y(有xn)hy(xn)h22y(n) y(xn)h(fxn,yn)h22y(n) n(xn,xn1).
在yn=y(xn)的前提下，f(xn,yn )=f(xn,y(xn))=y(xn).于是
化为显式，直接计算得
yn (1 )1ynh(x fn1,yn (0 )1), 然后再用 y n(1代)1 入(2.5)式，又有
yn (2 )1ynh(x fn1,yn (1 )1).
如此反复进行，得
y n ( k 1 1 ) y n h ( x n 1 ,f y n ( k ) 1 ) ( k 0 , 1 , )( . 2 . 6 )
函数f有关，当f的Lips.常数L比较小时，解对初值和
右端函数相对不敏感，可视为好条件. 若L较大则可
认为坏条件，即病态问题.
如果右端函数可导，由中值定理有
f(x,y1)f(x,y2)f( x y,)y1y2,在 y1,y2之 . 间
若假定
f (x, y) y
在域D内有界,
设
f ( x, y) y
L, 则
上页 下页
由于f(x, y)对y满足Lipschitz条件(1.3). 由(2.6)减(2.5)得
y n (k 1 1 ) y n 1 h f(x n 1 ,y n (k ) 1 ) f(x n 1 ,y n 1 ) hLyn(k)1yn1.
由此可知，只要hL<1，迭代法(2.6)就收敛到解yn+1. 关于后退欧拉方法的公式误差，从积分公式看到它 与欧拉法是相似的.
称为(隐式)后退的欧拉公式.它也可以通过利用均差 近似导数y (xn+1)，即
直接得到.y(x x n n 11 ) x y( nx n)y(x n 1)f(x n,y(x n 1)),
上页 下页
后退的欧拉公式与欧拉公式有着本质的区别, 后 者是关于yn+1的一个直接计算公式，这类公式称作是 显式的；前者公式的右端含有未知的yn+1，它实际上 是关于yn+1的一个函数方程,这类方程称作是隐式的.
定理2 设f在区域D (如定理1所定义) 上连续, 且 关于y满足利普希茨条件，设初值问题
y (x ) f(x ,y )y ,(x 0 ) s .
的解为y(x, s)，则
y (x ,s 1 ) y (x ,s 2 ) e L x x 0s 1 s 2.
上页 下页
这个定理表明解对初值依赖的敏感性，它与右端
似值
y
n
，称之为预测值，此预测值
1
y的n 1精度可能很差，
再用梯形公式(2.7)将它校正一次，即按(2.8)式迭代一
次得yn+1，这个结果称之为校正值.
上页 下页
这样建立的预测-校正系统通常称为改进的欧拉公式:
预测 校正
yn1ynh(fxn,yn), h
(2.9)
yn1yn2[f(xn,yn)f(xn1,yn1).]
上页 下页
9.2 简单的数值方法
9.2.1 欧拉法与后退欧拉法 我们知道，在xy平面上，微分方程(1.1)式的解
y=f(x)称作它的积分曲线，积分曲线上一点(x, y)的切 线斜率等于函数f(x, y)的值. 如果按f(x, y)在xy平面上 建立一个方向场，那么，积分曲线上每一点的切线 方向均与方向场在该点的方向相一致.
yn 1yn (k 11)h 2[f(xn 1,yn 1)f(xn 1,yn (k )1)],
于是有
yn1yn (k 11) h2L yn1yn (k)1,
使得
hL 1,
2
则当k→∞时有yn(k)1 yn1 , 这说明迭代过程(2.8)是收敛的.
上页 下页
9.2.3 改进的欧拉公式
笨，没有学问无颜见爹娘 ……” • “太阳当空照，花儿对我笑，小鸟说早早早……”
上页 下页
9.1 引 言
科学技术中很多问题都可用微分方程的定解问
题来描述，主要有初值问题与边值问题两大类，本章
只考虑初值问题. 常微分方程初值问题中最简单的例
子是人口模型，设某特定区域在t0时刻人口为y(t0)=y0
已知的，该区域的人口自然增长率为，人口增长与
上页 下页
初值问题的数值解法有个基本特点，他们都采 取“步进式”，即求解过程顺着节点排列的次序一步 一步地向前推进. 描述这类算法，只要给出用已知信 息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)的误差估计 及收敛性，还有递推公式的计算稳定性等问题.
基于上述几何解释，我们从初始点P0(x0, y0)出发, 先依方向场在该点的方向推进到x=x1上一点P1，然后 再从P1点依方向场在该点的方向推进到 x=x2 上一点 P2 , 循环前进做出一条折线P0 P1 P2.
上页 下页
一般地，设已做出该折线的顶点Pn,过Pn(xn, yn)依 方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个顶 点Pn,Pn+1的坐标有关系
我们看到，梯形方法虽然提高了精度，但其算法
复杂，在应用迭代公式(2.8)进行实际计算时，每迭代
一次，都要重新计算函数f(x, y )的值，而迭代又要反
复进行若干次，计算量很大，而且往往难以预测. 为
了控制计算量，通常只迭代一两次就转入下一步的计
算，这就简化了算法.
具体地说，我们先用欧拉公式求得一个初步的近
y1 y0h(y02yx00)10.11.1
y2 y1h(y12yx11)1.10.1(1.110..12)
上页 下页
依次计算下去，部分计算结果见下表.
xn 欧拉公式数值解yn 准确解y(xn)
0.2 1.191818
1.183216
0.4 1.358213
1.341641
0.6 1.508966
右端积分用左矩形公式hf(xn,y(xn))近似，再以yn代替 y(xn)，yn+1代替y(xn+1)也得到欧拉公式(2.1)，局部截 断误差也是(2.3).
如果右端积分用右矩形公式hf(xn+1,y(xn+1))近似， 则得到另一个公式
y n 1 y n h (x n f 1 ,y n 1 ), ( 2 .5 )