补充3：常微分方程数值解法

1数值计算方法第七章常微分方程数值解法2在工程和科学计算中，所建立的各种常微分方程的初值或边值问题，除很少几类的特殊方程能给出解析解，绝大多数的方程是很难甚至不可能给出解析解的，其主要原因在于积分工具的局限性。因此，人们转向用数值方法去解常微分方程，并获得相当大的成功，讨论和研究常微分方程的数值解法是有重要意义的。37.0基本概念1.一阶常微分方程的初值问题(7.0-1)注：若f在D={axb,|y|+}内连续，且满足Lip条件：L0，使|f(x，y1)–f(x，y2)|L|y1–y2|(7.0-2)则(7.0-1)的连续可微解y(x)在[a，b]上唯一存在。0)(),()),(,()(yaybaxxyxfxy42.初值问题的数值解称(7.0-1)的解y(x)在节点xi处的近似值yiy(xi)ax1x2...xn=b.为其数值解，方法称为数值方法。注：①考虑等距节点:xi=a+ih，h=(b–a)/n.②从初始条件y(a)=y0出发，依次逐个计算y1，y2，…，yn的值，称为步进法。两种：单步法、多步法。5③二阶常微分方程y''(x)=f(x，y(x)，y'(x))可设为一阶常微分方程组的初值问题：引进新的未知函数z(x)=y'(x)，则其初始条件为：称为一阶微分方程组的初值问题，方法类似。④边界问题，常用差分方法解。))(),(,()(')()('xzxyxfxzxzxy00')()(yazyay67.1初值问题数值解法的构造及其精度7.1.1构造方法对于(7.0-1)可借助Taylor展开（导数法）、差商法、积分法实现离散化来构造差分公式：1.Taylor展开：设yC[a，b]将y(xi+1)=y(xi+h)在xi处展开[xi，xi+1]y(xi+1)yi+hf(xi，yi)其中yiy(xi).称yi+1=yi+hf(xi，yi).i=0，1，2，...，n–1(7.1-1)为Euler求解公式，（Euler法）)(2)()()(21yhxyhxyxyiii)(2))(,()(2yhxyxhfxyiii72.用差商来表示：得差分方程：yi+1=yi+hf(xi，yi).即为Euler公式。若记yi+1=yi+hf(xi+1，yi+1).(7.1-2)称为向后Euler法。注：①Euler法为显式，向后Euler法为隐式——须解出yi+1.②可用迭代法yi+1(k+1)=yi+hf(xi+1，yi+1(k))k=0，1，2，…解得yi+1，其中yi+1(0)=yi+hf(xi，yi).).()()(1iiixyhxyxy),(1iiiiyxfhyy))(,()()()(1111iiiiixyxfxyhxyxy83.积分法:对(7.0-1)两边取积分得(7.1-3)取不同的数值积分可得不同的求解公式，如：①用矩形公式：左矩形y(xi+1)y(xi)+hf(xi，y(xi))Euler公式右矩形y(xi+1)y(xi)+hf(xi+1，y(xi+1))向后Euler公式1))(,()()(1iixxiidxxyxfxyxy))(,())(,())(,(111iixxiixyxhfxyxhfdxxyxfii9②用梯形公式：(7.1-4)称(7.1-4)为梯形公式隐式公式。显化：预估值：校正值：(7.1-5)称(7.1-5)为改进的Euler公式1311(,())[(,())(,())]'''(,())212iixiiiixhhfxyxdxfxyxfxyxfy))](,())(,([2)()(111iiiiiixyxfxyxfhxyxy)],(),([2111iiiiiiyxfyxfhyy)],(),([2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy104.几何意义Euler法折线法改进Euler法平均斜率折线法117.1.2截断误差与代数精度定义7.1-1:①称i=y(xi)–yi为数值解yi的（整体）截断误差。②若yk=y(xk)，k=0，1，2，…，i–1.由求解公式得数值解，则称为yi的局部截断误差。注：局部截断误差ei是指单步计算产生的误差，而（整体）截断误差i则考虑到每步误差对下一步的影响。)(~iixyyiiiyxye~)(12定义7.1-2若求解公式的（整体）截断误差为O(hp)，则称该方法是p阶方法，或是p阶精度。p阶公式定理7.1-1设数值解公式：yi+1=yi+h(xi，yi，h)中的函数(x，y，h)关于y满足Lipschitz条件：且其局部截断误差为hp+1阶，则其(整体)截断误差为hp阶，即该数值解公式为p阶方式。注：①局部截断误差较易估计,定理7.1-1表明：若ei=O(hp+1)则i=O(hp).②Euler局部截断误差为,所以一阶精度。向后Euler法也是一阶精度。③梯形公式为二阶精度。|~||),~,(),,(|yyLhyxhyx)()(2221hOyhei13例1：用Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：因为xi=1+0.1i，而f(x，y)=y+(1+x)y2，故f(xi，yi)=yi+(2+0.1i)yi2于是Euler计算公式为yi+1=yi+0.1[yi+(2+0.1i)yi2]，i=0，1，2，3，41)1(5.11),()1()()('2yxxyxxyxyxxy1)(14计算结果如下：xi1.01.11.21.31.4y(xi)-1.0000-0.9091-0.8333-0.7692-0.7143yi-0.9000-0.8199-0.7540-0.6986-0.6514注：Euler方法精度较低matlab代码：yi=-1;y=[yi];fori=1:5y2=yi+0.1*[yi+(2+0.1*i)*yi^2];y=[y,y2];yi=y2;end15例2：用改进Euler方法求解初值问题：取步长h=0.1，并与准确解比较解：xi=1+0.1i，于是改进Euler法的计算公式为i=0，1，2，3，4计算结果见P474表7.1-2注：改进Euler方法精度比Euler方法精度高5.0)1(5.11)),()((1)('2yxxyxyxxyxxxy1)(2(1)1(,)()10.1iiiiiiiyyfxyyyxi1111(1)0.110.1(1)(1)0.1210.110.1(1)iiiiiiiiiiyyyyiyyyyyyii1632(0)1EulerEuler,12dyxydxyyyx例设初值问题试分别用法和改进法求解并与精确解进行比较。17计算结果如下表所示：法：改进的法：上结果，此时计算解：取,...)2,1,0())1.0(2(1.0)2(1.0)(21Euler,...)2,1,0()2(1.0]1,0[,1.011212111nkyxkykyxykkkyynyxyyyEulerxhnnnnnnnnnnnnn18xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.0954450.21.1918181.1840971.1832160.31.2774381.2662011.2649110.41.3582131.3433601.3416410.51.4351331.4164021.4142140.61.5089661.4859561.48324019xEuler法y改进的Euler法y精确解01.0000001.0000001.0000000.11.0000001.0959091.095445…………0.71.5803381.5525141.5491930.81.6497831.6164751.6124520.91.7177791.6781661.6733201.01.7847701.7378671.732051207.2RungeKutta方法7.2.1构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：)()!1()(!...)(!2)(')()()()1(1)(21ppippiiiiiyphxyphxyhxhyxyhxyxy),(!...),('!2),()1(21iippiiiiiiyxfphyxfhyxhfyy1(1)111()()(1)!ppiiiheyxyyp21当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程：(7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Tayler公式构造，不实用，高阶导数f(i)不易计算。)()!1()(!...)(!2)(')()()()1(1)(21ppippiiiiiyphxyphxyhxhyxyhxyxy),(!...),('!2),()1(21iippiiiiiiyxfphyxfhyxhfyy1(1)111()()(1)!ppiiiheyxyyp227.2.2RungeKutta方法1.基本思想因为=y(xi)+hf(，y())=y(xi)+hK其中K=f(，y())称为y(x)在[xi，xi+1]上的平均斜率。若取K1=f(xi，y(xi))——Euler公式取K2=f(xi+1，y(xi+1))——向后Euler公式一阶精度取——梯形公式(改进)二阶精度猜想：若能多预测几个点的斜率，再取其加权平均作为K，可望得到较高精度的数值解，从而避免求f的高阶导数。1))(,()()(1iixxiidxxyxfxyxy121()2KKK232.RK公式(7.2-4)其中Kj为y=y(x)在xi+ajh(0aj1)处的斜率预测值。aj，bjs，cj为特定常数。111112),,(),(jssjsijijiipjjjiipjKbhyhaxfKyxfKKchyy243.常数的确定确定的原则是使精度尽可能高。以二阶为例：(7.2-5)(希望y(xi+1)–yi+1=O(hp)的阶数p尽可能高)首先：),(),()(12122122111hKbyhaxfKyxfKKcKchyyiiiiii)()(!21)(')()(321hOxyhxhyxyxyiiii25另一方面：由Taylor展开将K2在(xi，yi)处展开。K2=f(xi+a2h，yi+b21hK1)=f(xi，yi)+a2hfx'(xi，yi)+b21hK1fy'(xi，yi)+O(h2).代入(7.2-5)得：yi+1=yi+hc1f(xi，yi)+hc2f(xi，yi)+h2c2[a2fx'(xi，yi)+b21K1fy'(xi，yi)]+O(h3)=yi+h(c1+c2)f(xi，yi)+c2a2h2[fx'(xi，yi)+(b21/a2)f(xi，yi)fy'(xi，yi)]+O(h3)（希望）...)),()(),(),((000000yxfyyxxyxfyyxxf)()()(')(322221hOxyhacxycchyiii26希望：ei+1=y(xi