常微分方程的几种数值解法

1、常微分方程的几种数值解法数学与应用数学 肖振华 指导教师 张秀艳【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、RungeKutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用

2、范围有一个直观的感受。【关键词】 常微分方程 数值解法 MATLAB 误差分析【Abstract】 Many phenomena in nature and engineering can be attributed to the definite solution of the problem for differential equations. Among them, the ordinary differential equation solving is an important foundation for the content of the differential equat

3、ions. However, many of the differential equations are often difficult to obtain accurate analytical expression .At this time, the numerical solution provides a good idea. For the Numerical Solution of Ordinary Differential Equations in this article, we focuses on some commonly used numerical solutio

4、n, such as the Euler method, improved Euler method, Runge-Kutta method, Adams predictor corrector method as well as newer spectral methods. Through specific examples, combined with MATLAB solving and drawing, we initially know the solution process of general numerical solution of ordinary differenti

5、al equations . At the same time, according to the error analysis of various methods , everyone has an intuitive feel of the characteristics and scope of the various methods.【Keywords】Ordinary Differential Equations Numerical Solution MATLAB error analysis目录1 前言31.1 常微分方程概述31.2 常微分方程解与数值解法32 欧拉法和改进的欧

6、拉法42.1 欧拉法42.2 改进的欧拉方法52.3 算例53龙格-库塔法133.1 龙格-库塔法与泰勒展开133.2 龙格-库塔法公式与ode函数143.3算例174 阿达姆斯预估校正法214.1 阿达姆斯（Adams）公式214.2 预估校正方法234.3算例235 勒让德谱方法285.1 谱方法介绍285.2勒让德多项式与谱方法285.3 算例30参考文献391 前言1.1 常微分方程概述方程是一个在数学中非常熟悉的名词，在初等数学里，我们将我们要研究的问题作为一个或几个未知量，通过观察事物的规律，得出这些未知量与已知量之间的等式关系，这样就得到了一个简单的方程或方程组当然，这只是一个很

7、浅显粗略的定义。在数学上，物质的运动和变化规律是通过函数关系来表示的，在一些复杂的现象中，比如火箭的运动轨迹，我们要求的未知量就变成了满足特定条件的一个或一些未知函数（不再是简单的数值）。有时候，我们需要用到导数很微分的知识，即这些未知函数的导数与自变量满足某种关系，这种方程，我们称之为微分方程，即可以这样定义，凡是含有参数，未知函数和未知函数导数（或微分）的方程称之为微分方程。未知函数是一元函数的微分方程称为常微分方程（ODE），未知函数是多元函数的微分方程称为偏微分方程（PDE），我们要研究的，就是常微分方程。1.2 常微分方程解与数值解法微分方程的解，就是找出一个代入方程使之成为恒等式的

8、函数。.若微分方程的解中含有任意常数的个数与方程的阶数相同，且任意常数之间不能合并，则称此解为该方程的通解(或一般解).当通解中的各任意常数都取特定值时所得到的解，称为方程的特解。在实际问题中，这些函数往往还满足一些特定条件，这个称之为定解问题，如常微分方程的初值问题和边值问题，满足条件的解就是特解。在理想的状况下，我们解这个方程如果能够得出通解的解析表达式 ，那么代入数值很容易得到问题所要求的特解，但事实上，很多常微分方程得不到通解的解析表达式，或者表达式过于复杂。而且，常微分方程的特解是否存在，有几个解，牵涉到解的存在唯一性定理。在实际应用中，我们对解的精度要求往往并没有那么严格。在这种情

9、况下，数值解应运而生，通俗的讲，数值解就是所求函数在特定点的函数值，通过函数在特定点的取值，我们可以反过来观察这个函数的某些性质，所以，求常微分方程的数值解是很有意义的，事实上，解常微分方程很大一部分就是解数值解。 下面我们要讨论几种常见的常微分方程数值解法，目的是展示常微分方程数值解法的一般过程。以下章节，由于本人水平有限，只对数值解法的思路和过程进行讨论，而对于解的存在唯一性，方法精度的具体证明，不做过多讨论。一般而言，构造数值解法都是基于两个想法：一种是基于数值积分的构造方法，一种是基于泰勒展开的构造法。下面我们要讨论的是几种常见的数值解法，如欧拉法、改进的欧拉法、龙格库塔法、阿达姆斯预

