一阶微分方程的matlab数值解法

1、*大学毕业设计-PAGE - 2 -. z. - - .可修编 .毕业论文题 目一阶微分方程的matlab数值解法学 院数学科学学院专 业数学与应用数学班 级数学1303班 学 生邹健峰学 号7指导教师邢顺来讲师二一七年五月七日-PAGE . z.摘 要微分方程的数值解法是近现代数学家和科学家们研究的热点，微分方程的MATLAB数值解法可以帮助我们解决现实生活中的许多数学问题，提高计算机帮助我们解决问题的效率。本文主要讨论研究一阶微分方程的MATLAB数值解法中的三种Euler法和三种Runge-Kutta法，介绍以上六种方法的主要容，简单介绍了MATLAB的相关容和微分方程的开展简史。通过具

2、体的微分方程来研究以上算法的编程实现，通过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 diff

3、erential equation can help us solve many mathematical problems in real life,improve the efficiency of the puter 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,introd

4、uce the main contents of the above si* 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

5、 are e*plored 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.Keywords：First Order Differential Equation; Numerical

6、Solution; MATLAB; Euler Method; Runge-Kutta Method目录摘要 PAGEREF _Toc29420 IABSTRACT PAGEREF _Toc434 II1 前言 PAGEREF _Toc19134 11.1 引言 PAGEREF _Toc14582 11.2常微分方程的开展简史 PAGEREF _Toc21991 21.3 国外研究现状 PAGEREF _Toc15673 21.4 研究主要容及研究方案 PAGEREF _Toc11775 31.4.1 研究的主要容 PAGEREF _Toc1002 3 研究方案 PAGEREF _Toc297

7、35 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附录3

8、2-. z.1前言1.1引言James Bernoulli在1676年写给Newton的信中首次谈到了微分方程。大约到了十八世纪中期，微分方程才成为一门独立的学科，微分方程是研究和提醒自然规律不可或缺的重要工具。众所周知，在1846年数学家和天文学家共同发现的海王星就是微分方程的功绩。还有在1991年，科学考察队在Alps山脉发现了一个冰封保存完好的人类，然后通过仪器测得碳元素衰变的速率，建立微分方程数学模型，得出了这个人是生活在距今5000多年前的时代。Newton、Leibniz、Bernoulli family、Lagrange、Laplace等著名数学家们都对微分方程的开展做出了巨大的

9、奉献。微分方程的定义：如果知道自变量、未知函数以及函数的导数或微分组成的关系式我们就称为微分方程。常微分方程的定义：通过求解微分方程得出未知函数，并且自变量只有一个的微分方程我们就叫做常微分方程，自变量个数为两个或者两个以上的微分方程称为偏微分方程1。我们把含有n个互不相等的任意常数的解称为n阶常微分方程的通解。微分方程的特解是指满足微分方程初值条件的解1。我们在生活中会经常用到一些常微分方程模型，例如英国著名人口统计学家Malthus提出的著名的Malthus人口模型；荷兰著名生物学家Verhulst提出的logistic人口增长模型；传染病模型SI模型、SIS模型无免疫性、SIR强免疫性模

10、型；两生物种群生态模型；气象学家Lorenz提出的Lorenz方程用来预测天气变化；化学动力学模型2。可见微分方程与我们的生活有着密不可分的关系。MATLAB是分别取了Matri*和Laboratory两个单词的前三个字母组成的,是矩阵实验室的名称缩写。MATLAB是来自美国大学的计算机教师Cleve Moler创造的，并且由Cleve Moler和JohnLillte把其推向了市场，并得到了广阔用户的喜爱。MATLAB是数学科学和工程科学的高级计算机语言，用户可以用数学语言来编写程序，编程逻辑非常类似人们日常书写数学公式的逻辑，具有简单性、广泛性、开放性、效率高、价格低等特点。并且MATLA

