《数值分析》课程设计常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法常微分方程初值问题数值解的实现和分析—四阶Runge-kutta方法及预估-校正算法摘要求解常微分方程的初值问题，Euler方法，改进的Euler方法及梯形方法精度比较低，所以本文构造高精度单步的四级Runge-kutta方法及高精度的多步预估—校正算法及其Matlab编程来实现对常微分方程初值问题的求解，使在求解常微分方程时，对以前积分方法的收敛速度及精度都有了很高的提高。
关键词：Runge-kutta方法，Adams方法，预估—校正算法，Matlab目录1.前言 (1)2. 几个简单的数值积分法 (2)2.1Runge-kutta方法 (2)2.1.1 Runge-kutta方法的应用 (5)2.2预估—校正算法 (7)2.2.1 Adams数值积分方法简介及预估—校正算法 (7)2.2.2 预估—校正算法的应用 (12)3. 结果分析 (16)总结 (17)参考文献 (18)英文原文和中文翻译 (19)1英文原文 (19)2中文翻译 (20)1.前言常微分方程的初值问题是微分方程定解问题的一个典型代表，以下面的例子介绍常微分方程初值问题数值解的基本思想和原理。
例1.1 一重量垂直作用于弹簧所引起的震荡，当运动阻力与速度的平方成正比时，可借助如下二阶常微分方程描述若令和，则上述二阶常微分方程可化成等价的一阶常微分方程组类似于例1.1，对于m阶常微分方程其中。
若定义可得如下等价的一阶常微分方程组我们知道多数常微分方程主要靠数值解法。
所谓数值解法，就是寻求解在一系列离散节点上的近似值。
相邻两个节点之间的间距称为步长[1]。
2. 几个简单的数值积分法2.1 Runge-kutta方法Runge在1985年提出了一种基于Euler折线法的新的数值方法，此后这种新的数值方法又经过其同胞K.Heun和Kutta的努力[2]，发展完善成为后世所称的Runge-kutta 方法。
Runge重要发现的灵感主要来源于把Euler方法应用于初值问题(2.1.1) 和f不显含y时的初值问题(2.1.2) 之间的类比。
Runge观察到，当把Euler方法应用到(2.1.2)型初值问题时，得到差分方程事实上，这是在计算积分问题(2.1.3) 时采用了左矩形法则即用高为，宽为h的矩形代替了(2.1.3)式中的积分值。
类似的，当把Euler方法应用到初值问题(2.1.1)时，得到差分方程这是在计算积分问题(2.1.4) 时采用了左矩形法则显然(2.1.3)和(2.1.4)式右侧数值积分的精度决定着数值解的精度。
此时，Runge发现，这个矩形公式的逼近程度并不高。
如果在计算积分问题(2.1.4)时采用中点法则或梯形法则，则数值解的精度会更高。
采用中点法则和梯形法则分别代替(2.1.4)式中的积分会得到(2.1.5)(2.1.6) (2.1.5)和(2.1.6)是两个隐式的方程，它们在早期发现Runge-kutta方法过程中扮演了一个重要角色。
Runge的出发点是：分别用一个Euler步代替(2.1.5)中未知的和(2.1.6)式右侧未知的值，这样Runge得到如下的中点格式：和梯形格式：这两个方程具有很好的几何解释。
它们是有折线构成的，这些折线假定微分方程所确定的斜率在前面的点上已经计算出来了。
与他的后继者一样，Runge用Taylor展开和系数对比说明上述两个方法的阶为2. Runge意识到，用中点法和梯形法得出的两个方法的阶数还不够高，所以，他设想，如果使用比中点法和梯形法精度更高的Simpson法则得到的方法阶数会提高，如果用M和T分别表示用中点法和梯形法算得的数值积分。
Simpson法则可以写成众所周知的形式S=M+(T-M)/3。
Taylor展开表明，如果f依赖于y，则这个表达形式只是2阶的，接着Runge发现，通过重复使用Euler步骤对梯形法则做适当的调整，会使格式=M+(-M)/3成为3阶方法，他还把他的方法及其Taylor 展开式拓展到微分方程组[3]。
五年后，Heun在其1900年的文章中评论Runge的方法时说，Runge获得的上述方法是归纳性的而且是令人费解的，他主使用更具一般性的Gauss求积公式其中于是可以把一般的Gauss格式扩充为(2.1.7)把(2.1.7)式的右端进行二元Taylor展开并与的Taylor展开式的对应的系数比较，适当选取参数使方法具有尽可能高的精度。
上述算法公式中的系数希望如此确定，使得(2.1.7)的Taylor展开式中所有h幂次不超过p的那些项与中的相应项的系数相等[2]。
假定p=4并取s=3，用类似的推导，可建立下述常用的显型四级Runge-kutta方法：截断误差为，当右端函数f不依赖于y时，上述公式可简化为2.1.1 Runge-kutta方法的应用用四级Runge-kutta方法求初值问题=y2 e-x，y(1)=1在区间[1,2]上的数值解(取h=0.1) (1) 先求出常微分方程的精确解，再用四阶Runge-kutta方法进行求解。
Y1=dsolve(‘Dy=y^2*exp(-x)’,’y(1)=1’,’x’)运行后的精确解为：Y1=1/(1/exp(x) - 1/exp(1) + 1)现给出常用的四阶Runge-kutta方法求解常微分方程的Matlab主程序如下：function [X,Y]=Rungek(funfcn,x0,b,y0,h)x=x0;y=y0;n=fix((b-x0)/h); %求步数i=1;X=zeros(n,1);Y=zeros(n,1);X(i)=x0;Y(i)=y0; %赋初值for i=2:nk1=feval(funfcn,x,y);k2=feval(funfcn,x+h/2,y+h*k1/2);k3=feval(funfcn,x+h/2,y+h*k2/2);k4=feval(funfcn,x+h,y+h*k3);y=y+h/6*(k1+2*k2+2*k3+k4);Y(i)=y;x=x+h;X(i)=x;end %按照龙格-库塔方法进行求解X,Y1=1./(1./exp(X) – 1./exp(1) + 1), %精确解plot(X,Y,'mh',X,Y1,'bo') %绘图gridlegend('用四阶龙哥库塔方法计算的数值解','精确解y(x)'), %图形说明wcha=abs(Y-Y1), %截断误差自定义函数：function z=funfcn(x,y)z=y.^2.*exp(-x);在matlab工作窗口输入如下程序：x0=1,b=2;y0=1;h=0.1;[X,Y]=Rungek(funfcn,x0,b,y0,h)运行后得：X =1.00001.10001.20001.30001.40001.50001.60001.70001.80001.9000Y1 =1.00001.03631.07141.10541.1380 1.1692 1.1990 1.2273 1.2540 1.2793 wcha = 1.0e-007 * 0 0.0111 0.0290 0.0518 0.0777 0.1054 0.1338 0.1620 0.18960.215911.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9211.051.11.151.21.251.31.35用四阶龙哥库塔方法计算的数值解精确解y=(x)2.2 预估—校正算法2.2.1 Adams 数值积分方法简介及预估—校正算法 构造高阶精度格式的一个重要途径是，将方程改写成如下的积分形式(2.2.1.1) 然后以适当的差值多项式近似替代上式中的被积函数，这时后面将要介绍的Adams差值方法的出发点和思想。
选取k+1个插值结点：。
假定y(x)在上述结点上的近似值已被求出，记为，并令，它被视为在x=上的近似值。
利用上述结点和近似值作k次插值多项式用它近似(2.2.1.1)中的被积函数，得到近似公式(2.2.1.2) 将由此公式确定的当做的近似值。
在上述插值中，被插值点x位于包含所有插值结点的最小区间[]的外部，故称(2)为外插公式。
容易看出，当k=0时，(2.2.1.2)就是向前Euler公式。
为了得到(2.2.1.2)的具体计算公式，令x=，并将在区间写成Newton向后插值公式的形式其中代表j阶向前差分算子。
引进记号则有(2.2.1.3)将表达式(2.2.1.3)代入(2.2.1.2)得到(2.2.1.4) 其中系数(2.2.1.5) 它们与k和n无关。
例如，由于所以外插方法(2.2.1.4)和(2.2.1.5)的截断误差为(2.2.1.6)已知差分可用函数值的线性组合表示，故又可将外插公式(2.2.1.4)改写成(2.2.1.7) 其中系数已造表，应用时可查。
Adams外插公式表K 公式0 ,123仍然从积分关系式(2.2.1.1)出发，将其中的被积函数改用以为插值结点的k次Lagrange插值多项式近似。
这里，由于被插值点x位于包含所有插值结点的最小区间[]的部，故称相应的积分公式为插公式[1]。
用类似的方法可得到Adams插公式表如下Adams插公式表K 公式123在求解常微分方程初值问题的数值积分法中，还有一类高阶精度方法是多步方法。
一般来说，隐式多步方法（如Adams插法）比相应的显式多步方法精度高，并有较好的稳定性质。
然而，隐式方法计算的每一步需要求解非线性方程，所以它的计算比显式多步方法要复杂得多。
为了克服隐式多步方法计算方面的困难，人们经常隐式多步法和显式多步法结合起来使用，即预估—校正的算法。
下面介绍预估和校正的具体做法和计算公式：隐式k步方法可写成如下形式：(2.2.1.8) 其中，可认为是已知的，只是关于的非线性方程。
通常采用迭代法求解，计算公式如下：(2.2.1.9)容易证明，当h适当小并且恰当地选取初值时，上述迭代过程是收敛的。
显而易见，隐式方法(2.2.1.8)的每一步的计算量由(2.2.1.9)的迭代次数决定，因此选取好的初始值是十分重要的。
一种非常自然的做法就是取为用相应的k 步外插法计算的值，即令(2.2.1.10) 由(2.2.1.10)和(2.2.1.9)构成的算法被称为预估—校正算法，其中(2.2.1.10)为预估公式(记为P)，(2.2.1.9)为校正公式(记为C)。
由于这里的预估值比较精确，因此校正中所需的迭代次数一般不会很多。
如果出现迭代收敛慢的情况，应缩小步长在进行计算。
预估—校正算法的计算步骤：首先用预估公式(2.2.1.10)得到，然后计算，紧接着使用校正公式(2.2.1.9)求出，这样完成了一步校正。