常微分方程数值解法

1、.螂袃节莀袅聿膈荿薄袂肄莈蚇肇莃蒇蝿袀艿蒆袁肅膅蒅蚁袈膁蒄螃膄肇蒄袆羇莅蒃薅膂芁蒂蚈羅膇蒁螀膀肃薀袂羃莂蕿薂螆芈薈螄羁芄薈袆袄膀薇薆肀肆薆蚈袃莄薅螁肈芀蚄袃袁膆蚃薃肆肂蚂蚅衿莁蚂袇肅莇蚁羀羇芃蚀虿膃腿芆螂羆肅芆袄膁莄芅薄羄艿莄蚆腿膅莃螈羂肁莂羁螅蒀莁蚀肁莆莀螂袃节莀袅聿膈荿薄袂肄莈蚇肇莃蒇蝿袀艿蒆袁肅膅蒅蚁袈膁蒄螃膄肇蒄袆羇莅蒃薅膂芁蒂蚈羅膇蒁螀膀肃薀袂羃莂蕿薂螆芈薈螄羁芄薈袆袄膀薇薆肀肆薆蚈袃莄薅螁肈芀蚄袃袁膆蚃薃肆肂蚂蚅衿莁蚂袇肅莇蚁羀羇芃蚀虿膃腿芆螂羆肅芆袄膁莄芅薄羄艿莄蚆腿膅莃螈羂肁莂羁螅蒀莁蚀肁莆莀螂袃节莀袅聿膈荿薄袂肄莈蚇肇莃蒇蝿袀艿蒆袁肅膅蒅蚁袈膁蒄螃膄肇蒄袆羇莅蒃薅膂芁蒂蚈羅膇

2、蒁螀膀肃薀袂羃莂蕿薂螆芈薈螄羁芄薈袆袄膀薇薆肀肆薆蚈袃莄薅螁肈芀蚄袃袁膆蚃薃肆肂蚂蚅衿莁蚂袇肅莇蚁羀羇芃蚀虿膃腿芆螂羆肅芆袄膁莄芅薄羄艿莄蚆腿膅莃螈羂肁莂羁螅蒀莁蚀肁莆莀螂袃节莀袅聿膈荿薄袂肄莈蚇肇莃蒇蝿袀艿蒆袁肅膅蒅蚁袈膁蒄螃膄肇蒄袆羇莅蒃薅膂芁蒂蚈羅膇蒁螀膀肃薀袂羃莂蕿薂螆芈薈螄羁芄薈袆袄膀薇薆 第7章 常微分方程数值解法7.0 基本概念1. 一阶常微分方程的初值问题 (7.0-1)注：若f在D = a £ x £ b , |y|<+¥内连续，且满足Lip条件：$L ³0，使|f (x y1) f (x，y2)| £ L|y1 y2|