11、B工具箱储存的MATLAB函数库里面有大量的函数程序，可以解决许多的数学和工程科学问题。从2010版开场MATLAB的功能越来越多，编译器可以把M文件转换为C语言文件或者其他几种常用的计算机文件，并给用户提供源程序代码进展自我开发，结合Math Works公司提供的数据库，用户可以利用MATLAB来开发出功能强大的程序3,4,5。微分方程的MATLAB数值解法可以帮助我们解决许多现实生活中的数学问题，提高计算机帮助我们解决问题的效率，如今微分方程的MATLAB数值解法已经开展得比较成熟，有许多关于数学领域的MATLAB书籍可供我们选择阅读，很好的解决了我们的学习需求和使用问题。MATLAB应用

12、领域也非常广泛，比方在数学建模和数学实验领域，我们可以通过直接调用MATLAB里面的微分方程函数文件，在操作界面输入初始数据，还把抽象复杂的数学问题变成简单形象的二维或者三维图形，来帮助我们分析问题，找出解题方法，从而解决我们生活中负责的数学问题，表达了数学源自生活又联系生活的主要特征。除此之外，在物理学、工程科学等科学领域中MATLAB也具有重要地位6,7。1.2常微分方程的开展简史常微分方程的开展可以简单的分为以下几个阶段：“求通解阶段；“求定解阶段；“求所有解阶段；“求特殊解阶段。在常微分方程的“求通解阶段，人们希望能够用初等函数来表示解，Newton就希望能够用别离变量法来求解一阶微分

13、方程问题，但Euler却尝试用积分因子来求解一阶微分方程问题。而Bernoulli在研究微分方程的过程中提出了著名的Bernoulli方程，Riccati也在研究过程中提出了Riccati方程，他们的发现对微分方程未来的开展打下了根底，起到了积极的作用。直到Liouville证明了Riccati方程没有一般的初等解，还有Cauchy提出了微分方程的初值问题，常微分方程进入了“求定解阶段。直到19世纪末，由于天体力学的研究需求，常微分方程从“求定解阶段进入了“求所有解的阶段。首先是Poincar提出了定性理论，接着是Hilbert提出了关于极限环个数的问题，促进了微分方程定性理论的开展，紧接着L

14、yapunov提出了运动稳定性理论，广泛应用在物理和工程学中。而在动力系统方面，Birkhoff、Arnold、Smale都做出了奉献。到了20世纪70年代计算机技术的开展带动了常微分方程的开展，从“求所有解阶段进入到了“求特定解阶段，发现了一些特定解和方程，比方混沌解、孤立子、奇异吸引子、分形。关于常微分方程的研究还在继续，如今研究的领域更加广泛，更多的微分方程模型被提出和应用起来了1,8,9,10。 1.3国外研究现状关于微分方程数值解法的研究，在没有MATLAB以及类似的计算机软件问世之前，数学家和科学家们只能研究和解决一些低阶的微分方程问题，并且有些方程的解题过程过于繁琐和复杂。但是M

15、ATLAB的出现促进了数学科学研究的开展，对数学和相关领域的开展做出了巨大奉献。微分方程的MATLAB数值解法也迅速开展起来，国外现在主要研究通过MATLAB等计算机软件来研究：一些高阶微分方程的解的有关问题；一些复杂的微分方程是否存在定解问题；利用MATLAB来计算一些著名的数学模型问题，比方Volterra被捕食 捕食模型、竞争模型、共生模型、商品供需模型等等；研究混沌解、孤立子。在以上几方面的研究上都取得了不小的进步，使常微分方程的体系更加的完善了。在国，2009年中国数学会成功举办了第四届全国青年常微分方程学术会议，许多国知名学者和学术研究者都到场参加学术交流，这一项活动促进了我国常微

16、分方程领域的开展，为青年学者提供了珍贵经历，这次会议也受到了国外数学爱好者的高度关注。在书籍方面，我国也发行了许多关于微分方程数值分析和微分方程数值解法方面的书籍，都有详细介绍微分方程关于MATLAB的数值解法，比方林成森编著的第二版?数值计算方法?里面就有详细介绍常微分方程初值问题的数值解法，包括有：离散变量法和离散误差法；单步法；单步法的相容性、收敛性和稳定性；多步法；差分方程的简介；线性步法的相容性、收敛性和数值稳定性；常微分方程组和高阶微分方程的数值解法等几个方面，成为了21世纪高等院校的教材。周品、何正风等编著的?MATLAB数值分析?作为了理工科各专业本科生、研究生以及应用MATL

