常微分方程常用数值解法

1、第一章 绪论1.1 引言常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努

2、力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此

3、同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。1.2 常微分方程的概念1.常微分方程的定义含有未知量的等式称为方程，它表达了未知量所必须满足的某些条件。一般说来，凡含有自变量、未知函数以及未知函数的导数或

4、微分的方程称为微分方程。如果微分方程中的未知函数只依赖于一个自变量，则称为常微分方程；如果未知函数依赖于两个或多个的自变量，并且在方程中出现偏导数，则称为偏微分方程。在在一个微分方程中，未知函数最高阶导数的阶数，称为方程的阶数。如果一个微分方程关于未知函数及其各阶导数都是线性的，则称它为线性微分方程，否则称之为非线性微分方程。本论文主要介绍常微分方程，也简称微分方程。 以为未知函数，为自变量的一阶常微分方程的一般形式可表示为微分方程， (1.2.1) 将（1.2.1）中解出，则得到方程 (1.2.2)或 (1.2.3)也称（1.2.1）为一阶隐式微分方程，（1.2.2）为一阶显式微分方程，（1

5、.2.3）为一阶微分方程的微分形式。 阶隐式方程的一般形式为 (1.2.4) 阶显式方程的一般形式 (1.2.5)方程（1.2.4）中，如果函数对未知函数和它的各阶导数都是一次的，则称其为线性常微分方程，否则，称其为非线性微分方程。以为未知函数，为自变量的阶线性微分方程具有如下形式： (1.2.6) 2.常微分方程的解 设函数是定义在区间上的阶可微导数。如果把代入方程（1.2.4）后能使其成为恒等式，即 则称是微分方程（1.2.4）在区间上的一个解。例如，是微分方程在的一个解.是微分方程在区间的一个解.如果关系式决定的隐函数是方程（1.2.4）的解，则我们称是（1.2.4）的隐式解。例如一阶微

6、分方程 有隐式解 我们把含有个相互独立的任意常数的解 称为阶微分方程（1.2.4）的通解。在通解之中当一组任意常数确定时，所得到确定的解称为特解。例如，是二阶线性方程的通解，而都是其特解，其中是任意常数。一般地，方程的特解可由其通解中任意常数取确定的常数导出，且方程的通解不一定表示方程的所有解。3. 常微分方程初值问题为了确定微分方程的一个特解，我们可以给出这个微分方程的所满足的定解条件，常见的定解条件是初始条件，即方程（1.2.4）在某一点所满足的条件： （1.2.7）微分方程（1.2.4）连同初始条件（1.2.7）一起称为初始值问题。 一阶常微分方程初值问题是求解函数，且满足 （1.2.8

7、）其中为已知函数，为给定的初值。定理1.2.1 假设在矩形区域内连续，且关于变元Lipschitz连续，即存在正常数，使得对任意及，成立不等式 (1.2.9)其中常数称为Lipschitz常数,则处置问题存在惟一解，连续依赖于初值由常微分方程的基本理论，我们有：换句话说，若是如下初值问题 (1.2.10)的解，则 因此问题（1.2.8）是适定的。定理1.2.2 当假定函数满足定理1.2.1中的条件时，也即微分方程（1.2.8）的解是适定的。 综上，高阶常微分方程初值问题一般为 (1.2.11)其中为给定值。引进新的变量函数 (1.2.12)则初值问题（1.2.11）化成了一阶常微分方程组初值问

8、题 (1.2.13)通过求解（1.2.13）得到（1.2.11）的解1.3 常微分方程常用数值解法的思路一般说来，对于某些典型微分方程可以通过初等积分法求解，而对于复杂微分方程，需要得到解在若干个点上的近似值或者便于计算的近似表达式即可。本论文讲研究常微分方程常用数值解法及其matlab程序设计。对于一阶常微分方程的初值问题为 (1.3.1)求微分方程初值问题（1.3.1）的数值解，就是求解函在一系列离散点上精确值的近似值。由于计算的复杂及繁冗，在计算常微分方程初值问题的数值解是一般使用Matlab软件进行编程计算，按如下步骤：（1）引入点列，其中，称为步长。为了便于使用计算机进行编程计算，一

