《微分方程数值解法》【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。
其中，常微分方程求解是微分方程的重要基础内容。
但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。
，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。
同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。
【关键词】 常微分方程 数值解法 MA TLAB 误差分析引言在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。
2 欧拉法和改进的欧拉法2.1 欧拉法2.1.1 欧拉法介绍首先，我们考虑如下的一阶常微分方程初值问题 ⎩⎨⎧==00)(),('y x y y x f y（2--1）事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。
欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。
设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--1）可得：))(()(',n n n x y x f x y =，（2--2）又由差商的定义可得：hx y x y x y n n n )()(('1-≈+）（2--3） 所以有 ))(,()()(1n n n n x y x hf x y x y +≈+ （2--4）用)(k x y 的近似值k y )1,(+=n n k 代入（2--4），则有计算1+n y 的欧拉公式))(,(1n n n n x y x hf y y +=+ （2--5）2.1.2欧拉法误差分析从欧拉公式中可以看出，右端的n y 都是近似的，所以用它计算出来的1+n y 会有累计误差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为n y 是精确的前提下来估计11)(++-n n y x y ，记为1+n ε，泰勒展开有)()()(''2)(')()(1321++<<+++=n n n n n x x h O y h x hy x y x y ξξ（2--6）联立（2--5），（2--6）即得1+n ε=)(''22ξy h +)(3h O =)(2h O ，根据数值算法精度的定义，如果一个数值方法的局部截断误差1+n ε=)(1+p h O 则称这个算法具有P 阶精度，所以，欧拉方法具有一阶精度或者称欧拉方法为一阶方法。
2.2 改进的欧拉方法2.2.1 改进的欧拉法介绍用数值积分离散化问题（1），两边做积分有： dx x y x f x y x y n nx x n n ⎰+=-+1))(,()()(1（2--7）对右端积分使用梯形积分公式可得：[]))(,())(,(2))(,(111+++≈⎰+n n n n x x x y x f x y x f hdx x y x f n n（2--8）同欧拉方法，用)(k x y 的近似值k y )1,(+=n n k 代入（2--7），联立（2--8）得到改进的欧拉方法计算公式： ),(),((2111+++++=n n n n n n y x f y x f hy y （2--9）一般来说，如果求解1+n y 时，算法右端不包含1+n y ，称为显性计算公式，如果包含，则求解时还需要解方程，这种称为隐式计算公式。
显然公式（2--9）是一个隐式计算公式，事实上，改进的欧拉方法是用欧拉方法先求一个预测值1+n y ，再用这个预测值来计算1+n y ，即：⎪⎩⎪⎨⎧++=+=++++),(),((2),(1111n n n n n n n n n n y x f y x f hy y y x hf y y （2--10）2.2.2改进的欧拉法误差分析同欧拉法误差分析类似，用泰勒展开容易知道改进的欧拉方法具有二阶精度，证明略。
2.3 算例2.3.1（一阶常微分方程）求解初值问题⎪⎩⎪⎨⎧∈-==]1,0[21)0('x yx y y y解析：在MATLAB 中求解这个方程y=dsolve('Dy=y-2*x/y','y(0)=1','x') 得y =(2*x+1)^(1/2)它的解析解为x y 21+=，下面我们分别用欧拉方法和改进的欧拉方法来求其数值解。
欧拉方法：创建M 文件euler1.m,内容如下： function [x,y]=euler1(fun,x0,xfinal,y0,n) if nargin<5,n=50; endh=(xfinal-x0)/n; x(1)=x0;y(1)=y0;for i=1:nx(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+--')；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:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,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('数值解','精确解')，将图片保存为eulerprove1.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 取值增大时，可以发现改进的欧拉方法相比欧拉方法有更高的精度。
2.3.2（高阶微分方程）对于二阶常微分方程）（]1,0[1)0(2)0('''∈⎩⎨⎧==-=x x x xx ，求数值解 解析：先算出其解析解，在MATLAB 中输入：y=dsolve('D2x=-x','x(0)=1','Dx(0)=2')得到解为：Y=2*sin(t)+cos(t)，前面已经分别给出过欧拉方法和改进的欧拉方法的算例跟误差比较，这里我们就用精度更高的改进欧拉法进行数值求解。
改进的欧拉方法：先换元，令y x ='，则原方程可以转化为])1,0[(2)0(1)0(''∈⎪⎩⎪⎨⎧===-=x y x yx xy ，现在，二阶常微分方程转化为了一个一阶常微分方程组，同2.3.2的方法，建立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)+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,图形如下：上面，我们已经通过例子看出，改进的欧拉法相比于欧拉法，在每一个节点处的误差值更下，下面，我们来讨论节点的多少（步长大小）对误差的影响，创建erro3r.m文件，内容如下：function [N,Y]=error3(n0,nfinal)N(1)=n0,m=fix((nfinal-n0)/4);for i=1:mN(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]=error3(4,100)，返回节点个数值和Y值，Y代表在N个节点时，数值解与精确解差的绝对值的最大值的对数（10为底）。