资源描述：
《常微分方程的数值解法(III)》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。
1、对一阶常微分方程的初值问题，其一般形式是(1)在下面的讨论中，假定f(x,y)连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数L，使得则初值问题(1)的解必定存在且唯一。常微分方程的数值解法常微分方程的数值解法所谓数值解法，就是要求问题(1)的解在若干点：处的近似值yi(i=0,1,2…n)的方法，yi称为问题(1)的数值解。相邻两个节点的间距称为步长，步长可以相等，也可以不等。本章总是假定hn为定数，称为定步长，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散节点的
2、数值解。一、欧拉（Euler）法1.1Euler公式欧拉（Euler）方法是解初值问题的最简单的数值方法。采用差分方法解初值问题在的近似解即采用向前差商代替得：则Euler法的计算格式i=0,1,…,n(2)还可用数值积分法和泰勒展开法推导Euler格式。以数值积分为例进行推导。将方程的两端在区间上积分得，(3)用左矩形方法计算积分项代入(3)式,并用yi近似代替式中y(xi)即可得到向前欧拉（Euler）公式由于数值积分的矩形方法精度很低，所以欧拉（Euler）公式当然很粗糙。DIMENSIONX(
3、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)END1.2梯形公式为了提高精度,对方
4、程的两端在区间上积分得，改用梯形方法计算其积分项，即(4)代入(7.4)式,并用近似代替式中即可得到梯形公式(5)由于数值积分的梯形公式比矩形公式的精度高，因此梯形公式（5）比欧拉公式(2)的精度高。1.3改进的欧拉公式显式欧拉公式计算工作量小,但精度低。梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。先用欧拉公式(2)求出一个初步的近似值,称为预测值,它的精度不高,再用梯形公式(5)对它校正一次,即迭代一次,求得yi+1,称为校正值,这
5、种预测—校正方法称为改进的欧拉公式：（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)(
6、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设为节点上的近似解，则有改进的Euler格式为二、龙格-库塔（Runge-Kutta）法2.1龙格-库塔(Runge-Kutta)法的基本思想Euler公式可改写成则yi+1的表达式y(xi+1)与的Taylor展开式的前两项完全相同,即局部截断误差为。改进的Euler公式又可改写成（7）为yi+1在xi处的二阶Taylor多项式，为二阶方法，
7、其截断误差为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)和(1
8、1)式前三项相同：以上方程组有无穷多组解，如取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)