10、估校正法以及基于勒让德多项式的谱方法。这些方法既有基于数值积分的构造方法，如欧拉法，也有基于泰勒展开的构造法，如龙格库塔法。下面的章节，我们会一一介绍这些方法的思路，并通过一些简单的算例，让大家了解这些解法的一般过程。为了让大家看到各种数值解法的精确性和有效性，我们可能会尽量选取一些通解存在的算例来求数值解，这样做并不意味着我们实际中求数值解时要先求出精确解，而只是出于对方法对比的需要，望不致引起混淆。2 欧拉法和改进的欧拉法2.1 欧拉法 2.1.1 欧拉法介绍首先，我们考虑如下的一阶常微分方程初值问题 （2-1）事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将看做向量，（2-1

11、）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2-1）中的导数项用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。设在中取等距节点，因为在节点点上，由（2-1）可得： ， （2-2） 又由差商的定义可得： （2-3） 所以有 （2-4）用的近似值代入（2-4），则有计算的欧拉公式 （2-5）2.1.2欧拉法误差分析从欧拉公式中可以看出，右端的都是近似的，所以用它计算出来的会有累计误差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为是精确的前提下来估计，记为，

12、泰勒展开有 （2-6）联立（2-5），（2-6）即得=+=，根据数值算法精度的定义，如果一个数值方法的局部截断误差=则称这个算法具有P阶精度，所以，欧拉方法具有一阶精度或者称欧拉方法为一阶方法。2.2 改进的欧拉方法2.2.1 改进的欧拉法介绍用数值积分离散化问题（1），两边做积分有： （2-7）对右端积分使用梯形积分公式可得： （2-8）同欧拉方法，用的近似值代入（2-7），联立（2-8）得到改进的欧拉方法计算公式： （2-9）一般来说，如果求解时，算法右端不包含，称为显性计算公式，如果包含，则求解时还需要解方程，这种称为隐式计算公式。显然公式（2-9）是一个隐式计算公式，事实上，改进的欧拉