17、AB相关科技人员学习MATLAB数值分析、建模、仿真的教材或者参考书。现如今，由于信息时代的飞速开展，数值计算方面的计算机软件越来越多，功能也越来越全面，所以微分方程的数值解法受到了更多人的重视，开展前景变得更好了，这对于学习微分方程的数值解法的人来说是非常好的消息。1.4研究的主要容与方法1.4.1 研究的主要容本文针对一阶常微分方程的数值解法中的MATLAB数值解法进展算法介绍和编程实现，利用MATLAB对一阶微分方程中的Euler法、隐式Euler法、改进Euler法以及Runge-Kutta法进展精度分析并比较，提醒以上几种方法的优缺点，编程实现以上方法的数值解法，并通过具体事例来验证

18、程序的完整性、稳定性、以及缺乏。1.4.2研究方法一阶微分方程的MATLAB数值解法的编程实现需要通过Euler法、隐式Euler法、改进Euler法以及二到四阶Runge-Kutta法等方法的算法逻辑，针对不同一阶微分方程来编写MATLAB程序，然后再反复调试程序，找出其中的缺乏并改进，然后保存程序并做记录，以便下次使用时直接调用。1.5研究难点1微分方程的模型多样，研究方向多，不同的模型有不同的数值分析和数值解法；2微分方程的MATLAB数值解法需要符合实际问题，限制条件比较多；3微分方程的类型多样，MATLAB的函数库里只有局部微分方程的调用函数，不能满足研究的需求，需要自定义调用函数；

19、4微分方程的MATLAB数值解法需要反复调试，得出满意结果花费时间较长。2预备知识2.1显式Euler法显示Euler法也被人们简称为Euler法，是解决微分方程初值问题最简单的数值方法。初值问题的解代表通过点的一条称为微分方程的积分曲线。积分曲线上每一点的切线斜率等于函数在这点的值。Euler方法的求解过程是从初始点出发,作积分曲线在点的切线其斜率为，与直线相交与点，得到作为的近似值，过点，以为斜率的切线方程为当时，得这样就获得了点的坐标3,11。不断重复上述过程，就可以得到一系列的点。对已求得的点以为斜率做直线，当时，便可得到方程取，这样从逐一算出所对应的数值解，同样也获得了一条曲线作为原

20、始曲线的近似曲线，因此又称Euler法为折线法。在计算时，一般取常数，则显式Euler公式为.2.12.2隐式Euler法假设我们用近似等式来替代显式Euler公式右端的积分项，就得到再用替代，替代，我们就可以得到这就是隐式Euler公式。后退的Euler法也被称为隐式Euler法。而我们经常用到的梯形公式。对方程的两端在区间上积分，得要通过这个积分关系式得的近似值，只要近似的算出其中的积分项。假设使用梯形方法计算积分项再代入中，得将式子中的分别用代替，有.2.2这种差分格式就称为梯形格式，梯形格式上是显示Euler公式与隐式Euler公式的算术平均。梯形方法是隐式单步法，可以用迭代法求解，与

21、后退的Euler法一样，仍用Euler法进展初值迭代，所以梯形公式的迭代公式为2.32.3改进Euler法显示Euler公式计算工作量小，但精度不高，假设对计算结果的精度要求高，就必须使用就是精度更高的计算方法。微分方程初值问题的梯形公式相比显示Euler公式，精度提高了，但是用迭代的算法来求解计算工作量大。所以可以改进思路：综合表示Euler公式和梯形公式便可以得到改进的Euler公式。先用显示Euler公式求出一个初步的近似值，记为，称之为预测值，它的精度不高，再用梯形公式对它进展一次矫正，也就是迭代一次，求出，称为矫正值12,13。这种预测矫正的方法我们就叫做改进的Euler法.2.4它

