数值分析第九章常微分方程数值解法

1第九章常微分方程数值解法§1Euler方法§2改进的Euler方法§3Runge-Kutta法§4线性多步法2近似求解3一阶常微分方程初值问题为保证上述问题解的存在唯一性,对函数f

做如下假设:f关于y的偏导数有界,即存在常数L>0,使得本章中的L均指此数.常微分方程解的存在唯一性4数值解法:

就是寻求解y(x)在一系列离散点上的近似值相邻两个节点的间距称为步长.本章均假设此时节点常微分方程数值解法5

初值问题的各种差分方法都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.只要给出从已知信息计算的递推公式，此类计算格式统称为差分格式.常微分方程差分格式6数值求解一阶常微分方程初值问题难点:如何离散y

?常见离散方法差商近似导数数值积分方法

Taylor展开方法7结果有设用y(xn)的近似值yn代入上式右端，记所求结果为yn+1，这样导出的计算公式将方程y=f(x,y)中的y在xn点用向前差商近似差商近似导数(Euler公式)Euler公式,单步法,显式公式8Euler公式计算总结设的精确解为将区间[a,b]N等分,分点xn=x0+nh(n=1,2,…,N)步长用Euler公式求得yn(n=1,2,…,N.)9结果有由此得差分格式计算公式将方程y=f(x,y)中的y在xn+1点用向后差商近似向后Euler公式,单步法,隐式公式差商近似导数(向后Euler公式)10设将方程y=f(x,y)的两端从xn

到xn+1

求积分,得用不同的数值积分方法近似上式右端积分,可以得到计算y(xn+1)的不同的差分格式.用左矩形公式计算积分项再离散化，得Euler公式数值积分方法(Euler公式)11用右矩形公式计算积分项再离散化，得向后Euler公式数值积分方法(向后Euler公式)12用梯形法计算积分项再离散化，即可得如下计算公式梯形公式,单步法,隐式公式数值积分方法(梯形公式)13Taylor展开方法(Euler公式)由此亦可得Euler公式将y(xn+1)在xn

点作Taylor展开,得14§1Euler方法用差分格式来近似微分方程初值问题15例用Euler公式求解初值问题解取步长h=0.1,计算结果见下表计算公式为16y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321xn0.10.20.30.40.50.60.70.80.91.0yn1.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848两者比较可以看出Euler法的精度较差.17hhhStraightlineapproximationEuler法的几何意义x0x1x2x3y018向后Euler法这是一个隐式公式.已知yn,上式是关于yn+1的非线性方程,不能从中直接求出yn+1,需要用非线性方程求根的方法来求解,通常用迭代法.19已知yn求yn+1的迭代法则当hL<1时,迭代函数向后Euler法计算量大,但数值稳定性好!20Euler法的局部截断误差局部截断误差:用Euler法计算一步所产生的误差,即假设第n步的计算是精确的yn=y(xn),y(xn+1)与yn+1的差,即这种误差称为局部截断误差,简称截断误差.21当yn=y(xn)时,局部截断误差主项Euler法的局部截断误差22向后Euler法的局部截断误差定义其局部截断误差为向后Euler法的计算公式23向后Euler法的局部截断误差局部截断误差主项24在不考虑舍入误差的情况下,称y(xn+1)与yn+1之差为整体截断误差,记为令则Euler法的整体截断误差25Euler法的整体截断误差26注意e0=0,故27所以这表明,Euler方法的整体截断误差与h同阶.当h0时,eN

0.如果某种数值方法的局部截断误差为O(hp+1),则称它的精度是p阶的,或称之为

p

阶方法.由此可知Euler方法仅为一阶方法.28梯形公式§2改进的Euler法从直观上看,用梯形公式计算数值积分要比矩形公式精度高.相应地求解常微分方程的梯形公式的精度也要比Euler公式精度高.29梯形公式的局部截断误差局部截断误差主项故梯形公式是二阶方法30则当hL<2时,梯形公式是一个隐式公式.已知yn,用迭代法计算yn+1用梯形法计算,从yn到yn+1要进行迭代计算,计算量较大.实际计算常把Euler法和梯形法结合起来,这就是改进的Euler法.31先用Euler法求得一个初步的近似值，记为，称之为预测值，然后用它替代梯形法右端的yn+1

再直接计算fn+1,得到校正值

yn+1,这样建立的预测－校正系统称为改进的Euler法：预测校正改进的Euler法32实践表明,改进Euler法的精度明显高于Euler法.利用Taylor展开可以证明,改进Euler的局部截断误差为O(h3),故它是二阶方法.它有下列平均化形式：33例用改进Euler法求解解取步长h=0.1,计算结果见下表改进Euler法公式为34y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321xn0.10.20.30.40.50.60.70.80.91.0yn1.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379同Euler法的计算结果比较,改进Euler法明显改善了精度.35改进Euler法的几何意义(

K1

+

K2)/2K2K1xnxn+136考察差商根据微分中值定理，存在点，利用所给方程y=f(x,y)

