常微分方程的数值解法.ppt课件

1、对一阶常微分方程的初值问题，其一般形式是对一阶常微分方程的初值问题，其一般形式是 (1) (1)在下面的讨论中，假定在下面的讨论中，假定f(x,y)f(x,y)连续，且关于连续，且关于y y满足满足李普希兹李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数L L，使得，使得则初值问题则初值问题(1)(1)的解必定存在且唯一。的解必定存在且唯一。0)(),(yaybxayxfyyyLyxfyxf),(),(常微分方程的数值解法常微分方程的数值解法常微分方程的数值解法常微分方程的数值解法 所谓数值解法，就是要求问题所谓数值解法，就是要求问题(1)的解在若的解在若干

2、点：干点：处的近似值处的近似值yi(i=0,1,2n)的方法，的方法，yi称为问称为问题题(1)的数值解。相邻两个节点的间距的数值解。相邻两个节点的间距 称为步长，步长可以相等，也可以不等。本称为步长，步长可以相等，也可以不等。本章总是假定章总是假定hn为定数，称为定步长，这时节为定数，称为定步长，这时节点可表示为点可表示为 数值解法需要把连续性的问题加以离散化，数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。从而求出离散节点的数值解。 bxxxxann110iinxxh1niihxxi, 2 , 1 ,0一、一、 欧拉欧拉Euler法法1.1 Euler公式公式 欧拉欧拉Eu

3、ler方法是解初值问题的最简方法是解初值问题的最简单的数值方法。采用差分方法解初值问题单的数值方法。采用差分方法解初值问题在在 的近似解的近似解00)(),(yxyyxfybxxxxaxnni110,hxyxyii)()(1y)(,()()()(,()()(11iiiiiiiixyxhfxyxyxyxfhxyxy),()()(1iiiiyxhfxyxy即即采用向前差商采用向前差商 替代替代 得：得：则则EulerEuler法的计算格式法的计算格式 )(),(001xyyyxhfyyiiii i=0,1,n (2 ) 还可用数值积分法和泰勒展开法推导还可用数值积分法和泰勒展开法推导EulerEu

4、ler格式。以格式。以数值积分为例进行推导。数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得， ),(yxfy 1,iixx11),(iiiixxxxdxyxfdxy11)(,)(),()()(1iiiixxixxiidxxyxfxydxyxfxyxy(3) 用左矩形方法计算积分项用左矩形方法计算积分项 )(,)()(,11iiiixxxyxfxxdxxyxfii代入代入(3)式式,并用并用yi近似代替式中近似代替式中y(xi)即可得到向前即可得到向前欧拉欧拉Euler公式公式 ),(1iiiiyxhfyy 由于数值积分的矩形方法精度很低，由于数值积分的矩形方法

5、精度很低，所以欧拉所以欧拉EulerEuler公式当然很粗糙。公式当然很粗糙。 DIMENSION X(0:10),Y(0:10)DOUBLE PRECISION X,Y,H,X0,Y0,N1F(X,Y)=X-Yopen(2,file=aout.txt)H=1.0/10N=10X(0)=0.0Y(0)=1.0DO I=0,N-1 Y(I+1)=Y(I)+F(X(I),Y(I)*H X(I+1)=X(I)+HENDDO10FORMAT(1X,2D15.6/) WRITE(*,10) (X(I),Y(I),I=1,10)WRITE(2,10) (X(I),Y(I),I=1,10)END),(1ii

6、iiyxhfyy1.2 梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 ),(yxfy 1,iixx1)(,)()(1iixxiidxxyxfxyxy)(,()(,(2)(,1111iiiiiixxxyxfxyxfxxdxxyxfii(4 ) 代入代入(7.4)(7.4)式式, ,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 ),(),(2111iiiiiiyxfyxfhyy(5 ) 由于数值积分的梯形公式比矩形公式的精度高，由于数值积分的梯形公式比矩形公式的精度

7、高，因此梯形公式因此梯形公式5 5比欧拉公式比欧拉公式(2)(2)的精度高。的精度高。 1.3 改进的欧拉公式改进的欧拉公式 显式欧拉公式计算工作量小显式欧拉公式计算工作量小,但精度低。梯形公但精度低。梯形公式虽提高了精度式虽提高了精度,但为隐式公式但为隐式公式,需用迭代法求解需用迭代法求解,计算计算工作量大。综合欧拉公式和梯形公式便可得到改进工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。的欧拉公式。 先用欧拉公式先用欧拉公式(2)求出一个初步的近似值求出一个初步的近似值,称为预测值称为预测值, 它的精度不高它的精度不高, 再用梯形公式再用梯形公式(5)对它校对它校正一次正一次,即迭

8、代一次即迭代一次,求得求得yi+1,称为校正值称为校正值, 这种预这种预测测校正方法称为改进的欧拉公式：校正方法称为改进的欧拉公式：1iy),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy1,2 , 1 , 0ni（6 6） 预测预测 校正校正 DIMENSION X(0:10),Y(0:10)DOUBLE PRECISION X,Y,H,X0,Y0,N1F(X,Y)=Y-2*X/YH=1.0/10N=10X(0)=0.0Y(0)=1.0DO I=0,N-1 X(I+1)=(I+1)*H Y(I+1)=Y(I)+H*F(X(I),Y(I) Y(I+1)=Y(I)+(

9、F(X(I),Y(I)+F(X(I+1),Y(I+1)*H/2ENDDO10FORMAT(1X,2D15.6/) WRITE(*,10) (X(I),Y(I),I=1,10)ENDY(I+1)=Y(I)+(F(X(I),Y(I)+F(X(I+1),Y(I)+H*F(X(I),Y(I)*H/2),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy11),(iiiiipyxhfyy),(),(2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy111),(iiiiqpxhfy2/ )(111iiiqpy设设 为节点上的近似解，为节点上的近似解，则有改进的则有改

10、进的EulerEuler格式为格式为 iiizyiihxx,);, 3 , 2 , 1(0二、二、 龙格龙格- -库塔库塔Runge-KuttaRunge-Kutta法法2.1 2.1 龙格龙格- -库塔库塔(Runge-Kutta)(Runge-Kutta)法的基本思想法的基本思想 Euler Euler公式可改写成公式可改写成 ),(111iiiiyxfKhKyy则则yi+1yi+1的表达式的表达式y(xi+1)y(xi+1)与的与的TaylorTaylor展开式的前两展开式的前两项完全相同项完全相同, ,即局部截断误差为即局部截断误差为 。 改进的改进的EulerEuler公式又可改写成

11、公式又可改写成 )(2hO),(),()(21121211hKyxfKyxfKKKhyyiiiiii（7 7） 为为yi+1在在xi处的二阶处的二阶Taylor多项式，为二阶方法，其截多项式，为二阶方法，其截断误差为断误差为O(h3)2.1 2.1 高阶龙格高阶龙格- -库塔库塔(Runge-Kutta)(Runge-Kutta)法的构造法的构造),3,2(),(),(11111piKbhyhaxfKyxfKKchyyjilijnininnpiiinn一般地，一般地，RKRK方法设近似公式为：方法设近似公式为：（8 8） 其中其中ai，bij，ci都是参数，确定它们的原则是使近似公式在都是参数

12、，确定它们的原则是使近似公式在(xn，yn)处的处的Taylor展开式与展开式与y(x)在在xn处的处的Taylor展开式的前展开式的前面的项尽可能多地重合，这样就使近似公式有尽可能高地精面的项尽可能多地重合，这样就使近似公式有尽可能高地精度。以度。以p=2为例：为例：),(),()(12122122111KhbyhaxfKyxfKKcKchyynnnnnn（9 9） yn+1在在(xn，yn)处的处的Taylor展开式为：展开式为：),(),()(12122122111KhbyhaxfKyxfKKcKchyynnnnnn（9 9） )(),(),(),(),()()(),(),(),(),(

13、),(),(,(),(32212211321221212211hOhyxfyxfbyxfacyxfcchyhOyxfyxfhbyxhfayxfcyxfchyyxfhbyhaxfcyxfchyynnnnynnxnnnnnnnynnxnnnnnnnnnnnnnyn+1在在xn处的处的Taylor展开式为：展开式为：)(2),(),(),(),()()(2)()()()(323 21hOyxfyxfyxfhyxhfxyhOxyhxhyxyxynnnnynnxnnnnnnn（1010） （1111） 要使近似公式要使近似公式(8)的局部截断误差为的局部截断误差为O(h3)，则应要求，则应要求(10)和

14、和(11)式前三项相同：式前三项相同：2/12/112122221bcaccc以上方程组有无穷多组解，如取以上方程组有无穷多组解，如取c1=c2=1/2，a2=b21=1，近似，近似公式公式(8)即为改进的即为改进的Euler公式：公式：),(),(2/)(11211hKyhxfKyxfKKKhyynninnnnDIMENSION X(0:10),Y(0:10)DOUBLE PRECISION X,Y,H,X0,Y0,N1,k1,k2F(X,Y)=Y-2*X/YH=1.0/10N=10X(0)=0.0Y(0)=1.0DO I=0,N-1 X(I+1)=(I+1)*H K1=F(X(I),Y(I

15、) K2=F(X(I+1),Y(I)+H*K1) Y(I+1)=Y(I)+H*(K1+K2)/2.0ENDDO10FORMAT(1X,2D15.6/) WRITE(*,10) (X(I),Y(I),I=1,10)END),(),()(21121211hKyxfKyxfKKKhyyiiiiii2.2 四阶龙格四阶龙格库塔法库塔法 四阶经典龙格四阶经典龙格库塔公式。库塔公式。 ),()2,()2,(),()22(631422131212143211hkyxfkkhyxfkkhyxfkyxfkkkkkhyyiiiiiiiiii（1212） DIMENSION X(0:10),Y(0:10)DOUBL

16、E PRECISION X,Y,H,X0,Y0,N1,K1,K2,K3,K4F(X,Y)=X-YH=0.2N=5X(0)=0.0Y(0)=1.0DO I=0,N-1 X(I+1)=(I+1)*H K1=F(X(I),Y(I) K2=F(X(I)+0.5*H,Y(I)+0.5*H*K1) K3=F(X(I)+0.5*H,Y(I)+0.5*H*K2) K4=F(X(I)+H,Y(I)+H*K3) Y(I+1)=Y(I)+H*(K1+2*K2+2*K3+K4)/6ENDDO10FORMAT(1X,2D15.6/) WRITE(*,10) (X(I),Y(I),I=1,10)END),()2,()2,

17、(),()22(631422131212143211hkyxfkkhyxfkkhyxfkyxfkkkkkhyyiiiiiiiiii三、三、 一阶常微分方程组一阶常微分方程组 基于一阶常微分方程初值问题的各种数值解法基于一阶常微分方程初值问题的各种数值解法可类似得到一阶常微分方程组的各种解法，而高阶可类似得到一阶常微分方程组的各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。常微分方程可转化为一阶常微分方程组来求解。 3.1 一阶常微分方程组的初值问题一阶常微分方程组的初值问题pyxhfyyiiii),(1qpxhfyii),(12/ )(1qpyiSUBROUTINE GELR1(T,

18、Y,M,H,N,Z,F,D)DIMENSION Y(M),Z(M,N),D(M)DOUBLE PRECISION Y,Z,D,T,H,XDO I=1,M Z(I,1)=Y(I) ENDDODO J=2,N X=T+(J-2)*H !x(i) CALL F(X,Y,M,D) !求求D DO I=1,M Y(I)=Z(I,J-1)+H*D(I) ENDDO X=T+(J-1)*H !x(i+1) CALL F(X,Y,M,D) !求求D DO I=1,M D(I)=Z(I,J-1)+H*D(I) ENDDO DO I=1,M Y(I)=(Y(I)+D(I)/2.0 Z(I,J)=Y(I) ENDDOENDDOEND11),(iiiiipyxhfyy111),(iiiiqyxhfy2/ )(111iiiqpy SUBROUTINE F(T,Y,M,D)DIMENSION Y(M),D(M)DOUBLE PREC