韦联旺毕业汇报总结.doc

,编号： 毕业论文 常微分方程初值问题数值解法的计算机实现数学与计算科学学院韦联旺0400710424数学与计算科学学院袁朝晖副教授课 题： 院 （系）：专 业：学生姓名：学 号： 指导教师单位： 姓 名：职 称：理论研究 实验研究 工程设计 工程技术研究 软件开发题目类型： 2008年 06月 20 日,摘 要本文研究了常微分方程初值问题的数值解法及其计算机实现问题。给出了其理论描述和误差分析以及数值例子。我们讨论了常微分方程初值问题的一些常见数值方法，包括欧拉方法，休恩方法，泰勒方法，龙格库塔方法和预估校正方法，分析了各种方法的步长和截断误差的关系，给出了它们的精确度和稳定性分析，并对这些数值方法的优劣点进行了归纳、分析和比较。通过编写MATLAB和MATHEMATICA程序给出了相关数值例子的实验结果，进而分析实验结果进一步掌握步长和截断误差的关系，以及步长一定时，截断误差和步数的关系。特别是在分析预估校正方法的过程中，我们还得到的前4个已知值（包括）和预估校正法的精确度的关系。 最后，我们对这些数值方法做了进一步的分析比较，指出了不同情况下如何选择合适的数值算法来求解常微分方程初值问题。关键词：常微分方程；初值问题；数值解法；步长；截断误差；精确度.Abstract In this paper,we consider the numerical methor for initial value problem of ordinary differential equetions and how to complement to the computer,give the description of its theory and the analysis of error and numerical examples. We discusse the numerical methor for initial value problem of ordinary differen-tial equetions,including Euler method,Heun method,Taylor method,Runge-Hutta method, and Predictor-corrector method.We analyze the relationship between the step and the truncation error of various methods.we also show the precision of all methods and difference between them.Then we give some examples of the numerical methor for initial value problem of ordinary differential equetions.By editing corresponding MATLAB or MATHEMATICA programs,we give a way to more easily understand the process of solving ordinary differential equations.We also further grasp the relationship between the step and the truncation error and the relationship between truncation error and the number of step when the steps leght is a constant through analyzing the output of the programmings.Especially in the analyzing Predictor-corrector method,we obtain that the first four known values (including )of are related to the precision.Finally,we give the advantages and disadvantages of each method.Key words:ordinary differential equations;numerical methor for initial value problem of ordinary differential equetions;MATLAB；MATHEMATICA;steps leght;the truncation error;the precision目 录引言 11绪论 .21.1微分方程导论 21.2初值问题 . .31.3几何解析 .31.4 步长与误差 . . .42 常见的数值解法及其优劣点 . 42.1 欧拉方法 . 42.1.1欧拉方法的步长与误差 .52.1.2欧拉方法在MATLAB中的执行步骤52.2 休恩方法62.2.1 休恩方法的步长与误差 . 62.1.2休恩方法在MATLAB中的执行步骤72.3泰勒级数法.72.3.1泰勒定理 72.3.2四阶泰勒方法在MATLAB中的执行步骤 .92.4龙格库塔法 .92.4.1龙格库塔法的介绍 .92.4.2龙格库塔法在MATLAB中的执行步骤 .112.5预估校正法112.5.1 Milne-Simpon方法 .122.5.2误差估计与校正 .122.5.3正确的步长 .132.6 数值方法的收敛性分析 132.7 数值方法的稳定性分析 143数值解法的实用举例 153.1 欧拉法三种方法的比较 153.2各种方法的MATHEMATIC数值求解 .173.3 Milne-Simpon方法的一些思考 194结论 19谢辞 21参考文献 .22附录1 .23附录2 .26引言常微分方程诞生于运用数学分析方法解决物理与力学问题的过程中，人们通常认为常微分方程的开端工作是由意大利科学家伽利略(Galileo,15641642) 完成的。17世纪欧洲的建筑师们在建筑教堂和房屋时,需要考虑垂直梁与水平梁在外力作用下的变形，以及当外力撤销时梁的恢复程度，也就是梁的弹性问题。当时的建筑师用经验来处理这些问题。伽利略从数学角度对梁的性态进行了研究，将成果记录在关于两门新科学的对话一书中, 这些研究成果成为常微分方程的开端。从17世纪末开始，对天体问题、摆的运动及弹性理论等问题的数学刻画引出一系列常微分方程。微分方程是在解决实际问题的过程中产生的，微分方程的研究又促进实际问题的解决，同时也促进其他学科的发展。微分方程在物理、工程、力学、天文学、生物学、医学、经济学等诸多领域都有重要作用。如电子计算机与无线电装置的计算问题可归为微分方程求解；弹道计算与飞机飞行中的稳定性研究可归为微分方程的求解；化学反应中稳定性的研究也可归为微分方程求解等等。在天文学上，一般星体都是通过观察得到的，而海王星的发现却是个罕见的例外。牛顿研究天体运动的微分方程,从理论上得到行星运动的规律，而这些规律原来只是由开普勒通过观测归纳出的。而在1846 年,法国巴黎天文台的勒威耶(Leverrier,18111877)在对这个微分方程进行数值分析计算的基础上，预言太阳系中还有第八颗行星的存在，并计算出了第八颗行星的位置，这之后人们按照他的计算结果通过观察才找到海王星。这一事实既推动了天文学的发展,也促进了微分方程的发展。目前，常微分方程的实际背景广、应用性强的特点已受到广泛关注。许多国外教材和国内新版教材已在书中明确强调这一点,并在教材中编入实际应用的例子，希望通过大量的实际问题突出数学的应用，引导学生以常微分方程的形式建立数学模型解决各种实际问题。然而要给出一般方程解的解析表达式是十分困难的，而且往往从解析解得到的数值解也不容易。比如，求解一阶常微分方程初值问题只要少数十分简单的微分方程才能求得其精确解，即使求出解，也往往由于复杂或在解的表达式中有等初等函数值的计算，得到的仍不是精确值，多数情形只能用数值方法求其近似解。欧拉法和龙格库塔法等是求解常微分方程初值问题比较常用的方法，但在实际的应用中，这些求解方法有很多困难，因此借助计算机解决这个问题就显得比较方便。目前最常用的数学软件有MATLAB和Mathematica，借助这些软件来求解常微分初值问题,给出相应的计算机程序，方便在实际中应用，更好的服务于经济发展，也利于提高自己的计算机实际应用能力。如文献1和文献2中，给出了个别方法的计算机实现过程，但他们都没有给出相应的程序，切也没有给出数据分析，不利于掌握数值解法的精神。1 常微分方程初值问题及其常用数值解法的相关理论1.1 微分方程导论方程 (1.1)是一个微分方程，因为它包含“未知函数” 的导数，由于只有独立变量出现在式（1.1）的右端项中，因此的不定积分是方程的一个解。可由积分公式求解： （1.2）其中为积分常数。式（1.2）中的所有函数都是方程（1.1）的解，因为他们都能满足构成的曲线族，如下图。 积分方法可用于求解式（1.2）中函数的显式公式。在这样的解中有1个自由度，即积分变量。通常改变的值可以向下或向上“移动解曲线”， 可以找到过任意需要的曲线 。然而世界的奥妙极少表现为显式的公式，通常只能考查一个变量的变化如何影响另一个变量。将这种方法法医为数学模型，得到的就是包含未知变化率以及字变量和或应变量的方程考虑一个冷却物体的温度。可以想象，温度的变化率与该物体和周围介质的温差相关。经验规律能够验证这一猜想。牛顿冷却定律说明，温度变化直接与温差成正比。如果是周围介质的温度，而为物体在时刻的温度，则 （1.3）其中为正常数。负号是因为当物体的温度高于周围介质温度时，为负值。 如果时刻时的物体温度已知，则程之为初始条件，并将该信息包含在问题描述中。通常需要求解 （1.4） 可用变量分离技术来求解，得 （1.5） 选择不同的将得到不同的解曲线，而且没有把一条曲线变换为另一条曲线的简单方法。初始值决定了最终的解。由式（1.5）可以看出，当增大时。物体温度趋近于房间温度。当，则物体变热而不是变冷。1.2初值问题定义1.1初值问题4 (1.6)在区间上的解是一个可微函数,使得 , , (1.7)注意解曲线必须过初始点.1.3几何解析 在矩形区域中的每个点处，解曲线的斜率都可以有隐式公式得到。因此，整个矩形区域中的都可以计算出来,而且每个值边式与过点的解曲线相切的直线的斜率。 斜率场是指示该区域中斜率的图，用来显示解曲线如何被斜率约束。要沿一条解曲线运动，必须从初始点开始，检查斜率场，以确定沿哪个方向前进。然后在水平方向由到前进一小步,同时在垂直方向上移动大约的距离，使得结果有正确的斜率。解曲线上的下一个点为。重复该过程，沿该曲线运动，由于该过程只进行有限步，该方法将产生一个近视解。定义1.21给定矩形，设在上连续。对任意,，如果存在常数具有性质 (1.8) 则称在上的变量满足利普希茨条件。常数称为的利普希茨常数。定理1.14设定义在区域上,如果存在一个常数。使得 , (1.9) 则在区域上的变量满足利普希茨条件, 利普希茨条件常数为。定理1.2 (存在性和唯一性)12设在区域上连续，如果在区域的变量满足利普希茨条件,并且,则初值问题(6) ，在某个子区间上有唯一解。1.4 步长与误差 这里引入的求初值问题近似解的方法称为差分方法或离散变量法。它求离散点集上的近似解，这个离散点集称为网格（或格网）。基本的单步长方法形如，其中的称为增量函数。 任何求解初值问题的离散变量法都有误差源：离散误差和舍入误差。定义1.3(离散误差) 13设是初值问题的唯一解，是离散近似解集。 全局离散误差定义为, 其中 (1.10)它是唯一解与离散方法得到的解之间的差。 局部离散误差定义为 ,其中 (1.11)它是从到这一步计算的误差。2 常见的数值解法及其优劣点 所谓的微分方程数值解法，其基本思想是把求解区间和方程离散化，求出方程的解在一系列离散点上的近似值。因此，不同的离散方式就产生不同的数值解法。常用的数值解法有欧拉方法、泰勒方法、龙格-库塔方法等。在下问我们将介绍一些常见数值解法，以启发读者多解常微分方程的认识。2.1 欧拉方法不是所以的初值问题都有显式解，而且通常不可能找到解的公式。例如不存在，这种”闭形式表示”的解。因此，对于科学和工程目的，需要有计算近似解的方法。如果需要解具有多位有效数字，则需要更多的技术量和复杂的算法。第一种方法称为欧拉方法，它用来体现先进方法包含的概念，特别适合于快速编程。由于他随着步数增加而产生的积累误差较大，因此用途有限。然而研究它很重要，因为它的误差分析更易懂。欧拉法包含以下两种模式：前进欧拉法和后退欧拉法。前进欧拉法：为求解，良态初值问题的区间。将区间划分为的距子区间，并选择网格点 其中 （2.1）值称为步长。然后近似求解 在上， （2.2）设，连续，利用泰勒定理将在处展开，对每个值，存在一个和之间的值，使得 （2.3）将和代入等式（3），得到的表示： （2.4）如果步长足够小，则可以忽略2次项（包含的项），得到 （2.5）这就是欧拉近似。重复该过程，就能得到近似解曲线的一个点序列。前进欧拉方法的一般步骤是， ，其中 （2.6）后退欧拉法：后退欧拉法是基于后退微分近似并表示成 （2.7）后退欧拉法的精度和前进欧拉法的基本相同。2.1.1步长与误差在分别推导前进欧拉方法的公式(2.6)和后退欧拉方法的公式（2.7）时,每步忽略的项为，如果这是每步的误差，那么在区间的末端进行步之后,积累误差将是 （2.8）可能有其他的误差,但这一估计站了主要部分。所以欧拉方法的精确度可以这样计算： 设是初值问题（2.2）的解，如果且是有欧拉方法计算的近似值序列，则 （2.9） （2.10）区间末端的误差称为最终全局误差： （2.11）由此可看出，欧拉方法的步长减小，则可期望最终全局误差减小,为整数。2.1.2 欧拉法在 MATLAB 中的执行步骤： 定义 输入初始值 输入步长和步数 计算和: for i=1:n x=x+h y=y+hf(x,y) 或 y=y+hf(x+h,y+hf(x,y) end输出 x和 y结束2.2 休恩方法休恩（Heun）其实也是欧拉方法的改进型，因此也叫做欧拉改进法。该方法引入一种新的思路，来构造求解上的初值问题， （2.12）要得到解，可以用微积分基本定理，在上对积分得 （2.13）其中的不定积分为待定求函数。对求解方程（2.13） （2.14）然后可用数值积分方法逼近（2.14）中的定积分，如果采用步长的梯形公式，则结果为 （2.15）注意公式（2.15）的右端包含了待定值。要继续求解，需要的一个估计值，欧拉方法的解能够满足这一目的，将它代入（2.15）后，得到求解的公式，称为休恩（Heun）方法： （2.16）重复这个过程，得到逼近解曲线的一系列点，在每一步中都用欧拉方法作为预报，然后用梯形公式进行校正，得到最终的值。休恩方法的一般步骤为： (2.17 )2.2.1步长和误差 积分公式（2.14）中梯形公式的误差项为 （2.18）如果每步中的误差仅由式（2.18）给出，则步后休恩方法的累积误差将是 （2.19）所以休恩方法的精确度可以这样计算： 设是初值问题（2.12）的解，如果且是有欧拉方法计算的近似值序列，则 （2.20） （2.21）其中，特别地，区间终点处的最终全局误差满足， （2.22）这就是说，步长如果减小为原来的 （为整数），则可期望最终全局误差将降至大约 。2.2.2 休恩法在 MATLAB 中的执行步骤：定义输入初始值输入步长和步数计算和: for i=1:nx=x+hy=y+(h/2)(f(x,y)+f(x+h, y+hf(x,y) end输出 x和 y结束2.3 泰勒级数法泰勒级数法有着广泛的应用，并且是比较求解初值问题的各种不同数值方法的标准，它可设计为具有任意指定的精度。2.3.1泰勒定理 设，且在不动点处有次泰勒级数展开： （2.23）其中， （2.24）表示函数关于的次全导数。求导公式可以递归得计算： （2.25）并且一般有， （2.26）其中为导数算子 区间上的初值的近似数值解可由各子区间上的公式（2.23）来推导。次泰勒方法的一般步骤为： （2.27）其中在个步有，。 次泰勒方法的最终全局误差是阶的，因此可选择所需大小的，使得误差足够小。如果固定，则理论上可以推导出步长，使之满足任意想要的最终全局误差。然而在实际运算中，通常用和计算两个近似结果集，然后比较其结果。次泰勒方法的精度：设为初值问题的解,如果，并且是休恩方法产生的一个近似值序列，则 （2.28） 特别得，区间重点处的最终全局误差满足， （2.29）这就是说，步长如果减小为原来的，则可期望最终全局误差将降至大约，为整数。2.3.2 四阶泰勒法在 MATLAB 中的执行步骤：定义输入初始值输入步长和步数计算和: for i=1:nD=df(x,y)y=y+h(D(1)+h(D(2)/2+h*(D(3)/6+h*D(4)/24)end输出 x和 y结束2.4龙格库塔方法泰勒方法的优点是最终全局误差的阶为，并且可以通过选择较大的来得到较小的误差。然而泰勒方法的缺点是，需要先确定，并且要计算高阶导数，他们可能十分复杂。每个龙格库塔达到都由一个合适的泰勒方法推导而来，使得其全局误差为。一种折中方法是每步进行若干次函数求值，从而省去高阶导数计算。这种方法可构造任意阶精度的近似公式。最常用的是的龙格库塔方法，他适用于一般的应用，因为它非常的精确、稳定，且易于变成。但 5阶的 Runge - Kutta法每步要计算 6个函数值, 6阶方法要计算 7个函数值。而阶 () 方法， 则至少要计算个函数值。许多专家声称，没有必要使用更高阶的方法，因为提高的精度与增加的计算量相抵消。如果需要更高的精度，则应该使用更小的步长或某种自适应方法。2.4.1 四阶龙格库塔方法的介绍4阶龙格库塔方法可模拟的泰勒方法的精度。他基于如下对的计算： （2.30）其中，和形如， （2.31）通过与阶的泰勒级数方法的系数匹配，使得局部误差为，龙格和库塔得出了如下方程组： 该方程组有11个方程和13个未再的量，必须补充两个条件才能求解。最有用的选择是 （2.32）其余变量的解为： （2.33）将式（2.33）和式（2.34）中的值代入式（2.31）和式（2.30）。得到标准的阶龙格库塔方法，其描述如下。自初始点开始，利用 （2.34）生成近似值序列，其中 ， （2.35）龙格库塔方法的精度：设为初值问题的解，如果,并且是4阶龙格库塔方法产生的一个近似值序列，则 （2.36） （2.37）特别的，在区间的末端，最终全局误差满足 （2.38）这就是说，步长如果减小为原来的，则可期望最终全局误差将降至大约，为整数。2.4.2 龙格库塔法在 MATLAB 中的执行步骤：定义输入初始值输入步长和步数计算和: for i=1:nk1=hf(x,y)k2=hf(x+h/2,y+hk1/2)k3=hf(x+h/2,y+hk2/2)k4= hf(x+h,y+hk3)y=y+h(k1+2k2+2k3+k4)/6;end输出 x和 y结束2.5 预估校正方法欧拉方法、休恩方法、泰勒方法以及龙格库塔方法都称为单步方法，因为他们只利用前一个点的信息计算下一个点，即计算时只使用了初始点。一般地，只有用来计算。当计算下一个点，就可以利用几个已计算出的点来计算下一个点。 多步法的一个优点是，可以确定它的局部截断误差，并可以包含一个校正项，用于在每一步计算中改变解的精确度。该方法还可以确定步长是否小到能得到的精确值，同时又大到能够免除不必要的和费时的计算。使用预估子和校正子的组合在每一步只与要进行两次函数求值。 对于多步法，常用的预估校正公式有Milne-Hamning方法、Milne-Simpon方法、4阶隐预估校正公式等。下面，我们向大家介绍Milne-Simpon方法。2.5.1 Milne-Simpon方法 Milne-Simpon方法的预估子基于在区间上对的积分： （2.40）预报子使用的基于点，和的拉格朗日多项式逼近，在区间上对他积分，得到米尔恩预估子： （2.41）校正子的推导类似。此时值已知，基于点，和新点构造的新点的拉格朗日多项式，然后在区间上对该多项式积分，结果得Simpon公式： （2.42）2.5.2误差估计于校正 计算预估子和校正子的数值积分公式的误差项都是的，公式（2.41）和（2.42）的局部截断误差为 （预估子的截断误差） （2.43） （校正子的截断误差） （2.44）设足够小，使得在区间上近乎为常数，则可消去式（2.43）和式（2.44）中的5阶导数项，结果为 （2.45）公式（2.45）给出的预估子误差聚集基于两个计算值和，而没有使用高阶导数，可用它来改进预估值。假设每步中预估和校正值的差缓慢变化，则在式（2.45）中可用和分解替代和，得到如下的修正为： （2.46）在校正过程中用该修正值替代，公式（2.44）变为： （2.47）因此，改进（修正）的Milne-Simpon方法为 （2.48） （2.49） （2.50）由（2.48）和（2.49）知预估校正方法是休恩方法的一个改进。2.5.3正确的步长用预估校正方法在大区间上求解初值问题时，有时也出现问题。如果，而步长过大，则预计估校正方法可能不稳定。根据经验，当小误差递减传播时，结果稳定；当小误差递增传播时候结果不稳定。当在区间上使用的步长太大时，会出现不稳定，有时表现为计算解的振荡性。采用较小的步长可使振幅减小。如果使用步长控制，则Milne-Simpon方法应该使用如下的误差估计： （2.51）这是一类不动点迭代过程。可以证明，其的步长应该满足以下条件 （2.52）2.6 数值方法的收敛性分析 据泰勒公式有 （2.53） 称为数值方法： （2.54）的局部离散误差。假如是在无舍入误差的情形用式（2.53）计算得初值问题 的近似解，则称为数值（2.54）的整体离散误差。 要想近似格式收敛，则要使。由泰勒展式知,若单步法 （2.55）是 p阶的，则其局部离散误差为。从而，其中为一常数。 由 （2.56）以及 （2.57）式 (2.56) 减去式 (2.57) 得：据定理 1.2假设条件, 有。从而得到关系。如假设，则当时。从而得知单步法 ( 2.55) 收敛。对于其余方法，也可以类似得到其收敛性 。2.7 数值方法的稳定性分析稳定性的定义很多。在许多实用的数值方法中，通常总是取有限的固定步长进行计算。而在计算过程中某步误差（例如舍入误差）的传播能否受到控制是与步长的大小有关的，这种稳定性称为绝对稳定性。下面主要讨论数值方法的绝对稳定性问题。必须注意，收敛性与稳定性是两个不同的概念。收敛性是反映数值方法本身的截断误差对计算结果的影响，而稳定性是反映某一计算步中出现的误差（如舍入误差）对计算结果的影响。绝对稳定性是和步长密切相关的，对某一步长是稳定的数值计算公式，对另一种大一点的步长可能就不稳定。显然，只有既收敛而又稳定的数值计算公式才能在实际计算中使用。在考虑一个数值计算方法的绝对稳定性时，强烈依赖于微分方程的右端函数，这就给讨论带来的困难。为了讨论方便，通常只限与特定的方程 （2.58）来讨论绝对稳定性问题。着这个特定问题中，为了保证数值计算方法绝对稳，对参数和步长都要加以一定的限制。通常称和的允许范围为相应方法的绝对稳定区。下面以4阶龙格库塔公式为例来讨论绝对稳定区。对于方程（2.58），由4阶龙格库塔公式的古典形式可以推出 如果在上有一个扰动，且假设其后的全部计算都是精确的，则在处所引起的误差为 因此，为了使处的扰动不超过，应满足如下条件 在平面上，由不等式（2）所确定的区域就是4阶龙格库塔公式的绝对稳定区，其中 为绝对稳定区域的边界。如果绝对稳定区的边界为无限，则称相应的数值计算方法为无条件稳定的3 数值解法的实用举例 常见的通用数学软件包中，Matlab以数值计算见长，而MATHEMATICA以符号运算、公式推导见长。虽然这两种软件都有内置的函数来解常微分方程，但内置函数对我们掌握常微分方程近似数值求解过程的以及理论的分析没有什么帮助，因此我将在下文中列举一些常用方法的MATLAB或者MATHEMATICA程序。无论是MATLAB程序，还是MATHEMATICA程序，他们都有很强的灵活性，即只要改动步长或者迭代方式，就可以的到不同的结果。此外这些程序的运行结果还可以用来做理论分析。3.1 欧拉法三种方法的比较欧拉法在常微分方程近似数值求解,特别是理论分析中占有很重要的地位，它是最早的解决一阶常微分方程初值问题的一种数值解法 , 虽然它不够精确, 很少被采用 , 但它在某种程度上反映了数值方法的基本思想和特征。可以推广到一阶微分方程组及高阶微分方程 , 这是一种通用性强适用面广的简单方法 , 为深入研究问题奠定了基本条件。例1：求在区间上的数值解。容易得出本例的解析解：，下面使用自编的三个欧拉法计算程序计算其近似解 ,并与解析解比较 ,判断各种欧拉法精度的高低。在 MATLAB的 Editor中调用欧拉方法的三个函数，实行附录1中的MATLAB程序1，输出如下图1：