《计算方法第八章》ppt课件

1、第八章第八章 常微分方程初值常微分方程初值问题数值解法问题数值解法本章主要研讨常微分方程初值问题的数值求解：本章主要研讨常微分方程初值问题的数值求解：通常，假设函数通常，假设函数 f 关于第二个变量满足李普希茨条件关于第二个变量满足李普希茨条件L条条件，即为存在常数件，即为存在常数 L 0，使得，使得0( )( , ( )( )y xf x y xaxby ay1212( ,)( ,)f x yf x yL yy第一节第一节 普通概念普通概念1.1 欧拉法及其简单改良欧拉法及其简单改良0( )( , ( )( )y xf x y xaxby ay方法：选择适当的节点，用差分近似微分，将方程离散

2、化，方法：选择适当的节点，用差分近似微分，将方程离散化，从而求在这些节点上的函数值。从而求在这些节点上的函数值。012111NNnnnnnaxxxxxbhxxxx称为到的步长（通常取为常数h）1()()(, ()(,)nnnnnny xy xf xy xf xyh欧拉欧拉方法方法100(,),()nnnnyyhf xyy xy例子：例子：(01)(0)1yyxy 10( , )(1)1nnnnf x yyyyhyh yy xye精确解为下面我们分别取步长为下面我们分别取步长为0.1与与0.01进展计算，进展计算，计算结果显示在下面的图中。计算结果显示在下面的图中。步长为步长为0.1的计算结果。

3、的计算结果。步长为步长为0.01的计算结果的计算结果0.01 0.99005 0.99 0.1 0.90484 0.90438 0.2 0.81873 0.81791 0.3 0.74082 0.7397 0.4 0.67032 0.66897 0.41 0.66365 0.66228 0.59 0.55433 0.55268 0.6 0.54881 0.54716 0.9 0.40657 0.40473 0.91 0.40252 0.40068 0.99 0.37158 0.36973 1 0.36788 0.36603 DOUBLE PRECISION h,y(0：100) OPEN(20

4、,FILE=OUTPUT.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0do 10 i=1,100 y(i)=y(i-1)*(1.0-h)write(20,*) i*h,y(i)10 continueEND精度更高的欧拉公式：精度更高的欧拉公式：方法：选择计算导数值方法：选择计算导数值的精度更好的差分格式。的精度更好的差分格式。11()2()()2nnnnnyyy xhy xhy xhh3222()( )/3!( )()()()212nnny x hyhyy xhy xO hh11002(,),()nnnnyyhf xyy xy欧拉中点公式欧拉中点公式利用中点公式求解

5、微分方程时，有一个问题，就是计算时需利用中点公式求解微分方程时，有一个问题，就是计算时需求两个迭代初值！求两个迭代初值！对于这个问题，我们可以先用欧拉公式，经过给定的初值计对于这个问题，我们可以先用欧拉公式，经过给定的初值计算出第一个点的值，然后在利用这两个第一和第二个点算出第一个点的值，然后在利用这两个第一和第二个点的值进展计算，直到计算出全部节点上的值。的值进展计算，直到计算出全部节点上的值。下面，我们用中点公式求解刚刚的例子，计算的步长取下面，我们用中点公式求解刚刚的例子，计算的步长取0.01，可以看到，计算的精度比较高。可以看到，计算的精度比较高。此时，计算公式为此时，计算公式为11(

6、 , )2nnnf x yyyyhy 100(1),1yh yy计计算算结结果果0.02 0.9802 0.98020.1 0.90484 0.904840.2 0.81873 0.818740.3 0.74082 0.740840.4 0.67032 0.670350.5 0.60653 0.606560.6 0.54881 0.548850.7 0.49659 0.496630.8 0.44933 0.449380.9 0.40657 0.406631.2 欧拉方法的其他改良欧拉方法的其他改良微分方程数值解的关键在于对导数的处置，可以用差分来近似微分方程数值解的关键在于对导数的处置，可以用

7、差分来近似导数，也可以经过积分，将导数项化掉。导数，也可以经过积分，将导数项化掉。对于方程：对于方程：首先，作出划分首先，作出划分设曾经求出第设曾经求出第 n 个节点的函数值个节点的函数值 ，在区间，在区间 上对上对方程两边积分方程两边积分容易看出，要求第容易看出，要求第 n+1 个节点的函数值，关键在于选择适当个节点的函数值，关键在于选择适当的积分公式计算积分！的积分公式计算积分！( )( , ( )y xf x y x0121NNaxxxxxb1,nnxxny11()()( , ( )nnxnnxy xy xf x y x dx1如选择下矩形公式，那么如选择下矩形公式，那么这正是前面的欧拉

8、公式。这正是前面的欧拉公式。1(,)nnnnyyf xy h111(,)nnnnyyf xyh111 (,)(,)/2nnnnnnyyf xyf xyh2如选择上矩形公式，那么如选择上矩形公式，那么这是所谓的后退欧拉公式。这是所谓的后退欧拉公式。3如选择梯形公式，那么如选择梯形公式，那么这是所谓的欧拉梯形公式。这是所谓的欧拉梯形公式。直接利用曾经求得的知节点上的值计算未知节点上的函数值直接利用曾经求得的知节点上的值计算未知节点上的函数值的算法称为显式法。的算法称为显式法。例如：欧拉公式、欧拉中点公式例如：欧拉公式、欧拉中点公式计算未知节点上的函数值时，用到了未知节点上的函数值，计算未知节点上的

9、函数值时，用到了未知节点上的函数值，这种算法称为隐式法。这种算法称为隐式法。例如：后退欧拉法、欧拉梯形公式例如：后退欧拉法、欧拉梯形公式显然，利用隐式法求微分方程的数值解是，需求从表达式中显然，利用隐式法求微分方程的数值解是，需求从表达式中反解未知节点上的函数值。反解未知节点上的函数值。1.3 隐式法的详细计算：隐式法的详细计算：例如欧拉梯形公式例如欧拉梯形公式用迭代法计算用迭代法计算n+1步的值。步的值。1简单迭代简单迭代收敛的条件：收敛的条件：11111 (,)(,)/2(,)/2(,)/2nnnnnnnnnnnyyf xyf xyhyf xyhf xyh(0)1(,)nnnnyyhf x

10、y(1)( )111(,)/2(,)/2kknnnnnnyyf xyhf xyh/21Lh2牛顿迭代牛顿迭代必需指出，在真正计算中，常用的是简单迭代，而且只迭代一必需指出，在真正计算中，常用的是简单迭代，而且只迭代一步，迭代初值步，迭代初值 称为预测，迭代称为校正，这样组成的称为预测，迭代称为校正，这样组成的一组计算公式称为预测校正公式。一组计算公式称为预测校正公式。( )( )(1)( )11111( )11(,)/2(,)/21(,)/2kkkknnnnnnnnkynnyyf xyhf xyhyyfxyh(0)1ny(0)1(0)111(,) (,)(,)/2nnnnnnnnnnyyhf

11、xyyyf xyf xyh预测校正公式也称为改良的欧拉法，将上面的组合公式改预测校正公式也称为改良的欧拉法，将上面的组合公式改写为：写为：留意到留意到 ，将上式进一步改写为：，将上式进一步改写为：这是我们最终运用的计算格式。这是我们最终运用的计算格式。11 (,)(,(,)/2nnnnnnnnyyf xyf xyhf xyh1121211()2(,)(,)nnnnnnyyKKKhf xyKhf xh yK1nnxxh例子：例子：取步长为取步长为0.1计算，结果如图。计算，结果如图。(01)(0)1yyxy 11212111()2(,)(,)()nnnnnnnnyyKKKhf xyhyKhf x

12、h yKh yK 图：图： DOUBLE PRECISION h,y(0:10),ak1,ak2 OPEN(20,FILE=OUTPUT1.DAT,STATUS=UNKNOWN)h=1.0/10y(0)=1.0do 10 i=1,10 ak1=-h*y(i-1) ak2=-h*(y(i-1)+ak1) y(i)=y(i-1)+(ak1+ak2)/2.010 continue do 20 i=0,10 write(20,*) i*h,y(i),exp(-i*h)20 continueEND同理，对于后退欧拉公式同理，对于后退欧拉公式有预测校正公式有预测校正公式或改写为：或改写为：111(,)nn

13、nnyyf xyh(0)1(0)111(,)(,)nnnnnnnnyyhf xyyyf xyh11(,)(,)nnnnnnKhf xyyyf xyK h用此法解前面的例子用此法解前面的例子步长步长0.1步长步长0.011.4 误差估计误差估计定义：利用第定义：利用第n个节点的函数值，经过数值方法计算第个节点的函数值，经过数值方法计算第n+1个个节点的函数近似值，所引起的误差，称为节点的函数近似值，所引起的误差，称为n+1个节点上的部个节点上的部分截断误差。分截断误差。我们记我们记 为为n+1个节点上函数的准确值，个节点上函数的准确值， 为数值为数值解，解，那么部分截断误差为：那么部分截断误差为

14、：如部分截断误差为如部分截断误差为 ，称为具有，称为具有p阶部分截断误差。阶部分截断误差。1ny1()ny x111()nnny xy()pO h欧拉方法的误差分析：欧拉方法的误差分析：12()()()()()()( )/2()()( ) /2nnnnnnnny xy xy xhy xhhy xy x hyhy xy xyhh1()()()( )nnny xy xy xO hh121()()()( )(, ()()()(, ()()nnnnnnnnny xy xy xO hf xy xhy xy xhf xy xO h省略余项得公式：省略余项得公式：1()(, ()nnnnyy xhf xy

15、x2111()()nnny xyO h完全类似的可以得到完全类似的可以得到后退欧拉公式的部分截断误差为：后退欧拉公式的部分截断误差为：欧拉中点公式的部分截断误差为：欧拉中点公式的部分截断误差为：欧拉梯形公式的部分截断误差为：欧拉梯形公式的部分截断误差为：2111()()nnny xyO h3111()()nnny xyO h3111()()nnny xyO h定义：由初值引起的第定义：由初值引起的第 n 个节点的近似值与准确值之间的误差个节点的近似值与准确值之间的误差称为第称为第 n 个节点整体误差。个节点整体误差。定理：设下面求解微分方程的数值计算方法部分截断误差为定理：设下面求解微分方程的

16、数值计算方法部分截断误差为p+1阶，且函数阶，且函数 关于关于 y 满足利普希茨条件，满足利普希茨条件，同时初值是准确的，那么整体截断误差为同时初值是准确的，那么整体截断误差为p阶。阶。欧拉公式、后退欧拉公式的整体截断误差为欧拉公式、后退欧拉公式的整体截断误差为 1 阶。阶。欧拉中点公式、欧拉梯形公式的整体截断误差为欧拉中点公式、欧拉梯形公式的整体截断误差为 2 阶。阶。1(, )nnnnyyhxy h( , , )x y h( , , )( , , )x y hx y hL yy微分方程数值解法的进一步改良。再回到恒等式微分方程数值解法的进一步改良。再回到恒等式假设取假设取 作为节点，将被积

17、函数用插值多作为节点，将被积函数用插值多项式来近似，用插值多项式带到积分中去求出积分，那么可以项式来近似，用插值多项式带到积分中去求出积分，那么可以得得到所谓的亚当斯到所谓的亚当斯(Adams)显式公式显式公式部分截断误差：部分截断误差：11()( )( , ( )iixiixy xy xf x y x dx1,iii kx xx1011()iiiiki khyyb fb fb fA1(1) ( )kkkiR yB hy类似地，假设取类似地，假设取 作为作为节点，可得亚当斯节点，可得亚当斯(Adams)隐式公式隐式公式部分截断误差：部分截断误差：1,i iiii kxx xx*10112()i

18、iiiiki khyyb fb fb fb fA*2(2)1 ( )kkkiR yBhy进一步，假设我们将恒等式中的积分区间改为进一步，假设我们将恒等式中的积分区间改为 ，并在此区间上用辛普森公式，可得，并在此区间上用辛普森公式，可得1,i iixx1111()()( , ( )iixiixy xy xf x y x dx111143iiiiihyyfff1.5 绝对稳定性绝对稳定性一个常微分方程数值解法称为绝对稳定，是指在某一步如第一个常微分方程数值解法称为绝对稳定，是指在某一步如第一步产生的误差如计算机的存储误差，在计算中会逐一步产生的误差如计算机的存储误差，在计算中会逐渐减小。渐减小。称

19、方程称方程 为实验方程，设计算步长为为实验方程，设计算步长为 h ，设在计算开，设在计算开场时产生误差存储误差，此误差在以后会逐渐减弱，我场时产生误差存储误差，此误差在以后会逐渐减弱，我们称该算法相对于们称该算法相对于 是绝对稳定的，是绝对稳定的， 的全体的全体称为该算法的绝对稳定域。称为该算法的绝对稳定域。 yy hhh欧拉法的绝对稳定区域欧拉法的绝对稳定区域 后退欧拉法的绝对稳定区域后退欧拉法的绝对稳定区域11h11h1.6 部分截断误差的适用估计部分截断误差的适用估计1用两种阶数一样的算法求解，计算出用两种阶数一样的算法求解，计算出n+1步的近似值，步的近似值，从而得到部分截断误差估计。

20、从而得到部分截断误差估计。2用同样的公式，用不同步长计算出用同样的公式，用不同步长计算出n+1步的近似值，从步的近似值，从而得到部分截断误差估计。而得到部分截断误差估计。1.7 隐式法隐式法隐式法具有较好的绝对稳定性！隐式法具有较好的绝对稳定性！只不过在运用隐式法的时候，需求进展迭代，或者运用预测只不过在运用隐式法的时候，需求进展迭代，或者运用预测校正计算格式。校正计算格式。第二节第二节 泰勒级数法与龙格库塔法泰勒级数法与龙格库塔法对于方程：对于方程：取计算步长为取计算步长为 h ，那么，那么 ，将函数进展泰勒展，将函数进展泰勒展开开如函数如函数 y( x )有有p+1阶导数，容易得到阶导数，

21、容易得到p阶泰勒级数展开法：阶泰勒级数展开法：0( )( , ( )( )y xf x y xaxby ay1nnxxh21()()()()()/2!nnnnny xy xhy xhy xh yx( )21(1)12!( )(, )(1)!ppnnnnnppnyyyyhyhhpyE x hhp公式中的导数用下面公式计算：公式中的导数用下面公式计算：例子：例子：2(,)(,)(,)(,)2(,)(,)(,)nnnnxnnynnnnxxnnxynnnyynnnynnnyf xyyfxyfxyyyfxyfxyyfxyyfxyy ( )cos04(0)1y xxxycos,sin,cosnnnnnny

22、xyxyx 231cossin/2cos/6nnnnnyyhxhxhx步长步长0.1 步长步长0.01龙格库塔法：龙格库塔法：对于常微分方程的数值解法，一个关键在于选择精度高的算法对于常微分方程的数值解法，一个关键在于选择精度高的算法计算下面公式中的积分。计算下面公式中的积分。要高精度的计算积分，常用的方法是适当添加计算节点，无妨要高精度的计算积分，常用的方法是适当添加计算节点，无妨设用设用 m 个节点计算上面积分，节点为个节点计算上面积分，节点为那么积分为那么积分为11()()( , ( )nnxnnxy xy xf x y x dx,1,2,inixxhim11( , ( )(, ()nn

23、mxininixif x y x dxhf xh y xh将积分改写为：将积分改写为：那么得公式：那么得公式：这样的公式称为显式龙格库塔公式。这样的公式称为显式龙格库塔公式。11( , ( )nnmxiixif x y x dxK11mnniiiyyK111(,)(,),2,3,nniininijjjKhf xyKhf xh yKim确定二阶确定二阶 RK 法法:系数满足：系数满足：此为四个未知数的三个方程，恣意取此为四个未知数的三个方程，恣意取 ，得，得11122122211(,)(,)nnnnnnyyKKKhf xyKhf xh yK12222121110,0,022 21222111,2

24、 特别，取特别，取 ，得到，得到通常也称为变形欧拉法，也常写为通常也称为变形欧拉法，也常写为它具有二阶精度，也称为二阶二级它具有二阶精度，也称为二阶二级 RK 方法。方法。1/212121(,)(/2,/2)nnnnnnyyKKhf xyKhf xhyK1(/2,(,)/2)nnnnnnyyhf xhyhf xy例子例子( )cos04(0)1y xxxy三阶三级库塔法三阶三级库塔法部分截断误差为部分截断误差为4阶，整体截断误差为阶，整体截断误差为3阶阶11231123121(4)6(,)(,)22(,2)nnnnnnnnyyKKKKhf xyKhKhf xyKhf xh yKK最常用的四级四

25、阶公式，称为龙格库塔公式：最常用的四级四阶公式，称为龙格库塔公式：部分截断误差为部分截断误差为5阶，整体截断误差为阶，整体截断误差为4阶。阶。1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf xyKhKhf xyKhKhf xyKhf xh yK隐式龙格库塔法：隐式龙格库塔法：111221122(,)1,2,nnmmininiiimmyyKKKKhf xh yKKKim常用的有二阶常用的有二阶RK法：法：111111(,)22nnnnyyKKhf xh yK例子：例子：用隐式二阶用隐式二阶RK法：法：(01)(0)1yyxy 11111

26、0( , ),(/2)1/21/21nnnnnnhf x yyKh yKKyhhyyKyyhy 210(1/2)1nnyyhhy用显式二阶用显式二阶RK法：法： DOUBLE PRECISION h,y(0:100),yy(0:100)double precision ak OPEN(20,FILE=OUTPUT5.DAT,STATUS=UNKNOWN)h=1.0/100y(0)=1.0yy(0)=1.0do 10 i=1,100 y(i)=y(i-1)-(h/(1-h/2)*y(i-1) yy(i)=yy(i-1)*(1-h+h*h/2)10 continue do 20 i=0,100 w

27、rite(20,100) i*h,y(i),yy(i),exp(-i*h)100 format (1x,f4.2,f8.5,f8.5,f8.5)20 continueENDh=0.1h=0.01半隐式龙格库塔法：半隐式龙格库塔法：111221122,1111,11(,)(,)1,2,nnmmininiii iiiininii iiyyKKKKhf xh yKKKg K hfxh yKKim1122,11111,11(,)1(,)ininiii iiininii iiKhf xh yKKKhg fxh yKK最常用的二级三阶半隐式龙格库塔公式：最常用的二级三阶半隐式龙格库塔公式：11211121

28、10.413154321.4131543261 (1)(,)(,)661 (1)(,)6(,)0.17378667nnynnnnynnnnyyKKKhhfxyf xyKhhfxh yKf xh yK微分方程组微分方程组一阶方程组一阶方程组11122212121010202000( )( ,( ),( ),( )( )( ,( ),( ),( )( )( ,( ),( ),( )(),(),()mmmmmmmy xf x y xyxyxyxfx y xyxyxyxfx y xyxyxy xyyxyyxy00( )( , ( )()y xf x y xy xy 1212010200(,) ,(,)

29、 ,(,)mmmyy yyffffyyyy龙格库塔公式龙格库塔公式1123411223431(22)6(,)(,)22(,)22(,)nnnnnnnnnnyyKKKKKhf xyKhKhf xyKhKhf xyKhf xh yK写成分量方式：写成分量方式：,1,12341121112121221222312411322331(22)6(,)(,)2222(,)2222(,)1,2,i ni niiiiiinnnmnmiinnnmnmiinnnmniinnnmnmyyKKKKKhf xyyyKKKhKhf xyyyKKKhKhf xyyyKhf xh yKyKyKim例子：例子：1122212(

30、01)(0)0,(0)1yyyyyxyy 1,11,1112131411121121121212221312141132231(262)() ()22 ()22 ()nnnnnnnnnnyyKKKKKhyyKKKhyyKKKhyyKhyKyK2,12,212223242122122222232242231(262)()2()2()nnnnnnyyKKKKKhyKKh yKKh yKh yK DOUBLE PRECISION h,y1(0:100),y2(0:100)double precision ak1,ak2,ak3,ak4,bk1,bk2,bk3,bk4 OPEN(20,FILE=OUT

31、PUT6.DAT,STATUS=UNKNOWN)h=1.0/10y1(0)=0.0y2(0)=1.0do 10 i=1,10 ak1=h*(-y1(i-1)+y2(i-1) bk1=-h*y2(i-1) ak2=h*(-(y1(i-1)+ak1/2.0)+y2(i-1)+bk1/2.0) bk2=-h*(y2(i-1)+bk1/2.0) ak3=h*(-(y1(i-1)+ak2/2.0)+y2(i-1)+bk2/2.0) bk3=-h*(y2(i-1)+bk2/2.0) ak4=h*(-(y1(i-1)+ak3)+y2(i-1)+bk3) bk4=-h*(y2(i-1)+bk3) y1(i)=y1(i-1)+(ak1+2*ak2+2*ak3+ak4)/6.0 y2(i)=y2(i-1)+(bk1+2*bk2+2*bk3+bk4)/6.010 continue do 20 i=0,10 write(20,100) i*h,y1(i),i*h*exp(-i*h),y2(i),exp(-i*h)100 format (1x,f