22、也可以表示为嵌套模式.或者表示为以下的平均化形式2.52.4二阶Runge-Kutta法二阶Runge-Kutta的一般公式为2.6局部截断误差为.其中加权系数和选点调节系数都是特定系数。本质是用差商代替导数的情况下，把在区间上的假设干点的值加权平均，再作为斜率，作过点的直线交直线，求得下一点。适中选择是提高精度的方法14,15,16。我们常用的二阶Runge-Kutta公式有中点方法2.7和Heun方法2.82.5 三阶Runge-Kutta法三阶Runge-Kutta一般公式为2.9局部截断误差.常用的三阶Runge-Kutta公式有2.10和2.112.6 四阶Runge-Kutta法继

23、续以上的演算过程，可以推导出各种四阶Runge-Kutta公式，下面介绍一个被被称为经典Runge-Kutta公式，其公式的格式为2.12其截断误差为.四阶Runge-Kutta法要求所求的解具有比较好的光滑性，反之，使用四阶Runge-Kutta法求得的解精度比改进的Euler法求得的解精度低。2.7单步法的收敛性和稳定性定义1假设微分方程的右端函数在带形区域中连续，并且关于满足Lipschitz条件.假设对所有的，则称该单步法是收敛的.定义2如果存在正常数以及，使得对任意的初始出发值，*单步法的相应准确解以及，对所有的，恒有则称该单步法是稳定的.定义3对给定的微分方程和给定的步长h，如果由

24、单步法显式或隐式计算时有大小为的误差，即，而引起其后值的变化小于，则说该单步法是绝对稳定的.一般只限于典型的微分方程考虑数值方法的绝对稳定性我们仅限于为实数的情形.假设对所有，单步法都绝对稳定，则称为绝对稳定区间3.3微分方程数值解法的编程实现3.1 常微分方程的符号解法在我们开场着手编写程序前，我们先一起了解一下什么是算法，算法是计算机编程的基石。算法的定义是：算法是指解题方案的准确而完整的描述，是一系列解决问题的清晰指令，算法是一种用系统的方法描述如何解决问题的策略机制3,14,15。比方Euler算法向我们介绍了如何运用Euler法来求解问题，是一种解决问题的思路，只有先掌握了算法逻辑，

25、我们才可以根据具体问题需求来编写Euler法的MATLAB程序，或者编写其他编程软件的程序。而且根据实际问题的不同，我们可以灵活的使用不同的算法。我们接着来了解一下什么是常微分方程的符号解法。当一个微分方程能够求得它的解析解通解时，我们就可以通过MATLAB工具箱里面的函数求出该微分方程的准确解。但在实际生活中我们更需要求得微分方程的近似解或者特解，我们也可以通过MATLB求出它的数值解。一阶常微分方程First-order Ordinary Differential Equation, ODE写为其解为该解满足，并在初始条件下，可求得唯一解。在MATLAB中，常微分方程的符号解可使用dsol

26、ve函数求得，调用格式为dsolve(eqns)； %求出变量为t的微分方程的解析解dsolve(eqns,conds)； %求出变量为t的微分方程初值问题的解析解dsolve(eqns,Value)； %求出变量为*的微分方程的解析解dsolve(eqns,conds,Name,Value)； %求出变量为*的微分方程初值问题的解析解在MATLAB程序代码中，D表示对独立变量的微分，即，所以用户不能再自定义包含D的变量。微分方程在输入时，要输入成“Dy，要输入成“D2y假设初始条件的数目比被微量数目少，则结果中将包括不定常数等。系统默认“t为自变量，如果所求的方程中自变量是“*，必须要在语句

27、的末尾定义该方程中的自变量为“*，否则系统会把“t认为自变量，从而造成错误，所以要引起注意。我们通过两个例子来了解一下dsolve函数的使用方法：1、计算微分方程的解在MATLAB命令窗口输入clear allz=dsolve(Dy=-*,D*=y)运行程序输出结果z = *: 1*1 sym y: 1*1 sym z.* %查询*值ans =C2*cos(t) - C1*sin(t) z.y %查询y值ans = C1*cos(t) - C2*sin(t)2、计算微分方程的通解在命令窗口输入以下命令clear ally=dsolve(Dy=*2*y,*)输出结果y = C2*e*p(*2)所

