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

第七章常微分方程初值问题的数值解法1/13/20231一阶常微分方程初值问题是：

(7.1)(7.2)其中f(x,y)是已知的xoy平面上某个区域D上连续函数,式(7.1)是微分方程，有无穷多解,式(7.2)是确定解的初始条件。如果一元函数y(x)对一切a≤x≤b满足（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

如果存在正常数L>0,使得对一切x[a,b]及y1,y2[c,d]有

|f(x,y1)-f(x,y2)|L|y1-y2|则称f(x,y)关于y满足lipschitz(李普希滋）条件，其中L称为lipschitz常数。

1/13/20232关于初值问题解的存在,唯一及对初始条件的连续依赖性,有下列定理：定理7.1

设f(x,y)是在D={(x,y)|

a≤x≤b,y∈R}上的连续函数，其中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)|a≤x≤b,y∈R}上的连续函数，且f(x,y)关于y满足lipschitz条件,设y(x;s)是初值问题

y

=f(x,y)y(x0)=S的解,则当s取s1和s2时,有估计｜y(x;s1)-y(x;s2)｜≤EXP（L｜x-x0｜｜s1-s2｜)

定理7.2的意义：若初值问题(7.1),(7.2)中，初始值有一微小扰动，则解的扰动也是微小的，也即解连续依赖初始条件.1/13/202331.差商方法(化导数为差商的方法)

设初值问题(7.1)的准确解y(x)在节点xn之值为y(xn),记y(xn)的近似值为yn,则初值问题(7.1),(7.2)化为

差分方程问题(7.3)称为解一阶方程初值问题(7.1),(7.2)的欧拉(Euler)公式(欧拉方法)7.1欧拉方法(7.3)1/13/20234例7.1

用尤拉方法求初值问题的数值解(取h=0.1).

解因为f(x,y)=x-y2,x0=0,h=0.1,故由递推公式(4.3)得(7.4)由(7.4)计算所得的数值如下表n01234…xn00.10.20.30.4…yn000.010000.029990.05990…1/13/202352.数值积分方法对微分方程y´=f(x,y(x))在区间[x,x+k]上积分得特别当x=xn及k=h时有(7.5)对于(7.5),利用在区间[xn,xn+1]上的矩形积分公式再用yn代替y(xn),由此得到尤拉公式1/13/202363.隐式尤拉法和二步尤拉法在差商方法中,若改用向后差商[y(xn+1)-y(xn)]/h替代方程y’(xn+1)=f(xn+1,yn+1)中的导数项y’(xn+1),即可导出一个新的计算方法

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

y(xn+1)–yn+1称为局部截断误差。定义7.2：称一种数值方法的精度是p阶的,如果其局部误差为O(hp+1),其中p为正整数.

对于尤拉法,若yn=f(xn),则由(7.3)有而按台劳公式有因而有这说明尤拉法仅为一阶方法.同样可以证明隐式尤拉法的精度也是一阶的.1/13/20238二步尤拉法比尤拉法或隐式尤拉法具有更高的精度,由台劳展开式知设yn=y(xn),yn-1=Y(xn-1),将上式与二步尤拉法(7.7)相减得

y(xn+1)-yn+1=O(h3)因此,二步尤拉法的精度是二阶的. 1/13/202397.2改进的欧拉方法梯形公式将方程y’=f(x,y)的两端从xn到xn+1求积分得利用数值积分中的梯形公式并将其中的y(xn),y(xn+1)分别用yn,yn+1替代得计算公式此公式称为梯形公式.梯形公式实际上欧拉法与隐式欧拉法的算术平局平均,其精度是二阶的.(见P92)(7.8)1/13/202310例4.2用梯形公式求初值问题的解在x=0.01上的值y(0.01).解取h=0.01,从(4.8)得y1=y0+(h/2)(y0+y1),所以

另一方面,直接求解得y(x)=ex,由此有1/13/2023112.改进的尤拉法尤拉法是一个显式算法,计算量小,但精度低;梯形公式虽然提高了精度,但它是一种隐式算法,计算量大.综合使用这两种方法:先用尤拉法求得一个初步值,记作,称之为预报值;预报值的精度不高,但可用它替代(4.8)右端的yn+1再直接计算,得到校正值yn+1.这样建立的预报-校正系统为(4.9)(4.9)称为改进的尤拉法,或称作预估-校正法.

改进的尤拉法的两种等价形式:(4.10)(4.11)(4.10)称为嵌套形式,(4.11)称为平均形式,其中h为步长,n为步数,X0,y0为“老值”,即前一步的近似值,X1,y1为“新值”,即该步计算结果.1/13/202312例4.3

用改进的尤拉法求解初值问题在区间[0,1]上,步长h=0.1的数值解(用6位小数进行计算)解对于此题,改进的尤拉法的具体形式是取步长h=0.1,初始值y0=1,具体计算结果如下xn

ynxn

yn0.11.0959090.21.1840960.31.2662010.41.3433600.51.4164020.61.4859560.71.5525150.81.6164760.91.6781681.01.7378691/13/202313

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

的具体求法一方面y(xn+1)在xn处的台劳展开式为或其中f=f(xn,y(xn)),fx

(或fy)表示f关于第一(或二)个变量的偏导数在(xn,y(xn))处的值.另一方面由于k1,k2在xn处的台劳展开式为于是有若要求局部截断误差达到o(h3),则要求

c1+c2=1,c2a2=1/2,b21/a2=1(4.13)1/13/202315显然(4.13)无穷多个解,因此任意确定一组解,均将得到一个二阶龙格-库塔法的计算公式.例4.4若令a=1,则从(4.13)解得c1=c2=1/2,b21=1.因而得(4.14)(4.14)为改进的尤拉法.二级R-K方法是显示单步式,每前进一步需要计算两个函数值.由上面的讨论可知,适当选择四个参数c1,c2,a2,b21,可使每步计算两次函数值的二阶R-K方法达到二阶精度.能否在计算函数值次数不变的情况下,通过选择四个参数,使得二阶R-K方法的精度再提高呢?答案是否定的.无论四个参数怎样选择,都不能使公式(4.12)提高到三阶.这说明每一步计算两个函数值的二阶R-K方法最高阶为二阶