计算方法与实习常微分方程数值解法

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 ，求出

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

4、于函数f(x,y)在该点的值，即： y ( xk) = 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+

5、 f (x1,y1)(x-x1)与x=x2交点P2(x2,y2)的纵坐标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

6、， h=0.1 所以y(xi+1) yi+1 = yi+ h f (xi,yi) = yi +0.1( -2xiyi) = (1-0.2xi)yi 计算见P142 表7-2-1。,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) 这是局部截断误差(因为

7、认为yi = y(xi) 实际上可能yi y(xi)，从而有误差的积累，所以欧拉方法虽然简单，但精度不够，很少采用，只是说明这个方法，有助理解其他方法。,二、改进欧拉法 1. 梯形公式 y= f(x,y) 在xi , xi+1上的积分有 即 现被积函数中y (x)未知，用数值积分(矩形公式),用y(xi) yi ，则 y(xi+1) yi+1 = yi + h f (xi , yi) 若用梯形公式求积分，则 所以 实际上是取（ xi，yi）及（xi+1，yi+1）两点的平均斜率作为在（xi，yi）点的斜率代入欧拉公式中得到上式的。 现有yi+1在公式右边，如用欧拉公式求出y( xi+1)的一个

8、粗糙近似值yi+1（预测值），代入梯形公式中得yi+1（校正值）。,yi+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+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+

9、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 = yi + hf（xi+1，yp）= yi 0.2（ xi +0.1） yp yi+1= 1/2（ yp + yc ） 计算见P144 表7-2-2 2. 误差 y(xi+1)= y(xi+h)= y(xi) +h y( xi)+(h

10、2/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+ 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

11、”(xi)+O(h3) = y(xi)+hy(xi)+h/2 y”(xi)+O(h3) y(xi +1)yi+1 =O(h3 ) 局部截断误差O(hp+1)时，方法具有p阶精度。 因此欧拉公式有一阶精度，而改进欧拉公式有二阶精度。,3.改进欧拉法N-S图及程序,读入 x,y,读入数据 h,xn,输出x,y,xxn,yl=y+h f(x,y),x = x+h,输出x,y,结束,定义函数 f(x,y)=.,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 - x

12、i) 即 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 (xi +1 , yi +hk1) =0.5*k1 + k2 多取了一个点使精度得以提高。,如果设法在xi ， xi+1上多预报几个点的斜率值，然后将它们的加权平均值作为k* = f (x

13、i +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 ， yi ）而 k2 = f（ xi+l， yi+l）= f（ xi +l， yi + lhk1 ） yi+1= yi + hf（ 1 k1 + 2 k2 ） k1 = f（ xi ， yi ） k2

14、 = 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)+ 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

15、 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 + h k2 k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi + h/2 k1 ）,三、高阶龙格库塔公式 在xi ， xi+1上除xi ， xi+1

16、外再增加一点 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 + h（ 1 k1 + 2 k2 + 3 k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +lh， yi +lhk1 ） k3 = f（xi +mh， yi+mh

17、（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/6 1. 库塔公式 取l=1/2, m=1, 则1 =1/6, 2 =2/3, 3 =1/6, 1 =-1,2 =2 得库塔公式 yi+1= yi + h/6（ k1 + 4k2 + k3 ） k1 = f（ xi ， yi ） k2 = f（ xi +h/2， yi +h/2k1

18、 ） 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 ） k4 = f（ xi +h， yi+ hk3） 例3：用四阶龙格库塔方法求解 y=2xy 0 x1.8 y(0)=1取h=0.2,局部截断误差O(h5)，精度提高。,

19、解： 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表7-3-1，在表7-3-2中对比了几种方法的误差。表7-3-3说明再提高阶数没有多大意义。四阶是兼顾了精度和计算量的理想公式。,四、 四阶龙格-库塔公式的几

20、何意义。,xi,xi+h,A,B,E,C,D,R,S,Pi（ 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 = yi + 1/6*(k1+2k2+2k3 +k4 ),M,hk2,hk3,hk4,五、自动定步长四阶龙格库塔 设从xi出发 ，先以h为步长，用四阶龙格库塔公式求

21、得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) |，判断是否停止计算。,六、四阶龙格库塔法N-S图,读入x0, y0,读入数据 h,读入 xn,x0xn,k1, k2, k3, k4,y0= y0+ h/6(k1+2k2+2k3+k4),x0 = x0+h,结束,定义函数 f(

22、x, y),输出k1, k2, k3, k4 , x0 , y0,4 线形多步法,单步法中，求解yi+1时，实际已求出了y0，y1yi 。利用多个y值，期望得到较高精度的yi+1 y（xi+1）= y （xi） + xi+1 f（ x ， y （x） ） dx xi 用更精确的求积方法，可用更高次的插值多次来代f（x，y），选取的插值节点不同，则解法不同。 一 阿当姆斯内插公式及误差 除取xi ， xi+1两个节点外，其他节点取在xi ， xi+1时，f 的值未知，故选取xi ， xi+1以外的节点为xi-1 ， xi-2等。 现取xi-1 ， xi-2 ， xi ， xi+1四个插值节点，由

