第七章 常微分方程数值解法.

1、2022-3-1912022-3-192一阶常微分方程初值问题是：一阶常微分方程初值问题是： (7.1) (7.2)其中f(x,y)是已知的xoy平面上某个区域D上连续函数, 式(7.1)是微分方程，有无穷多解, 式(7.2)是确定解的初始条件。 如果一元函数y(x)对一切axb 满足（1） (x,y(x)D（2） y(xn)=yn（3） y存在，且y(x)=f(x,y(x)则称y(x)是初值问题(7.1),(7.2)在a,b上的解。定义1.1 如果存在正常数L0,使得对一切xa,b及y1 ,y2 c,d 有 |f(x,y1)-f(x,y2)| L|y1- y2|则称f(x,y)关于y满足li

2、pschitz(李普希滋）条件，其中L称为lipschitz常数。 )(),(00 xyyyxfy2022-3-193关于初值问题解的存在,唯一及对初始条件的连续依赖性,有下列定理：定理7.1 设f(x,y)是在D=(x,y)| axb,yR上的连续函数，其中a,b为有限实数,且f(x,y)关于y满足lipschitz条件，则对(x0,y0) D,初值问题(7.1),(7.2)在a,b存在唯一的连续函数微分解y(x).定理7.2 设f(x,y)是在D=(x,y)| axb, yR上的连续函数，且f(x,y)关于y满足lipschitz条件,设y(x;s)是初值问题 y =f(x,y) y(x0

3、)=S 的解,则当s取s1和s2时,有估计 y(x;s1)-y(x;s2) EXP（Lx-x0s1-s2) 定理7.2的意义：若初值问题(7.1),(7.2)中，初始值有一微小扰动，则解的扰动也是微小的，也即解连续依赖初始条件.2022-3-1941. 差商方法差商方法 (化导数为差商的方法化导数为差商的方法) 设初值问题(7.1)的准确解y(x)在节点xn之值为y(xn),记y(xn)的近似值为yn,则初值问题(7.1),(7.2)化为 差分方程问题hxyxyxxxyxyxynnnnnnn)()()()()(111)(,.),2 , 1 , 0(),(001xyynyxfhyynnnn(7.

4、3)称为解一阶方程初值问题(7.1),(7.2)的欧拉(Euler)公式(欧拉方法)(7.3)2022-3-195例7.1 用尤拉方法求初值问题0)0(2yyxy的数值解(取h=0.1). 解 因为f(x,y)=x-y2,x0=0,h=0.1,故由递推公式(4.3)得0,.)2 , 1 , 0(),( 1 . 0021ynyxyynnnn(7.4)由(7.4)计算所得的数值如下表n 0 1 2 3 4 xn 0 0.1 0.2 0.3 0.4 yn 0 0 0.01000 0.02999 0.05990 2022-3-1962. 数值积分方法数值积分方法对微分方程 y =f(x,y(x) 在区

5、间x,x+k上积分得特别当x=xn及k=h时有kxxdxxyxfxykxy)(,()()()(),(001xyyyxhfyynnnn1)(,()()(1kkxxkkdxxyxfxyxy(7.5)对于(7.5),利用在区间xn,xn+1上的矩形积分公式)(,()(,(1nnxxxyxhfdxxyxfnn再用yn代替y(xn),由此得到尤拉公式2022-3-1973. 隐式尤拉法和二步尤拉法在差商方法中,若改用向后差商y(xn+1)-y(xn)/h替代方程y(xn+1)=f(xn+1,yn+1)中的导数项y(xn+1),即可导出一个新的计算方法 )(),(00111xyyyxhfyynnnn(7.

6、6)此公式与尤拉法有本质的区别,后者是关于yn+1的一个直接的计算公式,这类公式称作是显式的;然而此公式的右端含有未知的yn+1,它实际上是关于yn+1的一个函数方程,因此称作隐式尤拉法.若改用中心差商y(xn+1)-y(xn-1)/2h替代方程y(xn)=f(xn,yn)中的导数项y(xn),便可得到下列公式)(),(20011xyyyxhfyynnnn(7.7)利用公式(7.7)计算yn+1时,需要调用前面两步的信息yn和yn-1因此称为二步尤拉法.2022-3-1984. 局部截断误差定义定义7.1 设yn是初始问题(7.1),(7.2)的精确解, 即yn=f(xn),则误差 y(xn+

7、1) yn+1称为局部截断误差。定义定义7.2：称一种数值方法的精度是p阶的,如果其局部误差为O(hp+1),其中p为正整数. 对于尤拉法,若yn=f(xn), 则由(7.3)有)()()(,()(1nnnnnnxhyxyxyxhfxyy而按台劳公式有1 21),()()()(2nnhnnnxxyxhyxyxy因而有1 211),()(2nnhnnxxyyxy这说明尤拉法仅为一阶方法.同样可以证明隐式尤拉法的精度也是一阶的.2022-3-199二步尤拉法比尤拉法或隐式尤拉法具有更高的精度,由台劳展开式知设 yn=y(xn),yn-1=Y(xn-1),将上式与二步尤拉法(7.7)相减得 y (x

8、 n+1)-yn+1=O(h3)因此,二步尤拉法的精度是二阶的.11 311),()(2)()(3nnhnnnxxyxhyxyxy2022-3-1910梯形公式将方程y=f(x,y)的两端从xn到xn+1求积分得1)(,()()(1nnxxnndxxyxfxyxy利用数值积分中的梯形公式)(,()(,()(,(1121nnnnxxhxyxfxyxfdxxyxfnn并将其中的y(xn),y(xn+1)分别用yn,yn+1替代得计算公式),(),(1121nnnnhnyxfyxfyyn此公式称为梯形公式.梯形公式实际上欧拉法与隐式欧拉法的算术平局平均,其精度是二阶的. (见P92)(7.8)202

9、2-3-1911例4.2 用梯形公式求初值问题的解在x=0.01上的值y(0.01).解 取h=0.01, 从(4.8)得y1=y0+(h/2)(y0+y1),所以 1)0(,yydxdy01005. 1000025. 0010025. 1)1 (.)1)(1 (422042201112222hhhhhyyyhh另一方面,直接求解得y(x)=ex,由此有01005.11.1)01.0(!2)01.0(!101.0!3)01.0(!2)01.0(!101.001.0232 ey2022-3-19122. 改进的尤拉法尤拉法是一个显式算法,计算量小,但精度低;梯形公式虽然提高了精度,但它是一种隐式

10、算法,计算量大.综合使用这两种方法:先用尤拉法求得一个初步值,记作 ,称之为预报值;预报值 的精度不高,但可用它替代(4.8)右端的yn+1再直接计算,得到校正值yn+1.这样建立的预报-校正系统为1ny1ny),(),(),(11211nnnnhnnnnnnyxfyxfyyyxhfyy校正预报(4.9)(4.9)称为改进的尤拉法,或称作预估-校正法. 改进的尤拉法的两种等价形式:),(,(),(21nnnnnnhnnyxhfyxfyxfyy)(),(),(2111cpnpnncnnnpyyyyxhfyyyxhfyy(4.10)(4.11)(4.10)称为嵌套形式,(4.11)称为平均形式,

11、其中h为步长,n为步数,X0,y0为“老值”,即前一步的近似值, X1,y1为“新值”,即该步计算结果.2022-3-1913例4.3 用改进的尤拉法求解初值问题1)0(2yyyxdxdy在区间0,1上,步长h=0.1的数值解(用6位小数进行计算)解 对于此题,改进的尤拉法的具体形式是)(,112122121nnnnnnyxnyxnhnnyxnnnyyyyyhyy取步长h=0.1,初始值y0=1,具体计算结果如下xn ynxn yn0.1 1.0959090.2 1.1840960.3 1.2662010.4 1.3433600.5 1.4164020.6 1.4859560.7 1.5525

12、150.8 1.6164760.9 1.6781681.0 1.737869 龙格-库塔法是采用在多个点计算f(x,y)的函数值,然后对这些函数值作线性组合,构造近似公式,再把近似公式和解的台劳展开式相比较,使前面的若干项吻合,从而使近似公式达到一定精度的数值计算公式.二阶龙格-库塔方法:二阶龙格-库塔方法是以计算二次f(x,y)为基础的,其公式为),(),(12122122111kbyhaxhfkyxhfkkckcyynnnnnn其中c1,c2,a2,b21是待定参数. 其选取原则是使近似公式的阶仅可能高,也就是使局部截断误差 达到最高阶,其中表示用精确值y(xn)替代公式中的yn计算得到的

13、近似值.111)(nnnyxyT1ny(4.12)2022-3-1915c1,c2,a2,b21 的具体求法一方面y(xn+1)在xn处的台劳展开式为)()( )( )()(3!212hoxyxhyxyxynhnnn或)()()()(3!212hofffhfxyxyyxhnn其中f=f(xn,y(xn),fx (或fy)表示f关于第一(或二)个变量的偏导数在(xn,y(xn)处的值. 另一方面由于k1,k2在xn处的台劳展开式为)()()(,()(,()(31212121221hofkbhfaxyxhfkbxyhaxhfkxhykyxnnnnn于是有)()()( )()(322211221ho

14、fffhacxhyccxyyyabxnnn若要求局部截断误差达到o(h3),则要求 c1+c2=1, c2a2=1/2, b21/a2=1 (4.13)2022-3-1916显然(4.13)无穷多个解, 因此任意确定一组解,均将得到一个二阶龙格-库塔法的计算公式.例4.4 若令a=1,则从(4.13)解得c1=c2=1/2, b21=1.因而得),(),(1212211211kyxhfkyxhfkkkyynnnnnn(4.14)(4.14)为改进的尤拉法.二级R-K方法是显示单步式,每前进一步需要计算两个函数值.由上面的讨论可知,适当选择四个参数c1,c2,a2,b21,可使每步计算两次函数值

15、的二阶R-K方法达到二阶精度.能否在计算函数值次数不变的情况下,通过选择四个参数,使得二阶R-K方法的精度再提高呢? 答案是否定的.无论四个参数怎样选择,都不能使公式(4.12)提高到三阶.这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶.若要获得更高阶得数值方法若要获得更高阶得数值方法, ,就必须增就必须增加计算函数加计算函数2022-3-1917高阶龙格-库塔方法:),(),(),(232131331212213322111kbkbyhaxhfkkbyhaxhfkyxhfkkckckcyynnnnnnnn其中c1,c2,c3,a2,a3,b21,b32是待定参数.(a) 三阶龙格-库塔方法(b) 标准四阶龙格-库塔公式(也称经典公式),(),(),(),()22(3422121312121214321611kyhxhfkkyhxhfkkyhxh