微分方程数值解.docx

1、微分方程数值解微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和 科学技术中的实际问题的研究中，常常需要求解微分方程.但往往只有少数较简单和 典型的微分方程可求出其解析解，在大多数情况下 ，只能用近似法求解，数值解法是 一类重要的近似方法.本文主要讨论一阶常微分方程的初值问题的数值解法， 探讨这些算法在处理来自生活实际问题中的应用，并结合 MATLA软件，动手编程予以解决.1微分方程的初值问题1.1预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题(1)dy f(x,y)dxy(xo) yo这里f x,y是矩形区域R : x xo

2、 a, y y。 b上的连续函数.对初值问题（1）需要考虑以下问题：方程是否一定有解呢？若有解，有多少个解 呢？下面给出相关的概念与定理.定义1 Lipschitz条件12:矩形区域R : XX。 a, y y| b上的连续函数f x, y若满足：存在常数L 0，使得不等式f x, % f x,y2 L % y?对所有x, y1 , x, y2 R都成立，则称f x, y在R上关于y满足Lipschitz条件定理1 解的存在唯一性定理 ：设f在区域D x, ya x b, y R上连续，关于y满足Lipschitz条件，则对任意的xo a,b ,yo R，常微分方程初值问题（1）当x a,b时

3、存在唯一的连续解y x该定理保证若一个函数f x,y关于y满足Lipschitz条件，它所对应的微分方程 的初值问题就有唯一解在解的存在唯一性得到保证的前提下，自然要考虑方程的求解问题.求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微 分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解， 也往往因形式过于复杂或计算量太大而不实用， 因此从实际问题中归结出来的微分方 程主要依靠数值解法.定义2 微分方程数值解：对初值问题（1）寻求数值解就是寻求解y x在一x1 x2 LXnXn 1 L系列离散节点上的近似解y,yi,y2丄，yn,yni,L，相邻两个节点的

4、间距 h Xn 1 Xn称为步长.在一 般情况下假定h h i 0,1,L为常数，这时节点为xn x0 nh, n 0,1,2丄.要求微分方程数值解，首先要建立数值算法，即对初值问题（ 1）中的方程离散化，建立求解数值解法的递推公式.一类是计算 ym时只用到前一点的值yn，称为单步法；另一类是用到yn 1前面k点的值yn, yn 1,L , yn k 1称为k步法.对初值问题（1）式的单步法可用一般形式表示为yn 1 yn h （Xn, yn, yn 1, h），其中多元函数 与f x, y有关，当 含有yn 1时，方法是隐式的；若 中不含yn 1，则为显式方法，所以显式单步法可表示为yn 1

5、 yn h （Xn,yn,h）. （2）设y x是初值问题（1 ）的准确解，称Tn 1 y xn 1 y xn h （xn, y xn ,h）为显式单步法（2）的局部截断误差.若存在最大正整数p,使显式单步法（2）式的局 部截断误差满足Tn1 y x h y x h x, y,h O hp 1 ，则称（2）式有p阶精度. 1.2几种常用的数值解法及其分析、比较1.2.1欧拉法与后退欧拉法1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式y（Xn1）h y（Xn）y（Xn）f（Xn,Yn）h将初值问题（1）离散化，则问题（1）可化为(3)yn 1 yn h f (Xn,yn),焉 X。 n h

6、,此方法称为欧拉法.欧拉方法的几何意义在数值计算公式中体现了出来.在xy平面上，一阶微分方程的解y y x称作它的积分曲线.积分曲线上一点x,y的切线斜率等于函数f x, y，按函数f x,y在xy平面上建立一个方向场，那么，积分曲线上每一点的切 线方向均与方向场在该点的方向相一致基于上述几何解释，从初始点Po（Xo,yo）出发，先依方向场在该点的方向上推进到x Xi上一点R，再从R依方向场的方向推进到x X2上一点F2，循环前进便作出一条折线RFP2L，因此欧拉方法又称为折线法.若初值yo已知，则由（3）式可逐步算出yi yo h f（xo, yo）,y2 yi h f（Xi,yJ,M为了分

7、析计算公式的精确度，通常可用泰勒展开将 y Xni在Xn处展开，则有容易看出，欧拉法（3）式具有一阶精度.如果（4）式右端积分用右矩形公式h f xn 1, y xn 1近似，则得到另一个公式(5)yn 1 yn hf Xn 1 , yn 1 ，称为后退欧拉法.值得一提的是：后退欧拉法与欧拉公式有着本质的区别，后者是关于 yn 1的直接计算公式，它是显式的，而(5)式的右端含有关于yn1的表达式，它是隐式的.在利 用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化 具体迭代过程如下：首先利用欧拉公式丫补铁yn h f(Xn,yn)给出迭代初值丫补铁，把它代入(5)式的 右端，使之转化为