28、以，通解为3.2Euler法的编程实现Euler法算法见附录13、用Euler法求解一阶常微分方程的数值解编写Euler法的MATLAB程序代码如下：function T= Euler(f,*0,y0,*n,N)% Euler.m函数为用Euler法求解微分方程% f为一阶常微分方程的一般表达式的右端函数；% *0,y0为初始条件% *n为取值围的一个端点；% N为区间的个数；% *为求解微分方程组的值*=zeros(1,N+1); %*为*n构成的向量y=zeros(1,N+1); %y为Yn构成的向量*(1)=*0;y(1)=y0;h=(*n-*0)/N;for n=1:N *(n+1)=

29、*(n)+h; y(n+1)=y(n)+h*feval(f,*(n),y(n);endT=*,yplot(*,y,b*:);legend(Eule求得的数值解);在MATLAB命令窗口输入clear all *=Euler (*,y) y-2*/y,0,1,1,10)运行程序，输出结果如下：表3.1 Euler法求解一阶常微分方程的数值解*yEuler0 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.71

30、781.0000 1.7848图3.13.3隐式Euler法的编程实现隐式Euler法算法见附录24、用隐式Euler法求解一阶常微分方程的数值解.编写隐式Euler法的MATLAB程序代码如下：function T=diEuler(f,*0,y0,*n,N)% diEuler.m为隐式Euler法求微分方程的数值解% f为一阶常微分方程的一般表达式的右端函数；% *0,y0为初始条件% *n为取值围的一个端点；% n为区间的个数；% *为求解微分方程组的值*=zeros(1,N+1); %*为*n构成的向量y=zeros(1,N+1); %y为Yn构成的向量*(1)=*0;y(1)=y0;h

31、=(*n-*0)/N;for n=1:N%用迭代法求y(n+1) *(n+1)=*(n)+h; z0=y(n)+h*feval(f,*(n),y(n); for k=1:3 z1=y(n)+h*feval(f,*(n+1),z0); if abs(z1-z0)1e-3; break; end z0=z1; end y(n+1)=z1;endT=*,yplot(*,y,b*:);legend(diEuler求得的数值解);在MATLAB命令窗口输入clear all*=diEuler(*,y) y-2*/y,0,1,1,10)运行程序，输出结果如下：表3.2 隐式Euler法求解一阶常微分方程的

32、数值解* 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法算法见附录35、用改进Euler法求解一阶常微分方程的数值解编写改进Euler法的MATLAB程序代码如下：function T=TranEuler(f,*0,y0,*n,N)% TranEuler.m函数为用改进Eu