得称K*=f(,y())为区间[xn,xn+1]上的平均斜率,这样只要对平均斜率K*提供一种算法，相应地便导出一种计算格式.§3Runge-Kutta法Runge-Kutta法的基本思想37梯形法：改进Euler法：其中Euler法：向后Euler法：38Runge－Kutta法的基本思想就是设法在[xn,xn+1]内多预报几个点的斜率值，然后把它们加权平均作为平均斜率，以期望构造出更高精度的计算格式.Runge-Kutta法的基本思想39考察区间[xn,xn+1]内一点xn+p=xn+ph,0<p1,用两个点xn,xn+p

的斜率K1,K2的加权平均代替平均斜率K*,得如下计算格式：其中有两个待定参数,p,适当选取它们的值，使上述格式有尽可能高的精度.二阶Runge-Kutta方法40K1K2xnxn+phxn+1

=xn+h(1)

K1

+

K2加权平均斜率41若p=1/2,则该格式是二阶的,故统称满足这一条件的一族公式为二阶Runge-Kutta方法.特别地,当p=1,=1/2

时,上述格式即为改进的Euler法，如果取p=1/2,=1

则上述格式称为中点法.42K1K2xnxn+1/2xn+1中点公式

K243由此可以构造三阶Runge-Kutta方法.为了进一步提高精度，可以考虑用三个点xn,xn+p,xn+q

的斜率值K1,K2,K3

加权平均得出平均斜率K*的近似值,其中三阶Runge-Kutta方法44经典的三阶Runge-Kutta方法45若用四个点的斜率值的加权平均来近似平均斜率,可以导出四阶Runge-Kutta方法.四阶Runge-Kutta方法46经典的四阶Runge-Kutta方法47xnxn+h/2xn+hK1K2K3K4经典的四阶Runge-Kutta方法48例设取步长h=0.2用四阶Runge-Kutta方法求解解49y(xn)1.18321.34161.48321.61251.7321xn0.20.40.60.81.0yn1.18321.34161.48321.61241.7320同Euler法和改进Euler法的计算结果比较,显然以四阶RK法的精度为高.四阶RK法每1步计算4次函数值,改进Euler法每1步只要计算2次函数值,但由于这里步长放大了2倍,两者所用计算量几乎相同.取步长h=0.2,计算结果见右表50Runge-Kutta方法的推导基于Taylor展开法,因而它要求解具有较好的光滑性.如果解的光滑性差,那么,使用四阶Runge-Kutta方法求得的数值解,其精度反而不如改进的Euler法.实际计算时,应当针对问题的具体特点选择合适的算法.51在逐步推进的求解过程中,计算yn+1之前事实上已经求出了一系列的近似值y0,y1,…,yn,如果充分利用前面多步的信息来预测yn+1,则可以期望会获得较高的精度.这就是线性多步法的基本思想.§4线性多步法52线性二步法的一般格式局部截断误差53一般的线性多步法公式可表示为其中i,i均为常数,fk=f(xk,yk)(k=n+1,n,…,nr).若1=0,则公式为显式,若10,则公式为隐式.定义:线性多步法在xn+1上的局部截断误差定义为54基于数值积分方法非常直观,能用它导出大部分的常用线性多步法,包括:Simpson公式

Adams显隐公式

Milne公式构造线性多步法的主要途径是基于数值积分方法和基于Taylor展开方法,前者可直接由方程两端积分后利用插值求积公式得到.55Simpson公式其中将y=f(x,y)的两端从xn1

到xn+1

求积分，得Simpson公式,两步法,隐式公式离散得其中F(x)=f(x,y(x)),用Simpson公式计算积分项56Simpson公式的截断误差故Simpson公式是隐式二步四阶方法57Adams显式公式右端积分中用过节点xn,xn1,

xn2,…,

xnr的F(x)的r次插值多项式Lr(x)代替F(x)求积分,再离散化得(r+1)阶Adams显式公式.一阶Adams显式公式就是Euler公式58右端积分中用过节点xn+1,xn,

xn1,…,

xnr+1的F(x)的r次插值多项式Lr(x)代替F(x)求积分,再离散化得(r+1)阶Adams隐式公式.Adams隐式公式一阶Adams隐式公式就是向后Euler公式59二阶Adams显式公式F(x)用过节点

xn,

xn1的线性插值函数L1(x)代替60再离散化，即可得如下二阶Adams显式公式截断误差61二阶Adams隐式公式:梯形公式再离散化，即可得如下计算公式F(x)用过节点

xn,

xn+1的线性插值函数代替62四阶Adams显式公式截断误差63四阶Adams隐式公式截断误差64Milne公式F(x)用过节点

xn,

xn1,

xn2的2次插值函数L2(x)代替再离散化，得Milne公式65Milne公式是四阶四步显式方法,截断误差为66基于Taylor展开方法能导出所有的常用线性多步法,下面以线性二步法为例来说明这种方法线性二步法的一般格式适当选取参数0,1,1,0,1,

