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

装订线毕业设计（论文）报告纸常微分方程的几种数值解法数学与应用数学 肖振华 指导教师 张秀艳【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、RungeKutta方法、Adams预估校正法以及勒让德谱方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。【关键词】 常微分方程 数值解法 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 equations. 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 solution, 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 differential 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 欧拉法和改进的欧拉法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 常微分方程概述方程是一个在数学中非常熟悉的名词，在初等数学里，我们将我们要研究的问题作为一个或几个未知量，通过观察事物的规律，得出这些未知量与已知量之间的等式关系，这样就得到了一个简单的方程或方程组当然，这只是一个很浅显粗略的定义。在数学上，物质的运动和变化规律是通过函数关系来表示的，在一些复杂的现象中，比如火箭的运动轨迹，我们要求的未知量就变成了满足特定条件的一个或一些未知函数（不再是简单的数值）。有时候，我们需要用到导数很微分的知识，即这些未知函数的导数与自变量满足某种关系，这种方程，我们称之为微分方程，即可以这样定义，凡是含有参数，未知函数和未知函数导数（或微分）的方程称之为微分方程。未知函数是一元函数的微分方程称为常微分方程（ODE），未知函数是多元函数的微分方程称为偏微分方程（PDE），我们要研究的，就是常微分方程。1.2 常微分方程解与数值解法微分方程的解，就是找出一个代入方程使之成为恒等式的函数。.若微分方程的解中含有任意常数的个数与方程的阶数相同，且任意常数之间不能合并，则称此解为该方程的通解(或一般解).当通解中的各任意常数都取特定值时所得到的解，称为方程的特解。在实际问题中，这些函数往往还满足一些特定条件，这个称之为定解问题，如常微分方程的初值问题和边值问题，满足条件的解就是特解。在理想的状况下，我们解这个方程如果能够得出通解的解析表达式 ，那么代入数值很容易得到问题所要求的特解，但事实上，很多常微分方程得不到通解的解析表达式，或者表达式过于复杂。而且，常微分方程的特解是否存在，有几个解，牵涉到解的存在唯一性定理。在实际应用中，我们对解的精度要求往往并没有那么严格。在这种情况下，数值解应运而生，通俗的讲，数值解就是所求函数在特定点的函数值，通过函数在特定点的取值，我们可以反过来观察这个函数的某些性质，所以，求常微分方程的数值解是很有意义的，事实上，解常微分方程很大一部分就是解数值解。 下面我们要讨论几种常见的常微分方程数值解法，目的是展示常微分方程数值解法的一般过程。以下章节，由于本人水平有限，只对数值解法的思路和过程进行讨论，而对于解的存在唯一性，方法精度的具体证明，不做过多讨论。一般而言，构造数值解法都是基于两个想法：一种是基于数值积分的构造方法，一种是基于泰勒展开的构造法。下面我们要讨论的是几种常见的数值解法，如欧拉法、改进的欧拉法、龙格库塔法、阿达姆斯预估校正法以及基于勒让德多项式的谱方法。这些方法既有基于数值积分的构造方法，如欧拉法，也有基于泰勒展开的构造法，如龙格库塔法。下面的章节，我们会一一介绍这些方法的思路，并通过一些简单的算例，让大家了解这些解法的一般过程。为了让大家看到各种数值解法的精确性和有效性，我们可能会尽量选取一些通解存在的算例来求数值解，这样做并不意味着我们实际中求数值解时要先求出精确解，而只是出于对方法对比的需要，望不致引起混淆。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中输入：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，内容如下: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 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; 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=euler2(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.fig图形如下：改进的欧拉方法：创建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 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,t(i),x(i),y(i); y1=y(i)+h*feval(f2,t(i),x(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),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*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，图形如下：从这个例子我们可以很直观的发现改进的欧拉方法相比欧拉方法具有更高的精度。注意，这里误差曲线的x轴变量为节点参数t，即对应的是固定步长下每一个节点处的误差值。2.3.3（高阶微分方程）对于二阶常微分方程，求数值解解析：先算出其解析解，在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,图形如下:从这个图的误差曲线，随着插值节点个数增多，误差也是越来越小的，所以，为了得到尽可能准确的解，使用改进的欧拉法时插值节点数应当越大越好。但是，随着节点个数的增多，计算量也会随之增大，所以，在实际应用时，应当结合具体要求灵活选取。注：这三个例子，一个是一阶常微分方程，一个是二元常微分方程组，第三个是二阶常微分方程，都是初值问题求数值解。通过这三个例子可以看出：用欧拉方法和改进的欧拉方法求数值解的一般思路是，通过换元，将高阶常微分方程转换为一阶常微分方程组，再用2.3.2的方法求解这个方程组。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算例3.3.1（三阶三段R-K法MATLAB编程解常微分方程组初值问题）计算初值问题： 其中解析：先求其解析解，在MATLAB输入：x,y,z=dsolve(Dx=t*x+Dy-t2,Dy=t*y+3*t,Dz=