常微分方程的数值解法2010

对一阶常微分方程的初值问题，其一般形式是(1)在下面的讨论中，假定f(x,y)连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使得则初值问题(1)的解必定存在且唯一。,常微分方程的数值解法,常微分方程的数值解法所谓数值解法，就是要求问题(1)的解在若干点：处的近似值yi(i=0,1,2n)的方法，yi称为问题(1)的数值解。相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定hn为定长，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。,一、欧拉（Euler）法1.1Euler公式欧拉（Euler）方法是解初值问题的最简单的数值方法。采用差分方法解初值问题在的近似解,即,采用向前差商代替得：,则Euler法的计算格式,i=0,1,n(2),还可用数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。将方程的两端在区间上积分得，,(3),用左矩形方法计算积分项,代入(3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式,由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。,DIMENSIONX(0:10),Y(0:10)DOUBLEPRECISIONX,Y,H,X0,Y0,N1F(X,Y)=X-Yopen(2,file=aout.txt)H=1.0/10N=10X(0)=0.0Y(0)=1.0DOI=0,N-1Y(I+1)=Y(I)+F(X(I),Y(I)*HX(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,1.2梯形公式为了提高精度,对方程的两端在区间上积分得，改用梯形方法计算其积分项，即,(4),代入(4)式,并用近似代替式中即可得到梯形公式,(5),由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（5）比欧拉公式(2)的精度高。,1.3改进的欧拉公式显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。先用欧拉公式(2)求出一个初步的近似值,称为预测值,它的精度不高,再用梯形公式(5)对它校正一次,即迭代一次,求得yi+1,称为校正值,这种预测校正方法称为改进的欧拉公式：,（6）,预测,校正,DIMENSIONX(0:10),Y(0:10)DOUBLEPRECISIONX,Y,H,X0,Y0,N1F(X,Y)=Y-2*X/YH=1.0/10N=10X(0)=0.0Y(0)=1.0DOI=0,N-1X(I+1)=(I+1)*HY(I+1)=Y(I)+H*F(X(I),Y(I)Y(I+1)=Y(I)+(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)END,Y(I+1)=Y(I)+(F(X(I),Y(I)+F(X(I+1),Y(I)+H*F(X(I),Y(I)*H/2,设为节点上的近似解，则有改进的Euler格式为,二、龙格-库塔（Runge-Kutta）法2.1龙格-库塔(Runge-Kutta)法的基本思想Euler公式可改写成,则yi+1的表达式y(xi+1)与的Taylor展开式的前两项完全相同,即局部截断误差为。改进的Euler公式又可改写成,（7）,为yi+1在xi处的二阶Taylor多项式，为二阶方法，其截断误差为O(h3),2.1高阶龙格-库塔(Runge-Kutta)法的构造,一般地，RK方法设近似公式为：,（8）,其中ai，bij，ci都是参数，确定它们的原则是使近似公式在(xn，yn)处的Taylor展开式与y(x)在xn处的Taylor展开式的前面的项尽可能多地重合，这样就使近似公式有尽可能高地精度。以p=2为例：,（9）,yn+1在(xn，yn)处的Taylor展开式为：,（9）,yn+1在xn处的Taylor展开式为：,（10）,（11）,要使近似公式(8)的局部截断误差为O(h3)，则应要求(10)和(11)式前三项相同：,以上方程组有无穷多组解，如取c1=c2=1/2，a2=b21=1，近似公式(8)即为改进的Euler公式：,DIMENSIONX(0:10),Y(0:10)DOUBLEPRECISIONX,Y,H,X0,Y0,N1,k1,k2F(X,Y)=Y-2*X/YH=1.0/10N=10X(0)=0.0Y(0)=1.0DOI=0,N-1X(I+1)=(I+1)*HK1=F(X(I),Y(I)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,2.2四阶龙格库塔法,四阶经典龙格库塔公式。,（12）,DIMENSIONX(0:10),Y(0:10)DOUBLEPRECISIONX,Y,H,X0,Y0,N1,K1,K2,K3,K4F(X,Y)=X-YH=0.2N=5X(0)=0.0Y(0)=1.0DOI=0,N-1X(I+1)=(I+1)*HK1=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,三、一阶常微分方程组基于一阶常微分方程初值问题的各种数值解法可类似得到一阶常微分方程组的各种解法，而高阶常微分方程可转化为一阶常微分方程组来求解。3.1一阶常微分方程组的初值问题,SUBROUTINEGELR1(T,Y,M,H,N,Z,F,D)DIMENSIONY(M),Z(M,N),D(M)DOUBLEPRECISIONY,Z,D,T,H,XDOI=1,MZ(I,1)=Y(I)ENDDODOJ=2,NX=T+(J-2)*H!x(i)CALLF(X,Y,M,D)!求DDOI=1,MY(I)=Z(I,J-1)+H*D(I)ENDDOX=T+(J-1)*H!x(i+1)CALLF(X,Y,M,D)!求DDOI=1,MD(I)=Z(I,J-1)+H*D(I)ENDDODOI=1,MY(I)=(Y(I)+D(I)/2.0Z(I,J)=Y(I)ENDDOENDDOEND,SUBROUTINEF(T,Y,M,D)DIMENSIONY(M),D(M)DOUBLEPRECISIONY,D,TD(1)=Y(2)D(2)=-Y(1)D(3)=-Y(3)END！D存放方程的右端函数值,EXTERNALFDIMENSIONY(3),Z(3,11),D(3)DOUBLEPRECISIONY,Z,D,T,H,XT=0.0Y(1)=-1.0Y(2)=0.0Y(3)=1.0H=0.01M=3N=11！积分点CALLGELR1(T,Y,M,H,N,Z,F,D)WRITE(*,*)DOI=1,NX=(I-1)*HWRITE(*,20)XWRITE(*,30)(Z(J,I),J=1,M)WRITE(*,*)enddo20FORMAT(1X,T=,F5.2)30FORMAT(1X,Y(1)=,D13.6,3X,Y(2)=,D13.6,&3X,Y(3)=,D13.6)END,与考试相关的复习指导,1.高斯直接消去法2.列选主元高斯消去法3.高斯约当消去法,三、解线性方程组的迭代法,1.雅可比（Jacobi）迭代法2.高斯-塞德尔（Gauss-Seidel）迭代法3.超松弛迭代法（SOR方