13、方法是用欧拉方法先求一个预测值，再用这个预测值来计算，即： （2-10）2.2.2改进的欧拉法误差分析同欧拉法误差分析类似，用泰勒展开容易知道改进的欧拉方法具有二阶精度，证明略。2.3 算例2.3.1（一阶常微分方程）求解初值问题解析：在MATLAB中求解这个方程y=dsolve('Dy=y-2*x/y','y(0)=1','x') 得y =(2*x+1)(1/2) 它的解析解为，下面我们分别用欧拉方法和改进的欧拉方法来求其数值解。欧拉方法：创建M文件euler1.m,内容如下：function x,y=euler1(fun,x0,xfinal,

14、y0,n)if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(fun,x(i),y(i);End再定义函数方程组中的函数f1，创建f1.m文件，内容如下：function f=f1(x,y)f=y-2*x/y在MATLAB中输入：x,y=euler1('f1',0,1,1,20)输出f,x,y的值，将数值解跟精确解画图表示，输入：plot(x,y,'r*-',x,sqrt(1+2*x),'g+-')；

15、 xlabel('x');ylabel('y');title('y=y-2x/y');legend('数值解','精确解')得到图形，保存为euler1.fig，图形如下：改进的欧拉方法：创建M文件eulerprove1.m,内容如下：function x,y=eulerprove1(fun,x0,xfinal,y0,n)if nargin<5,n=50;endh=(xfinal-x0)/n;y(1)=y0,x(1)=x0;for i=1:n x(i+1)=x(i)+h; y1=y(i)+h*feval(fu

16、n,x(i),y(i)/2; y2=h*feval(fun,x(i+1),y1)/2; y(i+1)=y1+y2End在MATLAB的command窗口输入: x,y1=eulerprove1('f1',0,1,1,20)返回f,x1,y1的值，作图，输入：plot(x,y1,'r*-',x,sqrt(1+2*x),'g+-')； xlabel('x');ylabel('y');title('y=y-2x/y');legend('数值解','精确解')，将图片保存为e

17、ulerprove1.fig,图形如下：为了便于比较两种方法的误差，将两者的误差作到同一个图上，继续输入：plot(x,abs(y-sqrt(1+2*x),'y*-',x,abs(y1-sqrt(1+2*x),'g+-')；xlabel('x');ylabel('y');title('误差曲线');legend('欧拉方法','改进的欧拉方法')将图片保存为error1.fig，图形如下：从该图形来看，改进的欧拉方法与欧拉方法误差接近，欧拉方法误差稍微大些，将x的取值扩宽，n取值增大

18、时，可以发现改进的欧拉方法相比欧拉方法有更高的精度。2.3.2（一阶常微分方程组）初值问题： ，先求出精确解，再用数值方法求解。解析：在MATLAB中输入：x,y=dsolve('Dx=x+4*y-exp(t)','Dy=x+y+2*exp(t)','x(0)=4','y(0)=5/4','t')得：x =4*exp(3*t)+2*exp(-t)-2*exp(t)，y =2*exp(3*t)-exp(-t)+1/4*exp(t)，下面继续用欧拉方法和改进的欧拉方法求数值解。欧拉方法：创建M文件euler2.m，内容

19、如下:function t,x,y=euler2(t0,tfinal,x0,y0,n)f1=inline('x+4*y-exp(t)');f2=inline('x+y+2*exp(t)')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;for i=1:n t(i+1)=t(i)+h; x(i+1)=x(i)+h*feval(f1,t(i),x(i),y(i); y(i+1)=y(i)+h*feval(f2,t(i),x(i),y(i);End在command窗口输入t,x,y=eule

20、r2(0,1,4,5/4,100)，得到t，x,y的值，作图，输入命令如下： plot(t,x,'r*-',t,y,'g*-',t,4*exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,2*exp(3*t)-exp(-t)+1/4*exp(t),'y+-')；legend('x数值解','y数值解','x精确解','y精确解')； xlabel('t'),ylabel('x(y)')，将图形保存为euler2.fi

21、g图形如下：改进的欧拉方法：创建M文件eulerprove2.m文件，内容如下：function t,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline('x+4*y-exp(t)');f2=inline('x+y+2*exp(t)')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;for i=1:n t(i+1)=t(i)+h; x1=x(i)+h*feval(f1,t(i),x(i),y(i); y1=y(i)+h*feval(f2,t(i),x

22、(i),y(i); x2=x(i)+h*feval(f1,t(i+1),x1,y1); y2=y(i)+h*feval(f2,t(i+1),x1,y1); x(i+1)=(x1+x2)/2; y(i+1)=(y1+y2)/2;End在command窗口输入t,x1,y1=euler2(0,1,4,5/4,100)，得到t，x,y的值，作图，输入命令如下： plot(t,x1,'r*-',t,y1,'g*-',t,4*exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,2*exp(3*t)-exp(-t)+1/4*exp(t),

23、'y+-')；legend('x数值解','y数值解','x精确解','y精确解')； xlabel('t'),ylabel('x(y)')，将图形保存为eulerprove2.fig，图形如下：我们再来比较改进的欧拉方法和欧拉方法的误差。输入： plot(t,abs(x-4*exp(3*t)+2*exp(-t)-2*exp(t),'r*-',t,abs(x1-4*exp(3*t)+2*exp(-t)-2*exp(t),'g*-',t,abs(y-4*

24、exp(3*t)+2*exp(-t)-2*exp(t),'b+-',t,abs(y1-2*exp(3*t)-exp(-t)+1/4*exp(t),'y+-')，legend('x欧拉方法','x改进的欧拉方法','y欧拉方法','y改进的欧拉方法')；xlabel('t');ylabel('x(y)')；xlabel('t');ylabel('x(y)')；title('误差曲线')，将图形保存为error2.fig，图

25、形如下：从这个例子我们可以很直观的发现改进的欧拉方法相比欧拉方法具有更高的精度。注意，这里误差曲线的x轴变量为节点参数t，即对应的是固定步长下每一个节点处的误差值。（高阶微分方程）对于二阶常微分方程，求数值解解析：先算出其解析解，在MATLAB中输入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解为：Y=2*sin(t)+cos(t)，前面已经分别给出过欧拉方法和改进的欧拉方法的算例跟误差比较，这里我们就用精度更高的改进欧拉法进行数值求解。改进的欧拉方法：先换元，令，则原方程可以转化为，现在，二阶常微分方程转化为

26、了一个一阶常微分方程组，同的方法，建立M文件eulerprove3.m,内容如下：function t,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline('y');f2=inline('-x')if nargin<5,n=50;endh=(tfinal-t0)/n;t(1)=t0,x(1)=x0;y(1)=y0;for i=1:n t(i+1)=t(i)+h; x1=x(i)+h*feval(f1,y(i); y1=y(i)+h*feval(f2,x(i); x2=x(i)+h*feval(f1,y1); y2=y(i

27、)+h*feval(f2,x1); x(i+1)=(x1+x2)/2; y(i+1)=(y1+y2)/2;end在command窗口输入t,x,y=eulerprove3(0,1,1,2,10)，得到t,x,y的值，其中x就是我们要求的数值解，作图,输入：plot(t,x,'r*-',t,2*sin(t)+cos(t),'b+-')；xlabel('x');ylabel('y）；legend('数值解','精确解')，将图形保存为eulerprove3.fig,图形如下：上面，我们已经通过例子看出，改进的欧

28、拉法相比于欧拉法，在每一个节点处的误差值更下，下面，我们来讨论节点的多少（步长大小）对误差的影响，创建erro3r.m文件，内容如下：function N,Y=error3(n0,nfinal)N(1)=n0,m=fix(nfinal-n0)/4);for i=1:m N(i+1)=N(i)+4; t,x1,y1=eulerprove3(0,1,1,2,N(i); Y(i)=log10(max(x1-2*sin(t)-cos(t); t,x2,y1=eulerprove3(0,1,1,2,N(i); Y(i+1)=log10(max(x2-2*sin(t)-cos(t);end输入N,Y=er

29、ror3(4,100)，返回节点个数值和Y值，Y代表在N个节点时，数值解与精确解差的绝对值的最大值的对数（10为底）。:plot(N,Y,'d-.'),xlabel('N'),ylabel('log_1_0max|error|) title('误差曲线'),将图片保存为error3.fig,图形如下:从这个图的误差曲线，随着插值节点个数增多，误差也是越来越小的，所以，为了得到尽可能准确的解，使用改进的欧拉法时插值节点数应当越大越好。但是，随着节点个数的增多，计算量也会随之增大，所以，在实际应用时，应当结合具体要求灵活选取。注：这三个例子，

30、一个是一阶常微分方程，一个是二元常微分方程组，第三个是二阶常微分方程，都是初值问题求数值解。通过这三个例子可以看出：用欧拉方法和改进的欧拉方法求数值解的一般思路是，通过换元，将高阶常微分方程转换为一阶常微分方程组，再用的方法求解这个方程组。3龙格-库塔法3.1 龙格-库塔法与泰勒展开泰勒（Taylor）展开首先，考虑考虑微分方程 （3-1）引入算子符号， （3-2）于是有 （3-3）设不带括号的及其各阶微商都在处取值，于是有泰勒展开式： （3-4）3.1.2 龙格-库塔（R-K）方法介绍龙格-库塔方法的基本思路是想办法计算在某些点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算

31、公式；再把近似的计算公式和解的泰勒展开式相比较，使得前面的若干项相吻合，从而达到较高的精度，一般的显示R-K方法的形式如下： （3-5）当式（3-5）的局部截断误差达到。龙格库塔方法是应用最广的求解常微分方程初值问题的单步法，下面我们给出几个推导公式。3.2 龙格-库塔法公式与ode函数3.2.1 二阶二段R-K法公式推导由式子（3-5），二阶二段的R-K方法可以写成 （3-6）不带括号的及其各阶微商都在处取值，由泰勒展开式（3-4）及截断误差的定义可知： （3-7）由于上式是二阶二段的R-K法，即必须有，所以 （3-8）式（3-8）有四个未知元，三个方程，故有无穷组解。若令，由式子（3-9）

32、可解得于是这个公式称为预估-校正公式。常用的二阶公式还有中点公式（取），这里不再列举。3.2.2 其他常见公式和ode函数从二阶二段R-K公式的推导可以看出，二段方法每一步需要计算两次函数值，而它也只能达到2阶精度，如果我们提高函数的计算次数，就可以得到精度更高的计算公式，高阶的R-K法的推导与二阶方法的推导完全相同，只是随着阶数的增高计算量会逐渐增大。下面给出几个常用的计算公式。（1） 库塔公式（3阶3段）（2）标准四阶R-K公式上述三阶公式具有三阶精度，四阶公式具有四阶精度，要验证这一点，我们观察如下泰勒展开： (3-9)其中落在平面上连接的线段上。只需要按照式（3.9）将k泰勒展开，代入

33、截断误差的计算公式中即可证明。(详细证明见参考文献2P253)由于R-K法是基于泰勒展开的数值方法，所以，如果解具有较好的光滑性，则能够得到较好的计算精度，反之，则可能四阶R-K法的数值精度还不如低阶的数值方法，这个需要根据具体情况自行选择数值方法。一般而言，R-K法在求解范围较大而精度要求较高时是比较好的方法，四阶的R-K法也可以用作阿达姆斯预估校正法的前几步计算，以保持四阶精度（下一节阿达姆斯预估校正法中也提到了这一点）。在MATLAB中，专门提供了求解微分方程的ode函数，如ode45，ode23，ode113，ode15s，ode23s等等。ode求解器分为变步长和固定步长两种求解模式

34、，其中，ode45是采用R-K法的变步长求解器，和它一样的还有ode23。目前，ode45是使用最多的求解函数，主要用于求解非刚性常微分方程。（如果求解时长时间没有响应，则方程是刚性的，可以换用ode23求解），ode函数的调用方式大都类似，以ode45为例，常用的语法格式有：T,Y = ode45(odefun,tspan,y0) ；T,Y = ode45(odefun,tspan,y0,options) ；sol = ode45(odefun,t0tf,y0.) ，其中odefun是函数句柄，可以是函数文件名或内联函数名，tspan是求解区间 t0 tf 或者一系列散点t0,t1,.,tf

35、，y0 是初始值向量，T 返回列向量的时间点，Y 返回对应T的求解列向量 ，options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 。3.3算例（三阶三段R-K法MATLAB编程解常微分方程组初值问题）计算初值问题： 其中解析：先求其解析解，在MATLAB输入：x,y,z=dsolve('Dx=t*x+Dy-t2','Dy=t*y+3*t','Dz=t*z-Dy+6*t3','x(0)=1','y(0)=2','z(0)=3','t')，解得：x=5/2

36、*exp(1/2*t2)*t2+t-1/2*exp(1/2*t2)*pi(1/2)*2(1/2)*erf(1/2*2(1/2)*t)+exp(1/2*t2);y =-3+5*exp(1/2*t2)；Z=-5/2*exp(1/2*t2)*t2-6*t2-12+15*exp(1/2*t2).下面我们用三阶三段的龙格-库塔公式来求数值解，建立名为rk33.m的M文件，内容如下：function t,x,y,z=rk33(t0,tfinal,x0,y0,z0,n)f1=inline('t*x+t*y+3*t-t2','t','x','y'

37、);f2=inline('t*y+3*t','t','y');f3=inline('t*z-t*y-3*t+6*t3','t','z','y');if nargin<6,n=50;endt(1)=t0,x(1)=1,y(1)=2,z(1)=3;h=(tfinal-t0)/n;for i=1:n t(i+1)=t(i)+h; x1=h*feval(f1,t(i),x(i),y(i); x2=h*feval(f1,t(i)+h/2,x(i)+x1/2,y(i)+h/2); x3=h

38、*feval(f1,t(i)+h,x(i)-x1+2*x2,y(i)+h); x(i+1)=x(i)+(x1+4*x2+x3)/6; y1=h*feval(f2,t(i),y(i); y2=h*feval(f2,t(i)+h/2,y(i)+y1/2); y3=h*feval(f2,t(i)+h,y(i)-y1+2*y2); y(i+1)=y(i)+(y1+4*y2+y3)/6; z1=h*feval(f3,t(i),z(i),y(i); z2=h*feval(f3,t(i)+h/2,z(i)+z1/2,y(i)+h/2); z3=h*feval(f3,t(i)+h,z(i)-z1+2*z2,y

39、(i)+h); z(i+1)=z(i)+(z1+4*z2+z3)/6;End在command窗口输入t,x,y,z=rk33(0,0.75,1,2,3,20)，返回t,x,y,z的值，我们分别将x,y,z的解析解跟数值解作在一张图上，则得到三幅图形，分别保存为rk33x.fig, rk33y.fig, rk33z.fig.作图命令如下：（以y曲线为例，另两幅类似） plot(t,-3+5*exp(1/2*t.2),'r*-',t,y,'g+-')；legend('精确解','数值解')；xlabel('t');yl

40、abel('y') ;title('函数y曲线图')。图形分别如下：下面我们来分析R-K法的误差跟数字节点个数（步长）之间的关系，任取上述三条曲线中的一条进行分析，以y曲线为例，建立errorrk.m文件，内容如下：function N,Y=errorrk(n0,nfinal)N(1)=n0,m=fix(nfinal-n0)/2);for i=1:m N(i+1)=N(i)+2; t,x,y,z=rk33(0,0.75,1,2,3,N(i); Y(i)=log10(max(abs(y+3-5*exp(1/2*t.2); t,x1,y1,z1=rk33(0,0.7

41、5,1,2,3,N(i+1); Y(i+1)=log10(max(abs(y1+3-5*exp(1/2*t.2);End在窗口输入：N,Y=errorrk(5,40);plot(N,Y,'+-.');xlabel('N');ylabel('log_1_0max|error|');title('y的误差曲线')，将图片保存为rk33error.fig，图形如下：从误差图形可以看出，R-K法随着数值节点的增多（步长减小），精度是不断提高的。（用自带ODE45求解高阶常微分方程组）将常微分方程组初值问题转化为一阶常微分方程组并用ode4

42、5求解解析：首先我们尝试用MATLAB求出该方程组的解析解，在command窗口输入：x,y=dsolve('D2x+x*(x2+y2)(-3/2)=0','D2y+y*(x2+y2)(-3/2)=0','x(0)=0.5','y(0)=0.25','Dx(0)=0.75','Dy(0)=1')，运行大概五分钟，提示：Warning: Explicit solution could not be found.x = empty sym ，y =，产生该问题的原因可能在于该非线性方程无法用solve求

43、解，或该方程的解析解不存在。下面，我们用数值方法求解该方程组，调用的是函数库中的ode45函数。首先，我们可以做变量替换，令，则该方程组可以化为一个一阶常微分方程组初值问题：现在，可以直接调用ode45函数进行数值计算。先建立一个名为fun45.m的M文件来定义方程组，内容如下：function y=fun45(t,x)y(1)=-x(3)*(x(3)2+x(4)2)(-3/2);y(2)=-x(4)*(x(3)2+x(4)2)(-3/2);y(3)=x(1);y(4)=x(2);y=y(:)在command窗口输入： x0=0.75,1,0.5,0.25;t,x=ode45('fun

44、45',0,1,x0),输出x的值，x是一个四列的矩阵，四列分别代表a,b,x,y在区间0,1上的值，可以读出 ：x(1)=0.5499, y(1)=0.7407，或者作图：plot(t,x)；legend('a','b','x','y')，xlabel('t')将图片保存为fun45.fig.图形如下：通过Data Cursor一样可以读出x(1)=0.5499, y(1)=0.7407。注：这两个例子一个是自己编的三阶三段的R-K法的一般固定步长程序，一个是系统自带的ode45函数求解，事实上MATLA

45、B提供了各种R-K法的算法函数，不仅使用方便，而且步长能够根据精度需要自动调整，这一节，我们重点需要了解ode函数的使用方法。4 阿达姆斯预估校正法4.1 阿达姆斯（Adams）公式 在第二节讲述改进的欧拉方法时，我们先通过积分得到了表达式（2-7），移项后如下: （4-1）在改进的欧拉方法中，我们是使用梯形公式来计算上式中的积分项，事实上，采用不同的数值方法计算它就能导出不同的计算方法，普遍采用的方法是用的插值多项式来代替，用来代替，这种做法也是线性多步法的一般思路，计算公式如下： （4-2）取等距节点，设已经求得这些点上的数值解，从而有，其中。 令，构造r次牛顿向后插值多项式如下： （4-

46、3）整理可以得到递推式如下（具体推导过程参见参考文献2）： (4-4)其中。下表给出了的值 012313/2-1/223/12-16/125/1255/24-59/2437/24-9/24当时，(4-4)就是欧拉公式。当时，(4-4)可以写为： (4-5)式(4-5)就称为阿达姆斯-巴士福德公式。用牛顿插值多项式的余项估算可证明该公式是一个4阶方法（证明略，详见参考文献2）若取个节点为，对应的函数为，则类似的可推得： （4-6）其中下表给出了的值. 012311/21/25/128/12-1/129/2419/24-5/241/24当时，式（4-6）为 （4-7）称为梯形公式（即改进的欧拉方法

47、）当时，式（4-6）为 （4-8）用牛顿插值多项式的余项可以上述公式是阶方法（证明略，详见参考文献2），即改进的欧拉方法是二阶方法（这与第二节证明相符），（4-8）是一个4阶方法，这个式子称为阿达姆斯-莫尔德公式。 从上述公式的形式可以看出，要求，欧拉方法与改进的欧拉方法只需要知道即可，这种公式称为单步法，第三节的R-K方法也是一个单步法，而用阿达姆斯公式计算，需要知道，这种方法称为多步法。4.2 预估校正方法 所谓预估校正方法，在数值算法中，主要就是用于隐式计算公式的计算的，可以看出式（4-8）就是一个隐式计算公式。回顾第二节的改进欧拉方法。为了能够使改进的欧拉方法在迭代时不需要每一步都解方

48、程，我们用欧拉方法先求一个预测值，再用这个预测值来计算，即：这其实就是一个预估校正算法。同样的，我们可以用公式（4-5）来预估，再用公式（4-8）来校正，则可以得到一个一个预估校正方法如下： （4-9）其中，（4-9）也是一个4阶算法。同4阶R-K算法相比较，阿达姆斯预估校正算法精度相同，但是迭代公式更为简单，计算量减少较多，但是，阿达姆斯预估算法无法自动进行，其开始的几个初值需要通过其他的方法计算得出（如泰勒展开），一般来说，多步法的前面几步通常是通过单步法来计算的，为了保证阿达姆斯预估算法的精度，我们可以用龙格-库塔法来计算前面几步。4.3算例4.3.1 （常微分方程初值问题）用Adams

49、预估校正算法在上求解初值问题： 将结果与ode45及改进的欧拉法求解结果想比较。解析：先求出解析解，输入：dsolve('Dy=x*y2','y(0)=1','x')，解得：y =-2/(x2-2)。下面用阿达姆斯预估校正法求数值解。创建M文件adams1.m，内容如下：function x,y=adams1(fun,x0,xfinal,y0,n)h=(xfinal-x0)/n,x(1)=x0,y(1)=y0;for i=1:3 x(i+1)=x(i)+h; y1=h*feval(fun,x(i),y(i); y2=h*feval(fun,x(i

50、)+h/2,y(i)+y1/2); y3=h*feval(fun,x(i)+h/2,y(i)+y2/2); y4=h*feval(fun,x(i)+h,y(i)+y3); y(i+1)=y(i)+1/6*(y1+2*y2+2*y3+y4);endfor i=4:n x(i+1)=x(i)+h; y5=y(i)+h/24*(55*feval(fun,x(i),y(i)-59*feval(fun,x(i-1),y(i-1)+37*feval(fun,x(i-2),y(i-2)-9*feval(fun,x(i-3),y(i-3); y(i+1)=y(i)+h/24*(9*feval(fun,x(i+

51、1),y5)+19*feval(fun,x(i),y(i)-5*feval(fun,x(i-1),y(i-1)+feval(fun,x(i-2),y(i-2);End（前几步使用的是4阶RK标准公示，下同）定义函数adams1f如下：function y=f1(x,y);y=x*y2.在command窗口输入：x,y=adams1('adams1f',0,1,1,40)，返回x,y和函数的值，作图，命令如下： plot(x,y,'r*-',x,-2./(x.2-2),'g+-');legend('数值解','精确解'

52、;);title('y=x*y2')，得到图形如下，保存为adams1a.fig调用ode45函数，输入x1,y1=ode45(''adams1f',0,1,1),得到用RK法的数值解（40个节点），现在来比较二者之间的误差，x2,y2=eulerprove1('adams1f',0,1,1,40)输入： subplot(311);plot(x,abs(y+2./(x.2-2);title('阿达姆斯法误差曲线');xlabel('x');ylabel('|error|')；subplot(

53、312);plot(x1,abs(y1+2./(x1.2-2);title('R-K（ode45）误差曲线');xlabel('x');ylabel('|error|');subplot(313);plot(x2,abs(y2+2./(x2.2-2);title('改进的欧拉法误差曲线');xlabel('x');ylabel('|error|')得到如下误差曲线图（都是40个数值节点），保存为adams1b.fig.从误差曲线可以直观的看出，阿达姆斯预估校正法与4阶RK法精度相当，相比改进的欧拉法

54、具有明显的精度提升，事实上，阿达姆斯法与ode45都是四阶方法，改进的欧拉法是二阶方法，在相同数目的取值节点情况下，前两者精度更高，这个图形也证明了这一观点。（常微分方程组初值问题）初值问题： ，先求出精确解，再用数值方法求解。解析：在MATLAB中输入：y,z=dsolve('Dy=y+4*z-exp(x)','Dz=y+z+2*exp(x)','y(0)=2','z(0)=3','x')得：y =-3/4*exp(-x)+19/4*exp(3*x)-2*exp(x)；z =3/8*exp(-x)+19/8*ex

55、p(3*x)+1/4*exp(x)。下面我们用阿达姆斯预估校正法来求其数值解，建立文件adams2.m，内容如下：function x,y,z=adams2(x0,xfinal,y0,z0,n)f1=inline('y+4*z-exp(x)','x','z','y');f2=inline('y+z+2*exp(x)','x','y','z');x(1)=x0,y(1)=y0,z(1)=z0;h=(xfinal-x0)/n;for i=1:3 x(i+1)=x(i)+h

56、; y1=h*feval(f1,x(i),z(i),y(i); y2=h*feval(f1,x(i)+h/2,z(i)+h/2,y(i)+y1/2); y3=h*feval(f1,x(i)+h/2,z(i)+h/2,y(i)+y2/2); y4=h*feval(f1,x(i)+h,z(i)+h,y(i)+y3); y(i+1)=y(i)+1/6*(y1+2*y2+2*y3+y4); z1=h*feval(f2,x(i),y(i),z(i); z2=h*feval(f2,x(i)+h/2,y(i)+h/2,z(i)+z1/2); z3=h*feval(f2,x(i)+h/2,y(i)+h/2,z(i)+z2/2); z4=h*feval(f2,x(i)+h,y(i)+h,z(i)+z3); z(i+1)=z(i)+1/6*(z1+2*z2+2*z3+z4);endfor i=4:n x(i+1)=x(i)+h; y5=y(i)+h/24*(55*feval(f1,x(i),z(i),y(i)-59*feval(f1,x(i-1),z(i-1),y(i-1)+37*feval(f1,x(i-2),z