第六章 常微分方程的数值解法-2012

1、第六章第六章 常微分方程的数值解法常微分方程的数值解法 对于圆柱形状容器壁上的容积刻度对于圆柱形状容器壁上的容积刻度,可以利用圆柱体可以利用圆柱体体积公式体积公式HDV22其中直径D为常数.由于体积V与相对于容器底部的任意高度H的函数关系明确,因此在容器上可以方便地标出容器刻度,而对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?下表是经过测量得到部分容器高度与直径的关系下表是经过测量得到部分容器高度与直径的关系. .H 0 0.2 0.4 0.6 0.8 1.0D 0 0.11 0.26 0.56 1.04 1.17xDVd41d2其中x表示高度,直径D是高度x的函数,记为D

2、(x),因此得到如下微分方程初值问题0)0()(41dd2VxDxV问题：函数问题：函数D(x)无解析表达式，无法求出其解析解。无解析表达式，无法求出其解析解。xDVd41d2 包含自变量、未知函数及未知函数的导数或微包含自变量、未知函数及未知函数的导数或微分的方程称为分的方程称为微分方程微分方程。 在微分方程中在微分方程中, , 自变量的个数只有一个自变量的个数只有一个, , 称为称为常微分方程常微分方程。 自变量的个数为两个或两个以上的微分方程叫自变量的个数为两个或两个以上的微分方程叫偏微分方程偏微分方程。 微分方程中出现的未知函数最高阶导数的阶数微分方程中出现的未知函数最高阶导数的阶数称

3、为微分方程的称为微分方程的阶阶数。如果未知函数数。如果未知函数y y及其各阶导及其各阶导数数 都是一次的都是一次的, ,则称它是则称它是线性线性的的, ,否则称为否则称为非线性非线性的。的。 )(,nyyy 常微分方程常微分方程(未知函数是一元函数未知函数是一元函数) 偏微分方程偏微分方程( (未知函数是多元函数未知函数是多元函数) ) vmcgdtdv 22xuxuutu 0yx222 ,2)( 的的解解xxy cxxy 2)(.,才能确定方程的解件方程加上适当的定解条边值条件初始值条件定解条件)2()1(:同一个微分方程，同一个微分方程，具有不同的初始条件具有不同的初始条件一阶一阶常微分方

4、程的常微分方程的初值问题初值问题0)(,),(ddyaybaxyxfxy解的存在唯一性解的存在唯一性（“常微分方程常微分方程”理论）：只要理论）：只要 f (x, y) 在在a, b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与，即存在与 x, y 无关的常数无关的常数 L 使使对任意定义在对任意定义在 a, b 上的上的 y1(x) 和和 y2(x) 都成立，则上述方程都成立，则上述方程存存在唯一解在唯一解。1212| ( , )( ,)|f x yf x yL yy0022220222( ,)()11221(0)11cos1cossin(0)11

5、1sinxyfx yyxyxydyxdxycxcyyyxcyyxyxdyxdxxcyyyyx 本 章 重 点 讨 论 一 阶 常 微 分 方 程 初 值 问 题dydxdydx当当x=0时时,y=1,可得可得c=1特解特解当当x=0时时,y=1,可得可得c=-1特解特解两边积分两边积分通解通解常微分方程求解的基本方法：常微分方程求解的基本方法：可分离变量法可分离变量法常系数齐次线性方程的解法常系数齐次线性方程的解法常系数非齐次线性方程的解法常系数非齐次线性方程的解法1)0( 1 , 0,dd2yxexyx 欧拉（Euler）方法是解初值问题的最简单的数值方法。初值问题的解y=y(x)代表通过点

6、 的一条称之为微分方程的积分曲线。积分曲线上每一点 的切线的斜率 等于函数 在这点的值。 ),(00yx),(yx)(xy),(yxf00)(),(yxyyxfy1 欧拉欧拉(Euler)法法Euler法的求解过程：法的求解过程：)(,(0000 xxyxfyy直线方程100010(,)()yyf xyxx交点P1 Pi+1 Pn y=y(x) P1 Pi P0 x0 x1 xi xi+1 xn P1 从初始点P0出发,作斜率为 的直线，与x=x1直线相交于P1点0( )y x00(,)f xy)(,(1111xxyxfyy)(,(121112xxyxfyy交点P2 Pi+1 Pn y=y(x

7、) P1 Pi P0 x0 x1 xi xi+1 xn P1P2 过点P1,作斜率为 的直线，与x=x2直线相交于P2点1( )y x11(,)f xy Pi+1 Pn y= y(x) P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 重复以上过程，就可获得一系列的点： P1,P1,Pn。对已求得点 ，以 为斜率作直线 (,)nnnP xy()ny x(,)nnf xy( ,)()nnnnyyf x yx x当 时,得1nxx11(,)()nnnnnnyyf xyxx 从图形上看,就获得了一条近似于曲线y=y(x) 的折线 。 Pi+ 1 Pn y= y(x) P

8、1 Pi Pn Pi+ 1 P0 x0 x1 xi xi+ 1 xn Pi P1 这样,从x0逐个算出对应的数值解 nxxx,21nyyy,21nPPPP3211nnxxh100(,)()nnnnyyhf xyyy x n=0,1,11(,)()nnnnnnyyf xyxx欧拉欧拉(Euler)法法00y)y(x );y,x(fydxdy x0 x1x2x30欧拉折线法欧拉折线法:, 0 x1 y(0)1, dyxydx例 用欧拉方法求下列方程的数值解 25. 1)15 . 0)(5 . 0(1)0 . 1(0 . 1)1(0.5)(01(0.5)11120001yhxyyyyhxyyy解：E

9、uler公式为 ),(1nnnnyxhfyy 当h=0.5时, 1 , 0 n nnnyhxy 25. 0 5 . 0 hh分分别别取取1)1(0.25)(0125. 0(0.25)0001 yxyyy), 1 , 0( 1 nyhxyynnnn当h=0.25时0625. 1)125. 0)(25. 0(1 25. 0)5 . 0(1112 yxyyy191347. 1)0625. 15 . 0)(25. 0(0625. 1 25. 0)75. 0(2223 yxyyy39600. 1)191347. 175. 0)(25. 0(191347. 1 )191347. 1 ,75. 0()75.

10、 0(25. 0)0 . 1(3334 hfyyxyyy; (0)1dyyxyydx00.50.751.010.25h = 0.5h = 0.25224) /x1(: y解析解为解析解为欧拉方法的收敛性欧拉方法的收敛性11212 () Taylor, (,),()()()()()2 ()(, ()() 2nnnnnnnnnnnnnny xxxxhy xy xhy xhy xyhy xhf xy xy将在点展开1(), () (, () nnnnny xyy xh f xy x假定已知准确值 利用欧拉公式，定义11112() ()()(, () (). 2nnnnnnnnRy xyy xy xh

11、f xy xhy2221 max |( )|,|()|(). 22a x bnnMy xhhRyMO h 令则局部截断误差局部截断误差定义定义2: 若给定方法的局部截断误差满足若给定方法的局部截断误差满足则称该方法是则称该方法是 P 阶的阶的，或称为，或称为具有具有 P 阶精度阶精度。11|(),pnRO hEuler公式是一阶收敛的！公式是一阶收敛的！2312()()hnnRy xO h21(). 2nnhRy x还可用以下方法推导还可用以下方法推导EulerEuler格式：格式： 数值微分数值微分 数值积分法数值积分法 对微分方程的离散，可以有多种思路，但最基本的想法是“以直代曲”。若用向

12、前差商近似导数若用向前差商近似导数)(,()()()()(1nnnnnnxyxhfxyxyhxyxy )(),(01ayyyxhfyynnnnEuler公式公式hxyxyxynnn)()()(1 )(,()()(1nnnnxyxfhxyxy )(,()(nnnxyxfxy （1）用数值微分方法）用数值微分方法若用向后差商近似导数，即若用向后差商近似导数，即)(,()()(111 nnnnxyxhfxyxy )(),(0111ayyyxhfyynnnnhxyxyxynnn)()()(11 后退后退Euler公式公式)(,()()(111 nnnnxyxfhxyxy隐式隐式Euler公式公式 隐式

13、隐式欧拉法的局部截断误差：欧拉法的局部截断误差：1111111112111, ,2nnnnynnnnnnnnnnnfxyfxy xfxyy xyy xhfxy xyxyxhyxyx由微分中值定理，得 其中 在，之间；又111()(,) nnnnyy xhf xy由向后欧拉公式， 定义 111132, 2nnynnnnnnyy xhfxyy xhhyxh yxyx23126nnnnnhhy xy xhyxyxyx而111111232311(), 231,23nnnynnnnnynnnnRy xyhfxy xyhhyxyxhhhfxRyxyx 从而即2211121,1,(1)1ynynynhfxh

14、fxhfxxxx 11 2322111231221 R1,23 3,226 R()2nynynnnnynnnnnhhhfxhfxy xyxhhy xfxy xy xhy xo h隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：111()nnnRy xy 232()()hny xO h即隐式欧拉公式也是一阶收敛的。即隐式欧拉公式也是一阶收敛的。1(,) 0, 1,.nnnnyyh f xyn比较欧拉显式公式和隐式公式及其局部截断误差231112()()()hnnnnTy xyy xO h显式公式显式公式111(,)nnnnyyh f xy隐式公式隐式公式231112()()()hnnnnTy

15、 xyy xO h梯形公式梯形公式 显、隐式两种算法的显、隐式两种算法的平均平均111 (,)(,)2nnnnnnhyyf xyf xy注：注：的确有局部截断误差的确有局部截断误差 ， 即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。但阶精度，比欧拉方法有了进步。但注意到该公式是注意到该公式是隐式隐式公式，计算时不得不用到迭代公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。法，其迭代收敛性与欧拉公式相似。3111()()nnnRy xyO h（2）用数值积分方法）用数值积分方法 1)(,()()(1nnxxnndxxyxfxyxy 11)(,()(nnnnxxxxdxxy

16、xfdxxy1111( , )(, (), ( , )(, (),()()nnnnnnnnf x yf xy xf x yf xy xy xyy xy分别用左矩形和右矩形公式，即 代替上式右端的积分,并注意 ,分别得到1111 (,) (,) nnnnnnnnyyh f xyyyh f xy，。向前欧拉公式和向后欧拉公式：)(,()(nnnxyxfxy 若对积分用梯形公式，则得若对积分用梯形公式，则得 )(,()(,(2)()(111 nnnnnnxyxfxyxfhxyxy )(),(),(20111ayyyxfyxfhyynnnnnn梯形公式梯形公式公式公式局部截断误差局部截断误差精精度度显

17、显隐隐稳稳定定性性步数步数欧拉欧拉显式显式公式公式1 1阶阶显显差差单步单步欧拉隐欧拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形梯形公式公式2 2阶阶隐隐好好单步单步3(3)3nhyx2(2)2nhyx2(2)2nhyx欧拉法小结欧拉法小结8.3 改进欧拉(Euler)方法8.3 改进欧拉(Euler)方法2 改进欧拉改进欧拉(Euler)方法方法显式欧拉公式：显式欧拉公式：计算工作量小，但精度低。梯形公式：梯形公式：提高了精度，但为隐式公式，需用迭代法求解,计算工作量大。改进的欧拉公式：改进的欧拉公式：综合欧拉公式和梯形公式 先用欧拉公式求出一个初步的近似值 ，称为预测值，它的精度不高，再

18、用梯形公式对它校正一次，即迭代一次，求得yn+1，称为校正值， 这种预测-校正方法称为改进的欧拉公式：改进的欧拉公式：1ny1111( ,) ( ,)(,) 2nnnnnnnnnnyyhf x yhyyf x yf xy预测校正称为称为Euler公式与梯形公式的公式与梯形公式的预测预测校正系统校正系统。 改进的欧拉公式的精度为改进的欧拉公式的精度为二阶二阶。这是一种一步。这是一种一步显式格式显式格式，它可以表示为嵌套形式，它可以表示为嵌套形式。11(,)(,(,)2nnnnnnnnhyyf xyf xyhf xy1111( ,) ( ,)(,) 2nnnnnnnnnnyyhf x yhyyf

19、x yf xy )(21),(),(11qpnpnnqnnnpyyyyxhfyyyxhfyy实际计算时，常改写成以下形式1111( ,) ( ,)(,) 2nnnnnnnnnnyyhf x yhyyf x yf xy 校校正正预预测测 ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyy:Euler解 改进的法计算公式 nnnnnnnnnnnyhhyyhyyyhyyhyyy2)(2 211例：写出如下初值问题的改近Euler 法的近似解的表达式(0)1yyy )(21),(),(11qpnpnnqnnnpyyyyxhfyyyxhfyy.5, 2 . 0,4 . 12

20、 . 1 , 1)1(,sin:2位位小小数数点点后后至至少少保保留留取取步步长长处处的的近近似似值值和和的的解解在在方方法法计计算算初初值值问问题题用用改改进进的的例例 hyxyyyEuler方方法法改改进进的的解解Euler: 校校正正预预测测 ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyyxyyyxfsin),(2 2211(sin)(sin)1()2pnnnnqnppnnpqyyhyyxyyhyyxyyy0(1)1,0.2yyh2211(sin)(sin)1()2pnnnnqnppnnpqyyhyyxyyhyyxyyy00,1 (1)1nyy10,0.

21、631706 0.799272 (1.2)0.71549pqnyyyy时21,0.476964 0.575259(1.4)0.52611pqnyyyy时例：考察初值问题例：考察初值问题 在区间在区间0, 0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。 1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.0000

22、2.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 72 2 单步法的稳定性单步法的稳定性定义：定义：若某算法在计算过程中任一步产生的误差在以后的计算中都逐步衰减，则称该算法是绝对稳定的。一般分析某算法的稳定性时，为简单起见，只考虑模型方程或试验方程：yy , . yy 为一复常数引进试验方程：*111,nnneyy1*1

23、nnyy为由公式得到的准确值(无舍入误差)，为由公式实际得到的值(有舍入误差).记数值误差为：向前欧拉公式的稳定性向前欧拉公式的稳定性试验方程的欧拉公式若每步计算有舍入误差，则 1 (,), 0,1,2,nnnnyyh f xyn1 (,) (1)nnnnnnyyh f xyyh y*1 (2)nnnyyh y11(1)(2): (1), |1|1 nnnnneeh eehhe当 稳定。当为复数时当为实数时 |1| 1 Re()11hh 表示 复平面上以为中心，以 为半径的圆,此圆称为向前欧拉方法的绝对稳定区域。|1| 1 20.hh 绝对稳定区域为带状。3 3 龙格龙格- -库塔法库塔法),

24、(),(,()()()()(11nnnnnnnnnxxyhfxyhyxyxy微分中值定理.,)()(,(1上的平均斜率在区间为曲线nnnnxxxyyf0)(),(yaybxayxfy初值问题)(,()(,(nnnnxyxfyfEuler公式)(,()(,(11nnnnxyxfyf后退Euler公式2/)( ),( ),(2111121kkhyyhkyxfkyxfknnnnnn改进Euler公式步长一定是步长一定是h 吗？吗？斜率斜率一定取一定取k1 k2 的的平均值平均值吗？吗？.,1公公式式可可能能构构造造出出精精度度更更高高的的代代替替平平均均斜斜率率用用他他们们的的加加权权平平均均值值上

25、上多多预预测测几几个个点点的的斜斜率率在在区区间间 nnxx例如取例如取 m 个点的斜率构造如下形式的公式个点的斜率构造如下形式的公式 )(),(),(),(),(2211111,1123213133121221mmnnmmmmnmnmnnnnnnkckckchyyhkbhkbyhaxfkhkbhkbyhaxfkhkbyhaxfkyxfk,的的常常数数无无关关是是与与其其中中nfcbaiiji该公式称为该公式称为龙格库塔龙格库塔(Runge-Kutta)公式公式,简称简称R-K公式公式.求解：求解：只需将公式的局部截断误差在只需将公式的局部截断误差在xn点进行点进行Taylor展开展开,令其令

26、其前面尽可能多的项为前面尽可能多的项为0, 便可导出便可导出ai，bij，ci所满足的方程组所满足的方程组,即可即可从中求出这些系数从中求出这些系数. 以二阶R-K公式为例，说明R-K公式的构造过程)(),(),(22111121221kckchyyhkbyhaxfkyxfknnnnnn.,21212待定其中ccba其局部截断误差为：)()()(2211111kckchyxyyxyRnnnn点进行泰勒展开在nx),(),()(1212211hkbyhaxfhcyxfhcyxynnnnnn)(! 32)(4321hOyhyhyhyxynnnnn 点进行泰勒展开在nnxxy)(1fyxfynnn)

27、,(其中f f fxyyfxfyxfdxd ydxd yyxnnnn),()(2)(2yxyyyxyxxyxnnff f f ff ff ff f fdxd ydxd y泰勒展开在nnnxhkbyhaxf),(1212)(),(22121212hOhffbhfafhkbyhaxfyxnn)()21()21()1 (3221222221hOhffbch facfhccRyx代入局部截断误差公式要使R关于h的阶数尽量高，应满足021021012122221bcaccc021021012122221bcaccc无穷组解若取1,2121221bacc)(21),(),(211121kkhyyhkyhx

28、fkyxfknnnnnn)(),(),(22111121221kckchyyhkbyhaxfkyxfknnnnnn改进改进Euler公式公式021021012122221bcaccc无穷组解若取21, 1, 021221bacc21121)21,21(),(hkyyhkyhxfkyxfknnnnnn)(),(),(22111121221kckchyyhkbyhaxfkyxfknnnnnn中点公式中点公式)22( ),( ),( ),( ),(43216134222312221kkkkyyhkyhxfkkyxfkkyxfkyxfkhnnnnhnhnhnhnnn 常用的标准四阶常用的标准四阶RK公

29、式公式（经典（经典R-K方法）方法）最常用的四阶标准最常用的四阶标准RK公式（经典公式（经典R-K方法）为：方法）为：例例 用四阶标准用四阶标准R-K公式解初值问题公式解初值问题 1)0(10 ,2yxyxyy解：取解：取 h=0.2 ，四阶标准，四阶标准R-K法的具体格式如下法的具体格式如下: )22(6 )(2)(),(2)2(2)2()2,2(2)2(2)2()2,2( /2),(432113334222311121kkkkhyyhkyhxhkyhkyhxfkkhyhxkhykhyhxfkkhyhxkhykhyhxfkyxyyxfknnnnnnnnnnnnnnnnnnnnnn10y5 ,

30、 1 , 0 , 2 . 0nnxn已知已知 183229. 1)22(6 84324. 033849. 0181728. 1)(2)(90864. 018318. 0091818. 12)2(2)2( 91818. 01 . 12 . 01 . 12)2(2- )2( 1012-4321013003042002031001020001kkkkhyyhkyhxhkykkhyhxkhykkhyhxkhykyxyk同理可计算得同理可计算得， , ,5432yyyy具体结果见下表具体结果见下表 至少具有四位有效数字至少具有四位有效数字. 4 一阶方程组和高阶方程的数值解法一阶方程组和高阶方程的数值解

31、法bxaayayyyxfdxdyyyxfdxdymmmmmm , )()(),(),(111111写成向量的形式：)(),(aYYxFdxdYm个未知函数，m个初始条件各种方法都可以直接运用过来。bxazazyayzyxgdxdzzyxfdxdy , )()(),(),(00Euler公式),(),(11nnnnnnnnnnzyxhgzzzyxhfyy以两个方程的方程组为例R-K公式112341(22)6nnnnyyhKKKKzz)2,2(12KhYhxFK1(1)(2)112(1)(2)11(,)(,)(,)222(,)222nnnnnnnnnnnnf xyzKg xyzhhhf xyKzKKhhhg xyKzK bxazazyayzyxgdxdzzyxfdxdy , )()(),(),(00)(),(aYYxFdxdY(1)(1)223(1)(1)22(1)(1)334(1)(1)33(,)222(,)222(,)(,)nnnnnnnnnnnnhhhf xyKzKKhhhg xyKzKf xh yhKzhKKg xh yhKzhK)2,2(23KhYhxFK),(34