数值分析71常微分方程初值问题数值解引言

1、第第7章章 常微分方程初值问题数值解法常微分方程初值问题数值解法 7.1 引言引言 7.2 欧拉方法与改进的欧拉方法欧拉方法与改进的欧拉方法 7.3 龙格龙格- -库塔方法库塔方法7.1 引引 言言本章将考察本章将考察一阶方程的初值问题一阶方程的初值问题 )2 . 1(.)()1 . 1(),(00yxyyxfy 我们知道，只要我们知道，只要 f(x, y) 适当光滑适当光滑譬如关于譬如关于 y 满满足足利普希茨利普希茨(Lipschitz)条件条件理论上就可以保证初值问题的解理论上就可以保证初值问题的解 yf(x) 存在并且唯一存在并且唯一.),(),(yyLyxfyxf 虽然求解常微分方程

2、有各种各样的虽然求解常微分方程有各种各样的解析方法解析方法，但，但解析方法只能用来求解一些特殊类型的方程，实际问解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠题中归结出来的微分方程主要靠数值解法数值解法. 所谓所谓数值解法数值解法, 就是寻求解就是寻求解y(x)在一系列离散节点在一系列离散节点 121nnxxxx上的近似值上的近似值 y1, y2,yn, yn+1,. 相邻两个节点的间距相邻两个节点的间距hn=xn+1- -xn 称为称为步长步长. 假定假定 hi=h (i=1,2,)为为定数定数, 这时这时节点为节点为 xi=x0+ih (i=0,1,2,) (等

3、距节点等距节点). 初值问题的初值问题的数值解法数值解法有个有个基本特点基本特点: 都采取都采取“步步进式进式”，即求解过程顺着节点排列的次序一步一步地，即求解过程顺着节点排列的次序一步一步地向前推进向前推进. 首先，要对微分方程离散化，建立求解数值解的首先，要对微分方程离散化，建立求解数值解的递推公式递推公式. 一类是计算一类是计算 yn+1 时只用到前一点的值时只用到前一点的值 yn，称，称为为单步法单步法. 另一类是用到另一类是用到 yn+1前面前面 k 点的值点的值 yn, yn-1, yn-k+1，称为，称为 k步法步法. 其次，要研究公式的其次，要研究公式的局部截断误差局部截断误差

4、和和阶阶，数值解，数值解 yn 与精确解与精确解 y(xn) 的的误差估计误差估计及及收敛性收敛性，还有递推公式的还有递推公式的计算稳定性计算稳定性等问题等问题.7.2 欧拉方法与改进的欧拉方法欧拉方法与改进的欧拉方法7.2.1 欧拉法与后退欧拉法欧拉法与后退欧拉法 在在xoy平面上，微分方程平面上，微分方程(1.1)式式的解的解y=f(x)称作称作它的它的积分曲线积分曲线，积分曲线积分曲线上一点上一点(x, y)的切线斜率等的切线斜率等于函数于函数f(x, y)的值的值. 如果按如果按 f(x, y) 在在xoy平面上建立一平面上建立一个个方向场方向场，那么，那么，积分曲线积分曲线上每一点的

5、切线方向均上每一点的切线方向均与方向场在该点的方向相一致与方向场在该点的方向相一致.基于上述几何解释，我们从初始点基于上述几何解释，我们从初始点P0(x0, y0)出发出发,先依方向场在该点的方向推进到先依方向场在该点的方向推进到x=x1上一点上一点P1，然后，然后再从再从P1点依方向场在该点的方向推进到点依方向场在该点的方向推进到 x=x2 上一点上一点P2 , 循环前进做出一条循环前进做出一条折线折线P0 P1 P2. 一般地，设已做出该折线的顶点一般地，设已做出该折线的顶点Pn, ,过过Pn(xn, yn)依方向场的方向再推进到依方向场的方向再推进到Pn+1(xn+1, yn+1)，显然

6、两个，显然两个顶点顶点Pn, ,Pn+1的坐标有关系的坐标有关系斜斜率率),(111nnnnnnnnyxfhyyxxyy ) 1 .2 (),(1nnnnyxhfyy 这就是著名的这就是著名的( (显式显式) )欧拉欧拉( (Euler) )公式公式. . 若初值若初值y0已已知，则依公式知，则依公式(2.1)可逐次逐步算出各点数值解可逐次逐步算出各点数值解. .即即),(0001yxhfyy ),(1112yxhfyy 例例1 用欧拉公式求解初值问题用欧拉公式求解初值问题)2 . 2(. 1)0(),10(2 yxyxyy 解解 取步长取步长h=0.1，欧拉公式的具体形式为，欧拉公式的具体形

7、式为)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy依次计算下去，依次计算下去，部分计算结果部分计算结果见下表见下表. xy21 与准确解与准确解 相比，可看出欧拉公式的计算结相比，可看出欧拉公式的计算结果精度很差果精度很差. xn 欧拉公式数值解欧拉公式数值解yn准确解准确解y(xn) 误差误差 0.2 0.4 0.6 0.8 1.0 1.191818 1.3