8、显式，直接计算得 y(1) yn h f (xn 1, y(0).如此反复进行，得ynki1) yn h f (Xn i,ynk1) k 0,1,L ，则得到后退欧拉法的迭代公式yn0)1 yn h f (Xn,yn)(k 1) (k)、yn 1 yn h f(Xn 1, yn 1)可以看出，后退欧拉法具有一阶精度，且计算比较麻烦 1.2.2梯形方法为得到比欧拉法精确度高的计算公式，在等式(4)式右端积分中若用梯形求积 公式近似，并用yn代替y Xn，yn 1代替y Xn 1，则得Yn 1 Yn 2 f Xn,Yn f Xn 1, Yn 1 ， (6)称其为梯形方法梯形方法与后退欧拉法一样，都

9、是隐式单步法，可用迭代法求解，其迭代公式为 yn0)1 yn h f(Xn，yn)(k 1) h (k) . (7)yn 1 yn - f Xn, % f X. 1, % 1为了分析梯形公式的收敛性，将(6)与(7)式相减，得yn 1 ynk11) h f Xn 1, yn 1 f ，k 0,1,2,L2因为f X,y满足Lipschitz条件，于是有yn 1 丫二1 吐yn 1 ynk1，其中L为f X,y2关于y的Lipschitz常数.如果选取h充分小，使得匹1，则当k 时有2ynk11) yn1，这说明迭代过程(7)式是收敛的.容易推导得出梯形法(7)式是二阶方法.经分析，梯形方法虽然

10、提高了精度，但是以增加计算量为代价的 从上述的迭代公式可以看出，每迭代一次都要重新计算 f x,y的值，而且迭代又要进行若干次，计算相当的复杂为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法 1.2.3改进的欧拉方法由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方 法是只计算一两次就转入下一步的计算，先用欧拉公式（ 3）求得一个初步的近似解yni，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常 称为改进的欧拉公式.具体公式如下hyn 1 yn - f Xn, yn f X.仆 y. hf （8）2改进的欧拉法与梯形法一样，是二阶方法.1.2

11、.4 Ru nge-Kutta 方法由前面讨论可知，从（4）式Xn 1y Xn 1 y Xn X f X, y X dxXn可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高， 它必然要增加求积积点，为此将（4）式的右端用求积公式表示为Xn 1 rf x, y x dx h f Xn ih, y Xn ih ， （9）Xn i 1一般来说，点数r越多，精度越高，上式右端相当于增量函数 x,y,h，为得到便于计算的显式方法，将公式（9）表示为：yn 1 yn h Xn,yn,h , （ 10）其中rXn,yn,h CjKii 1K1 f Xn , yn （ 11）i 1Ki f

12、 Xn ih,yn h jKj ,i 2,L rj 1这里c, i, j均为常数.c为加权因子，Ki为第i段斜率，共有r段.我们把（10）和(11)称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用的一个公式：K4 f Xn h, yn hK3Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算 ym，只用到前面一步的值yn即可，因此每步的步长可以独立取定 常用的Runge-Kutta方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长 h可取得大些，求解区间上的

13、总步数可以少些.但Runge-Kutta方法也有些缺点，比如四阶Runge-Kutta方法每算一步需要四次计算f x, y的值，计算量较大(对于复杂的f x, y而言).2数值方法的应用实例网例1对于初值问题 y 10y，分别用欧拉法、改进的欧拉法，梯形法求y 1的y 0 1近似值.解：易得该方程的解析解y x e10x，y 1 4.5400e-005,为比较，将按不同数值计算方法所得结果列表如下：表1三种不同方法的数值结果h欧拉法改进的欧拉法梯形法0.2-1110.109.7656E-0045.4994E-0050.012.6561 E-0054.6223 E-0054.5026 E -00

14、50.0014.3717 E-0054.5408 E-0054.5396 E -0050.00014.5173 E-0054.5400 E-0054.5400 E -0050.020.0150.010.0050欧拉法误差改进欧拉法误差 梯形法误差图1三种不同方法数值解与精确解的误差曲线从表1中可以看出：当h 0.2时，三种方法均不稳定，计算结果严重偏离精确值； h 0.1时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当 h 0.01时，三种 方法均稳定，但精确度有区别可以看出，h越小，计算结果越好，要想计算结果充 分接近于解析解还须取较小的h值.图1反映的步长h 0.01时，三种数值方法的

15、所得数值解与解析解在0,1区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进 的欧拉法和梯形法精确度都明显高于欧拉法.例2用欧拉法、改进的欧拉法和 Runge-Kutta法求解初值问题dydx0,12xy ,X y并比较三种方法的结果.解：方程为n 1的伯努利方程，可求得解析解为现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果:图2 ( a)步长为0.2时R-K法和解析解比较图2 ( b)步长为0.2时改进的Euler法和解析解比较图2 ( c)步长为0.2时欧拉法和解析解比较上图2(a)，(b)，(c)描述的是步长为0.2时，用欧拉法、改进的欧拉法， Runge-Kutta法求解方程所得的数值解与解析解之间的对比图由图可知， Runge-Kutta法所得数值解曲线和解析解曲线吻合的很好，改进的欧拉法和欧拉法随着计