常微分方程数值解法

1、第一节 问题的提出含有一个或多个导数及其函数的方程式称为微分方程，在工程中常遇到求解微分方程的问题。很多微分方程的解不能用初等函数来表示，有时即使能够用解析式表示其解，但计算量太大而不实用(表达式过于复杂)。需要用数值方法来求解，一般只要求得到若干个点上的近似值或者解的简单的近似表达式(精度要求满足即可)。 常微分方程（一个自变量）微分方程 偏微分方程（一个以上自变量）一一. 定解问题定解问题实际中求解常微分方程的所谓定解问题有两类：初值问题和边值问题定解指已知因变量和/或其导数在某些点上是已知的(约束条件)。1. 边解问题约束条件为已知，在自变量的任一非初值上，已知函数值和/或其导数值，如

2、y ” = f (x,y,y) 求解y y(a)=，y(b)=常常可以将边解问题转化为初值问题求解。2. 初值问题约束条件为在自变量的初值上已知函数值，如：y= f (x,y) dy/dx = f (x,y) xx0 y(x0) = y0 y(x0) = y0 求解y(x)，以满足上述两式，即在 ax0 x1 xnb上的y(xi)的近似值 yi (i=0,1,2,n)。通常取等距节点，即h= xi+1-xi ，有 xi = x0+ih(i=0,1,2,n)初值问题的数值解法特点：按节点顺序依次推进，由已知的y0 , y1 , , yi ，求出yi+1，这可以通过递推公式得到。 单步法：利用前一

3、个单步的信息 (一个点)，在y=f(x) 上找下 一点yi， 有欧拉法， 龙格库格法。 初值问题的常见解法 预测校正法：多步法，利用一个以 上的前点信息求f(x)上 的下一个yi，常用迭 代法，如改进欧拉法， 阿当姆斯法。第二节 欧拉法y= f (x,y) 求解计算y(xi) (i=1,2,n)的近似值yiy(x0) = y0 一、欧拉折线法一、欧拉折线法1. 解一阶方程初值问题的几何意义 y= f (x,y) (x,y)R，R=axb，-y+也即解y = y(x)所表示的积分曲线 y = f(x,y) dx+c上，每一点(x，y)的切线斜率等于函数f(x,y)在该点的值，即： y ( xk)

4、 = f (xk，yk)y(x0)=y0，表示积分曲线从P0（x0，y0）点出发且在P0（x0，y0）的切线斜率为f（x0，y0），故设想积分曲线y = y (x)在x = x0 附近可用切线近似代替。2. 欧拉折线法在 （ x0，y0 ）点上的切线方程为y=y0+ f (x0 , y0)(x-x0)设方程与x = x1交点P1（x1，y1）纵坐标y1取为y(x1)的近似值，则有y(x1) y1 = y0+ f (x0,y0)(x1-x0) = y0+ h f (x0,y0)同理 ：在(x1，y1 )上的切线方程 y = y1+ f (x1,y1)(x-x1)与x=x2交点P2(x2,y2)的

5、纵坐标y2取为y(x2)的近似值有y(x2) y2 = y1+ f (x1,y1)(x2-x1) = y1+ h f (x1,y1)上述的h = x2x1= x1x0过Pi(xi,yi)作斜率为f(xi,yi)的切线方程与x = xi+1的交点，有欧拉公式 y(xi+1) yi+1 = yi+ h f (xi，yi) (i=0,1,n-1) xi = x0 +i h (i=1,2,n)例1：求解初值问题y=-2xy 0 x1.8 y(0)=1取h=0.1解： f(x,y) =-2xy， x0=0， y(x0)=1， h=0.1 所以y(xi+1) yi+1 = yi+ h f (xi,yi)

