1、第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法第第6章章 常微分方程数值解法常微分方程数值解法 1 引言引言 2 欧拉法和改进的欧拉法欧拉法和改进的欧拉法3 龙格龙格-库塔法库塔法4 阿当姆斯方法阿当姆斯方法第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法1 引言引言 在高等数学里我们已经接触过常微分方程，对于一在高等数学里我们已经接触过常微分方程，对于一些典型的常微分方程，有求解析解的基本方法，但多数些典型的常微分方程，有求解析解的基本方法，但多数情况下，遇到的问题比较复杂，此时，只能利用近似方情况下，遇到的问题比较复杂，此时，只能利用近似方法求解，一般有两种近
2、似方法法求解，一般有两种近似方法。近似解析方法近似解析方法数值方法数值方法第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法实际求解的常微分方程，大多是实际求解的常微分方程，大多是定解问题定解问题满足指定满足指定条件的特解条件的特解 初值问题初值问题 边值问题边值问题 本章讨论常微分方程，数值解的最简单问题本章讨论常微分方程，数值解的最简单问题 一阶方程一阶方程初值问题初值问题，即函数，即函数f(x)满足下列微分方程和初值条件：满足下列微分方程和初值条件：在几何问题是在几何问题是(6-1)表现为一簇曲线，称表现为一簇曲线，称(6-1)的积分曲线，的积分曲线，初值问题初值问题(6-1
3、)(6-2)就是要求一条过就是要求一条过(x0,y0)的积分曲线的积分曲线00(,)()dyfxydxyxy(61)(62)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法方程的精确解方程的精确解y(x)称称为积分曲线。为积分曲线。方程是否有解，解是否唯一？方程是否有解，解是否唯一？第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法定理定理1 对初值问题对初值问题(6-1)(6-2)，若，若f(x,y)在区域在区域G=axb ,|y|内连续，且关于内连续，且关于y满足李普希兹条件，即存在常数满足李普希兹条件，即存在常数L，使，使|f(x,y1)-f(x,y2)|L|y
4、1-y2|（6-3）对对G中任意两个中任意两个y1,y2均成立，其中均成立，其中L是与是与x,y无关的常数，无关的常数，则初值问题则初值问题(6-1)(6-2)在在(a,b)内存在唯一解，且解是连续内存在唯一解，且解是连续可微的。可微的。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法设设f(x，y)在带形区域在带形区域R：axb，-y+上为上为x，y的连续函数，且对任意的的连续函数，且对任意的y满足李普希茨满足李普希茨(Libusize)条件条件 f(x，y1)-f(x，y2)Ly1-y2 其中其中(x，y1)、(x，y2)R，L为正常数。为正常数。在求初值问题在求初值问题(6
5、-1)(6-2)的数值解时，我们通常采用的数值解时，我们通常采用离离散化方法散化方法，求在自变量，求在自变量x的离散点的离散点a=x0 x1x2xn=b 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 上的上的准确解准确解y(x)的近似值的近似值 y0，y1，y2，yn常取离散点常取离散点x0，x1，x2，xn为等距，即为等距，即 x i+1-xi=h，i=0，1，2，n-1 h称为步长。图称为步长。图6.1表示为初值问题表示为初值问题(61)(62)在在n+1个离散点上的准确解个离散点上的准确解y(x)的近似值。的近似值。第第6 6章章 常微分方程数值解法常微分方程数值解法
6、计算方法图 6.1第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法数值解法的重点不在于求准确解（即解析解），而是直接数值解法的重点不在于求准确解（即解析解），而是直接求求一系列点上的近似解一系列点上的近似解。求解过程顺着节点排列的顺序一步步向前推进求解过程顺着节点排列的顺序一步步向前推进，也即按递推公式由也即按递推公式由y0，y1.yi推推yi+1，下面各种方法的，下面各种方法的实质是建立实质是建立递推公式递推公式。初值问题（初值问题（6.1）(6.2)的数值解法的基本特点是：的数值解法的基本特点是：第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法2 欧拉法和改进的
7、欧拉法欧拉法和改进的欧拉法 一、一、欧拉方法欧拉方法1.基本思想基本思想 区间区间a，b上给定上给定n+1个点个点x0，x1，x2，xn),()()()(x),(1iiiiiiyxfxyhxyxyyxfdxdy用差商的导数考虑在节点),()()(1iiiiyxhfxyxy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 再用再用yi近似地代替近似地代替y(xi)，则初值问题，则初值问题(6-1)(6-2)就化为就化为从从x0出发根据初值问题，出发根据初值问题，y(x0)=y0 再利用上式得再利用上式得y(x1)y1=y(0)+hf(x0，y0)，再以再以y1作为作为y(x1)的近
8、似值，代入上式求的近似值，代入上式求y2.yn y(x2)y2=y1+hf(x1,y1).y(xi+1)yi+1=yi+hf(xi,yi)i=0,1,(6.4)称为解初值问题的称为解初值问题的欧拉方法欧拉方法 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法2几何意义几何意义 欧拉公式有很明显的几何意义。我们知道初值问题欧拉公式有很明显的几何意义。我们知道初值问题(6.1)中的微分方程的解是中的微分方程的解是xoy平面上的一簇积分曲线平面上的一簇积分曲线这簇积分曲线上任意点这簇积分曲线上任意点(x，y)的斜率为的斜率为f(x，y)，而初值问题而初值问题(6.1)(6.2)的解是过
9、点的解是过点(x0，y0)的一条特定的积的一条特定的积分曲线。分曲线。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法Ox0 x1x2xn(x0,y0)(x1,y1)(x2,y2)(xn,yn)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法y y(x)过点过点P0(x0,y0)，从，从P0出发以出发以f(x0,y0)为斜率做一直线与直线为斜率做一直线与直线x=x1交于点交于点p1(x1,y1),显然有显然有:y1=y0+hf(x0,y0),再从再从p1出发出发,以以f(x1,y1)为斜率做一直线推进到为斜率做一直线推进到x=x2上一点上一点p2(x2,y2),依此
10、类推依此类推,这样得到解曲线的一条近似曲线这样得到解曲线的一条近似曲线,它就是折它就是折线线p0p1p2所以欧拉方法又叫所以欧拉方法又叫欧拉折线法欧拉折线法 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法欧拉法是用欧拉法是用yi通过通过 yi+1=y i+hf(x i,y i)i=0,1,求求yi+1,这样利用这样利用y0y1y2计算计算yi+1用前一步的用前一步的y i单步法单步法计算计算yi+1用前几步的用前几步的yi-n多步法多步法 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法例例1：用欧拉法求解方程：用欧拉法求解方程 1)0(2yxyy0 x1.2 h
11、=0.2解解:欧拉法的具体形式为欧拉法的具体形式为:yi+1=yi+hf(xi,yi)=(1-0.4xi)yi 所以所以:y1=y0+h f(x0,y0)y2=y1+h f(x1,y1)=(1-0.4x1)y1 =0.920000 =(1-0.4x0)y0=1第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法xiyi y(xi)(xi)0110.210.9607890.40.9200000.8521440.60.7728000.6976760.80.5873220.5277921.00.3993830.3678791.20.2396300.236938可见欧拉法的精度是很差的可见欧
12、拉法的精度是很差的 所求值用下表列出所求值用下表列出,并与精确值对比并与精确值对比 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法二、欧拉方法的误差分析二、欧拉方法的误差分析定义定义1(p146)对于初值问题，当假设对于初值问题，当假设yi是准确的时，用某是准确的时，用某种方法求种方法求yi+1时所产生的截断误差称为该方法的时所产生的截断误差称为该方法的局局部截断误差部截断误差。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法我们来看在第我们来看在第i+1步使用欧拉方法所得步使用欧拉方法所得yi+1的局部截断误差的局部截断误差y(xi+1)-yi+1 假定假定yi
13、是准确的是准确的,即即yi=y(xi)由由y(xi+1)=y(xi+h)，应用泰勒展开，应用泰勒展开y(xi+1)=y(xi+h)=y(xi)+hy(xi)+y()/2*h2而由欧拉公式算出而由欧拉公式算出yi+1=yi+hf(xi,yi)=y(xi)+hf(xi,y(xi)=y(xi)+hy(xi)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法y(xi+1)=y(xi+h)=y(xi)+hy(xi)+y()/2*h2yi+1=yi+hf(xi,yi)=y(xi)+hf(xi,y(xi)=y(xi)+hy(xi)两式相减得两式相减得 y(xi+1)-yi+1=(h2/2)*y(
14、)=O(h()=O(h2 2)即欧拉方法所得即欧拉方法所得y yi+1i+1的局部截断误差为的局部截断误差为O(hO(h2 2)注意注意:只计算了一步，事实上每一步都有可能产生误差，只计算了一步，事实上每一步都有可能产生误差，有时误差会原来越大，有时又会得到很好的控制，因此还有时误差会原来越大，有时又会得到很好的控制，因此还要考虑整体截断误差。要考虑整体截断误差。)(22yh称为局部截断误差的主项称为局部截断误差的主项第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法定义定义2(p147)设设y yi i是用某种方法计算初值问题是用某种方法计算初值问题(6-1)(6-2)(6-1)
15、(6-2)在在x xi i点的近点的近似解，而似解，而y(xy(xi i)是它的精确解，则称是它的精确解，则称为该方法的整体截断误差，也称为该方法的精度。为该方法的整体截断误差，也称为该方法的精度。补：若某方法的局部截断误差为补：若某方法的局部截断误差为O(hp+1)，则该方法的精度，则该方法的精度为为p阶的。阶的。iiiyxy)(欧拉方法的精度为欧拉方法的精度为O(h)，一阶的，一阶的第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法三、改进的欧拉法三、改进的欧拉法 欧拉法虽然形式简单，计算方便，但比较粗糙，欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当精度也低。特别
16、当y=y(x)的曲线曲率较大时，欧拉法的效的曲线曲率较大时，欧拉法的效果更差。果更差。为了构造较高精度的数值解法，对初值问题为了构造较高精度的数值解法，对初值问题再做分析再做分析00)(),(yxyyxfy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法对对y=f（x，y）等式两边在（）等式两边在（xi，xi+1）上取积分）上取积分 11)(,(iiiixxxxdxxyxfdxy1)(,()()(1iixxiidxxyxfxyxyxixi+1(xi,yi)(xi+1,yi+1)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法在上图也可利用在上图也可利用),(1iii
17、iyxhfyy矩形公式精度不高，再次说明欧拉精度低矩形公式精度不高，再次说明欧拉精度低后退欧拉公式右矩形公式),(111iiiiyxhfyy这样就利用数值积分公式计算这样就利用数值积分公式计算y(xi+1)的近似值。如的近似值。如果用左矩阵计算右面的积分：果用左矩阵计算右面的积分：)8.6()(,()()(11iixxiidxxyxfxyxy显式欧拉公式),(1iiiiyxhfyy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法为了得到更精确的方法我们可使用梯形公式为了得到更精确的方法我们可使用梯形公式),(),(2),(111iiiixxyxfyxfhdxyxfii,.2,1,
18、0),(),(2111iyxfyxfhyyiiiiii此时(6-9)改进的欧拉公式，又称为改进的欧拉公式，又称为梯形公式梯形公式xixi+1(xi,yi)(xi+1,yi+1)这样得到的点列仍这样得到的点列仍为一折线，只是用为一折线，只是用平均斜率来代替原平均斜率来代替原来一点处的斜率。来一点处的斜率。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法不难发现不难发现:欧拉公式欧拉公式 yi+1=yi+hf(x i,yi)是关于是关于yi+1的显式的显式，只要已知，只要已知yi，经一次计算可立即得到，经一次计算可立即得到
19、yi+1的值；的值；而而改进欧拉公式改进欧拉公式,.2,1,0),(),(2111iyxfyxfhyyiiiiii中的中的yi+1以隐式给出以隐式给出，且，且yi+1含在函数含在函数f(xi+1，yi+1)中，中，所以所以梯形法是梯形法是隐式单步法隐式单步法，第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 一般来说，这是一个非线性方程（除非一般来说，这是一个非线性方程（除非f对对y是是线性的），可用我们前面讲过的非线性方程的各线性的），可用我们前面讲过的非线性方程的各种方法求解，比如用迭代法种方法求解，比如用迭代法 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法
20、),(),(2),()(11)1(1)0(1kiiiiikiiiiiyxfyxfhyyyxhfyy可证明其收敛，局部截断误差可证明其收敛，局部截断误差O(h3)k=0,1,2,具体做法是：先用欧拉公式具体做法是：先用欧拉公式(64)求出一个求出一个y(0)i+1作作为初始近似值，然后再用改进的欧拉公式为初始近似值，然后再用改进的欧拉公式(69)进进行迭代，即行迭代，即直到满足(1)()11(1)11kkiikiiyyyy取再转到下一步计算第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法例例2:证明解常微分方程初值问题的梯形方法精度是二阶的证明解常微分方程初值问题的梯形方法精度是二
21、阶的证明证明:),(),(211iiiiihyxfyxfhyT令T=y(xi+1)-Th=y(x i+h)-T h)(2)(2)(.!3)(!2)()()(T32hiiiiiiixyhxyhxyhxyhxyxhyxy则.!3 !2)()()()()4(32yhyhxhyxyxytailorxxxyiihiihi展开处在将第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法所以所以:)()(12)()(4)(!343433hoxyhhoxyhxyhT局部截断误差为局部截断误差为O(h3)所以精度为所以精度为2阶的阶的第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法四、预报
22、四、预报-校正方法校正方法 我们看到梯形法虽然提高了精度，但其算法复杂，我们看到梯形法虽然提高了精度，但其算法复杂，每算一点，都需进行反复迭代，为了控制计算量，通常每算一点，都需进行反复迭代，为了控制计算量，通常只迭代一两次就转入下一步计算，以简化算法。只迭代一两次就转入下一步计算，以简化算法。具体地说，我们先用欧拉公式求一个初步的近似值具体地说，我们先用欧拉公式求一个初步的近似值),(1iiiiyxhfyy 称为预测值，其精度可能很差，再用梯形公式将其称为预测值，其精度可能很差，再用梯形公式将其校正一次校正一次 1iy),(),(2111iiiiiiyxfyxfhyy第第6 6章章 常微分方
23、程数值解法常微分方程数值解法 计算方法预报预报-校正公式校正公式),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy(6-11)在实际计算时，还常常将式在实际计算时，还常常将式(611)写成下列形式：写成下列形式：(6-12),.2,1,0)(2),(),(*2*111*2*1ikkhyyhkyhxfkyxfkiiiiii第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法例例3 用预报用预报-校正公式求解初值问题校正公式求解初值问题 2(0)1,01,0.1xyyyyxh解：由解：由预报-校正公式有公式有)2()2(2211111iiiiiiiiiiiii
24、yxyyxyhyyyxyhyyh=0.1,x0=0,x1=0.1,x2=0.2,.,x9=0.9,x10=1 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 1.1)101(1.01)2(00001yxyhyyi=0095909.1)2()2(211100001yxyyxyhyyi=1191818.1)2(11112yxyhyy184097.1)2()2(222211112yxyyxyhyyi=2最后最后,得下表得下表 其解析解为其解析解为 21yx第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 表 62 第第6 6章章 常微分方程数值解法常微分方程数值解法 计
25、算方法预报预报-校正公式的截断误差为校正公式的截断误差为O(h3)（证明，（证明，p150)预报预报-校正公式的整体截断误差为校正公式的整体截断误差为O(h2)预报预报-校正公式是校正公式是单步显式单步显式第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法练习练习1:用用Euler法解法解1)0()6.00(2yxxyydxdy取取 h=0.2 计算到计算到x=0.6第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法练习练习2:用改进用改进Euler法（梯形公式）解法（梯形公式）解2)1()21(38yxydxdy取取 h=0.2 小数点后至少保留小数点后至少保留5位位
26、第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法练习练习3:用预报用预报校正方法求校正方法求 1)1(0sin2yxyyy取步长取步长 h=0.2 计算计算y(1.2)及及y(1.4)的近似值，小数点后保留的近似值，小数点后保留5位位第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法3 龙格库塔法龙格库塔法 一、基本思想一、基本思想 对初值问题对初值问题00)(),(yxyyxfdxdy(axb)目的是求目的是求y(x)在一些给定点上的值在一些给定点上的值y(x1)，y(x2)，y(xn)研究差商研究差商微分中值定理微分中值定理第第6 6章章 常微分方程数值解法常微分
27、方程数值解法 计算方法)(,()()()()(1hxyhxhfxyhxhyxyxyiiiiii得到一种微分方程的数值计算公式，得到一种微分方程的数值计算公式，关键在求关键在求)(,(hxyhxfii称为在区间称为在区间xi,xi+1上的上的平均斜率平均斜率，记为记为)(,(*hxyhxfKii*1)()(hKxyxyii则K*第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法若用这种方法研究欧拉公式，可以发现欧拉公式仅取一个若用这种方法研究欧拉公式，可以发现欧拉公式仅取一个点点(xi,yi)的斜率的斜率f(xi,yi)代替代替k*，梯形公式却是利用，梯形公式却是利用xi与与xi+1
28、两两个点的斜率值个点的斜率值),(),(11*2*1iiiiyxfkyxfk然后取值然后取值)(21*2*1*kkk由此也可知，梯形公式精度要高些。由此也可知，梯形公式精度要高些。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法它也启发我们它也启发我们若设法在若设法在xi,xi+1上多找几个点的斜率值，上多找几个点的斜率值，然后将它们然后将它们加权平均加权平均作为作为K*的近似值，则有可能构造出的近似值，则有可能构造出精度更高的计算公式精度更高的计算公式，这就是，这就是Runge-Kutta方法的基本思方法的基本思想。想。即即为权系数其中jjijirjjhyhxfK),(1*第第
29、6 6章章 常微分方程数值解法常微分方程数值解法 计算方法于是得到一般公式于是得到一般公式),(11hyhxfhyyjijirjjii其中其中j，j，j为为待定系数待定系数(6-14)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法二、二、Runge-Kutta公式公式 我们以两点我们以两点xi，xi+p=xi+h为例，说明为例，说明(6-14)式中系数式中系数的求法的求法),(),(12122111kyhxhfkyxhfkkkyyiiiiii其中其中1，2，为待定参数为待定参数(6-15)(ixhy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法由二元函数由二元函
30、数Taylor展开式有展开式有)()(),(),()()(),(),(),(32212hOxyyxfyxfhxhyhOyxfkyxhfyxfhkiiiyiixiiiyiixii代入式代入式(6-15)，得到，得到yi+1的展开式的展开式)()(),(),()()(3222211hOxyyxfyxfhxyhyyiiiyiixiii(6-16)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法利用利用),(),()(),(yyxfyxfxyyxfyyx得到得到)()(),(),(2)()()()(2)()()(32321hOxyyxfyxfhxhyxyhOxyhxhyxyxyiiiyi
31、ixiiiiii第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法把它与式把它与式(6-16)比较，使它们比较，使它们前三项前三项相等，相等，得得212112221则则y(xi+1)-yi+1=O(h3)(6-17)满足满足(6-17)式的解有无穷多，对应于式的解有无穷多，对应于(6-15)的公式均具的公式均具有二阶精度，统称为有二阶精度，统称为二阶龙格二阶龙格-库塔公式库塔公式。第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法特别地，取特别地，取则则(6-15)式正好是式正好是预报预报-校正公式校正公式。如取如取)21,21(),(12121kyhxhfkyxhfk
32、kyyiiiiii(6-18)1,21,1211,1,021则有则有第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法三阶龙格三阶龙格-库塔公式为库塔公式为)2,()2/,2/(),()4(612131213211kkyhxhfkkyhxhfkyxhfkkkkyyiiiiiiii局部截断误差为局部截断误差为O(h4)(6-19)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法四阶龙格四阶龙格-库塔公式为库塔公式为),()2/,2/()2/,2/(),()22(61342312143211kyhxhfkkyhxhfkkyhxhfkyxhfkkkkkyyiiiiiiiii
33、i(6-20)此公式又称为此公式又称为经典的龙格经典的龙格-库塔公式库塔公式，局部截断误差为，局部截断误差为O(h5)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法例例 4 用经典的龙格用经典的龙格-库塔法计算库塔法计算 取步长取步长h=0.2解：由x0=0，y0=1，h=0.2，利用公式(6-20)可计算出1)0(9.0,0,/2yxyxyy16965.0),(18173.0)2/,2/(18364.0)1.1,1.0(2.0)2/,2/(2.0)121(2.0),(3004200310020001kyhxhfkkyhxhfkfkyhxhfkxyxhfk从而从而18323.1
34、)22(61432101kkkkyy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法表表6-3ixiyi00110.21.1832320.41.34166730.61.48328140.81.61251351.01.732140第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法补充补充-单步法的收敛性和稳定性单步法的收敛性和稳定性计算过程中舍入误差总是存在的，以计算过程中舍入误差总是存在的，以Euler法为例。法为例。假设，由于舍入误差的影响，实际得到的是假设，由于舍入误差的影响，实际得到的是.iiiyy),(1iiiiyxhfyy),(1iiiiyxhfyy两式相减
35、两式相减,可得可得:),(),(1iiiiiiyxfyxfh第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法由微分中值定理得由微分中值定理得:),(),(1iiiiiiyxfyxfh由此可知由此可知:只有当只有当1|),(1|nyxhf误差不增误差不增,算法稳定算法稳定第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法1.单步法的收敛性单步法的收敛性补充定义补充定义1:设设y(x)是初值问题是初值问题(6.2)的精确解的精确解,对单步法对单步法),(1hyxgyyiiii若有若有:ixyyiih),(lim0则称单步法收敛则称单步法收敛第第6 6章章 常微分方程数值解
36、法常微分方程数值解法 计算方法补充定义补充定义2：用单步法解模型方程用单步法解模型方程 得到解满足稳定性方程得到解满足稳定性方程 若若 ,就称此方法是就称此方法是绝对稳定绝对稳定的的 在复平面上,所有满足 所围成的区域称为方所围成的区域称为方法的绝对稳定区域法的绝对稳定区域.0)Re(,yynnyhEy)(11|)(|hE1|)(|hE的交称为绝对稳定区间在实轴为实数时，此区域与（当h)0第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法例例5：求：求Euler法的绝对稳定区间法的绝对稳定区间 yyxfy),(解：),(1iiiiyxhfyyiiiiyhEyhyhy)()1(hh
37、1)E（所以021|1|)E|hhh得（由第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法6.4 阿达姆斯方法阿达姆斯方法 我们已经知道，初值问题我们已经知道，初值问题(61)等价于积分方程，等价于积分方程，即即1),()()(1iixxiidxyxfxyxy11),(iiiixxxxdxyxfdxy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 对积分式分别采用对积分式分别采用矩形矩形和和梯形面积公式梯形面积公式 可得到欧拉公式和改进欧拉公式可得到欧拉公式和改进欧拉公式(梯形法梯形法)，截断误差分别为截断误差分别为O(h2)和和O(h3)。若追溯到数值积分原理，
38、梯形法则若追溯到数值积分原理，梯形法则 实际上是用节点实际上是用节点 x i和和x i+1 的线性插值函数代替的线性插值函数代替f(x，y)而得到的。而得到的。为此，我们自然可以想到，若用为此，我们自然可以想到，若用更高次的插值多项式来更高次的插值多项式来代替代替f(x，y)，则所得公式的精度会更高，则所得公式的精度会更高。-这就是线性多步法的这就是线性多步法的起源思想起源思想。(624)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 对于对于线性多步法线性多步法 是要利用前面已经算出的若干个值是要利用前面已经算出的若干个值 yi-k，yi-1，yi来求yi+1。现用k次多项式
39、Pk(x)来代替f(x，y)有:f(x,y)=P k(x)+R k(x)111()()()()iiiixxiikkxxy xy xP x dxRx dx1),()()(1iixxiidxyxfxyxy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法 舍去余项舍去余项 并设并设yi=y(xi)，而，而yi+1为为y(xi+1)的近似值，于是可得的近似值，于是可得到线性多步法的计算公式到线性多步法的计算公式1()iixkkxRRx dx11()iixiikxyyP k dx第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法一、一、阿达姆斯阿达姆斯(Adams)显式显式（外
40、推）公式（外推）公式 设设xi=x0+ih,i=0,1,2为等距分布，步长为为等距分布，步长为h,选选 x i-3 ,x i-2,x i-1,x i 为插值节点为插值节点 做做f(x,y)的三次拉格朗日的三次拉格朗日插值多项式插值多项式L3(x)代入代入得：称为称为Adams 四步显式公式，也称四步显式公式，也称Adams外推公式外推公式P201(6-22)1),()()(1iixxiidxyxfxyxy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法局部截断误差 为：=O(h5)第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法二、二、阿达姆斯阿达姆斯(Adams)
41、隐式隐式（内插）公式（内插）公式选选 x i-2,x i-1,x i,x i+1 为插值节点为插值节点 做做f(x,y)的三次拉格朗日的三次拉格朗日插值多项式插值多项式L3(x)代入代入得得P202 (6-23)称为称为Adams 三步隐式公式，也称三步隐式公式，也称Adams内插公式内插公式1),()()(1iixxiidxyxfxyxy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法局部截断误差：5(5)319()720Rh y 第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法三三.阿达姆斯预估校正公式阿达姆斯预估校正公式 我们常把阿达姆斯显式及隐式联立使用，即
42、构造所谓阿达姆斯预估校正公式第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法方法方法 公式显式 局部截断误差 步数Euler yi+1=yi+hf(xi,yi)O(h2)单步后退Euler yi+1=yi+hf(xi+1,yi+1)O(h2)单步梯形 O(h3)单步预报-校正 O(h3)单步四阶龙格库P198 （6-20)O(h5)单步Adams外插外插 O(h5)多步（4)Adams内插内插 O(h5)多步（3)几种方法的比较几种方法的比较),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy第第6 6章章 常微分方程数值解法常微分方程数值解法 计算方法练习3：解解 的单步法为：的单步法为：1、求其局部截断误差主项，并回答是几阶的、求其局部截断误差主项，并回答是几阶的2、证明该方法是无条件稳定的、证明该方法是无条件稳定的 00(,)()dyf x ydxy xy),(2),(3111iiiiiiyxfyxfhyy