33、ler法求解微分方程% f为一阶常微分方程的一般表达式的右端函数；% *0,y0为初始条件% *n为取值围的一个端点；% N为区间的个数；% *为求解微分方程组的值*=zeros(1,N+1); %*为*n构成的向量y=zeros(1,N+1); %y为Yn构成的向量*(1)=*0;y(1)=y0;h=(*n-*0)/N;for n=1:N *(n+1)=*(n)+h; z0=y(n)+h*feval(f,*(n),y(n); y(n+1)=y(n)+h/2*(feval(f,*(n),y(n)+feval(f,*(n+1),z0);endT=*,yplot(*,y,b*:);legend(E

34、uler求得的数值解);clear all*=TranEuler(*,y) y-2*/y,0,1,1,10)运行程序，输出结果如下：表3.3 改进Euler法求解一阶常微分方程的数值解*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的编程过程以及计

35、算结果，但是以上三种方法不止一种MATLAB程序，也许有更简单的方法在等着读者去发现。通过输出的数据看来，每种算法得出的结果各不一样，则我们终究应该选择哪一种算法更适宜呢？让我们把以上三种方法的输出结果放在一起作为比照，相信你很快会得你的结论。首先我们先编写一个算法比较程序：clear all%bijiao.m一阶微分方程数值解比照程序%输入导函数，*0，y0，*n，NS=(*,y) y-2*/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*-.)

36、hold onplot(R1(:,1),R1(:,2),g+-)hold onplot(R2(:,1),R2(:,2),b*:)legend(Euler法求得的解,diEuler法求得的解,TranEuler法求得的解)运行程序输出结果如下：表3.4 比照Euler法，idEuler法以及TranEuler法的数值解* 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

37、 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。通过比照可以说

38、明，当步长取得较大时，Euler法的精度不够高。当步长变小时，Euler精度会有所提高。隐式Euler与Euler法情况相似，而改进Euler法的精度比前两种方法的精度都要高，所以当时，隐式Euler与Euler法的数值解曲线都在向改进Euler法的曲线贴合，误差不断减小；当时，这三种方法之间的误差可以忽略不计。所以我们在编程计算一阶微分方程的时候，当步长当时，我们应该选用改进Euler法来求数值解；当时，三种方法我们可以任选一种使用。另外当步长取定以后，步数越多，误差也会越大。3.5二阶Runge-Kutta法的编程实现二阶Runge-Kutta法算法见附录46、用二阶Runge-Kutta

39、法来求解一阶常微分方程的数值解编写二阶Runge-Kutta法的MATLAB程序代码如下：function R=RK2(f,*0,y0,*n,N)% RK2.m函数为用二阶Runge-Kutta法求微分方程的解% f为微分方程% *0,y0为左右端点% *n为给定的初始值% N为给定迭代步长；% R为求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;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(n+1)=y

40、(n)+h*(k1+3*k2)/4;endR=T,y编写二阶Heun法函数图像的MATLAB程序代码如下：S=(*,y) *-y+1;t1,y1=ode23(S,0:0.1:1,1);R=RK2(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b* -.)legend(二阶Heun方法求得的解,ode23求得的解)hold offK=t1,y1将这两个M文件分别命名为RK2.m和RK2_Heun.m放在一个文件夹中，然后直接运行RK2_Heun.m程序，输出结果如下：表3.5二阶Runge-Kutta法求解一阶常微分方程的数值解* y(RK2

41、) y(ode23)0 1.0000 1.00000.1000 1.0000 1.00030.2000 1.0003 1.00250.3000 1.0009 1.00840.4000 1.0021 1.01940.5000 1.0041 1.03690.6000 1.0071 1.06240.7000 1.0112 1.09680.8000 1.0167 1.14130.9000 1.0238 1.19681.0000 1.0325 1.2642图3.53.6三阶Runge-Kutta法的编程实现三阶Runge-Kutta法算法见附录57、用三阶Runge-Kutta法来求解一阶常微分方程的数

42、值解编写三阶Runge-Kutta法的MATLAB程序代码如下：function R=RK3(f,*0,y0,*n,N)% RK3.m函数为用三阶Runge-Kutta法求微分方程的解% f为微分方程% *0,y0为左右端点% *n为给定的初始值% N为给定迭代步长；% R为求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+3*h/4,y(n

43、)+3*h/4*k2); y(n+1)=y(n)+(2*k1+3*k2+4*k3)/9;endR=T,y编写三阶Runge-Kutta法函数图像的MATLAB程序代码如下：S=(*,y) *2-y+1;t1,y1=ode23(S,0:0.1:1,1)R=RK3(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b*-.)legend(三阶RK法求得的解,ode23求得的解)hold off将这两个M文件分别命名为RK3.m和RK3_1.m放在一个文件夹中，然后直接运行RK3_1.m程序，输出结果如下：表3.6 三阶Runge-Kutta法求解一

44、阶常微分方程的数值解*yRK3 yode230.1000 1.0003 1.00030.4000 1.0199 1.01940.6000 1.0641 1.06240.7000 1.0994 1.09680.8000 1.1450 1.14130.9000 1.2018 1.19681.0000 1.2707 1.2642图3.63.7四阶Runge-Kutta法的编程实现四阶Runge-Kutta法算法见附录68、用四阶Runge-Kutta法来求解一阶常微分方程的数值解编写四阶Runge-Kutta法的MATLAB程序代码如下：function R=RK4(f,*0,y0,*n,N)% R