使其局部截断误差达到最高阶(使公式的阶尽可能高)局部截断误差:假设第n步以前的计算是精确的,yi=y(xi)(in),用多步法计算一步所产生的误差,即y(xn+1)与yn+1的差.67为符号简单记,把简记为由Taylor展开式得局部截断误差6869要使局部截断误差的阶高,Rn+1的前面几项应尽可能多为0,因只有5个未知数,可令前5项为0.70其中fi=f(xi,yi),i=n

1,n,n+1.Simpson公式两步四阶隐式方法局部截断误差71

答案:该公式为(r+1)阶Adams显式公式适当选取0,1,…,r

使下列公式的阶尽可能高适当选取0,1使公式的阶尽可能高.例

r=1的情形局部截断误差72得二阶Adams公式要使阶尽可能高,须令局部截断误差73

答案:该公式为(r+1)阶Adams隐式公式.适当选取1,0,1,…,r1

使下列公式的阶尽可能高74

Hamming公式四阶三步隐式公式适当选取0,2,1,0,1,

使下列公式的阶尽可能高

答案:所求公式为截断误差75例分别用四阶Adams显式和隐式公式求解初值问题的数值解,取h=0.1.利用线性多步法求解初值问题,必须先用其他方法算出开始几个点的近似值.一般可用同阶的单步法.从以上例子看到同阶的Adams方法,隐式方法要比显式方法误差小,这可以从两种方法的局部截断误差主项的系数大小来加以解释.76一般地,同阶的隐式法比显式法精确,而且数值稳定性好.但在隐式公式中,通常很难解出yn+1,需用迭代法求解,这样又增加了计算量.因此实际计算时,很少单独使用显式公式或隐式公式,而是将它们联合使用:先用显式公式求出y(xn+1)的预测值,再用隐式公式对预测值进行校正.77仿照改进的Euler格式的构造方法，将显式和隐式两种Adams格式相匹配，可构成下列Adams预测－校正系统：预测校正Adams预测-校正系统78预测校正Milne-Hamming预测-校正系统以上两种预测-校正公式均为四阶公式,其起步值通常用四阶RK公式计算.79§5单步法的相容性、收敛性与稳定性单步法的一般形式当含有yn+1时,方法是隐式的,若不含yn+1则为显式方法.称(x,y,h)为增量函数显式单步法的一般形式为80例改进Euler法的增量函数(x,y,h)=f(x,y)例

Euler法的增量函数单步法的局部截断误差为p(1)阶精度的单步法是指81数值解法的基本思想是,通过某种离散化手段将微分方程转化为差分方程.转化是否合理?即(1)当h0时,差分方程是否能无限逼近微分方程?(相容性问题)(2)差分方程的解yn是否趋于微分方程的准确解y(xn)?(收敛性问题)82单步法可以看成是由下述方程近似而来即固定xn

,

当h0时上式左边极限为y(xn),如果近似是合理的,右边应有83定义若增量函数满足(x,y,0)=f(x,y)则称单步法与常微分方程初值问题是相容的.本章介绍的所有显式单步法(RK方法)均是相容的84定义若一种数值方法对于任意固定的xn=x0+nh,当h0(同时n

）时数值解yn趋于准确解y(xn).则称该方法是收敛的.判断一种数值方法的收敛性问题本质上就是估计该方法整体截断误差的问题.例

Euler方法的整体截断误差为当h0时,en0,故Euler法是收敛的.收敛性85定理设单步法yn+1=yn+h(xn,yn,h)

具有p(1)阶精度,其增量函数(x,y,h)关于y满足Lipschitz条件又设初值是准确的,即y0=y(x0),则其整体截断误差为y(xn)yn=O(hp)该定理意味着,满足定理条件的单步法是收敛的只要微分方程y=f(x,y)中的函数f关于y满足Lipschitz条件,则本章介绍的所有显式单步法(RK方法)均满足定理条件,故所有RK方法均收敛.86关于收敛性的讨论有个前提，即必须假定差分方法的每一步计算都是准确的(没有舍入误差).然而实际计算中往往由于有舍入误差等原因而产生扰动，而这些扰动有可能“淹没”真解,所以还要考虑稳定性问题.稳定性问题87稳定性比较复杂.为简化讨论,仅考察下列模型方程定义若一种数值方法在节点值yn上大小为的扰动,引起以后各节点值ym(m>n)上产生的偏差值均不超过,则称该方法是绝对稳定的.用某种数值方法求解上述模型方程，对固定的h,使得该方法绝对稳定的h和的全体称为该方法的绝对稳定区域.88Euler方法的绝对稳定区域模型问题计算格式设实际计算yn时产生误差n,它使节点值yn+1产生误差n+1,则从而为使误差不扩大,需89向后Euler方法的绝对稳定区域模型问题计算格式设实际计算yn时产生误差n,它使节点值yn+1产生误差n+1,则向后Euler方法的绝对稳定区域90梯形公式的绝对稳定区域模型问题计算格式91类似于Euler方法的讨论,得梯形法的绝对稳定区域为即梯形公式的绝对稳定区域92RK方法的绝对稳定区域二阶RK方法的绝对稳定区域对于模型问题改进Euler法93当是实数时,2.51h0当是实数时,2.78h0二阶RK