一阶微分方程的MATLAB数值解法.doc

精选优质文档-倾情为你奉上 毕业论文题 目 一阶微分方程的matlab数值解法 学 院 数学科学学院 专 业 数学与应用数学 班 级 数学1303班 学 生 邹健峰 学 号 7 指导教师 邢顺来讲师 二一七年五月七日专心-专注-专业摘 要微分方程的数值解法是近现代数学家和科学家们研究的热点，微分方程的MATLAB数值解法可以帮助我们解决现实生活中的许多数学问题，提高计算机帮助我们解决问题的效率。本文主要讨论研究一阶微分方程的MATLAB数值解法中的三种Euler法和三种Runge-Kutta法，介绍以上六种方法的主要内容，简单介绍了MATLAB的相关内容和微分方程的发展简史。通过具体的微分方程来研究以上算法的编程实现，通过MATBLAB求解具体的一阶微分方程来探究以上方法的优缺点，通过图表来分析得出结论：改进Euler法和四阶Runge-Kutta法的阶的精度较高，具有较好的算法稳定性。关键词：一阶微分方程；数值解法；MATLAB；Euler法；Runge-Kutta法ABSTRACTThe numerical solution of differential equations is a hot spot for modern mathematicians and scientists, the MATLAB numerical solution of the differential equation can help us solve many mathematical problems in real life, improve the efficiency of the computer to help us solve the problem. This paper mainly discusses three kinds of Euler methods and three Runge-Kutta methods in MATLAB numerical solution of first order differential equation, introduce the main contents of the above six methods, a brief introduction to the development of MATLAB related content and differential equations of a brief history. The programming of the above algorithm is studied by the concrete differential equation, he advantages and disadvantages of the above methods are explored through MATBLAB solving the specific first-order differential equation, through the chart to analyze the conclusions: The improved Euler method and the fourth order Runge-Kutta method have higher accuracy，has good algorithm stability.Key words：First Order Differential Equation; Numerical Solution; MATLAB; Euler Method; Runge-Kutta Method目 录摘要IABSTRACTII1 前言11.1 引言11.2常微分方程的发展简史21.3 国内外研究现状21.4 研究主要内容及研究方案31.4.1 研究的主要内容31.4.2 研究方案31.5 研究难点32 预备知识42.1 显式Euler法42.2 隐式Euler法52.3 改进Euler法52.4 二阶Runge-Kutta法62.5 三阶Runge-Kutta法72.6 四阶Runge-Kutta法72.7 单步法的收敛性和稳定性93 微分方程数值解法的编程实现113.1 常微分方程的符号解法113.2 显式Euler法的编程实现123.3 隐式Euler法的编程实现143.4 改进Euler法的编程实现173.5 二阶Runge-Kutta法的编程实现203.6 三阶Runge-Kutta法的编程实现223.7 四阶Runge-Kutta法的编程实现24结论29参考文献30致谢31附录321 前言1.1 引言James Bernoulli在1676年写给Newton的信中首次谈到了微分方程。大约到了十八世纪中期，微分方程才成为一门独立的学科，微分方程是研究和揭示自然规律不可或缺的重要工具。众所周知，在1846年数学家和天文学家共同发现的海王星就是微分方程的功劳。还有在1991年，科学考察队在Alps山脉发现了一个冰封保存完好的人类，然后通过仪器测得碳元素衰变的速率，建立微分方程数学模型，得出了这个人是生活在距今5000多年前的时代。Newton、Leibniz、Bernoulli family、Lagrange、Laplace等著名数学家们都对微分方程的发展做出了巨大的贡献。微分方程的定义：如果知道自变量、未知函数以及函数的导数（或微分）组成的关系式我们就称为微分方程。常微分方程的定义：通过求解微分方程得出未知函数，并且自变量只有一个的微分方程我们就叫做常微分方程，自变量个数为两个或者两个以上的微分方程称为偏微分方程1。我们把含有n个互不相等的任意常数的解称为n阶常微分方程的通解。微分方程的特解是指满足微分方程初值条件的解1。我们在生活中会经常用到一些常微分方程模型，例如英国著名人口统计学家Malthus提出的著名的Malthus人口模型；荷兰著名生物学家Verhulst提出的logistic人口增长模型；传染病模型SI模型、SIS模型（无免疫性）、SIR（强免疫性）模型；两生物种群生态模型；气象学家Lorenz提出的Lorenz方程用来预测天气变化；化学动力学模型2。可见微分方程与我们的生活有着密不可分的关系。MATLAB是分别取了Matrix和Laboratory两个单词的前三个字母组成的,是矩阵实验室的名称缩写。MATLAB是来自美国大学的计算机老师Cleve Moler创造的，并且由Cleve Moler和John Lillte把其推向了市场，并得到了广大用户的喜爱。MATLAB是数学科学和工程科学的高级计算机语言，用户可以用数学语言来编写程序，编程逻辑非常类似人们日常书写数学公式的逻辑，具有简单性、广泛性、开放性、效率高、价格低等特点。并且MATLAB工具箱储存的MATLAB函数库里面有大量的函数程序，可以解决许多的数学和工程科学问题。从2010版开始MATLAB的功能越来越多，编译器可以把M文件转换为C语言文件或者其他几种常用的计算机文件，并给用户提供源程序代码进行自我开发，结合Math Works公司提供的数据库，用户可以利用MATLAB来开发出功能强大的程序3,4,5。微分方程的MATLAB数值解法可以帮助我们解决许多现实生活中的数学问题，提高计算机帮助我们解决问题的效率，如今微分方程的MATLAB数值解法已经发展得比较成熟，有许多关于数学领域的MATLAB书籍可供我们选择阅读，很好的解决了我们的学习需求和使用问题。MATLAB应用领域也非常广泛，比如在数学建模和数学实验领域，我们可以通过直接调用MATLAB里面的微分方程函数文件，在操作界面输入初始数据，还把抽象复杂的数学问题变成简单形象的二维或者三维图形，来帮助我们分析问题，找出解题方法，从而解决我们生活中负责的数学问题，体现了数学源自生活又联系生活的主要特征。除此之外，在物理学、工程科学等科学领域中MATLAB也具有重要地位 6,7。1.2 常微分方程的发展简史常微分方程的发展可以简单的分为以下几个阶段：“求通解”阶段；“求定解”阶段；“求所有解”阶段；“求特殊解”阶段。在常微分方程的“求通解”阶段，人们希望能够用初等函数来表示解，Newton就希望能够用分离变量法来求解一阶微分方程问题，但Euler却尝试用积分因子来求解一阶微分方程问题。而Bernoulli在研究微分方程的过程中提出了著名的Bernoulli方程，Riccati也在研究过程中提出了Riccati方程，他们的发现对微分方程未来的发展打下了基础，起到了积极的作用。直到Liouville证明了Riccati方程没有一般的初等解，还有Cauchy提出了微分方程的初值问题，常微分方程进入了“求定解”阶段。直到19世纪末，由于天体力学的研究需求，常微分方程从“求定解”阶段进入了“求所有解”的阶段。首先是Poincaré提出了定性理论，接着是Hilbert提出了关于极限环个数的问题，促进了微分方程定性理论的发展，紧接着Lyapunov提出了运动稳定性理论，广泛应用在物理和工程学中。而在动力系统方面，Birkhoff、Arnold、Smale都做出了贡献。到了20世纪70年代计算机技术的发展带动了常微分方程的发展，从“求所有解”阶段进入到了“求特定解”阶段，发现了一些特定解和方程，比如混沌（解）、孤立子、奇异吸引子、分形。关于常微分方程的研究还在继续，如今研究的领域更加广泛，更多的微分方程模型被提出和应用起来了1,8,9,10。 1.3 国内外研究现状关于微分方程数值解法的研究，在没有MATLAB以及类似的计算机软件问世之前，数学家和科学家们只能研究和解决一些低阶的微分方程问题，并且有些方程的解题过程过于繁琐和复杂。但是MATLAB的出现促进了数学科学研究的发展，对数学和相关领域的发展做出了巨大贡献。微分方程的MATLAB数值解法也迅速发展起来，国内外现在主要研究通过MATLAB等计算机软件来研究：一些高阶微分方程的解的有关问题；一些复杂的微分方程是否存在定解问题；利用MATLAB来计算一些著名的数学模型问题，比如Volterra被捕食 捕食模型、竞争模型、共生模型、商品供需模型等等；研究混沌（解）、孤立子。在以上几方面的研究上都取得了不小的进步，使常微分方程的体系更加的完善了。在国内，2009年中国数学会成功举办了第四届全国青年常微分方程学术会议，许多国内知名学者和学术研究者都到场参加学术交流，这一项活动促进了我国常微分方程领域的发展，为青年学者提供了宝贵经验，这次会议也受到了国外数学爱好者的高度关注。在书籍方面，我国也发行了许多关于微分方程数值分析和微分方程数值解法方面的书籍，都有详细介绍微分方程关于MATLAB的数值解法，比如林成森编著的第二版数值计算方法里面就有详细介绍常微分方程初值问题的数值解法，包括有：离散变量法和离散误差法；单步法；单步法的相容性、收敛性和稳定性；多步法；差分方程的简介；线性步法的相容性、收敛性和数值稳定性；常微分方程组和高阶微分方程的数值解法等几个方面，成为了21世纪高等院校的教材。周品、何正风等编著的MATLAB数值分析作为了理工科各专业本科生、研究生以及应用MATLAB相关科技人员学习MATLAB数值分析、建模、仿真的教材或者参考书。现如今，由于信息时代的飞速发展，数值计算方面的计算机软件越来越多，功能也越来越全面，所以微分方程的数值解法受到了更多人的重视，发展前景变得更好了，这对于学习微分方程的数值解法的人来说是非常好的消息。1.4 研究的主要内容与方法1.4.1 研究的主要内容本文针对一阶常微分方程的数值解法中的MATLAB数值解法进行算法介绍和编程实现，利用MATLAB对一阶微分方程中的Euler法、隐式Euler法、改进Euler法以及Runge-Kutta法进行精度分析并比较，揭示以上几种方法的优缺点，编程实现以上方法的数值解法，并通过具体事例来验证程序的完整性、稳定性、以及不足。1.4.2 研究方法一阶微分方程的MATLAB数值解法的编程实现需要通过Euler法、隐式Euler法、改进Euler法以及二到四阶Runge-Kutta法等方法的算法逻辑，针对不同一阶微分方程来编写MATLAB程序，然后再反复调试程序，找出其中的不足并改进，然后保存程序并做记录，以便下次使用时直接调用。1.5 研究难点（1）微分方程的模型多样，研究方向多，不同的模型有不同的数值分析和数值解法；（2）微分方程的MATLAB数值解法需要符合实际问题，限制条件比较多；（3）微分方程的类型多样，MATLAB的函数库里只有部分微分方程的调用函数，不能满足研究的需求，需要自定义调用函数；（4）微分方程的MATLAB数值解法需要反复调试，得出满意结果花费时间较长。2 预备知识2.1 显式Euler法显示Euler法也被人们简称为Euler法，是解决微分方程初值问题最简单的数值方法。初值问题的解代表通过点的一条称为微分方程的积分曲线。积分曲线上每一点的切线斜率等于函数在这点的值。Euler方法的求解过程是从初始点出发,作积分曲线在点的切线（其斜率为），与直线相交与点，得到作为的近似值，过点，以为斜率的切线方程为当时，得这样就获得了点的坐标3,11。不断重复上述过程，就可以得到一系列的点。对已求得的点以为斜率做直线，当时，便可得到方程取，这样从逐一算出所对应的数值解，同样也获得了一条曲线作为原始曲线的近似曲线，因此又称Euler法为折线法。在计算时，一般取（常数），则显式Euler公式为 . （2.1）2.2 隐式Euler法假设我们用近似等式来替代显式Euler公式右端的积分项，就得到再用替代，替代，我们就可以得到这就是隐式Euler公式。后退的Euler法也被称为隐式Euler法。而我们经常用到的梯形公式。对方程的两端在区间上积分，得要通过这个积分关系式得的近似值，只要近似的算出其中的积分项。若使用梯形方法计算积分项再代入中，得将式子中的分别用代替，有. （2.2）这种差分格式就称为梯形格式，梯形格式上是显示Euler公式与隐式Euler公式的算术平均。梯形方法是隐式单步法，可以用迭代法求解，与后退的Euler法一样，仍用Euler法进行初值迭代，所以梯形公式的迭代公式为 （2.3）2.3 改进Euler法显示Euler公式计算工作量小，但精度不高，若对计算结果的精度要求高，就必须使用就是精度更高的计算方法。微分方程初值问题的梯形公式相比显示Euler公式，精度提高了，但是用迭代的算法来求解计算工作量大。所以可以改进思路：综合表示Euler公式和梯形公式便可以得到改进的Euler公式。先用显示Euler公式求出一个初步的近似值，记为，称之为预测值，它的精度不高，再用梯形公式对它进行一次矫正，也就是迭代一次，求出，称为矫正值12,13。这种预测矫正的方法我们就叫做改进的Euler法. （2.4）它也可以表示为嵌套模式.或者表示为下列的平均化形式 （2.5） 2.4 二阶Runge-Kutta法二阶Runge-Kutta的一般公式为 （2.6）局部截断误差为.其中加权系数和选点调节系数都是特定系数。本质是用差商代替导数的情况下，把在区间上的若干点的值加权平均，再作为斜率，作过点的直线交直线，求得下一点。适当选择是提高精度的方法14,15,16。我们常用的二阶Runge-Kutta公式有 中点方法 （2.7）和Heun方法 （2.8）2.5 三阶Runge-Kutta法三阶Runge-Kutta一般公式为 （2.9）局部截断误差.常用的三阶Runge-Kutta公式有 （2.10）和 （2.11） 2.6 四阶Runge-Kutta法继续以上的演算过程，可以推导出各种四阶Runge-Kutta公式，下面介绍一个被被称为经典Runge-Kutta公式，其公式的格式为 （2.12）其截断误差为. 四阶Runge-Kutta法要求所求的解具有比较好的光滑性，反之，使用四阶Runge-Kutta法求得的解精度比改进的Euler法求得的解精度低。2.7 单步法的收敛性和稳定性定义1 假设微分方程的右端函数在带形区域中连续，并且关于满足Lipschitz条件.若对所有的， 则称该单步法是收敛的.定义2 如果存在正常数以及，使得对任意的初始出发值，某单步法的相应精确解以及，对所有的，恒有则称该单步法是稳定的.定义3 对给定的微分方程和给定的步长h，如果由单步法（显式或隐式）计算时有大小为的误差，即，而引起其后值的变化小于，则说该单步法是绝对稳定的.一般只限于典型的微分方程考虑数值方法的绝对稳定性（我们仅限于为实数的情形）.若对所有，单步法都绝对稳定，则称为绝对稳定区间3.3微分方程数值解法的编程实现3.1 常微分方程的符号解法在我们开始着手编写程序前，我们先一起了解一下什么是算法，算法是计算机编程的基石。算法的定义是：算法是指解题方案的准确而完整的描述，是一系列解决问题的清晰指令，算法是一种用系统的方法描述如何解决问题的策略机制3,14,15。比如Euler算法向我们介绍了如何运用Euler法来求解问题，是一种解决问题的思路，只有先掌握了算法逻辑，我们才可以根据具体问题需求来编写Euler法的MATLAB程序，或者编写其他编程软件的程序。而且根据实际问题的不同，我们可以灵活的使用不同的算法。我们接着来了解一下什么是常微分方程的符号解法。当一个微分方程能够求得它的解析解（通解）时，我们就可以通过MATLAB工具箱里面的函数求出该微分方程的精确解。但在实际生活中我们更需要求得微分方程的近似解或者特解，我们也可以通过MATLB求出它的数值解。一阶常微分方程（First-order Ordinary Differential Equation, ODE）写为其解为该解满足，并在初始条件下，可求得唯一解。在MATLAB中，常微分方程的符号解可使用dsolve函数求得，调用格式为dsolve(eqns)； %求出变量为t的微分方程的解析解dsolve(eqns,conds)； %求出变量为t的微分方程初值问题的解析解dsolve(eqns,Value)； %求出变量为x的微分方程的解析解dsolve(eqns,conds,Name,Value)； %求出变量为x的微分方程初值问题的解析解在MATLAB程序代码中，D表示对独立变量的微分，即，所以用户不能再自定义包含D的变量。微分方程在输入时，要输入成“Dy”，要输入成“D2y”若初始条件的数目比被微量数目少，那么结果中将包括不定常数等。系统默认“t”为自变量，如果所求的方程中自变量是“x”，必须要在语句的末尾定义该方程中的自变量为“x”，否则系统会把“t”认为自变量，从而造成错误，所以要引起注意。 我们通过两个例子来了解一下dsolve函数的使用方法：1、计算微分方程的解在MATLAB命令窗口输入clear allz=dsolve('Dy=-x','Dx=y')运行程序输出结果z = x: 1x1 sym y: 1x1 sym>> z.x %查询x值ans =C2*cos(t) - C1*sin(t)>> z.y %查询y值ans = C1*cos(t) - C2*sin(t)2、计算微分方程的通解在命令窗口输入以下命令clear ally=dsolve('Dy=x*2*y','x')输出结果y = C2*exp(x2)所以，通解为3.2 Euler法的编程实现（Euler法算法见附录1）3、用Euler法求解一阶常微分方程的数值解编写Euler法的MATLAB程序代码如下：function T= Euler(f,x0,y0,xn,N)% Euler.m函数为用Euler法求解微分方程% f为一阶常微分方程的一般表达式的右端函数；% x0,y0为初始条件% xn为取值范围的一个端点；% N为区间的个数；% x为求解微分方程组的值x=zeros(1,N+1); %x为Xn构成的向量y=zeros(1,N+1); %y为Yn构成的向量x(1)=x0;y(1)=y0;h=(xn-x0)/N;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(f,x(n),y(n);endT=x',y'plot(x,y,'b*:');legend('Eule求得的数值解');在MATLAB命令窗口输入clear all>> x=Euler (x,y) y-2*x/y,0,1,1,10)运行程序，输出结果如下：表3.1 Euler法求解一阶常微分方程的数值解 x y（Euler）0 1.00000.1000 1.10000.2000 1.19180.3000 1.27740.4000 1.35820.5000 1.43510.6000 1.50900.7000 1.58030.8000 1.64980.9000 1.71781.0000 1.7848图 3.13.3 隐式Euler法的编程实现（隐式Euler法算法见附录2）4、用隐式Euler法求解一阶常微分方程的数值解.编写隐式Euler法的MATLAB程序代码如下：function T=diEuler(f,x0,y0,xn,N)% diEuler.m为隐式Euler法求微分方程的数值解% f为一阶常微分方程的一般表达式的右端函数；% x0,y0为初始条件% xn为取值范围的一个端点；% n为区间的个数；% x为求解微分方程组的值x=zeros(1,N+1); %x为Xn构成的向量y=zeros(1,N+1); %y为Yn构成的向量x(1)=x0;y(1)=y0;h=(xn-x0)/N;for n=1:N%用迭代法求y(n+1) x(n+1)=x(n)+h; z0=y(n)+h*feval(f,x(n),y(n); for k=1:3 z1=y(n)+h*feval(f,x(n+1),z0); if abs(z1-z0)<1e-3; break; end z0=z1; end y(n+1)=z1;endT=x',y'plot(x,y,'b*:');legend('diEuler求得的数值解');在MATLAB命令窗口输入clear allx=diEuler(x,y) y-2*x/y,0,1,1,10)运行程序，输出结果如下：表3.2 隐式Euler法求解一阶常微分方程的数值解x y(diEuler)0 1.0000 0.1000 1.0909 0.2000 1.1743 0.3000 1.2517 0.4000 1.3237 0.5000 1.3910 0.6000 1.4540 0.7000 1.5128 0.8000 1.5676 0.9000 1.6183 1.0000 1.6647图 3.23.4 改进Euler法的编程实现（改进Euler法算法见附录3）5、用改进Euler法求解一阶常微分方程的数值解编写改进Euler法的MATLAB程序代码如下：function T=TranEuler(f,x0,y0,xn,N)% TranEuler.m函数为用改进Euler法求解微分方程% f为一阶常微分方程的一般表达式的右端函数；% x0,y0为初始条件% xn为取值范围的一个端点；% N为区间的个数；% x为求解微分方程组的值x=zeros(1,N+1); %x为Xn构成的向量y=zeros(1,N+1); %y为Yn构成的向量x(1)=x0;y(1)=y0;h=(xn-x0)/N;for n=1:N x(n+1)=x(n)+h; z0=y(n)+h*feval(f,x(n),y(n); y(n+1)=y(n)+h/2*(feval(f,x(n),y(n)+feval(f,x(n+1),z0);endT=x',y'plot(x,y,'b*:');legend('Euler求得的数值解');clear allx=TranEuler(x,y) y-2*x/y,0,1,1,10)运行程序，输出结果如下： 表3.3 改进Euler法求解一阶常微分方程的数值解x y0 1.0000 0.1000 1.0959 0.2000 1.1841 0.3000 1.2662 0.4000 1.3434 0.5000 1.4164 0.6000 1.4860 0.7000 1.5525 0.8000 1.6165 0.9000 1.6782 1.0000 1.7379图 3.3我们已经介绍完了关于微分方程的Euler法、隐式Euler法、改进Euler法的MATLAB的编程过程以及计算结果，但是以上三种方法不止一种MATLAB程序，也许有更简单的方法在等着读者去发现。通过输出的数据看来，每种算法得出的结果各不相同，那么我们究竟应该选择哪一种算法更合适呢？让我们把以上三种方法的输出结果放在一起作为对比，相信你很快会得你的结论。首先我们先编写一个算法比较程序：clear all%bijiao.m一阶微分方程数值解对比程序%输入导函数，x0，y0，xn，NS=(x,y) y-2*x/y;R=Euler(S,0,1,1,10);R1=diEuler(S,0,1,1,10);R2=TranEuler(S,0,1,1,10);hold onplot(R(:,1)',R(:,2)','r*-.')hold onplot(R1(:,1)',R1(:,2)','g+-')hold onplot(R2(:,1)',R2(:,2)','bx:')legend('Euler法求得的解','diEuler法求得的解','TranEuler法求得的解')运行程序输出结果如下：表3.4 对比Euler法，idEuler法以及TranEuler法的数值解x y(Euler) y(diEuler) y(TranEuler) 0 1.0000 1.0000 1.0000 0.1000 1.1000 1.0909 1.0959 0.2000 1.1918 1.1743 1.1841 0.3000 1.2774 1.2517 1.2662 0.4000 1.3582 1.3237 1.3434 0.5000 1.4351 1.3910 1.4164 0.6000 1.5090 1.4540 1.4860 0.7000 1.5803 1.5128 1.5525 0.8000 1.6498 1.5676 1.6165 0.9000 1.7178 1.6183 1.6782 1.0000 1.7848 1.6647 1.7379图 3.4我们可以直观的观察到，在步长时，Euler法与改进Euler法的曲线贴合较近，隐式Euler法与改进Euler法的曲线贴合较远。让我们看看当步长分别为；时的曲线贴合情况（曲线图像见附录7）。通过对比可以说明，当步长取得较大时，Euler法的精度不够高。当步长变小时，Euler精度会有所提高。隐式Euler与Euler法情况相似，而改进Euler法的精度比前两种方法的精度都要高，所以当时，隐式Euler与Euler法的数值解曲线都在向改进Euler法的曲线贴合，误差不断减小；当时，这三种方法之间的误差可以忽略不计。所以我们在编程计算一阶微分方程的时候，当步长当时，我们应该选用改进Euler法来求数值解；当时，三种方法我们可以任选一种使用。另外当步长取定以后，步数越多，误差也会越大。3.5 二阶Runge-Kutta法的编程实现（二阶Runge-Kutta法算法见附录4）6、用二阶Runge-Kutta法来求解一阶常微分方程的数值解编写二阶Runge-Kutta法的MATLAB程序代码如下：function R=RK2(f,x0,y0,xn,N)% RK2.m函数为用二阶Runge-Kutta法求微分方程的解% f为微分方程% x0,y0为左右端点% xn为给定的初始值% N为给定迭代步长；% R为求微分方程的解h=(y0-x0)/N;x=zeros(1,N+1);y=zeros(1,N+1);T=x0:h:y0;y(1)=xn;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+2*h/3,y(n)+2*h*k1/3); y