3、 (7.0-2)则(7.0-1)的连续可微解y(x)在a，b上唯一存在。2. 初值问题的数值解称(7.0-1)的解y(x)在节点xi处的近似值yi » y(xi) a < x1 <x2 < . < xn = b.为其数值解，方法称为数值方法。注： 考虑等距节点: xi = a + ih，h = (b a)/n. 从初始条件y(a) = y0出发，依次逐个计算y1，y2，yn的值，称为步进法。两种：单步法、多步法。 二阶常微分方程y''(x) = f (x，y(x)，y'(x)可设为一阶常微分方程组的初值问题：引进新的未知函数z(x) =

4、 y'(x)，则其初始条件为：称为一阶微分方程组的初值问题，方法类似。 边界问题，常用差分方法解。7.1 初值问题数值解法的构造及其精度7.1.1 构造方法对于(7.0-1)可借助Taylor展开（导数法）、差商法、积分法实现离散化来构造求积公式：1. 设y Î Ca，b将y(xi+1) = y(xi+h)在xi处展开 x Î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求解公式，（E

5、uler法）2. 用差商来表示：得差分方程：Þ 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).3. 对(7.0-1)两边取积分得 (7.1-3)取不同的数值积分可得不同的求解公式，为： 用矩形公式：Þ y(xi+

6、1) » y(xi) + hf (xi，y(xi) Þ Euler 公式 y(xi+1) » y(xi) + hf (xi+1，y(xi+1) Þ 向后Euler 公式 用梯形公式：Þ Þ (7.1-4)称(7.1-4)为梯形公式¾¾隐式公式。显化：预估值： 校正值：.4. 几何意义Euler法¾¾折线法改进Euler法¾¾平均斜率折线法例1：例2： P473, P4747.1.2 截断误差与代数精度定义7.1-1 称 ei = y(xi) yi 为数值解yi的（整体）截断误

7、差。 若yk = y(xk)，k = 0，1，2，i 1. 由求解公式得数值解，则称为yi的局部截断误差。注：局部截断误差是指单步计算产生的误差，而（整体）截断误差则考虑到每步误差对下一步的影响。定义7.1-2 若求解公式的（整体）截断误差为O(h p)则称该方法是p阶方法，或是p阶精度。<p阶公式>定理7.1-1 设数值解公式：yi+1 = yi + hj(xi，yi，h)中的函数j(x，y，h)关于y满足Lipschitz条件：，且其局部截断误差为hp+1阶，则其（整体）截断误差为hp阶，即该数值解公式为p阶方式。注： 局部截断误差较易估计 定理7.1-1表明：若ei = O(

8、hp+1) 则ei = O(hp). Euler局部截断误差为 所以一阶精度。向后Euler法也是一阶精度。 梯形公式为二阶精度。例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.1yi + (2 + 0.1i)yi2，i = 0，1，2，3，4计算结果见P473表7.1-1注：Euler方法精度较低例2：用改进Euler方法求解初值问题：取步长h = 0.1，并与准确解比较解：xi

9、 = 1 + 0.1i，于是改进Euler法的计算公式为i = 0，1，2，3，4计算结果见P474表7.1-2注：改进Euler方法精度比Euler方法精度高7.2 Runge¾Kutta方法 7.2.1 构造高阶单步法的直接方法由Taylor公式：当h充分小时，略去Taylor公式余项，并以yi、yi+1分别代替y(xi)、y(xi+1)，得到差分方程： (7.2-1)其局部截断误差为：即(7.2-1)为p阶方式，上述方式称为Taylor方式。注：利用Taylor公式构造，不实用，高阶导数f (i)不易计算。7.2.2 Runge¾Kutta方法1. 基本思想因为 =

10、y(xi) + hf (x，y(x) = y(xi) + hKx其中Kx = f (x，y(x)称为y(x)在xi，xi+1上的平均斜率。若取 K1 = f (xi，y(xi) Euler公式 取 K2 = f (xi+1，y(xi+1) 向后Euler公式 一阶精度 取 梯形公式 二阶精度猜想：若能多预测几个点的斜率，再取其加权平均作为Kx，可望得到较高精度的数值解，从而避免求f的高阶导数。2. R¾K公式 (7.2-4)其中Kj为y = y(x)在xi + ajh (0 £ aj £ 1)处的斜率预测值。aj，bjs，cj为特定常数。3. 常数的确定确定的原则

11、是使精度尽可能高。以二阶为例： (7.2-5)(希望y(xi+1) yi+1 = O(hp)的阶数p尽可能高)一方面：另一方面：将K2在(xi，yi)处展开。K2 = f (xi，yi) + a2hf x'(xi，yi) + b21hK1 f y'(xi，yi) + O(h2).代入(7.2-5)得：yi+1 = yi + hc1 f (xi，yi) + hc2 f (xi，yi) + h2c2a2 f x'(xi，yi) + b21K1 f y'(xi，yi) + O(h3) = yi + h(c1 + c2) f (xi，yi) + c2a2h2f x

12、9;(xi，yi) + (b21/a12) f (xi，yi) f y'(xi，yi)+O(h3)（希望）希望：ei+1 = y(xi+1) yi+1 = O(h3). 则应：特例：a2 = 1 Þ c1 = c2 = 1/2，b21 = 1，得2阶R-K公式： ¾¾改进欧拉公式。c1 = 0 Þ c2 = 1，a2 = 1/2，b21 = 1/2，得： (7.2-7)称为中点公式。4. 最常用的R-K公式 标准4阶R-K公式 (7.2-8)算法：输入a，b，n，y0h=(b-a)n，x0 = afor i = 1, i<=n, i+K1

13、= f(x0, y0)K2 = f(x0+h/2, y0+h*K1/2)K3 = f(x0+h/2, y0+h*K2/2)K4 = f(x0+h, y0+h*K3)x0 = x0+hy0 = y0 + h*(K1+2*K2+2*K3+K4)/6输出x0，y0function Runge_Kutta4(a,b,h,y0)n=(b-a)/h;x0=a;for i=1:n K1=f(x0,y0) K2=f(x0+h/2,y0+h*K1/2) K3=f(x0+h/2,y0+h*K2/2) K4=f(x0+h,y0+h*K3) x0=x0+h y1=y0+h*(K1+2*K2+2*K3+K4)/6; y

14、0=y1end;end;function f=f(x,y) f=x+y;end;例1：（P478）用标准4阶R-K公式求：的数值解。取h = 0.2，并与标准解y = 2ex x 1比较。解：因为f(x，y) = x + y，从而由(7.2-8)得：xi+1K1K2K3K4yi+1y(xi+1)0.211.21.221.4441.2428001.2428060.41.4428001.6870801.7115081.9851021.5836361.5836490.61.9836362.2820002.3118362.6460032.0442132.0442380.82.6442133.00863

15、43.0450763.4532282.6510422.6510821.03.4510423.8961463.9406574.4391733.436503注：步长h的选择 (P479)使用数值解法求解初值问题(7.0-1)，选择步长h是一个重要问题。从每一步看，h小局部截断误差小，整体截断误差也就小；但从整个区间看，h小则节点多，这不仅使计算工作量增大，而且也使舍入误差的累积严重，导致数值不稳定，所以步长h的大小要适当。在实际计算中，如前所述，常常采用事后估计误差及自动选择步长的办法来保证节点xi处的数值解yi(i = 1，2，N)满足所要求的精度。设所用的一步法是p阶方法，其整体截断误差有渐近

16、公式y(xi+1) yi+1 = Mhp + O(hp+1)其中M是与h无关的常数。于是可用Richardson外推法提高数值解的精度。事实上，从节点xi出发，先以h为步长经一步计算求出数值解，则整体截断误差如上式所示为 (7.2-9)然后以h/2为步长，从节点xi出发经过两步计算求出数值解，则整体截断误差为 (7.2-10)用2p乘(7.2-10)式减去(7.2-9)式，得故有： (7.2-11) (7.2-12)或： (7.2-13)注：(1)(7.2-11)表明若取，则整体截断误差为hp+1，即精度提高一阶。(7.2-12)表明数值解的整体截断误差近似为(7.2-13)表明数值解的整体截

17、断误差近似为(2) 记，e为精度要求，则当时，有且即步长h / 2是合适的；当s > e 时，说明步长h / 2仍然偏大，须将步长减半，继续计算；当s < e / 2p时，说明已有，步长h / 2偏小，应取步长h。例2(P481) 用变步长的标准4阶R-K方法求初值问题：的数值解，要求精度为7.3 线性多步法希望避免求多个点上f (x，y)的值，并且充分利用前面几步的结果。一般形式：yi+1 = a0 yi + a1 yi-1 + + ap yi-p + h(b-1 y'i+1 +b0 y'i + + bp y'i-p) (7.3-1)其中y'k =

18、 f (xk，yk)，(k = i p，i p +1，i，i+1)而ak，bk为待定常数。注：(1) 若b1 = 0时，(7.3-1)为显式公式，否则为隐式公式。(2) 推导方法主要有：数值积分法和Taylor展开法。7.3.1 数值积分法取节点：xi，xi-1，xi-2，xi-3，作f的三次L-插值多项式，有记：f k = f(xk，yk) (k = i，i 1，i 2，i 3)xi-k = xi kh，x = xi + th.代入上式有：所以： Þ (7.3-2)其中f k = f(xk，yk) (k = i，i 1，i 2，i 3)称(7.3-2)为4阶Adams显式公式。其局

19、部截断误差 (7.3-3)若取节点：xi+1，xi，xi-1，xi-2，作f的3次L-插值多项式，可得4阶Adams隐式公式： (7.3-4)其局部截断误差： (7.3-5)注： 隐式公式的显化：（预测校正） (7.3-6) 并非所有线性多步法公式(7.3-1)都可用数值积分法得到，但都可用Taylor展开法得到。7.3.2 Taylor 展开法设yi-k = y(xi kh)，y'i-k = y'(xi kh)展开为：代入(7.3-1)得： (7.3-7)为使(7.3-1)有m阶精度，只须(7.3-7)的前m+1项与y(xi+1)的展式：对应相等，即有方程组。 (7.3-8)

20、此时有(7.3-9)特别取p = 3，m = 4有：令a0 = a1 = a2 = b-1 = 0，可得a3 = 1，b0 = 8/3，b1 = -4/3，b2 = 8/3，b3 = 0.代入(7.3-1)得Milne公式： (7.3-10)即其中： (7.3-11)令a1 = a3 = b2 = b3 = 0，得Hamming 公式： (7.3-12)其中： (7.3-13)注： 与单步法相比，多步法不须反复计算f在(xi，xi+1)上某些点处的值，工作量大大减少。 多步法的前几步须用同阶单步法求之。例1 用4阶Adams显式公式求初值问题：的数值解，取h = 0.05。解：先用标准4阶Ru

21、nge-Kutta公式求出此初值问题在x1 = 0.05，x2 = 0.1，x3 = 0.15处的数值解，然后用公式(7.3-2)求其余节点的数值解（p487） Milne显，Hamming隐，显化得Milne-Hamming公式： (7.3-14)7.4 预估校正系统直接预估校正格式：先用显式公式算出预估值，再用同阶隐式公式进行校正，没有充分利用局部截断误差的信息。利用误差补偿的办法，对预估值和校正值进行修正，可以使计算结果的精度更高一些。以4阶Adams为例：显式：局部截断误差 （7.3-3）隐式：局部截断误差 (7.3-5)设pi+1，ci+1分别表示xi+1处数值解的预估值和校正值，则

22、 (7.4-1) (7.4-2)两式相减： (7.4-3)代入(7.4-1)得 (7.4-4)即：以作y(xi+1)的预估值，精度可提高一阶。同理，将(7.4-3)代入(7.4-2)得 （7.4-5）即：以作为y(xi+1)的校正值，其精度可提高一阶。注：由于预估值ci+1还未算出，所以用前一步的ci，pi对pi+1进行修正，由此得： 预估 修正 校正 修正同理，对Milne公式（显）和Hamming公式（隐）可得带有误差补偿的预估校正系统公式： (7.4-7)注：在(7.4-6)和(7.4-7)中显然无c3，p3，可取p3 = c3 = 0.例1：取步长h = 0.1，用带误差补偿的预估校正公式(7.4-6)求解解：所需y1，y2，y3用标准4阶R-K公式计算，然后按(7.4-6)计算y4，y5. P4927.5 边值问题的差分法基本思想：运用数值微分将导数用离散点上函数值表示，从而将边值问题的微分方程和边界条件转化为只含有限个未知数的差分方程组，并将此差分方程组的解作为该边值问题的数值解。1. 二阶常微分方程的第一边值问题 (7.5-1)其中q(x)(³0)，f(x)在a，b)上连续，a，b为常数。设等距节点：xi = a + ih，i = 0，1，2，n，对其中内节点应用三点微分公式： (6.6-9)