23、于计算yi+1时， yi-2 ， yi-1 ， yi ，已求解出来，插值多项式,L3（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 ， y （ xi ） ）+ （x- xi+1 ）（x- xi）（x- xi-2）/（xi-1- xi+1 ）（ xi-1 - xi）（ xi-1 - xi-2） f （ xi-1

24、， 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 ） ） f（ x ， y （x） ） L3（x） 所以yi+1 y （ xi） + xi+1 L3（x）dx xi 等距节点时， xi+1 - xi = xi - xi-1 = xi+1 - xi-2 = h，设yi= f（ xi ， y （xi） ） 令x= xi+th 0t 1,yi+1 = yi + h/6 f（ xi+1 ， y （xi+1） ） 1t（t-1）（t+2）dt 0 -

25、h/2 f（ xi ， y （xi） ） 1（t-1）（t+1）（t+2）dt 0 +h/2 f（ xi-1 ， y （xi-1） ） 1（t-1）t（t+2）dt 0 -h/6 f（ xi-2， y （xi-2） ） 1（t-1）t（t+1）dt 0 =y（ xi ）+h/249 yi+1+19 yi-5 yi-1+ yi-2 阿当姆斯内插公式是四阶公式 y（ xi-1）- yi+1 =（-19/720）h5yi（5）+ =0（h5）,二 阿当姆斯外推公式及误差,内插公式中由于选取了xi+1 点为节点，从而使公式成为隐式，若取xi-3 ， xi-2 ， xi-1 ， xi 作为插值节点，即：

26、 i i L3（x）= (x- xj )/(xk-xj)f(xj -y(xj ) k=i-3 j=I-3jk 所以y(xj +1)= yi+1 = xi+1 L3（x）dx xi = yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 外推公式 y(xj +1)- yi+1 = （25/720）h5yi（5）+ =0（h5）,计算f（x，y）的次数，算yi+1 时 yi-1 ，yi-2， yi-3 已在计算yi时得到，故只要算yi就可，故无需算四次，故计算变化四阶龙格-库塔公式，而0（h5）相同，但需用其它方法先行计算y1 ， y2 ， y3 后才能用该法。依次计算y4

27、，y5 ，所以不具有自启动性，此外h要改变也很麻烦。 若将外推公式作预测式，内插公式作校正，可得： 阿当姆斯预测校正式 yj +1（0）= yi +h/2455 yi+59 yi-1+37 yi-2-9 yi-3 yj +1= yi +h/249f(xj +1 ， yj +1（0)+19 yi+5 yi-1- yi-2 ,例4： y=-2xy 0 x1.2 y(0)=1 取h=0.2 解：四阶龙格-库塔求出， y1， y2 ， y3 ，由 y0 ， y1 ， y2 ， y3阿当姆斯预测校正公式。计算见p159 表7-6，7-7,5 一阶方程组与高阶方程,一 一阶方程组 u=（x，u，v）, u(x 0)= u0 v=(x，u，v) , v(x 0)= v0 初值解 记y=（u，v）T，y0=（ u0 ， v0 ）T F=（ ， ）T y= F（ x ， y ） y(x 0)= y0 用龙格库塔解, yi+1= yi + h/6（ k1 + 2k2 +2 k3 + k4 ） k1 = f（ xi ， yi ） k2 = f（xi +h/2， yi + h/