45、K4.m函数为用四阶Runge-Kutta法求微分方程的解% f为微分方程% *0,y0为左右端点% *n为给定的初始值% N为给定迭代步长；% R为求微分方程的解h=(y0-*0)/N;*=zeros(1,N+1);y=zeros(1,N+1);T=*0:h:y0;y(1)=*n;for n=1:N k1=h*feval(f,T(n),y(n); k2=h*feval(f,T(n)+h/2,y(n)+k1/2); k3=h*feval(f,T(n)+h/2,y(n)+k2/2); k4=h*feval(f,T(n)+h,y(n)+k3); y(n+1)=y(n)+(k1+2*k2+2*k3+

46、k4)/6;endR=T,y编写四阶Runge-Kutta法函数图像的MATLAB程序代码如下：S=(*,y) *2-y+1;t1,y1=ode45(S,0:0.1:1,1)R=RK4(S,0,1,1,10)plot(R(:,1),R(:,2),r)hold onplot(t1,y1,b*-.)legend(四阶RK法求得的解,ode45求得的解)hold off将这两个M文件分别命名为RK4.m和RK4_1.m放在一个文件夹中，然后直接运行RK4_1.m程序，输出结果如下：表3.7 四阶Runge-Kutta法求解一阶微分方程的数值解* y(RK4) y(ode45) 0 1.0000 1.

47、0000 0.1000 1.0003 1.0003 0.2000 1.0025 1.0025 0.3000 1.0084 1.0084 0.4000 1.0194 1.0194 0.5000 1.0369 1.0369 0.6000 1.0624 1.0624 0.7000 1.0968 1.0968 0.8000 1.1413 1.1413 0.9000 1.1969 1.1969 1.0000 1.2642 1.2642图 3.7我们可以看到，用四阶Runge-Kutta法求一阶微分方程的数值解曲线和ode45求得的数值解曲线完全重合，说明在给定步长的情况下，四阶Runge-Kutta法阶

48、的精度到达了我们求一阶微分方程的精度要求，与准确解之间的误差非常小。接下来我们可把以上三种Runge-Kutta方法进展比较一下，看看对于同样的微分方程各数值方法的数值解之间的误差是多少，以帮助我们更好更恰当的选择数值方法。编写一阶微分方程数值解比照的MATLAB程序代码如下：clear all%duibi.m一阶微分方程数值解比照程序%输入导函数，*0，y0，*n，h，NS=(*,y) *-y+1;t1,y1=ode23(S,0:0.1:1,1);t2,y2=ode45(S,0:0.1:1,1)R=RK2(S,0,1,1,10);K=RK3(S,0,1,1,10);G=RK4(S,0,1,1

49、,10);plot(R(:,1),R(:,2),r,K(:,1),K(:,2),c,G(:,1),G(:,2),r)hold onplot(t1,y1,b* -,t2,y2,m+ :)legend(二阶Heun方法求得的解,三阶RK法求得的解,四阶RK法求得的解,ode23求得的解,ode45求得的解)hold off将程序命名为duibi.m，和RK2.m，RK3.m，RK4.m放在一起，运行duibi.m程序，输出函数图像结果请看图 3.7。表3.8一阶微分方程五种数值解比照表* y(RK2) y(RK3) y(RK4) y(ode23) y(ode45)0 1.0000 1.00001.

50、0000 1.0000 1.00000.1000 1.00001.00031.0003 1.0003 1.00030.2000 1.00031.00261.0025 1.0025 1.00250.3000 1.00091.0086 1.0084 1.0084 1.00840.4000 1.0021 1.0199 1.0194 1.0194 1.01940.5000 1.00411.0380 1.0369 1.0369 1.03690.6000 1.00711.0641 1.0624 1.0624 1.06240.7000 1.01121.0994 1.0968 1.0968 1.09680.8

51、000 1.01671.1450 1.1413 1.1413 1.14130.9000 1.0238 1.2018 1.1969 1.1968 1.19691.0000 1.0325 1.2707 1.2642 1.2642 1.2642图3.8我们通过观察图 3.7我们可以发现二阶Heun法的数值解曲线与其它数值解曲线偏离程度较大，而且是增大趋势，而其它数值解曲线相互之间偏离程度小，其中四阶Runge-Kutta法与ode45法的数值解曲线重合，与ode23法的数值解曲线几乎重合。在MATLAB中，函数ode23和ode45采用的是变步长的Runge-Kutta法，并且ode45的精度更高一

52、些。ode23适合的常微分方程问题类型为低精度容差或适度非刚性问题；ode45适合的常微分方程问题类型为非刚性方程。所以在精度误差的允许围，我们优先选择四阶Runge-Kutta法或者函数ode45来解决一阶微分方程的初值问题。则在编程实现的时候，我们该如何选择适当的Runge-Kutta法呢？我们同样通过改变步长的方法来研究该问题，分别取步长为；，观察曲线的贴合情况Runge-Kutta法曲线比照见附录8。当步长在变小的过程中，RK3与RK4的曲线越来越贴近，重合程度越来越大，并且与RK2偏离程度也越来越大，说明RK2的精度比较低的，当步长之后，RK2的精度可能会比改进的Euler法还要低，

53、所以我们在编程实现一阶微分方程的数值解并对步长严格控制的时候，我们应该优先选择三阶或者四阶的Runge-Kutta法，控制数值解与准确解之间的误差。结 论本文探究的一阶微分方程的MATLAB数值解法有Euler法、隐式Euler、改进Euler法、二阶Runge-Kutta法、三阶Runge-Kutta法以及四阶Runge-Kutta法。采用比照分析的方法探究了以上方法的精度问题，如何编程实现算法，以及如何选择程序求解一阶微分方程初值问题。在给定步长的微分方程初值问题中，得出以下结论：利用Euler法可以求解简单的微分方程初值问题，Euler法是最简单的单步法，MATLAB编程实现操作简单，很

54、好的表达了算法的决策思想，在步长充分小的时候，Euler法才可用，一般不用于实际计算，但是使用方便，是一阶精度的算法；隐式Euler法也被称为后退的Euler法，隐式Euler法的算法在稳定性方面表现出色，但是迭代次数比Euler法多，计算过程比较久，在特殊场合人们也会使用该算法来计算；改进Euler法相比前面的两种算法来说，算法的精度是二阶的，具有良好的稳定性，受步长的影响不大，被称为预测矫正系统，MATLAB编程实现过程简单，人们通常会选用该算法来求解一阶微分方程的初值问题；二阶Runge-Kutta法的算法精度是二阶的，二阶Runge-Kutta法可以推导出改进的Euler法，也可以换算

55、成文中有介绍的中点方法和Heun方法；三阶Runge-Kutta法的算法精度是三阶的，稳定性也比教好；四阶Runge-Kutta法的算法精度是四阶的，是以上的算法中阶的精度最高的算法，稳定性也以上的算法中最好的，四阶Runge-Kutta法的公式比较多，人们经常选择经典的四阶Runge-Kutta法来解决一阶微分方程的初值问题。Runge-Kutta法的编程实现复杂程度比Euler法的要难一些，需要考虑Runge-Kutta算法的步长的取值，并且由于迭代的次数较多，所以算法的计算量大；单步法的收敛性与稳定性：针对给定步长的一阶微分方程初值问题，三种Euler法和二到四阶Runge-Kutta法

56、都具有收敛性；在实际应用中，选择数值方法的步长时，我们应该考虑稳定性，步长的取值变化应该在绝对稳定区间，以保证数值解与准确解之间的误差。参 考 文 献1王高雄，朱四铭，周之铭，王寿松. 常微分方程M. ：高等教育，20062C.W.吉尔. 常微分方程初值问题的数值解法M. 费景高 译. ：科学，19783德丰. MATLAB数值分析M. 第二版.：机械工业，20124楼顺天，于卫，闫华梁. MATLAB程序设计语言M. ：电子科技大学，19975夏省祥，于正文. 常用数值算法及其MATLAB实现M. ：清华大学，20146白峰杉. 数值计算引论M.第二版. ：高等教育，20107白峰杉，庆扬，关治. 数值计算原理M.：清华大学，20008文林. 数学史概论M. ：高等教育，20029杜瑞芝. 数学史辞典M. ：教育，200010VictorJ.katz. 数学史通论M. 文林 译. ：高等教育，200411Gerald Recktenwald. 数值方法和MATLAB实现与应用M.伍卫国，万群，辉译.：机械工业，200412吴开腾，谭燕梅，