9、般取步长为定值，即， (2.2)（2）使用常微分方程常用数值解方法，即输入由计算出的递推公式。（3）利用（2）中的公式逐步求出近似解第二章 单步法所谓单步法是指这类方法在计算时，只用到前一步的值，然后逐步往下计算。单步法的精度不高，这是它的缺点。但我们只要给定初值就可以计算，所以还是比较方便的。以下我们介绍两种典型的单步法：Euler法和Runge-Kutta法。2.1 Euler法为了计算出初始值问题 （2.1.1） 在这一区间上若干个点上解的近似解，我们将此区间等分，令为（2.1.1）的解，即当时精确值的近似解。我们的任务是要用尽可能简单的方法求出与尽可能接近的值。Euler法又称Eule

10、r折线法，Euler折线法是数值解法最简单的一种。Euler折线法的基本思想是利用微分中值定理对（2.1.1）的解函数进行近似。由于是可微函数，故 其中，是介于和之间的一个值。当比较小时，和相差不大，故用代替上式中的，就得到了的近似值 当知道以后，再用类似方法来求在点的近似值 （2.1.2）这样就得到了（2.1.1）的解在点的近似解，从几何上看，Euler折线法就是在局部范围内用切线上的值去替代解曲线上的值。Euler公式在实际计算时较少采用，但由于它的结构简单，易于分析，在理论上具有非常重要的意义。2.1.1 向前Euler法Euler折线法的计算量小，但误差大，后来人们对Euler折线法进

11、行了改进：由方程（2.1.1），在节点处成立 （2.1.3）将在处利用Taylor展开式展开，得 其中将将式（2.1.3）代入上式得 （2.1.4） 在式（2.1.4）中略去高阶项 （2.1.5）并用和分别代替式（2.1.4）中的和，可得Euler公式（也称向前Euler公式： （2.1.6）2.1.2 向后Euler法向后Euler法和向前Euler法相差不大，只是将在处利用Taylor展开式展开，得 其中将将式（2.1.3）代入上式得 （2.1.7） 在式（2.1.7）中略去高阶项，并用和分别代替式（2.1.7）中的和，可得向后Euler公式： （2.1.8）其中，局部截断误差为 因此，向

12、后Euler公式同向前Euler公式，为一阶方法。 向前Euler公式（2.1.6）和向后Euler公式（2.1.8）相比，向前Euler公式关于是显式格式。而向后Euler公式中，隐含在方程中，称这样的格式为隐式格式。2.1.3 中点差分公式法若将展开式（2.1.4）和（2.1.7）多展开一项，得 其中将上述两式相减并整理得 其中略去高阶项，并用和代替上式中 和，即得中心差分公式 （2.1.9）其局部截断误差为 即中心差分公式为二阶方法。2.1.4 梯形公式如果将（2.1.1）中微分方程两端从到积分，则得 b 如果用梯形求积公式计算上式中的积分，则有 略去高阶项可得梯形公式： (2.1.10

13、)其中，局部截断误差为 显然，梯形公式为二阶隐式公式。2.1.5 改进Euler法一般来说，隐式格式比显示格式具有较好的数值稳定性。但是，由于在隐式格式里，往往满足的是一个非线性方程，因而需要使用迭代法得到的一个近似值，然后代入隐式格式作校正，并以这个校正值作为的近似值。印版地，可以选择适当的显示格式计算预测值，我们把这样的格式称为预估一校正公式。预估一校正公式往往既便于计算又具有较好的数值稳定性。如，对于梯形公式，可用向前Euler公式计算出的预测值，再用梯形公式修正得到： 上述公式为改进Euler公式。 其总体截断误差为在步长h取得足够小的时侯，我们用改进Euler方法解得数值与理论值非常

