常微分方程数值解法缩

1、 第六章常微分方程数值解法n本章研究常微分方程初值问题 (6.1.1) 的数值解法. 并且假定f(x, y)满足解的存在唯一性定理及相当光滑等条件.n初值问题(6.1.1)的精确解记为y(x). d( , ),d( )yf x yaxbxy a建立数值解法的基本思想n 本章的数值解法，它不是求(6.1.1)的解y(x)的解析表达式或近似表达式，而是通过某种离散化方法，将连续变量的初值问题转化为关于离散量的差分方程的初值问题来求一系列离散点上的解值y(xi) 的近似值yi. 利用计算机解微分方程主要使用数值方法. n取一系列点 x0, x1, , xn , y(x0)=y0 , y(x1)y1

2、, , y(xn)yn , y0 , y1 , , yn , 称为数值解用离散化方法建立求y(xn)的近似值yn的递推格式(差分方程)，求得yn h= xn xn-1称为步长本章都取定步长.n解初值问题(6.1.1)的数值解法，其特点都是采取步进式的方法，即求解过程顺着节点排列的次序一步一步向前推进.n这种数值解法分为两大类：(1)单步法：若求yn+1 ，只需利用它前一步的 信息yn ，则称这种方法为单步法。它由 y0出发，可求得y1 ，y2， y3 (2)多步法：若求yn+1 ，需利用它前面至少两个点的息，则称这种方法为多步法.数值解法研究的主要问题：(1)方法的推导：采用的离散化手段，精度

