微分方程数值解法

微分方程数值解法【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、RungeKutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。【关键词】 常微分方程 数值解法 MATLAB 误差分析引言 在我国高校，微分方程数值解法作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程近四十年来，微分方程数值解法不论在理论上还是在方法上都获得了很大的发展同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。2 欧拉法和改进的欧拉法2.1 欧拉法 2.1.1 欧拉法介绍首先，我们考虑如下的一阶常微分方程初值问题 （2-1）事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将看做向量，（2-1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2-1）中的导数项用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。设在中取等距节点，因为在节点点上，由（2-1）可得： ， （2-2） 又由差商的定义可得： （2-3） 所以有 （2-4）用的近似值代入（2-4），则有计算的欧拉公式 （2-5）2.1.2欧拉法误差分析从欧拉公式中可以看出，右端的都是近似的，所以用它计算出来的会有累计误差，累计误差比较复杂，为简化分析，我们考虑局部截断误差，即认为是精确的前提下来估计，记为，泰勒展开有 （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）是一个隐式计算公式，事实上，改进的欧拉方法是用欧拉方法先求一个预测值，再用这个预测值来计算，即： （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,y0,n)if nargin5,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+-)； 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 nargin5,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(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（高阶微分方程）对于二阶常微分方程，求数值解解析：先算出其解析解，在MATLAB中输入：y=dsolve(D2x=-x,x(0)=1,Dx(0)=2)得到解为：Y=2*sin(t)+cos(t)，前面已经分别给出过欧拉方法和改进的欧拉方法的算例跟误差比较，这里我们就用精度更高的改进欧拉法进行数值求解。改进的欧拉方法：先换元，令，则原方程可以转化为，现在，二阶常微分方程转化为了一个一阶常微分方程组，同2.3.2的方法，建立M文件eulerprove3.m,内容如下：function t,x,y=eulerprove2(t0,tfinal,x0,y0,n)f1=inline(y);f2=inline(-x)if nargin5,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: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=error3(4,100)，返回节点个数值和Y值，Y代表在N个节点时，数值解与精确解差的绝对值的最大值的对数（10为底）。:plot(N,Y,d-.),xlabel(N),ylabel(log_1_0max|error|) title(误差曲线),将图片保存为error3.fig,图形如下:从这个图的误差曲线，随着插值节点个数增多，误差也是越来越小的，所以，为了得到尽可能准确的解，使用改进的欧拉法时插值节点数应当越大越好。但是，随着节点个数的增多，计算量也会随之增大，所以，在实际应用时，应当结合具体要求灵活选取。3龙格-库塔法3.1 龙格-库塔法与泰勒展开3.1.1泰勒（Taylor）展开首先，考虑考虑微分方程 （3-1）引入算子符号， （3-2）于是有 （3-3）设不带括号的及其各阶微商都在处取值，于是有泰勒展开式： （3-4）3.1.2 龙格-库塔（R-K）方法介绍龙格-库塔方法的基本思路是想办法计算在某些点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式和解的泰勒展开式相比较，使得前面的若干项相吻合，从而达到较高的精度，一般的显示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）可解得于是 这个公式称为预估-校正公式。常用的二阶公式还有中点公式（取），这里不再列举。3.2.2 其他常见公式和ode函数从二阶二段R-K公式的推导可以看出，二段方法每一步需要计算两次函数值，而它也只能达到2阶精度，如果我们提高函数的计算次数，就可以得到精度更高的计算公式，高阶的R-K法的推导与二阶方法的推导完全相同，只是随着阶数的增高计算量会逐渐增大。下面给出几个常用的计算公式。（1） 库塔公式（3阶3段） （2）标准四阶R-K公式 上述三阶公式具有三阶精度，四阶公式具有四阶精度，要验证这一点，我们观察如下泰勒展开： (3-9)其中落在平面上连接的线段上。只需要按照式（3.9）将k泰勒展开，代入截断误差的计算公式中即可证明。(详细证明见参考文献2P253)由于R-K法是基于泰勒展开的数值方法，所以，如果解具有较好的光滑性，则能够得到较好的计算精度，反之，则可能四阶R-K法的数值精度还不如低阶的数值方法，这个需要根据具体情况自行选择数值方法。一般而言，R-K法在求解范围较大而精度要求较高时是比较好的方法，四阶的R-K法也可以用作阿达姆斯预估校正法的前几步计算，以保持四阶精度（下一节阿达姆斯预估校正法中也提到了这一点）。在MATLAB中，专门提供了求解微分方程的ode函数，如ode45，ode23，ode113，ode15s，ode23s等等。ode求解器分为变步长和固定步长两种求解模式，其中，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，y0 是初始值向量，T 返回列向量的时间点，Y 返回对应T的求解列向量 ，options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 。3.3算例用自带ODE45求解高阶常微分方程组将常微分方程组初值问题转化为一阶常微分方程组并用ode45求解解析：首先我们尝试用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求解，或该方程的解析解不存在。下面，我们用数值方法求解该方程组，调用的是函数库中的ode4