14、接近。2.2 Runge-Kutta法Runge-Kutta法是一种间接使用Taylor展开式来获取高阶方法的数值方法，也是一种特殊的单步法，具有相当实用价值。其基本做法是：首先，用函数在附近的个点上的函数值的线性组合来代替的导数。然后，通过Taylor展开式，确定其中的组合系数。这样既可以避免计算高阶导数，又可以获得较高的方法阶。设为一个正整数（代表使用的函数值的个数）， 和 为一些特定参数，则数值方法 (2.2.1)称为级Runge-Kutta公式。2.2.1 二阶Runge-Kutta法显然，公式（2.2.1）是显式的。如果对于该代数方程组的一组解，对应的方法阶为，则得到一个级阶显式Ru

15、nge-Kutta公式。当时，公式（2.2.1）变为 其局部截断误差为 经计算得参数方程组： （2.2.2） 1.若取则对应两二级二阶方法（又称变形Euler公式）： 2. 若取则对应两级二阶方法（又称两级Huen公式）： 2.2.2 三阶Runge-Kutta法1. 三级Kutta公式（三级三阶格式） 2. 三级Huen公式（三级三阶格式） 这两个常用方法在实际使用中能达到较低的精确要求。2.2.3 古典Runge-Kutta法古典Runge-Kutta公式是精度高、最实用的一种常微分方程常用数值解法。 2.2.4 一级隐式中点公式设为一个正整数（代表使用的函数值的个数）， 和 为一些特定参

16、数，则数值方法 (2.2.2) 称为级隐式Runge-Kutta公式。当时，有一级隐式中点公式为 2.2.5 二级隐式中点公式当时，有二级隐式中点公式为 相对显式公式而言，对于同样的，隐式格式不仅可获取较高的精度，而且具有较好的数值稳定性。第三章 多步法 所谓多步法指：这类方法在计算时，除了用到前一步的值，之外，还要用到 ，这前面步的值，这个算法的代表就是阿达姆斯(Adams)方法。 本章介绍线性多步法。线性多步法多使用于求解步骤较多的情况。一般地，步线性多步法公式可写成 （3.1）当时，式（3.1）为显式格式。否则为隐式格式。1. 显式格式首先，将常微分方程（2.1.1）两端在区间上求定积分

17、，得 如果在上面等式中的定积分用关于结点的数值积分近似，可得如下显式格式 其中 而为函数关于结点的Lagrange插值函数的基函数。 当时，四步Adams显式格式 其局部截断误差为 2.隐式格式如果在等式（3.1）中的定积分用关于结点的数值积分近似，可得如下隐式格式 其中 而为函数关于结点的Lagrange插值函数的基函数。当时，三步Adams隐式格式 其局部截断误差为 Adams显式与隐式法相比较而言：同阶数情况下，隐式格式的局部截断误差系数的绝对值比显式格式的小，隐式格式的稳定性相对显式格式的要大。 总结 本论文总结了单步法（向前Euler法，向后Euler法，中心差分法，梯形法，改进Eu

18、ler法，Runge-Kutta法）及线性多步法（Adams法）几种比较常用的常微分方程数值解法，并通过数值算例对同一问题的解的误差计算比较，总结各种算法的优缺点。一般而言，对Euler法而言，求解常微分方程初值问题数值时计算简单、程序简洁并且计算量小，但是，计算精度却很低。中心差分法及梯形法较好地提高计算精确度，同时又保持计算简单，程序简洁的优点。改进的Euler法比Euler法得到的计算结果要好很多。我们用改进Euler方法解得数值与理论值非常接近。对Runge-Kutta法来说，越高阶的方法计算精度越高，但相应的计算量也越大。古典Runge-Kutta公式是精度高、最实用的一种常微分方程

19、常用数值解法。而线性多步法多使用于求解步骤较多的情况。 我们要把常微分方程常用数值解法掌握并灵活运用并加以推广应用，来解决我们的实际生活中所遇到问题，为我国的发展服务，充分体现了常微分方程常用数值解法的重大意义。作为一名数学工作者，我们有责任也有义务继续探索常微分方程组的数值解法，并且在不断探索其解法的前提下，开发出更好更优的软件来解决实际算法问题。但是，其现有理论也还远远不能满足需要，有待于进一步研究，使这门学科的理论更加完善。 数值实例1.分别用改进Euler公式和古典Runge-Kutta公式解常微分方程初值问题取步长，计算结果保存于下表1中，精确解为。解：fun函数： function