6、= yi +0.1( -2xiyi) = (1-0.2xi)yi 3. 误差若y(xi) = yi，则y(xi+1) yi+1 0由泰勒展开式在xi上展开有y(xi+1)= y(xi + h)= y(xi) + y (xi) h+(h2/2!) y”()而yi+1 = yi+ h f (xi,yi) = y(xi) + h f(xi, y(xi) = y(xi)+h y(xi)所以y(xi+1) yi+1 = (h2/2) y”() = O(h2)这是局部截断误差(因为认为yi = y(xi)实际上可能yi y(xi)，从而有误差的积累，所以欧拉方法虽然简单，但精度不够，很少采用，只是说明这个

7、方法，有助理解其他方法。二、改进欧拉法二、改进欧拉法1. 梯形公式 y= f(x,y) 在xi , xi+1上的积分有即现被积函数中y (x)未知，用数值积分(矩形公式)11),(iiiixxxxdxyxfdxy1)(,()()(1iixxiidxxyxfxyxy11)(,()(,()()(1iiiixxxxiiidxxyxfydxxyxfxyxyhxyxfdxxyxfiixxii)(,()(,(1用y(xi) yi ，则y(xi+1) yi+1 = yi + h f (xi , yi)若用梯形公式求积分，则所以实际上是取（ xi，yi）及（xi+1，yi+1）两点的平均斜率作为在（xi，yi

8、）点的斜率代入欧拉公式中得到上式的。现有yi+1在公式右边，如用欧拉公式求出y( xi+1)的一个粗糙近似值yi+1（预测值），代入梯形公式中得yi+1（校正值）。)(,()(,(2)(,(111iiiixxxyxfxyxfhdxxyxfii)(,()(,(2)(1111iiiiiiixyxfxyxfhyyxyyi+1（校正值） 较yi+1 （预测值）准确，即有改进的欧拉公式（预测-校正公式） yi+1 = yi + hf (xi，yi) yi+1= yi + h/2f (xi , y(xi), f (xi+1 , y( xi+1)为便于编程，改为 yp= yi + h f（xi，yi） yi

9、+1= yi + h/2（k1 + k2） yc = yi + h f（xi+1，yp） 或 k1 = f（xi，yi） yi+1= 1/2（ yp + yc ） k2 = f（xi+1，yi + h k1 ）此外为提高精度，可迭代求yi+1，即有 yi+1（0） = yi + hf（xi，yi） yi+1（k+1） = yi + h/2f（xi，yi ）， f（xi+1，yi+1 （k） ）直到 yi+1（k+1） - yi+1（k） 时，取y（ xi+1） yi+1（k+1）例2 用改进欧拉法求例1中的初值问题解： yp= yi + hf（xi，yi）=（10.2 xi ） yi yc =

10、 yi + hf（xi+1，yp）= yi 0.2（ xi +0.1） yp yi+1= 1/2（ yp + yc ） 计算见P144 表7-2-22. 误差 y(xi+1)= y(xi+h)= y(xi) +h y( xi)+(h2/2!) y”(xi)+ k1= f（xi，yi）= f（xi，y （ xi ）= y（ xi）k2 = f（xi+1，yi + hk1 ）= f（xi +h ， y （ xi ）+ h k1 ）= f（xi，y （ xi ）+h f（xi，y （ xi ）/x + h k1 f（xi，y （ xi ）/y += y（ xi）+h f（xi，y （ xi ）/x+

11、 y（ xi） f（xi，y （ xi ）/y+O(h2)= y（ xi）+h y”（ xi）+ O(h2)所以：yi+1 = yi + h/2 (k1 + k2) = y(xi)+h/2y(xi)+ y(xi)+h y”(xi)+O(h3) = y(xi)+hy(xi)+h2/2 y”(xi)+O(h3)y(xi +1)yi+1 =O(h3 )局部截断误差O(hp+1)时，方法具有p阶精度。因此欧拉公式有一阶精度，而改进欧拉公式有二阶精度。3.改进欧拉法N-S图及程序读入 x,y读入数据 h,xn输出x,yxxnyl=y+h f(x,y)x = x+h输出x,y结束定义函数 f(x,y)=.

12、y=y+0.5*h*f(x,y)+f(x+h,yl)第三节 龙格库塔方法一、基本思想一、基本思想根据微分中值定理有：y(xi +h)= y(xi+1)y (xi)/h (0 1，h= xi+1 - xi)即f (xi +h, y(xi +h) = y(xi+1)yi/h 是平均斜率y(xi+1)= yi +hf (xi +h , y(xi +h) = yi +hk*记 k* = f (xi +h , y(xi +h)而计算k*的方法是欧拉公式中 k*= f (xi , yi )精度低改进欧拉公式中 k*=0.5*f (xi , yi )+ f (xi+1, yi+1) =0.5*k1+ f (

13、xi +1 , yi +hk1) =0.5*k1 + k2 多取了一个点使精度得以提高。如果设法在xi ， xi+1上多预报几个点的斜率值，然后将它们的加权平均值作为k* = f (xi +h , y(xi +h)的近似值，则有可能构造出更高精度的计算公式，这就是龙格库塔方法的基本思想。二二 、二阶龙格库塔公式、二阶龙格库塔公式在xi ， xi+1内有 xi+l = xi +lh（0l 1）取k*=1 k1 + 2 k2 (是xi, xi+1两个点的斜率值k1, k2加权平均值）则 yi+1 = yi +hk* = yi +h（1 k1 + 2 k2 ）用欧拉公式，取k1= f（ xi ， y

14、i ）而k2 = f（ xi+l， yi+l）= f（ xi +l， yi + lhk1 ） yi+1= yi + h（ 1 k1 + 2 k2 ） k1 = f（ xi ， yi ） k2 = f（ xi +l， yi + lhk1 ）现计算1 ， 2 及l。若要使上述公式具有二阶精度，即y(xi +1)yi+1 = O(h3)设y(xi) = yi ，展开y(xi +1)= y (xi)+h y(xi)+ h2/2 y”(xi)+O(h3) k1= f (xi , yi ) = y(xi)k2 = f (xi +l, yi+lhk1) = f(xi , yi)+ lhfx(xi , yi)

15、+ fy(xi , yi)+ f (xi , yi)+O(h2) = y(xi)+ l h y”(xi)+O(h2)yi+1= y(xi)+ h(1 +2) y(xi)+ h22 l y”(xi)+O(h3)具有二阶精度，只需 1 + 2 = 1 三个未知量，两个方程 2 l = 1/2（1）当l=1时，即xi +l =xi +1，解得1 = 2 = 1 /2 yi+1= yi + h/2（ k1 + k2 ） k1 = f（ xi ， yi ） k2 = f（ xi +l， yi + hk1 ）即为改进欧拉公式（2）取l=1/2，则1 =0，2 =1，有变形的欧拉公式 yi+1= yi +

16、h k2 k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi + h/2 k1 ）三、高阶龙格库塔公式三、高阶龙格库塔公式在xi ， xi+1上除xi ， xi+1外再增加一点xi +m = xi +mh （l m 1） 设xi +m 处的斜率为k3 ，则：k*= 1k1 +2 k2 +3 k3 yi+1= yi+ hk*= yi+ h (1k1 +2 k2 +3 k3)k1= f（ xi ， yi ） k2 = f（xi +lh， yi+lhk1） k3 = f（xi +m， yi+m）= f（xi +mh， yi+mh（1k1 + 2k2 ）得 yi+1= yi +

17、 h（ 1 k1 + 2 k2 + 3 k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +lh， yi +lhk1 ） k3 = f（xi +mh， yi+mh（1k1 + 2k2 ）利用泰勒展开方法，选取1 , 2 , 3 , l , m及 1 , 2使上述公式具有三阶精度O(h4)。得 1 + 2 =1 1 + 2 + 3 =1有7个未知数，5个方程， 2 l+ 3 m=1/2可得一族解，所得公式 2 l2+ 3 m2=1/3为三阶龙格库塔公式 3 lm 2 =1/61. 库塔公式取l=1/2, m=1, 则1 =1/6, 2 =2/3, 3 =1/6, 1 =-1,2

18、 =2得库塔公式 yi+1= yi + h/6（ k1 + 4k2 + k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi +h/2k1 ） k3 = f（xi +h， yi - hk1+ 2hk2 ）2. 四阶龙格库塔公式在xi ， xi+1上用四个点的ki (i=1,2,3,4)做k*的加权平均值，可得一族经典四阶龙格库塔公式。 yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ） k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi +h/2k1 ） k3 = f（ xi +h/2 ， yi+ h/2k2 ）

19、k4 = f（ xi +h， yi+ hk3）例3：用四阶龙格库塔方法求解y=2xy 0 x1.8 y(0)=1取h=0.2局部截断误差O(h5)，精度提高。解： yi+1= yi + 0.2/6（ k1 + 2k2 +2 k3 + k4 ） k1 = -2 xi yi k2 = f（xi +0.1， yi +0.1k1 ）=-2（xi +0.1）（yi +0.1k1 ） k3 = f（xi +0.1， yi+ 0.1k2 ）=-2（xi +0.1）（yi +0.1k2 ） k4 = f（ xi +0.2，yi+ 0.2k3）=-2（xi +0.2）（yi +0.2k3 ）计算结果见P150表

20、7-3-1，在表7-3-2中对比了几种方法的误差。表7-3-3说明再提高阶数没有多大意义。四阶是兼顾了精度和计算量的理想公式。四、四、 四阶龙格四阶龙格-库塔公式的几何意义。库塔公式的几何意义。xixi+hABECDR SPi（ xi， yi ）hk1斜率为f（ xi， yi ）的PA，AM=hki，PA中点R(xi+h/2,yi+h/2ki)的斜率作PB, BM=hk2 ，有PB中点S (xi+h/2,yi+h/2k2)的斜率作PC，CM=hk3，C点(xi+h,yi+hk3)的斜率作PD，DM=hk4，取EM=hk*=h/6*(k1+2k2+2k3+k4 )y(xi+h)= yi+h =

21、yi + 1/6*(k1+2k2+2k3 +k4 )Mhk2hk3hk4五、五、(步长的选择步长的选择)自动定步长四阶龙格库塔自动定步长四阶龙格库塔设从xi出发 ，先以h为步长，用四阶龙格库塔公式求得y(xi+1)的近似值yi+1(h) ，则y(xi+1)yi+1(h)ch5再以h/2为步长，两次计算后(先yi+h/2，再算yi+1)可得yi+1 y(xi+1)yi+1(h/2)c(h/2)5 = c/16* h5有 y(xi+1)yi+1(h/2) yi+1(h/2) yi+1(h) /15即利用事后误差估计法 | yi+1(h/2) yi+1(h) |，判断是否停止计算。六、四阶龙格库塔法

22、N-S图读入x0, y0读入数据 h读入 xnx0 xnk1, k2, k3, k4y0= y0+ h/6(k1+2k2+2k3+k4)x0 = x0+h结束定义函数 f(x, y)输出k1, k2, k3, k4 , x0 , y04 线形多步法线形多步法l单步法中，求解yi+1时，实际已求出了y0，y1yi 。利用多个y值，期望得到较高精度的yi+1ly（xi+1）= y （xi） + xi+1 f（ x ， y （x） ） dxl xil用更精确的求积方法，可用更高次的插值多次来代f（x，y），选取的插值节点不同，则解法不同。l一 阿当姆斯内插公式及误差l除取xi ， xi+1两个节点外

23、，其他节点取在xi ， xi+1时，f 的值未知，故选取xi ， xi+1以外的节点为xi-1 ， xi-2等。l现取xi-1 ， xi-2 ， xi ， xi+1四个插值节点，由于计算yi+1时， yi-2 ， yi-1 ， yi ，已求解出来，插值多项式lL3（x）=（x- xi ）（x- xi-1）（x- xi-2）/（xi+1- xi ）（ xi+1 - xi-1）（ xi+1 - xi-2） f（ xi+1 ， y （ xi+1 ） ）+ （x- xi+1 ）（x- xi-1）（x- xi-2）/（xi- xi+1 i ）（ xi - xi-1）（ xi - xi-2） f（ xi

24、， y （ xi ） ）+ （x- xi+1 ）（x- xi）（x- xi-2）/（xi-1- xi+1 ）（ xi-1 - xi）（ xi-1 - xi-2） f （ xi-1 ， y （ xi-1 ） ）+ （x- xi+1 ）（x- xi）（x- xi-1）/（xi-2- xi ）（ xi-2 - xi-1）（ xi-2 - xi+1） f （ xi-2 ， y （ xi-2 ） ）lf（ x ， y （x） ） L3（x）l所以yi+1 y （ xi） + xi+1 L3（x）dxl xil等距节点时， xi+1 - xi = xi - xi-1 = xi+1 - xi-2 = h，设

25、yi= f（ xi ， y （xi） ）l令x= xi+th 0t 1lyi+1 = yi + h/6 f（ xi+1 ， y （xi+1） ） 1t（t-1）（t+2）dt l 0l-h/2 f（ xi ， y （xi） ） 1（t-1）（t+1）（t+2）dtl 0l+h/2 f（ xi-1 ， y （xi-1） ） 1（t-1）t（t+2）dtl 0l-h/6 f（ xi-2， y （xi-2） ） 1（t-1）t（t+1）dtl 0l=y（ xi ）+h/249 yi+1+19 yi-5 yi-1+ yi-2 l阿当姆斯内插公式是四阶公式ly（ xi+1）- yi+1 =（-19/72

26、0）h5yi（5）+ =0（h5）l内插公式中由于选取了xi+1 点为节点，从而使公式成为隐式，若取xi-3 ， xi-2 ， xi-1 ， xi 作为插值节点，即：l i ilL3（x）= (x- xj )/(xk-xj)f(xj -y(xj )l k=i-3 j=I-3jkl所以y(xj +1)= yi+1 = xi+1 L3（x）dxl xil= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 l外推公式ly(xj +1)- yi+1 = （25/720）h5yi（5）+ =0（h5）l计算f（x，y）的次数，算yi+1 时 yi-1 ，yi-2， yi-3 已

27、在计算yi时得到，故只要算yi就可，故无需算四次，故计算变化四阶龙格-库塔公式，而0（h5）相同，但需用其它方法先行计算y1 ， y2 ， y3 后才能用该法。依次计算y4 ，y5 ，所以不具有自启动性，此外h要改变也很麻烦。l若将外推公式作预测式，内插公式作校正，可得：l 阿当姆斯预测校正式lyj +1（0）= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 lyj +1= yi +h/249f(xj +1 ， yj +1（0)+19 yi+5 yi-1- yi-2 l例4：ly=-2xy 0 x1.2l ly(0)=1l取h=0.2l解：四阶龙格-库塔求出， y1， y2 ， y3 ，由 y0 ， y1 ， y2 ， y3阿当姆斯预测校正公式。计算见p159 表7-6，7-7l一 一阶方程组lu=（x，u，v）, u(x 0)= u0l lv=(x，u，v) , v(x 0)= v0 初值解l l记y=（u，v）T，y0=（ u0 ， v0 ）T F=（ ， ）Tl l y= F（ x ， y ）l l y(x 0)= y0l用龙格库塔解l yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ）l k1 = f（ xi ， yi ）l k2 = f（xi +h/2， yi + h/2k1 ）l k3 = f