8、58213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.052719 欧拉公式具有明显的几何意义欧拉公式具有明显的几何意义, , 就是就是用折线近似用折线近似代替方程的解曲线代替方程的解曲线，因而常称公式，因而常称公式(2.1)为为欧拉折线法欧拉折线法. .( )yy xxynx1nxnp1np1np x 还可以通过几何直观来考察欧拉方法的精度还可以通过几何直观来考察欧拉方法的精度. .假假设设yn=y(xn), ,即顶

9、点即顶点Pn落在积分曲线落在积分曲线y=y(x)上，那么，上，那么，按欧拉方法做出的折线按欧拉方法做出的折线PnPn+1便是便是y=y(x)过点过点Pn的切线的切线. .从图形上看从图形上看, ,这这样定出的顶点样定出的顶点Pn+1偏离偏离了原来的积分曲线，可了原来的积分曲线，可见欧拉方法是见欧拉方法是相当粗糙相当粗糙的的. . 为了分析计算公式的精度，通常可用为了分析计算公式的精度，通常可用泰勒展开泰勒展开将将y(xn+1)在在xn处展开，则有处展开，则有).,()(2),()()(2)()()()(1221 nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy 在在yn=y

10、(xn)的前提下，的前提下，f(xn,yn )=f(xn,y(xn)=y ( (xn n) ). .于是于是可得欧拉法可得欧拉法(2.1)的的公式误差公式误差为为)3 . 2(),(2)(2)(2211nnnnxyhyhyxy 称为此方法的称为此方法的局部截断误差局部截断误差. . 如果对方程如果对方程(1.1)从从xn到到xn+1积分，得积分，得) 4 . 2(.)(,()()(11 nnxxnndttytfxyxy右端积分用右端积分用左矩形公式左矩形公式 hf(xn,y(xn) 近似，再以近似，再以 yn代代替替 y(xn)，yn+1代替代替 y(xn+1)也得到欧拉公式也得到欧拉公式(2

11、.1)，局部，局部截断误差也是截断误差也是(2.3).称为称为( (隐式隐式) )后退的欧拉公式后退的欧拉公式. . 如果右端积分用如果右端积分用右矩形公式右矩形公式 hf(xn+1,y(xn+1)近似，近似，则得到另一个公式则得到另一个公式) 5 . 2 (),(111 nnnnyxhfyy 后退的欧拉公式与欧拉公式有着后退的欧拉公式与欧拉公式有着本质的区别本质的区别, 后后者是关于者是关于yn+1的一个直接计算公式，这类公式称作是的一个直接计算公式，这类公式称作是显式的显式的；前者公式的；前者公式的右端含有未知的右端含有未知的yn+1，它实际上，它实际上是关于是关于yn+1的一个函数方程的

12、一个函数方程, ,这类方程称作是这类方程称作是隐式的隐式的. . 显式显式与与隐式隐式两类方法各有特点，考虑到数值稳两类方法各有特点，考虑到数值稳定性等其他因素，人们有时需要选用定性等其他因素，人们有时需要选用隐式隐式方法，但方法，但使用使用显式显式算法远比算法远比隐式隐式方便方便. .注：注：隐式方程通常用隐式方程通常用迭代法迭代法求解，而迭代过程的实求解，而迭代过程的实质是质是逐步逐步显式化显式化. . 设用欧拉公式设用欧拉公式),() 0(1nnnnyxhfyy 给出迭代初值给出迭代初值 ，用它代入，用它代入(2.5)式的式的右端，使之转右端，使之转化为显式，直接计算得化为显式，直接计算

13、得) 0(1 ny),() 0(11) 1 (1 nnnnyxhfyy然后再用然后再用 代入代入(2.5)式，又有式，又有) 1 (1 ny).,() 1 (11) 2(1 nnnnyxhfyy如此反复进行，得如此反复进行，得) 6 . 2 ()., 1 , 0(),()(11) 1(1 kyxhfyyknnnkn如果如果 f(x, y)对对y满足满足Lipschitz条件条件(1.3). 由由(2.6)减减(2.5)得得),(),(11)(111) 1(1 nnknnnknyxfyxfhyy.1)(1 nknyyhL由此可知，只要由此可知，只要hL1，迭代法，迭代法(2.6)就收敛到解就收敛

14、到解. .7.2.2 梯形方法梯形方法 为得到比为得到比欧拉法精度高的计算公式，在等欧拉法精度高的计算公式，在等式式(2.4) 右端积分用梯形求积公式近似右端积分用梯形求积公式近似, 并用并用yn代替代替y(xn), yn+1代替代替y(xn+1)，则得，则得 ) 7 . 2(,),(),(2111 nnnnnnyxfyxfhyy称为称为矩形方法矩形方法. 矩形方法矩形方法是是隐式单步法隐式单步法，用迭代法求解，同后，用迭代法求解，同后退的欧拉方法一样，仍用退的欧拉方法一样，仍用欧拉法欧拉法提供迭代初值，则提供迭代初值，则矩形迭代公式矩形迭代公式为为 ) 8 . 2()., 1 , 0(),(

15、),(2);,()(11) 1()0(11 kyxfyxfhyyyxhfyyknnnnnknnnnn为了分析迭代过程的收敛性为了分析迭代过程的收敛性, 将将(2.7)与与(2.8)相减相减, 得得),(),(2)(1111) 1(11knnnnknnyxfyxfhyy 于是有于是有,2)(11) 1(11knnknnyyhLyy 使得使得, 12 hL则当则当k时有时有 , 这说明迭代过程这说明迭代过程(2.8)是收敛是收敛的的.1)(1 nknyy7.2.3 改进的欧拉公式改进的欧拉公式我们看到，梯形方法虽然提高了精度，但其算法我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式复

16、杂，在应用迭代公式(2.8)进行实际计算时，每迭代进行实际计算时，每迭代一次，都要重新计算函数一次，都要重新计算函数f(x, y )的值，而迭代又要反的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测复进行若干次，计算量很大，而且往往难以预测. 为为了控制计算量，通常只迭代一两次就转入下一步的计了控制计算量，通常只迭代一两次就转入下一步的计算，这算，这就简化了算法就简化了算法. .具体地说，我们先用欧拉公式求得一个初步的近具体地说，我们先用欧拉公式求得一个初步的近似值似值 ，称之为，称之为预测值预测值，此预测值，此预测值 的精度可能很的精度可能很差，再用梯形公式差，再用梯形公式(2.

17、7)将它校正一次，即按将它校正一次，即按(2.8)式迭式迭代一次，这个结果代一次，这个结果称之为称之为校正值校正值. .1 ny1 ny这样建立的这样建立的预测预测校正系统校正系统通常称为通常称为改进的欧拉公式改进的欧拉公式:或表为下列或表为下列平均化形式平均化形式 ).(21),(),(11cpnpnncnnnpyyyyxhfyyyxhfyy( (2.9) ), ),(1nnnnyxhfyy 预测预测.),(),(2111 nnnnnnyxfyxfhyy校正校正例例2 用改进的欧拉法解用改进的欧拉法解例例1中的初值问题中的初值问题(2.2). . 解解 仍取步长仍取步长h=0.1, ,改进的

18、欧拉公式为改进的欧拉公式为 ).(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy部分计算结果部分计算结果见下表见下表 xnyn 误差误差 xnyn 误差误差00.20.4 1. 1.184096 1.34336000.0000880.0017190.60.81.01.4859561.6164761.7378690.0027160.0040240.05818同例同例1中的欧拉法的计算结果比较，明显改善了精度中的欧拉法的计算结果比较，明显改善了精度. 例例3 3 (两种方法的精度比较两种方法的精度比较)用欧拉公式和改进的欧拉公式解下述初值问题，并与用欧拉公式和改进

19、的欧拉公式解下述初值问题，并与其准确解其准确解y=e- -x+x进行比较进行比较. . 1) 0(6 . 00, 1yxxyy解解 取步长取步长h=0.1，xk=kh (k=0,1,6). .用两种方法用两种方法进行计算对应结果及绝对误差见下表进行计算对应结果及绝对误差见下表 xn 欧拉公式欧拉公式 改进欧拉公式改进欧拉公式 yn 误差误差 yn 误差误差 0.0 0.1 0.2 0.3 0.4 0.5 0.6 1.000000 1.000000 1.010000 1.029000 1.056100 1.090490 1.13144104.810-38.710-3 1.210-21.410-2

20、1.610-21.810-2 1 1.005000 1.019025 1.041218 1.070802 1.107076 1.14940401.610-42.910-44.010-44.810-45.510-4 5.910-47.2.3 单步法的局部截断误差与阶单步法的局部截断误差与阶 初值问题初值问题(1.1),(1.2)的的单步法单步法可用可用一般形式一般形式表示为表示为11( , ),(2.10)nnnnnyyhx y yh其中多元函数其中多元函数 与与f(x, y )有关，当有关，当 含有含有yn+1时，方法时，方法是是隐式隐式的，若不含的，若不含yn+1则为则为显式方法显式方法，所

21、以，所以显式单步显式单步法法可表示为可表示为 (x, y, h)称为称为增量函数增量函数，例如对欧拉法，例如对欧拉法(2.1)有有1( , , ),(2.11)nnnnyyh x y h. ),(),(yxfhyxnn 它的局部截断误差已由它的局部截断误差已由(2.3) 给出给出, 对一般显式单步对一般显式单步法则可如下定义法则可如下定义. . 定义定义1 设设y(x)是初值问题是初值问题(1.1),(1.2)的的准确解准确解, 称称11()( )( , ( ), ).(2.12)nnnnnTy xy xh x y xh为显式单步法为显式单步法( (2.11) )的的局部截断误差局部截断误差.

22、 . Tn+1之所以称为之所以称为局部的局部的，是假设在，是假设在 xn前各步前各步没没有误差有误差. .当当 yn= =y(xn) 时，计算一步，则有时，计算一步，则有.),(,()()(),()()(11111 nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy 所以，局部截断误差可理解为用方法所以，局部截断误差可理解为用方法(2.11)计算计算一步一步的误差，也即公式的误差，也即公式(2.11)中用准确解中用准确解y(x)代替数代替数值解产生的公式误差值解产生的公式误差. 根据定义根据定义, 显然显然欧拉法的局部欧拉法的局部截断误差截断误差为为).()(2)()()()(,()()(3211hOxyhxyhxyhxyxyxhfxyx