《常微分方程数值解及其MATLAB实现(共31页).doc》由会员分享，可在线阅读，更多相关《常微分方程数值解及其MATLAB实现(共31页).doc（31页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上茂 名 学 院毕 业 论 文题 目 常微分方程数值解及其MATLAB实现 英文并列题目Numerical Solution of Ordinary Differential Equations and MATLAB Implementation 学院 理学院 专业 数学与应用数学（师范）班级 数学05-1班 学生 李尧光 指导教师（职称） 李伟勋 （副教授） 完成时间 2009 年 1 月 15 日至 2009 年 6 月 10 日 系主任林全文批准日期茂 名 学 院 毕 业 论 文 任 务 书 数学 系 数学与应用数学（师范）专业 数学05-1班 学生 李尧光 一、
2、毕业论文课题 常微分方程数值解及其MATLAB实现 二、毕业论文工作自 2009 年 1 月 15 日起至 2009 年 6 月 10 日止三、毕业论文进行地点 茂名学院图书馆、学生宿舍 四、毕业论文的内容要求 1内容要求 常微分方程数值解及其MATLAB实现主要介绍常微分方程初值问题的数值解法，探讨了采用单步法求解常微分方程初值问题的数值解，并运用MATLAB进行编程求解。 2过程要求 （1）4月10日前按要求组织所需资料，完成提纲、摘要。 （2）5月26日前按要求查阅相关文献、撰写论文初稿。 （3）6月9日前按要求完成论文终稿、打印装订以及答辩前的准备工作。 （4）6月10日至14日为小组
3、答辩，15日至16日为公开答辩。 （5）论文的语言文字要规范，格式要按规定的格式进行编辑，字数在1.2万至1.5万之间，参考文献至少10项,至少1篇外文文献。 （6）在整个毕业论文写作过程中要经常与指导老师交流。 教研室负责人 江蓉 指导教师 李伟勋 （副教授） 接受设计论文任务开始执行日期 2009 年 1 月 16 日 学生签名 专心-专注-专业摘 要科学技术和工程中大量的问题都表达为常微分方程的形式，然而，只有少数十分简单的微分方程能够用初等方法求解，在很多情况下，我们只能求得解的近似值，因此，就有研究常微分方程数值解显得非常重要。本文介绍常微分方程初值问题的数值解法，探讨了采用单步法求
4、解常微分方程初值问题的数值解，并运用MATLAB进行编程求解，程序简洁、直观, 求解速度快, 方法实用性较强。关键词：微分方程；数值解法；MATLABAbstractScience and technology and engineering problems are mostly expressed as a form of ordinary differential equations. However, only a few very simple differential equations can be solved with elmentary method. In most ca
5、ses, we can only obtain approximate answers. Therefore, to take a research of differential equation numerical is obviously important. This thesis is about introduction of numerical method for elementary value problem of ordinary differential equations and to talk about applying the single step metho
6、d to solve the elementary value problem of ordinary differential equations, as well as applying MATLAB programming which is concise, visual, solving the problems quickly and practical.Key words: differential equations; numerical method of equations; MATLAB目录摘 要IAbstractII第一章 绪言11.1 关于微分方程11.2 关于MATL
7、AB21.3 文献综述3第二章 微分方程导论5第三章 常微分方程数值解法73.1 概念73.1.1 单步法73.1.2 单步法的离散误差与阶73.2 欧拉方法83.2休恩方法93.3泰勒级数法103.4龙格库塔方法12xcvfsdf3.4.1 的龙格库塔方法123.4.2 的龙格库塔方法153.4.3 龙格-库塔-费尔伯格方法16第四章 常微分方程数值解法的MATLAB实现204.1 欧拉方法204.2 休恩方法214.3 4次泰勒方法224.4 4阶龙格库塔方法234.5 龙格-库塔-费尔伯格方法24结束语26参考文献27致 谢28第一章 绪言1.1 关于微分方程常微分方程是数学的一个重要分
8、支，也是偏微分方程、变分法、控制论等数学分支的基础。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力的工具。在1718世纪，在力学、天文、物理和技术科学中，就已借助微分方程取得了巨大的成就。质点动力学和刚体动力学的问题很容易化微分方程的求解问题。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为微分方程的进一步提出了新的问题，促使人们对微分方程进行更深的研究，以便适应科学技术飞速发展的需要。微分方程的首要问题是如何求一个给定的方程的通解和特解。到目前为止，人们已经对许多微分方程得出了求解的一部方法。例如一阶微分方程中的变量可分离的方程、线性方程、恰当方程、以
9、及常系数二阶线性常系数方程组等。求一个方程的解，最自然的想法是用初等函数来表达方程的解，这是不容易做到的，例如对非常简单的方程，我们就无法用初等函数来表示它的解，但我们可以用初等函数积分形式来表达该微分方程的解，即将该方程的通解表达为微分方程的解也称为积分，求解过程称为对这个微分方程进行积分，习惯上也称为初等积分法。应该看到，能用初等积分法求解的微分方程为数很少，绝大部分的微分方程都无法求出通解。由于求通解存在很多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始值问题解的存在惟一性进行了研究。目前解答存在惟一性、延拓性、大范围的存在性以及解对初始值和参数的连续性和可微
10、性等理论问题都已发展成熟。从而推动人们从其他方面来研究微分。例如，采用将未知函数表示成一致收敛的级数形式，从而扩大了微分方程的可解领域。以后，人们引进新的特殊函数（非初等函数），例如椭圆函数、阿贝尔函数、贝塞尔函数、勒让得函数等来表达微分方程的解。使微分方程和函数论，特别是和复变函数论紧密地联系起来，产生了微分方程的解析理论。与此同时，人们开始采用各种近似方法来求微分方程的特解，例如用函数近似的逐次比较法、Taylor级数法、待定系数法，都可以在一个区间上求近似解。又如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程解的近似值。今年来随着电子计算机的飞
11、速发展，开发出了许多功能强大、使用方便的软件包，在微分方程的求解和应用之中发挥了巨大的作用，真正使微分方程在科学技术和经济发展中得到了充分的应用，解决了许多重大问题。1.2 关于MATLABMATLAB软件是由美国MathWorks公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB是英文MATrix LABoratory(矩阵实验室)的缩写。它的第1版(DoS版本1.0)发行于1984年，经过20余年的不断改进，它的功能正在不断的完善。MATLAB集中了日常数学处理中的各种功能，包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。在MATLAB环境下，用户可以集成地进行程序设
12、计、数值计算、图形绘制、输入输出、文件管理等各项操作。MATLAB提供了一个人机交互的数学系统环境，该系统的基本数据结构是矩阵,在生成矩阵对象时,不要求作明确的维数说明。与利用c语言或FORTRAN语言作数值计算的程序设计相比，利用MATLAB可以节省大量的编程时间。MATLAB系统最初是由CIeve Moler教授用FORTRAN语言设计的，现在的MATLAB程序是MathWorks公司用C语言开发的，第一版由steve Bangert主持开发编译解释程序，Steve Kleiman完成图形功能的设计，John Little和Cleve Moler主持开发了各类数学分分析的子模块，撰写用户指
13、南和大部分的M文件。自从第1版发行以来，已有众多的科技工作者加入到MATLAB的开发队伍中，并为形成今天约MATLAB系统做出了巨大的贡献。MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性，使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。在国际学术界，MATLAB已经被确认为准确、可靠的科学计算标准软件。在许多国际一流学术刊物上，（尤其是信息科学刊物），都可以看到MATLAB的应用。MTAL
14、AB系统由五个主要部分组成。(1)MATALB语言体系MATLAB是高层次的矩阵数组语言。具有条件控制、函数调用、数据结构、输入输出、面向对象等程序语言特性。利用它既可以进行小规模端程，完成算法设计和算法实验的基本任务，也可以进行大规模编程，开发复杂的应用程序。 (2)MATLAB工作环境 这是对MATLAB提供给用户使用的管理功能的总称包括管理工作空间中的变量据输入输出的方式和方法，以及开发、调试、管理M文件的各种工具。 (3)图形句相系统 这是MATLAB图形系统的基础，包括完成2D和3D数据图示、图像处理、动画生成、图形显示等功能的高层MATLAB命令，也包括用户对图形图像等对象进行特性
15、控制的低层MATLAB命令，以及开发GUI应用程序的各种工具。(4)MATLAB数学函数库这是对MATLAB使用的各种数学算法的总称包括各种初等函数的算法，也包括矩阵运算、矩阵分析等高层次数学算法。 (5)MATLAB应用程序接口(API) 这是MATLAB为用户提供的一个函数库，使得用户能够在MATLAB环境中使用c程序或FORTRAN程序，包括从MATLAB中调用于程序(动态链接)，读写MAT文件的功能。可以看出MATLAB是一个功能十分强大的系统，是集数值计算、图形管理、程序开发为一体的环境。除此之外，MATLAB还具有根强的功能扩展能力，与它的主系统一起，可以配备各种各样的工具箱，以完
16、成一些特定的任务。用户可以根据自己的工作任务，开发自己的工具箱。1.3 文献综述爱恩斯坦曾经说过，数学受到高度尊崇的另一个原因在于：恰恰是数学，给精密的自然科学提供了无可置疑的可靠保证，没有数学，它们无法达到这样的可靠程度。在自然科学的很多领域中，都会遇到常微分方程初值问题，研究常微分方程的数值解法，借助数学软件MATLAB对解法进行实现，对于生产实际问题有着十分重要的意义。在写这篇论文时我阅读了大量的相关文献，这使得我对常微分方程初值问题的数值解有了比较深刻的认识，懂得了运用MATLAB,以下是对这些文献作简要的分析。李庆杨的数值分析是我参考的第一篇文献，该书介绍了多种常微分方程初值问题的数
17、值解法，如：欧拉法、龙格-库塔方法等，对这些数值方法给出了理论基础。由此可见数学的严谨性，也启发本人写这篇论文的思路，掌握常微分方程初值问题的数值解法，必须了解常微分方程初值问题的数值解的概念，各种解法的构造、意义及其理论依据。周义仓的常微分方程及其应用将常微分方程理论、方法与应用有机结合，理论、方法和应用互相渗透与补充，体现出数学在解决实际问题中的巨大作用，并通过实际的例子，使得各种解法更加直观。而王兵团的数学实验基础侧重于介绍各种方法的构造过程，这有助于本人后面的程序编写。在利用MATLAB编写程序时，我也查阅不少的文献。胡亮剑的MATLAB数学实验系统的介绍了MATLAB的基础知识和操作
18、，加深了我对MATLAB的认识。邓黎的常微分方程数值解及其Matlab实现说到的利用内嵌程序和根据数值解法的思想编程两种实现方法，令我对本论文的MATLAB的编程有了基本的构思。特别是姜健飞的数值分析及其MATLAB实验详细的介绍了基于MATLAB命令的常微分方程初边值问题的数值解法，开阔了我的思维，指导我对常微分方程数值解法进行编程。所查阅的资料文献，为我写作常微分方程数值解及其MATLAB奠定了扎实的基础。但本论文虽然经过多次的修改，但由于水平有限，论文肯定会存在不足之处，甚至是手误。第二章 微分方程导论数学分析中所研究的函数，是反映客观现实世界运动过程中的量与量之间的一种关系。但在大量的
19、实际问题中遇到稍微复杂的一些运动过程时，反映运动规律的量与量之间的关系（即函数）往往不能直接写出来，却比较容易地建立这些变量和它们的导数（或微分）间的关系式，这种联系着自变量、未知函数及它的导数（微分）的关系式，数学上称之为微分方程，如果在微分方程中，自变量的个数只有一个，我们称这种微分方程为常微分方程。定义1 初值问题（initial value problem,简称L.V.P.） （1）在区间上的解是一个可微函数，使得 注意解曲线必须过初始点。几何解释 在矩形区域中的每一点处，解曲线的斜率都可由隐式公式得到。因此，整个矩形区域中的都可以计算出来，并且每个值表示与过点的解曲线相切的直线的斜率
20、。定义2 给定矩形，设在上连续。对任意，如果存在常数具有 则称在上的变量满足利普希茨（Lipschitz）条件。常数称为的利普希茨常数。定理1 设定义在区域上，如果存在一个常数,使得 则在区域上的变量满足利普希茨条件，利普希茨常数为。定理2（存在性和唯一性） 设在区域上上连续，如果在区域上的变量满足利普希茨条件，并且,则初值问题（1），在某个子区间上有唯一解。上述定理从理论上解决初值问题（1）的解的存在唯一性但是只有少数特殊类型的初值问题能得到解的解析表达式，而大量的实际问题得到的常微分方程初值问题只能用数值方法计算其数值解常微分方程初值问题的数值解是指：对自变量的变化区间,在其一系列离散节点
21、上，计算解在的值的近似值，。称为初值问题（1）的数值解计算初值问题的数值解的方法称为数值解法。下面的章节将对常微分方程初值问题的数值解法进行探讨, 并运用MATLAB进行编程实现常微分方程初值问题的数值解。第三章 常微分方程数值解法在自然科学的很多领域中，都会遇到常微分方程初值问题，然而，只有少数十分简单的微分方程能够用初等方法求解，生产实际问题中所遇到的常微分方程往往具有相当复杂的形式有很多不能把它的解用解析表达式表示出来，或者即使能表示出来，它的表达式也拄拄非常复杂而不切实用况且生产实际中提出的问题往往也只要在一定的范围内求出微分方程曲解在某些点上的数值，并不一定要求出它的解析表达式。因此
22、，就有研究常微分方程数值解的必要。3.1 概念由于本文研究的是用单步法求常微分方程的数值解，因此先介绍相关的概念。3.1.1 单步法在计算之时只用到的方法，其计算公式有：显式单步法计算公式 隐式单步法计算公式 这里引入的求初值问题近似解的方法称为差分方法或离散变量法。它求离散点集上的近似解，这个离散点集称为网格（或格网）。式中的称为增量函数。3.1.2 单步法的离散误差与阶设是初值问题的唯一解，是离散近似解集。全局离散误差定义为 其中 它是唯一解与离散方法得到的解之间的差。 局部离散误差定义为 其中 它是从到这一步计算的误差。设是初值问题的唯一解，是离散近似解集，若存在最大整数使显式单步法的局
23、部离散误差满足 则称该方法具有阶精度。若将上式展开写成 则称为局部离散误差主项。3.2 欧拉方法 欧拉方法是最简单的求微分方程数值解的方法，这种方法用来体现先进方法包含的概念。由于它随着步数增加而产生的累积误差较大，因此用途有限，然而研究它也很重要，因为它的误差分析更易懂。设为求解良态初值问题的区间。实际上，下面的过程不是要找到满足该初值问题的可微函数，而是要生成点集，并且将这些点作为近似解，即。如何构造“近似满足微方程”呢？首先为这些点选择横坐标，为方便起见，将区间划分为个等距子区间，并选择网格点 其中 （1）其中值称为步长。然后近似求解 在上， （2）设和连续，利用泰勒定理将在处展开，对每
24、个值，存在一个和之间的值，使得 （3）将和代入等式（3），得到的表示： （4）如果步长足够小，则可以忽略次项（包含的项），得到 （5）这就是欧拉近似。重复该过程，就能得到近似解曲线的一个点序列。欧拉方法的一般步骤是 ，其中 （6）欧拉方法的精度 设是初值问题（2）的解，如果且是由欧拉方法计算的近似值序列，则 （7）所以，欧拉方法具有阶精度。3.2休恩方法为得到比欧拉法精度高的计算公式，下面研究休恩（heun）方法。休恩方法引入一种新的思路，来构造求解上的初值问题 （1）要得到解，可以用微积分基本定理，在上对积分得 （2）其中的不定积分为待求函数，对求解方程（2），结果为 （3）然后可用数值积分
25、方法逼近（3）中的积分，如果采用步长为的梯形公式，则结果为 （4）注意公式（4）的右端包含了待定值。要继续求解，需要的一个估计值，欧拉方法的解能够满足这一目的，将它代入（4）后，得到求解的公式，称为休恩（heun）方法。 （5）重复这个过程，得到逼近解曲线的一系列点，在每一步中都用欧拉方法作为预报，然后用梯形公式进行校正，得到最终的值。休恩方法的一般步骤为， （6）注意休恩方法中微分和积分的作用，在点处画出解曲线的切线，用它来得预报点。积分公式（3）梯形公式的误差项为 （7）如果每步中的误差仅由式（7）给出，则在M步后休恩方法的累积误差将是 （8）休恩方法的精度 设的初值问题（1）的解，如果，
26、并且是休恩方法产生的一个近似值序列，则 （9）其中，。所以，休恩方法具有阶精度。3.3泰勒级数法泰勒级数法有着广泛的应用，并且是比较求解初值问题的各种不同数值方法的标准，它可设计为具有任意指定的精度。下面首先将泰勒定理用新的公式表示，使之适合于求解微分方程。泰勒定理 设，且在不动点处有处泰勒公式级数展开： (1)其中， (2)表示函数关于的次全导数。求导公式可以递归地计算： （3）并且一般有 （4）其中为导数算子 区间上初值问题的近似数值解可由各子区间上的公式（1）来推导。次泰勒方法的一般步骤为 （5）其中在各步有。次泰勒方法的最终全局误差是阶的，因此可选择所需大小的，使得误差足够小。如果固定
27、，则理论上可以推导出步长，使之满足任意想要的最终全局误差。然而在实际运算中，通常用和计算两个近似结果集，然后比较其结果。次泰勒方法的精度 设为初值问题的解。如果，而为次泰勒方法产生的近似序列，则 （6）所以，次泰勒方法具有阶精度。3.4龙格库塔方法泰勒方法的优点是最终全局误差的阶为，并且可以通过选择较大的来得到较小的误差。然而泰勒方法的缺点是，需要先确定,并且要计算高阶导数，他们可能十分复杂。每一个龙格库塔（Runge-Kutta）方法都由一个合适的泰勒方法推导而来，使得其最终全局误差为。一种折中方法是每步进行若干次函数求值，从而省去高阶导数计算。这种方法可构造任意阶精度的近似公式。最常用的是
28、的龙格库塔方法，它适用于一般的应用，因为它非常精确、稳定，且易于编程。许多专家声称，没有必要使用更高阶的方法，因为提高的精度与增加的计算量相抵消。如果需要更高的精度，则应该使用更小的步长或某种自适应方法。3.4.1 的龙格库塔方法公式： 其中 这里均为常数。上面的式子称为级显式龙格-库塔方法（Runge-Kutta），简称R-K方法。下面我们讨论4阶龙格-库塔方法（RK4）4阶龙格-库塔方法（RK4）可模拟的泰勒方法的精度。它基于如下对的计算： （1）其中和形如 （2）通过与阶段泰勒级数方法的系数匹配，使得局部误差为，龙格和库塔得出了如下方程组： （3）该方程组有11个方程和13个未知量，必须
29、补充两个条件才能求解。最有用的选择是 （4）其余的变量的解为 （5） 将式（4）和式（5）中的值代入式（2）和（1），得到标准的阶龙格-库塔方法，其描述如下。自初始点开始，利用 （6）生成近似值序列，其中 （7）步长为的辛普森方法的误差项为 （8）如果每一步都唯一误差入式（8）所示，那么步之后方法的累积误差为 （9）下面的定理阐述最终全局误差与步长的关系，它说明使用方法时需要多大的计算量。4阶龙格-库塔方法的精度 设定初值问题的解，如果，且为4阶龙格-库塔方法产生的近似解序列，则 （10）所以，4阶龙格-库塔方法具有阶精度。3.4.2 的龙格-库塔方法2阶龙格-库塔方法模拟的2阶泰勒级数方法的
30、精度。尽管该方法用起来不如,但它的证明更易于理解，更能显其中的原理。首先写出的泰勒级数公式： （12）其中为包含的3阶导数的常数，而级数中的其它项包含幂。式（12）中的和必须表示为及其偏导数的函数。由于 （13）可以用双变量函数的链式求导规则对式（13）关于求导，结果为 利用式（13）可以得到 （14）将式（13）和式（14）代入（12），得到的泰勒展开： （15）考虑阶的龙格库塔方法，它用两个函数值的线性组合来表示： （16）其中 （17）下面用两个独立变量的函数的泰勒多项式逼近来展开，得出的表达式： （18）其中包含的2阶偏导数。然后在式（16）中应用式（18），得到的表示： （19）比较
31、式（15）和式（19）中相似的项，得出如下结论： 表明 表明 表明因此，如果要求和满足关系 （20）则式（19）中的方法将与式（15）的泰勒方法有相同的精度。由于方程组（25）只有个等式，却有个未知量，它的解是不定的，因此允许选定其中的一个系数。这里举出两种特殊的选择为例。情况():选择，由此可得。将这些参数代入式（16）,得到： （21）当用该方法生成时，结果为休恩方法。情况():选择，由此可得。将这些参数代入式（16）,得到： （22）当用此方法生成是，称为改进的欧拉柯西方法。3.4.3 龙格-库塔-费尔伯格方法要保证初值问题的解的精确性，一种方法是分别由步长和进行两次求解，并比较较大步长
32、所对应的网格处的结果。但这样对较小的步长将需要大量计算，而且当结果不够好时必须重新计算。龙格-库塔-费尔伯格方法，试图解决这一问题。它用一个过程来确定是否使用了正确的步长。在每一步中，使用两个不同的求近似解的方法，并比较其结果。如果两个结果相近，则接受该近似；如果两个答案的差超出了指定的精度，则减小步长；如果答案超过了要求的有效位数，则增加步长。每一步要求使用下面个值： （23）然后用阶龙格库塔方法求出初值问题的一个近似解： （24）其中用了个函数值和，注意公式（24）中没有使用。用阶龙格库塔方法得到更好的解： （25）最佳步长可以通过当前步长值乘以标量得到，标量为 （26）其中为指定的误差控
33、制容差。第四章 常微分方程数值解法的MATLAB实现为了得到第三章各种常微分方程数值解法的MATLAB程序，这章将研究各种解法的MATLAB实现，把程序保存在命名为chengxu的文件里。4.1 欧拉方法（欧拉方法）通过计算 ，其中，求上的初值问题的近似解。、欧拉方法的算法（1） 输入函数、自变量区间点，初值和步数。（2） 计算步长。（3） 用欧拉公式求数值解。、欧拉方法的MATLAB程序function E=euler(f,a,b,ya,M)h=(b-a)/M;T=zeros(1,M+1);Y=zeros(1,1+M);T=a:h:b；Y(1)=ya;for j=1:M Y(j+1)=Y(j
34、)+h*feval(f,T(j),Y(J);endE=T,Y;保存在名为euler的M文件中。例、用欧拉方法求解微分方程。步长为，程序执行步。解：由得，在MATLAB输入：f=inline(t2-y,t,y); euler(f,0,2,1,40)运行结果：4.2 休恩方法（休恩方法）通过计算 ，求上的初值问题的近似解，其中。、休恩方法的算法（1） 输入函数、自变量区间点，初值和步数。（2） 计算步长。（3） 用休恩公式求数值解。、休恩方法的MATLAB程序function H=heun(f,a,b,ya,M)h=(b-a)/M;T=zeros(1,M+1);Y=zeros(1,M+1);T=a:h:b;Y(1)=ya;for j=1:M k1=feval(f,T(j),Y(j); k2=feval(f,T(j+1),Y(j)+h*k1); Y(j+1)=Y(j)+(h/2)*(k1+k2);endH=T,Y;保存在名为heun的M文件中。 例、用休恩方法求解微分方程。步长为，执行步。解：由得在MATLAB输入：f=inline(exp(-2*t)-2*y,t,y); heun(f,0,2,0.1,40) 运行结果， 4.3 4次泰勒方法4次泰勒方法 通过计算和，并在每一步中使用泰勒多项式，求上初值问题的近似解。、4次泰勒方法的算法（1） 输入函数df、自变量区间端点