计算机毕业设计（论文）微分方程数值解

1、山 东 英 才 学 院毕 业 论 文 设 计 论 文 题 目： 微分方程数值解 二级学院 ： 计算机电子信息工程学院 学科专业： 计算机及应用 学 号： 姓 名： 班 级： 指导教师： 论文提交时间： 山东英才学院教务处制 2011年 3 月 1 日毕业论文（设计）内容介绍论文（设计）题 目微分方程数值解选题时间2010.11.20完成时间2011.3.9论文（设计）字数33840关 键 词微分方程，周期解，边值问题论文（设计）题目的来源、理论和实践意义： 微分方程是数学科学联系实际问题的主要桥梁之一，它是含有未知函数及其导数的方程。常微分方程的求解是现代科学研究和工程技术中经常遇到的实际问题

2、，然而，从实际问趣中建立出来的微分方程往往具有非常复杂的形式，有些解析式难以计算，有些则根本不能用解析式来表达，所以利用数值解法叫求解实际问题就显得非常重要。论文（设计）的主要内容及创新点： 如果未知函数的自变量是一个，称为常微分方程;自变量多于一个，称为偏微分方程。在科学研究和工程计算中碰到的许多微分方程，根本不存在解析解，或者求解析解的代价很大，求解过程过于复杂，在这种情况下，我们只能借助于数值计算来求方程的数值解。附：论文（设计）本人签名： 年 月 日目 录 中文摘要5第一章 常微分方程的解6第一节 常微分方程的基本概念 6 第二节 常微分方程的12步骤 10第三节 偏导数的方程 14第

3、二章 递增方程的应用17第一节 递增数列17第二节 数列的极限20第三章 与积分有关的数列的极限问题24第一节 积分的应用24第二节 单调定性的松弛法26第三节 松弛算法法的证明33第四章 简单的单步法及基本概念36 第一节 解初值问题的梯形法36第二节 左矩形公式39 第三节 隐式 euler方法 40第4节 预估 校正euler方法 42参考文献45摘要：常微分方程的形成和发展是与力学、天文学、物理学及其他自然科学技术的发展互相促进和互相推动的。本文第一章讲述了常微分方程的发展历史，第二章介绍了一系列常微分方程的周期解和边值问题，说明了其研究现状，第三章举例说明了其在生态学和军事上的应用。

4、无论在数学研究还是在自然科学以及其他应用科学，常微分方程都显现出其重要的理论和应用价值。随着科学技术的发展和社会进步，常微分方程的理论和应用不断扩大和深入，其作用也越来越被人们所重视。在数学应用方面，它有着比通常导数更广泛的应用，对于导数不存在而对称导数存在的函数，我们就可以用对称导数研究此类函数的一些重要性质. 常微分方程研究的内容包括解的基本性质（如存在性、惟一性等）、解的解析表达式或近似的解析表达式、解的定性性质以及解的数值解法。常微分方程的形成和发展是与力学、天文学、物理学及其他自然科学技术的发展互相促进和互相推动的。数学的其他分支的新发展如复变函数、李群、组合拓扑学等都给常微分方程的

5、发展以深刻的影响。关键词：微分方程，周期解，边值问题； 一 常微分方程的解（一）：常微分方程的基本概念1常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求函数，满足其中是已知函数，是已知值。假设在区域上满足条件：（1）在上连续；（2）在上关于变量满足lipschitz条件：， （1.3） 其中常数称为lipschitz常数。我们简称条件（1）、（2）的基本条件。由常微分方程的基本理论，我们有：定理1 当在上满足基本条件时，一阶常微分方程初值问题（1.1）、（1.2）对任意给定存在唯一解在上连续可微。定义1 方程（1.1）、（1.2）的解称为适定的，若存在常数和，对任意满足条件及的和，

6、常微分方程初值问题 （1.4）存在唯一解，且适定问题的解连续依赖于（1.1）右端的和初值。由常微分方程的基本理论，还有：定理2 当在上满足基本条件时，微分方程（1.1）、（1.2）的解是适定的。我们在本章中假设在上满足基本条件，从而（1.1）、（1.2）的解存在且适定。一般的一阶常微分方程组初值问题是求解 （1.5）（15）的向量形式是 （1.5）其中记。类似于定理1和定理2，我们有：定理3 若映射满足条件(1) 在上是从到上的连续映射;(2) 在上关于满足lipschits条件;任意。则常微分方程组初值问题（1.5）存在的唯一的连续可微解而且解是适定的。高阶常微分方程初值问题一般为 （1.6

7、）其中是给定多元函数，为给定值。引进新的变量函数 （1.7.）则初值问题（1.6）化成了一阶常微分方程组初值问题通过求解（1.8）得到(1.6)的解。2.初值问题数值解基本概念初值问题的数值解法，是通过微分方程离散化而给出解在某些节点上的近似值。在上引入节点称为步长。在多数情况下，采用等步长，即。记（1.1）,(1.2)的为准确解为,记的近似值为,记为.。 求值问题数值解的方法是步进法，即在计算出后计算。数值的方法有单步与单步法之分。单步法在计算时只利用而多步法在计算时不仅要利用还要利用前面已算出的若干个。我们称要用到的多步法为步方法。单步法可以看作多步法，但两者有很大差别。步方法只能用于的计

8、算，要用其它的方法计算；而且在稳定性上单性法比的多步法容易分析；此外单步法容易改变步长。单步法和多步法又都有显式方法和稳式方法之分。单步显式法的计算公式可写成 （1.9）隐式单步法的计算公式可写成 （1.10）在(1.10)中右端项显含。从而（1.10）是的方程式，要通过解方程求出。 显式多步法计算公式为 （1.11） 而隐式多步法计算公式为 (1.12)右端项含。多步法中一类常用方法是线性多步法 （1.13）其中是独立于和的常数。时（1.13）是显式的，时是隐式的。（二）：常微分方程的12步骤我们常分12步，来说明微分方程的解法。包括一个自变量和它的未知函数以及未知函数的微商的等式。运动着的

9、物体的位置(x,y,z)是随时间t的变化而变化的。按牛顿运动定律,力等于质量乘加速度,它牵涉到x、y、z对于t的二阶微商, 这就可用常微分方程表示。例如，在大地上的自由落体运动，可以用下面的常微分方程来描述： (1)式中g是重力加速度，z是铅直位置。 一般说来,如果y是自变量x的函数，则y的常微分方程可以表达为 (2) 式中f是它所依赖的n+2个变量的函数,n为正整数。由自变量x的n个未知函数y1,y2,ym的m个常微分方程 (3) 所形成的一组方程称为常微分方程组，其中n1,n2,，nm为非负整数。如果一个常微分方程（组）关于所有未知函数及其各阶微商都是线性的，则称为线性常微分方程（组）；否

10、则，称为非线性常微分方程（组）。如果能由(2)解出最高阶微商，则得到 (4) 式中是它所依赖的n+1个自变量的函数。这种就最高阶微商解出的微分方程,称为正规型微分方程；而称(2)为隐微分方程。任一正规型微分方程(4)与微分方程组 是等价的，因此(4)总可以化成一个与之等价而形如(5) 的正规型方程组。对于微分方程组(3),也有上述相似的结果，即任一正规型微分方程组也可化为等价的而形如(5)的正规型方程组。 满足常微分方程的函数称为常微分方程的解，也就是说,对方程(2),如果有函数(x),在x轴的某区间i上有定义,具有从1阶到n阶的微商且满足 对所有xi,则称(t)为方程(2)在区间i上的解。当

11、前计算机的发展为常微分方程的应用及理论研究提供了非常有力的工具。质点动力学是这个阶段研究的问题的主要来源之一。例如牛顿建立了太阳系行星运动方程 (6) 并求出其通解的显式解析表达式。这里t是时间,r(t)=(x(t)，y(t)，z(t)是以太阳为原点的直角坐标系中行星的位置，g是万有引力常数，m是太阳质量。 这个阶段主要是求常微分方程的通解，亦即对于方程(2)求含有n个任意常数c1,c2,cn的形如 (7) 的解。对于方程组(3)，则要求含有n1+n2+nm个任意常数的解组。 这个阶段的成果有：g.w.莱布尼茨关于齐次方程和线性方程的通解；雅各布第一·伯努利提出并解决、现命名为伯努利

12、方程的特殊非线性方程 (8) l.欧拉等得到的常系数线性常微分方程的通解；以及利用变换xet将欧拉方程 (9) 化为常系数线性方程。 但是,求显式通解的可能性十分有限,上述努力经过一段时间便停滞下来。当时若干实际问题又迫切需要解决,类似地，1841年刘维尔证明了下述的黎卡提方程 (10) 有且只有在v为非负整数时才有“初等解”,亦即经过有限次的初等运算求得的显式解，从而结束了一般常微分方程求通解的企图。极限环，即孤立周期解，已成为许多实际问题的核心，例如在无线电技术中的范德坡方程 (11) 便是非线性方程产生孤立周期振荡的典型例子。奇点大范围分布的研究与组合拓扑学的研究在庞加莱的工作中是双生的

13、兄弟关系。例如工程控制论中火箭发动机中燃烧过程由于时滞现象而产生的带有时滞的常微分方程或称微分差分方程 (12) 这里不仅依赖于t时的(t)，而且还依赖于(t-)时的(t-)和(t-)。这样便产生了微分差分方程，以及更广义的泛函微分方程。又如由于空气中的湍流对于飞行中飞机机翼所引起的干扰对飞机运动的影响，这里常微分方程带有随机摄动项，如 式中x是随机输入,是反馈参数。这类问题产生了常微分方程和概率论之间的一个新分支随机微分方程。 （三）：偏导数的方程含有未知函数及其各阶偏导数的方程。如(余类此)uta2（uxxuyyuzz）0（1）其中uu（x，y，z，t）为未知函数 ，x ，y，z，t是自变

14、 量。18世纪 ，数学家们已开始用偏 微分方程来研究问题 。方程（1）便是用来描述热的传导规律的。1746年 ，j.ler.达朗贝尔给出了一维波动方程（两端固定的弦的振动问题）： 由于弦的两端固定，故在x0和xl处（l为弦的长度）应满足边界条件： u（0，t）0 u（l，t）0 t0（3）又当t0时的状态，即初始条件是 u（0，x）j（x） ut（0，x）（x）（4） 例 设,.证明:数列的极限存在并求出此极限.例1可以作如下推广:命题1 若,则数列的极限存在且为.证明 由知.由算术几何平均不等式知,假设,再次用算术几何平均不等式知,由数学归纳法知,对任意正整数均有,因而数列有界.又当时,故,

15、即数列单调递增.由数列的单调有界定理知存在,设为,对两边同时取极限得:,可解得或（舍去）.故.注 由命题1立得例1的极限存在且为.例 证明数列收敛,其中,并求极限.通过观察、猜想、分析可将例2推广为以下更一般的形式:命题2 若,定义,则数列存在极限且为.证明 由可知,当且仅当时取等号.设,则=,当且仅当时取等号.二 递增方程的应用(一)：递增数列由数学归纳法可知,对一切自然数均有成立.所以数列是单调递增且有上界的数列.（i）当时,同理可证数列是单调递减且有下界的数列.由（i）可知数列是单调有界数列,从而命题3得证.注 由以上第一章的命题可知数列是单调有界数列,则必收敛,设,对两边同时取极限的得

16、: 即 .所以是方程的一个正根,从而例3得证.例1 已知,.证明存在并求其值.例1可以作如下推广:命题1 若,则数列极限存在且为的正根.证明 由得.又,则.由递推关系知.因函数是递增函数,则由知 与 的符号相同.而 的符号又与 的符号相同,故依次下去便知最终与 的符号相同.而,即,所以 ,从而 ,于是便有,故数列是单调递增数列.又,假设当时都有 成立,则当 时,由数学归纳法知,对一切自然树都有 ,即数列有界.由数列的单调有界定理知数列必存在极限,设,对两边同时取极限的得 即 .所以数列收敛于方程的正根.例2用经典四阶r-k方法计算初值问题 步长取h=0.1及0.2,给出计算误差并分析其稳定性.

17、解 本题直接按r-k方法的公式计算.因精确解为，其计算误差如表所示.图2-1从计算结果看到，h=0.2时误差很大，这是由于在=-20,h=0.2时h=-4，而四阶r-k方法的绝对稳定区间为-2.785,0,故h=0.2时计算不稳定，误差很大.而h=0.1时=-2，其值在绝对稳定区间-2.785,0内，计算稳定，故结果是可靠的.(二)：数列的极限推论 若, ,则数列的极限存在且收敛于方程的一个正根,即 .注 利用该推论易知例4中的数列的极限存在且为3.1文4、5中的一些题也可由此推论直接得出.例1设,数列、分别定义为, ,证明.例1可以作如下推广:命题 设,数列、分别定义为, ,则仍有 成立.证

18、明 由已知 可得 ,.假设当 时均有,则当时有, ,由数学归纳法知,对任意自然数都有,成立.由算术几何平均不等式知,当且仅当时取等号.而当时, 即,即.故有.而当时,由于,所以就有, ,因此对任意自然数都有下式成立,所以数列、均为单调有界数列.故由数列的单调有界定理知、存在,分别设为、,对两边同时取极限得,可解得 ,即.注 例5是命题5的特殊形式,证明类似.通过以上这简单的五个例子很容易看出它们是各自推广后命题参量的特殊值,还有好多题都可直接根据这些命题很快得出其极限值.例2将例1的初值问题用修正的milne-hamming预测-校正公式计算及，初值，仍用已算出的精确解，即，给出计算结果及误差

19、 从结果看，此方法误差比四阶adams隐式法和四阶hamming方法小，这与理论分析一致. 讲解：线性多步法的局部截断误差定义为与单步法相似，可表示为三：与积分有关的数列的极限问题(一)：积分的应用本节主要是就一些与积分有一定联系的数列举了两个例子,判断它们的方法也主要是单调有界定理,进一步证明了单调有界定理在判断数列收敛即极限存在问题中的应用.例1证明 存在.证明 设,可知在上非负单调递减,所以,即,亦即,所以数列有界.又= =0,即,所以数列是单调递减的数列.由数列的单调有界定理知,数列的极限存在,也就是存在.例2设连续函数在上是正的单调递减的函数,且,证明数列收敛.证明 由假设及积分中值

20、定理有 =,其中.又在上是正的单调递减的,所以,从而,即数列是单调递减的数列.又=,其中,所以.又因为,所以,所以单调递减且有下界.所以由单调有界定理知存在,即数列收敛.(二)：单调定性的松弛法通过本文可以看到单调有界定理在求一些有特殊形式的极限时常常很有效.解边值问题的松弛法     方法：     public   class   d16r1f         void   sor(double   a,   double   b,   do

21、uble   c,       double   d,   double   e,   double   f,       double   u,   int   jmax,   double   rjac)         int   n,   j,   l;     double   maxits,   eps,

22、  zero,   half,   qtr,   one,   anormf,   omega,   aaa,   bbb;     double   resid,   anorm;                     maxits   =   1000;              

23、      eps   =   0.00001;                     zero   =   0.0;                     half   =   0.5;                 &#

24、160;   qtr   =   0.25;                     one   =   1.0;                     anormf   =   zero;                    

25、; for   (j   =   2;   j   <=   jmax   -   1;   j+)                             for   (l   =   2;   l   <=   jmax   -   1;   l+)  

26、                                   anormf   =   anormf   +   math.abs(fjl);                     omega   =   one;       &

27、#160;             for   (n   =   1;   n   <=   maxits;   n+)                                                 ano

28、rm   =   zero;                             for   (j   =   2;   j   <=   jmax   -   1;   j+)                      

29、60;                                         for   (l   =   2;   l   <=   jmax   -   1;   l+)                

30、60;                                                               if   (j   +   l)   %   2)   =(   n   %

31、  2)                                                                                  

32、60;             aaa   =   ajl   *   uj   +   1l   +   bjl   *   uj   -   1l;                                      

33、              bbb   =   cjl   *   ujl   +   1   +   djl   *   ujl   -   1;                                     

34、0;               resid   =   aaa   +   bbb   +   ejl   *   ujl   -   fjl;                                         &

35、#160;           anorm   =   anorm   +   math.abs(resid);                                                     ujl   =

36、60; ujl   -   omega   *   resid   /   ejl;                                                                

37、60;                                                                       if   (n   =   1   )   

38、0;                                 omega   =   one   /   (one   -   half   *   rjac   *   rjac);                    

39、60;       else                                     omega   =   one   /   (one   -   qtr   *   rjac   *   rjac*   omega);      

40、;                       if   (n   >   1)   &&   (anorm   <   eps   *   anormf)   )   return;                        

41、                 system.out.println(“   maxits   exceeded”);                     system.exit(1);                      例子：     import &#

42、160; java.text.*;     public   class   d16r1         public   static   void   main   (string   args)         /program   d16r1     /driver   for   routine   sor     int  

43、jmax,   i,   j,   midl;     double   pi,   rjac,   aaa;     double   a   =   new   double1212;     double   b   =   new   double1212;     double   c   =   new   double121

44、2;     double   d   =   new   double1212;     double   e   =   new   double1212;     double   f   =   new   double1212;     double   u   =   new   double1212;     d16r1f

45、  g   =   new   d16r1f();     decimalformat   form   =   new   decimalformat(“0.00”);     jmax   =   11;     pi   =   3.1415926;     for   (i   =   1;   i   <=   jmax;

46、   i+)                                 for   (j   =   1;   j   <=   jmax;   j+)                            

47、             aij   =   1.0;                                     bij   =   1.0;                       

48、0;             cij   =   1.0;                                     dij   =   1.0;                       &#

49、160;             eij   =   -4.0;                                     fij   =   0.0;                      

50、               uij   =   0.0;                                     midl   =   jmax   /   2   +   1;     fmidlmidl   =

51、60; 2.0;     rjac   =   math.cos(pi   /   jmax);     g.sor(a,   b,   c,   d,   e,   f,   u,   jmax,   rjac);     system.out.println();     system.out.println(“sor   solution:”);     sys

52、tem.out.println();     for   (i   =   1;   i   <=   jmax;   i+)                 for   (j   =   1;   j   <=   jmax;   j+)               

53、0;             system.out.print(form.format(uij)   +   “     “);             system.out.println();         system.out.println();     system.out.println(“test   that   sulotion 

54、0; satisfies   difference   eqns:”);     system.out.println();     for   (i   =   2;   i   <=   jmax   -   1;   i+)         for   (j   =   2;   j   <=   jmax   - &#

55、160; 1;   j+)         aaa   =   ui   +   1j   +   ui   -   1j   +   uij   +   1   +   uij   -   1;     fij   =   aaa   -   4.0   *   uij;    

56、60;   for   (j   =   2;   j   <=   jmax   -   1;   j+)     system.out.print(form.format(fij)   +   “     “);     system.out.println();                运行结果：   

57、0; sor   solution:         0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00         0.00     -0.02

58、    -0.04     -0.06     -0.08     -0.09     -0.08     -0.06     -0.04     -0.02     0.00         0.00     -0.04     -0.09     -0.13     -0.17   &

59、#160; -0.19     -0.17     -0.13     -0.09     -0.04     0.00         0.00     -0.06     -0.13     -0.20     -0.28     -0.32     -0.28     -0.20     -0

60、.13     -0.06     0.00         0.00     -0.08     -0.17     -0.28     -0.41     -0.55     -0.41     -0.28     -0.17     -0.08     0.00       

61、0; 0.00     -0.09     -0.19     -0.32     -0.55     -1.05     -0.55     -0.32     -0.19     -0.09     0.00         0.00     -0.08     -0.17     -0.28

62、    -0.41     -0.55     -0.41     -0.28     -0.17     -0.08     0.00         0.00     -0.06     -0.13     -0.20     -0.28     -0.32     -0.28   &

63、#160; -0.20     -0.13     -0.06     0.00         0.00     -0.04     -0.09     -0.13     -0.17     -0.19     -0.17     -0.13     -0.09     -0.04     0.

64、00         0.00     -0.02     -0.04     -0.06     -0.08     -0.09     -0.08     -0.06     -0.04     -0.02     0.00         0.00     0.00    

65、 0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00             test   that   sulotion   satisfies   difference   eqns:      

66、;   0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00         0.00     -0.00     0.00     -0.00     0.00     -0.00

67、60;   0.00     -0.00     0.00         0.00     0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00     0.00         0.00     -0.00 

68、0;   0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00         0.00     0.00     -0.00     0.00     2.00     0.00     -0.00     0.00

69、    0.00         0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00         0.00     0.00     -0.00     0.00 &

70、#160;   -0.00     0.00     -0.00     0.00     0.00         0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00     -0.00     0.00     &

71、#160;   0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00     0.00(三)：松弛算法法的证明证明 令，由已知条件，不妨设在上严格增加，于是对，必有，使得，且，即（1） 若，则结论成立. （2） ，则令由于-=所以在上连续，又.，故由连续函数介值定理知存在，使得，即.再记即得结论.(3)若，则令用与情形（2）类似的方法可证得存在，

72、使得. 即.再记，即得结论. 当（2）中的时，即得推论1 设函数在内单调连续，则对，使得 .令（2）式中的即得推论2 设函数在内单调连续，在内连续，,且不变号，则对，使得.推论3 若ca,b且，则.定理：设函数在上有阶连续导数，且函数在上有阶连续导数，且，若是由（1）式所确定的，则有四 简单的单步法及基本概念（一）：解初值问题的梯形方法若在(3.4)的积分中用梯形公式，则得(3.6)称为梯形方法.上述三个公式(3.2)，(3.5)及(3.6)都是由计算，这种只用前一步即可算出的公式称为单步法，其中(3.2)可由逐次求出的值，称为显式方法，而(3.5)及(3.6)右端含有当f对y非线性时它不能直

73、接求出，此时应把它看作一个方程，求解，这类方法称为稳式方法.此时可将(3.5)或(3.6)写成不动点形式的方程这里对式(3.5)有，对(7.2.6)则，g与无关，可构造迭代法(3.7)由于对y满足条件(3.2)，故有当或，迭代法(3.7)收敛到，因此只要步长h足够小，就可保证迭代(3.7)收敛.对后退euler法(3.5)，当时迭代收敛，对梯形法(3.6)，当时迭代序列收敛.例3.1用euler法、隐式euler法、梯形法解取h=0.1，计算到x=0.5，并与精确解比较.解 本题可直接用给出公式计算.由于，euler法的计算公式为n=0时，.其余n=1,2,3,4的计算结果见表7-1.对隐式euler法，计算公式为解出当n=0时，.其余n=1,2,3,4的计算结果见表7-1.表4-2的三种方法及精确解的计算结果 图4-2对梯形法，计算公式为证明 作辅助函数 一方面，由已知条件及法则有= （1）另一方面，又