3、准则.(2)收敛性：差分方程的解是否充分逼近初值问题的解.(3)稳定性：初始数据、计算过程中每步产生的误差对以后各步解的影响，这种误差传播是否在衰减.n具体的数值方法还应考虑(1)误差估计(2)解的起动方法(3)步长如何选取(4)隐式方法的如何计算建立数值解法的基本途径n常用的离散化方法 (1)Taylor展开 (2)化导数为差商 (3)数值积分6.1 单步法及基本概念Euler折线法n利用Taylor展开法将y(xn+1)在xn处Taylor展开得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,）称此方法为Euler折线法或矩形法1()()(, ()nnnny xy xhf

4、xy x21()()()()()2!nnnnnhy xy xhy xhy xyn 利用化导数为差商的方法得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,） (6.1.2)(,()()()(1nnnnnxyxfxyhxyxy1()()(, ()nnnny xy xhf xy xn利用数值积分的方法在xn,xn+1上对y(x)=f(x, y (x)积分得用左矩形求积公式计算定积分有 y(xn+1)y(xn)+hf(xn,y(xn)以此得差分方程 yn+1=yn+hf(xn,yn) （n=0,1,2,）1d)(,()()(1nnxxnnxxyxfxyxy梯形法n用不同的近似公式计算

5、定积分的值,就得到解初值问题的不同数值解法.n用梯形求积公式计算积分得 (6.1.3) 这个方法称为梯形法. 它是隐格式.n运用它常采用下面的迭代格式111 (,)(,)(0,1,2,)2nnnnnnhyyf xyf xyn(0)1(1)( )111(,) (,)(,)2(0,1,2,;0,1,2,)nnnnkknnnnnnyyhf xyhyyf xyf xykn改进的Euler方法n若梯形法只迭代一次,便得改进的Euler方法 (6.1.4)也可写成形式 (6.1.4)11(1)(0)11(,)(0,1,2,) (,)(,)2nnnnnnnnnnyyhf xynhyyf xyf xy1121

6、211122(,)(,)nnnnnnyyhkhkkf xykf xh yhk数值方法精度的衡量准则n定义定义 设yn=y(xn)，则称Tn=y(xn+1)-yn+1为方法的从xn到xn+1这一步的局部截断误差. n定义若差分方程对所有y(x)Mr都精确成立，而对于某个r+1次多项式不能精确成立，则称这个数值方法是r阶的. 等价定义等价定义 若数值方法的局部截断误差为O(hr+1)，则称这种方法为r阶方法，这里r为非负整数.n方法的阶数越高，逼近效果越好.n由Taylor展开式知 (6.1.5)对于Euler法Tn=y(xn+1)- y(xn)-hf(xn,y(xn)=y(xn+1)-y(xn)

7、-hy(xn) Euler法是1阶方法 )(! 2)( )()()(21nnnnnxyhxhyxyhxyxy2()2!nhyn对于梯形法由Taylor展开式知 (6.1.6)(6.1.5)与(6.1.6)相减并利用y(xn+1)=y (xn)+O(h)梯形法法是2阶方法.11()()()()(0,1,2,)2nnnnnhTy xy xy xy xn2111()()()()()2!nnnnnhy xy xhy xhy xy x311()()()()()2nnnnnhTy xy xy xy xO hn对于改进Euler法，设yn=y(xn)，利用Taylor展开式nK1=f(xn,y(xn)=y(

8、xn) 211(, ()(, ()(, ()(, ()(, ()(, ()()(, ()()()nnnnnnnnnnnnnnnnnKf xh y xhkf xy xhf xy xhkf xy xxyf xy xhf xy xy xf xy xxyy xhy x将其代入中Tn有改进Euler法是2阶方法.31121()()()()2nnnTy xy xkkO hRunge-Kutta方法nTaylor展开式n理论上讲，只要解y(x)充分光滑，通过保留Taylor展开式的若干项就可得到任意阶的近似公式，但计算y(x)的各阶导数很麻烦。可间接利用这种思想 )(,()(,()(,(! 2)(,()()

9、(! 2)( )()()(221nnnnnnnnnnnnnnxyxfyxyxfxyxfxhxyxhfxyxyhxhyxyhxyxynEuler法也可写成形式n其局部截断误差为O(h2) ，是一阶方法.每步计算f的值一次111(,)nnnnyyk hkf xynEuler预报-校正公式也可写成形式n局部截断误差为O(h3)，是二阶方法每步计算f的值二次n可以通过增加计算f的值的次数，提高公式的阶数（精度） 1121211122(,)(,)nnnnnnyyhkhkkf xykf xh yk hn以f在不同点上的函数值的线性组合来代替yn+1 yn，其中有一些可待定选取的待定参数，通过Taylor展

10、开确定这些待定参数使建立的数值方法按要求达到一定的阶数，这种思想就是Runge-Kutta方法的思想nRunge-Kutta法的一般形式其中Ri,ai,bij都是常数，a0=0，b1j=0， ( j=1,2,i-1)1111(,),1,2, .snniiiiininijjjyyhR Kkf xa h yhb Kisn二阶Runge-Kutta法的一般形式 其中R1, R2, a, b为待定常数 其局部截断误差为O(h3)，是二阶方法每步计算f的值二次11122121(,)(,)nnnnnnyyRhkR hkKf xyKf xah ybhkn设yn=y(xn) K1=f(xn,yn)=f(xn,

11、y(xn)把K2中f在(xn,y(xn)处Taylor展开21212(, ()(, ()(, ()(, ()()(, ()(, ()(, ()(, ()()nnnnnnnnnnnnnnnnKf xah y xbhKf xy xahf xy xbhKf xy xO hxyf xy xh af xy xbf xy xf xy xO hxyn再将K1,K2代入yn+1中，n将其与y(xn+1)泰勒展开式比较，要使y(xn+1)-yn+1=O(h3)，含h0, h1, h2的项相同即有11122122322()() (, ()(, ()(, ()(, ()()nnnnnnnnnnnyyRhKR hKy

12、 xh RRf xy xh aRf xy xbR f xy xf xy xO hxy个未知数，个方程满足条件的解不止一组取R1=R2=1/2,a=b=1,就是改进Euler法；取R1=0,R2=1,a=1/2,b=1/2的中点方法.122211212RRaRbRn四阶Runge-Kutta法的一般形式其中有13个待定常数局部截断误差为O(h)，是四阶方法每步计算f的值四次1112233441222113331132244411422433(,)(,)(,)(,)nnnnnnnnnnyyRhKR hKR hKR hKKf xyKf xa h yb hkKf xa h yb hkb hkKf xa

13、 h yb hkb hkb hkn设yn=y(xn)n把K2,K3,K4中f在(xn,y(xn)处泰勒展开后，将K1,K2,K3,K4代入yn+1中，再将yn+1按h的幂重新整理后与y(xn+1)泰勒展开式比较，要使y(xn+1)-yn+1=O(h5)，含h0, h1, h2 , h3, h4的项相同从而确定各参数n 24181)(121)(61)(4131211443322433422443323243234222433222433422433224343332324243232224433224342414323132124321RbbababaRaRbaababaRRbababaRRba

14、RaRaRaRaRaRaRaRaRabbbabbabaRRRR13个未知数，11个方程满足条件的解不止一组最常用的是n标准四阶Runge-Kutta公式 (6.1.7)112341122343(22)6(,)(,)22(,)22(,)nnnnnnnnnnhyyKKKKKf xyhKhKf xyhKhKf xyKf xh yhKn单步法的一般形式为 yn+1=yn+h(xn,yn,h) (6.1.8)其中多元函数依赖于f.单步法的收敛性nen= y(xn) yn称为整体截断误差. 收敛性就是讨论当x=xn固定且h=(xn-x0)/n0时 en 0的问题.n定义定义 如果对f(x, y)满足解的存

15、在唯一性条件的初值问题(6.1.1)，差分方程 (6.1.8)的解对每个确定xa,b满足则称单步法(6.1.8)是收敛的.0lim( ),nhnyy xnhxan定义定义 若增量函数(x,y,h)使 (x,y,0)=f(x,y)成立，则称单步法(6.1.8)与(6.1.1) 相容.n满足相容条件的方法至少是1阶的.n定理定理 设增量函数(x,y,h)在区域axb, -y, 0hh0中连续，并且对变量y满足Lipschitz条件, 即(x,y1,h)-(x,y2,h)Ly1-y2在这个前提下单步法(6.1.8)收敛的充分必要条件是相容性条件成立.单步法的绝对稳定性n当步长取定后，计算中的误差随着

16、步数的增加会不会积累到超出我们许可的范围，这就是稳定性问题.n单步法(6.1.8)应用于模型方程 y=y (0) 设得到的解为yn+1=R(h)yn当0时,实验方程的精确解y(x)=e(x-a)按模递减的,这要求 (6.1.8)的解yn是递减的，误差也是递减的，即要求满足R(h) 1n定义定义 单步法(6.1.8)应用于模型方程 y=y(0) ，若得到的解yn+1=R(h)yn满 足R(h) 1, 则称单步法(6.1.8)是绝对稳定 的.使R(h) 1成立的 =h所在区间称为 绝对稳定区间.n例 标准四阶Runge-Kutta公式 由R(h) 0时,实验方程的精确解y(x)=e(x-a)按模递

17、增的; (6.4.1)的解当n时无界；误差也是无界的.这种情况下，如果误差相对于真解是小的，就说方法是相对稳定的. 这时要求 ri(h)r0(h) ，i=1,p.n当0时,实验方程的精确解y(x)=e(x-a)按模递减的; 要求(6.4.1)的解当n时递减到0，误差也是递减的.这种情况下，就需要ri(h)1, i=0,1,p来保证.这时就称方法是绝对稳定的. 方法(6.2.1)的稳定性取决于特征多项式(r; h)零点 的性质，故又称(r; h)为(6.2.1)的稳定多项式. 稳定性的定义 记=hn定义 设方法(6.2.1)是收敛的. ri() (i=0,1,p)是稳定多项式(r; )=0的根,

18、 r0()是形如r0(h)=1+h+O(h2)的根.(1)若对任意的,R有 ri()r0() ，i=1,p.且当ri()=r0() 时, ri()是单根,则称方法(6.2.1)在,上是相对稳定的， ,称为此方法的相对稳定区间. (2)若对任意的(,)R有 ri()0时,相对稳定性重要； 当0时,绝对稳定重要，才有意义.n定义 若一个方法的绝对稳定区间是 (-,0), 则称该方法是A稳定的.n(1)收敛的方法的相对稳定区间不空。因为方法稳定的前提是收敛，故满足根条件和相容条件,0,.n(2)只有0的情况讨论绝对收敛性才有意义.n(3)从误差分析的角度看，绝对稳定的方法是理想的.n(3)相对稳定区

19、间和绝对稳定区间越大越好，可使尽可能大的一类微分方程是稳定的.例 讨论梯形方法的相对和绝对稳定区间.n解 梯形方法为对于试验方程y=y ，梯形方法成为即其稳定多形式为它只有一个根，记为且r0(0)=r0=1( ; )(1)(1)22hhr hr1(1)(1)22nnhhyy11()2nnnnhyyyy11()2nnnnhyyyy012( )12hr hh 由于由于其只有一个根，所以梯形方法的相对稳定区间是(-，+) 当 0时, r0() 1 当 0时, r0()1 这样，梯形方法的绝对稳定区间是(-，0).即梯形方法是A稳定的.例 讨论Simpson方法的相对和绝对稳定区间nSimpson方法

20、为利用该格式求解试验方程y=y 得其稳定多形式为它有两个根24( ; )(1)(1)333hhr hrhr1111(4)3nnnnnhyyyyy20021133( )(0,( )1),113hhr hhr hh2121133( )113hhr hh114(1)(1)333nnnhhyhyy考察当 0 时，上式小于等于1；当 0时,当 0时,n因此，Simpson方法不存在绝对稳定区间.2021211333|( )|1111133hhhr hhh2112211333| ( )|1111133hhhr hhhn步长的选取:稳定性与步长有关，步长的选取一定要保证方法是稳定的，即h属于稳定区间.n收敛

21、性是反映递推公式本身的整体截断误差对计算结果的影响；稳定性反映某一计算步骤中出现的误差对计算结果的影。6.5 预测-校正方法线性多步法(6.2.1)当b-10时，是隐格式估计出yn+1的一个初值y(0)n+1，利用迭代y(j)n+1是yn+1的第j次近似值.当hb-11时，此迭代收敛. 1110()pnnin iin iiyb hfa yhb f(1)( )11110(,)()pjjnnnin iin iiyb hf xya yhb fn用显式方法来作预测值y(0)n+1 ，用隐式方法迭代校正一次得y(1)n+1值，这种显式与隐式联合使用构成的方法称为预报-校正法.作为预报的显式公式称为预测式

22、，用于校正的隐式公式称为校正式.n例如，改进的Euler方法.11(1)(0)11(,)(0,1,2,) (,)(,)2nnnnnnnnnnyyhf xynhyyf xyf xyn预测公式与校正公式选取同阶方法，可使截断误差用预测值和校正值表示，再用截断误差来修正，得到提高精度的方法.n4阶显式Adams方法 局部截断误差为n4阶隐式Adams方法 局部截断误差为 1112919524nnnnnnhyyffff1123555937924nnnnnnhyyffff5(5)11251()()720nnnnTy xyh yx5(5)1119()()720nnnnTy xyh yx n4阶Adams预报-校正方法称为4阶Adams预报-校正格式(PECE模式).1(0)123(0)(0)111(1)(0)1112(1)(1)111:555937924:(,):919524:(,)nnnnnnnnnnnnnnnnnnhP yyffffEff xyhC yyffffEff xy(0)5(5)11(1)5(5)11251()()72019()()720nnnnnny xyh yxy xyh yx (0)(0