20、 f=fun(x,y); f=y-2*x/y;精确解： x=0:0.1:1; y=sqrt(1+2*x)y = Columns 1 through 9 1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 Columns 10 through 11 1.6733 1.73211. 应用改进Euler法，运用matlab软件，程序为：format long h=0.1;x=0:h:1.0;y=zeros(size(x);y(1)=1;y1=sqrt(1+2*x);N=length(x);for n=1:N-1 K1=feva

21、l(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); y(n+1)=y(n)+h*K1;endy2. 应用古典Runge-Kutta法，运用matlab软件，程序为：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y(1)=1;y1=sqrt(1+2*x);N=length(x);for n=1:N-1 K1=feval(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h

22、*K2); K4=feval(fun,x(n)+h,y(n)+h*K3); y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endy表1精确解改进Euler公式古典Runge-Kutta公式0.11.0954451150103321.1000000000000001.0954455316930940.21.183215956619923 1.1918181818181821.1832167455059930.31.2649110640673521.2774378337147221.2649122283403920.41.3416407864998741.3582125995

23、602891.3416423537503710.51.4142135623730951.4351329186577961.4142155778900850.61.4832396974191331.5089662535663321.4832422227719920.71.5491933384829671.5803382376552171.5491964523021430.81.6124515496597101.6497834310477111.6124553496589860.91.6733200530681511.7177793478600871.6733246590162561.01.732

24、0508075688771.7847708324979821.7320563651655652.分别用向前Euler公式、向后Euler公式、中心差分公式和梯形公式解常微分方程初值问题取步长。解：M文件： x=0:0.1:1; y=-exp(x);1.应用向前Euler法，运用matlab软件，程序为：h=0.1;x=0:h:1;y=-exp(x);y1=zeros(size(x);n=length(x);y1(1)=1;for i=1:n-1; y1(i+1)=(1+(h/2)*y1(i);enderror=abs(y-y1)2. 应用向后Euler法，运用matlab软件，程序为：h=0.

25、1;x=0:h:1;y=-exp(x);y2=zeros(size(x);n=length(x);y2(1)=1;for i=1:n-1; y2(i+1)=(1/(1-(h/2)*y2(i);enderror=abs(y-y2)3. 应用中心差分法，运用matlab软件，程序为：h=0.1;x=0:h:1;y=-exp(x);y3=zeros(size(x);n=length(x);y3(1)=1;y3(2)=-exp(h);for i=1:n-1; y3(i+2)=y3(i)+h*y3(i+1);enderror=abs(y-y3)4. 应用梯形法，运用matlab软件，程序为：h=0.1;

26、x=0:h:1;y=-exp(x);y4=zeros(size(x);n=length(x);y4(1)=1;for i=1:n-1; y4(i+2)=(8/(4-h)-1)*y4(i);enderror=abs(y-y4).3.分别用一级隐式中点公式、二级隐式中点公式解常微分方程初值问题取步长，计算结果保存于下表1中，精确解为。解：fun函数： function f=fun(x,y) f=y-2*x/y;1. 应用一级隐式中点公式法，运用matlab软件，程序为：format longx=0:0.1:1.0;y=sqrt(1+2*x);h=0.1;x=0:h:1.0;y1=zeros(siz

27、e(x);y1=sqrt(1+2*x);y1(1)=1;N=length(x);for n=1:N-1 error=1; K1=feval(fun,x(n),y1(n); while(error=1.0e-010) K11=feval(fun,x(n)+0.5*h,y1(n)+0.5*h*K1); error=abs(K11-K1); K1=K11; end y1(n+1)=y1(n)+h*K1;ende=abs(y-y1)2. 应用二级隐式中点公式法，运用matlab软件，程序为：format longx=0:0.1:1.0;y=sqrt(1+2*x);h=0.1;x=0:h:1.0;y2=

28、zeros(size(x);y2=sqrt(1+2*x);y2(1)=1;N=length(x);for n=1:N-1 error=1; K1=feval(fun,x(n),y2(n); K2=feval(fun,x(n),y2(n); while(error=1.0e-010)K11=feval(fun,x(n)+(3-sqrt(3)/6*h,y2(n)+1/4*h*K1+(3-2*sqrt(3)/12*h*K2);K22=feval(fun,x(n)+(3+sqrt(3)/6*h,y2(n)+1/4*h*K2+(3+2*sqrt(3)/12*h*K1); error=max(abs(K1

29、1-K1),abs(K22-K2); K1=K11; K2=K22; end y2(n+1)=y2(n)+0.5*h*(K1+K2);ende=abs(y-y2)4.分别用Adans显式公式和Adans隐式公式解常微分方程初值问题取步长。解：fun函数： function f=fun(x,y); f=-y+x.2;1.应用Adans显式法，运用matlab软件，程序为：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y1=x.2-2*x+4-3*exp(-x);y(1)=1;N=length(x);for n=1:N-1 K1=feval(fun,x(n)

30、,y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(fun,x(n)+h,y(n)+h*K3); y(n+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endfor n=5:N-1 M1=feval(fun,x(n-1),y(n-1); M2=feval(fun,x(n-2),y(n-2); M3=feval(fun,x(n-3),y(n-3); M4=feval(fun,x(n-4),y(n-4); y(n)=y(n-1)+h/24*(55

31、*M1-59*M2+37*M3-9*M4);endy=abs(y-y1)2.应用Adans隐式法，运用matlab软件，程序为：format longh=0.1;x=0:h:1.0;y=zeros(size(x);y1=x.2-2*x+4-3*exp(-x);y(1)=1;N=length(x);for n=1:N-1 K1=feval(fun,x(n),y(n); K2=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K1); K3=feval(fun,x(n)+0.5*h,y(n)+0.5*h*K2); K4=feval(fun,x(n)+h,y(n)+h*K3); y(n

32、+1)=y(n)+h/6*(K1+2*K2+2*K3+K4);endfor n=3:N-1 error=1 M2=feval(fun,x(n),y(n); M3=feval(fun,x(n-1),y(n-1); M4=feval(fun,x(n-2),y(n-2); temp1=y(n); While(error=1.0e-013) M1=feval(fun,x(n+1),temp1); tempt2=y(n)+h/24*(9*M1+19*M2-5*M3+M4); error=abs(temp1-temp2); temp1=temp2;endy(n+1)=temp2;endy=abs(y-y1

33、).参考文献1 戴嘉尊，邱建贤，微分方程数值解.南京：东南大学出版社，2002.2 冯康. 数值计算方法.北京:国防工业出版社，1978.3 周义仓，常微分方程及其应用.北京，科学出版社.2009.4 李瑞遐，何志庆.微分方程数值解法.上海：华东理工大学出版社.2005.5 徐树方，高平，张平文.数值线性代数.北京：北京大学出版社.2000.6 余徳浩，汤华中.微分方程数值解法.北京：高等教育出版社.2003.7 曾金平.微分计算方法.长沙：湖南大学出版社.2006.8 谭浩强，C程序设计.北京：清华大学出版社.2009.9 陆金甫，关治偏微分方程数值解法.北京：清华大学出版社.1987.10

34、 王高雄，周之铭，朱思铭，王寿松.常微分方程.北京：高等教育出版社. 1987.11 胡健伟，汤怀民.微分方程数值方法.北京:科学出版社.2000.12 李信真，车刚明，欧阳洁，封建湖.计算方法.西安:西北工业大学出版社. 2000.13 玛康.数值计算方法.杭州:浙江大学出版社.2003.14 李荣华，玛果忱.微分方程数值方法(第三版).北京:高等教育出版社.1996.15 刘会灯，朱飞.MATLAB编程基础与典型应用.北京:人民邮电出版社.2006.16 苏金明，阮沈勇.MATLAB实用教程（第二版）.北京:电子工业出版社.2008.17 刘卫国.MATLAB程序设计教程.北京:中国水利水电出版社.2005.18 马知恩，周义仓.常微分方程定性与稳定性方法.北京:科学出版社.2007.19 C.W.吉尔.常微分方程初值问题的数值解法,费景高译.北京:科学出版社. 1978.20 楼顺天，于卫，闰华梁.MATLAB程序设计语言西安: