微分方程数值解

..微分方程数值解及其应用绪论自然界中的许多事物的运动和变化规律都可以用微分方程来描述，因此对工程和科学技术中的实际问题的研究中,常常需要求解微分方程．但往往只有少数较简单和典型的微分方程可求出其解析解，在大多数情况下,只能用近似法求解，数值解法是一类重要的近似方法．本文主要讨论一阶常微分方程的初值问题的数值解法，探讨这些算法在处理来自生活实际问题中的应用，并结合MATLAB软件，动手编程予以解决．１微分方程的初值问题[1]1.1预备知识在对生活实际问题的研究中，通常需要考虑一阶微分方程的初值问题00(,)()dyfxydxyxy（1）这里,fxy是矩形区域R：00,xxayyb上的连续函数．对初值问题（1）需要考虑以下问题：方程是否一定有解呢?若有解，有多少个解呢?下面给出相关的概念与定理．定义1Lipschitz条件[1][2]：矩形区域R：00,xxayyb上的连续函数,fxy若满足：存在常数0L，使得不等式1212,,fxyfxyLyy对所有12,,,xyxyR都成立，则称,fxy在R上关于y满足Lipschitz条件.定理1解的存在唯一性定理[1][3]：设f在区域,,DxyaxbyR上连续，关于y满足Lipschitz条件，则对任意的00,,xabyR，常微分方程初值问题（1）当,xab时存在唯一的连续解yx．该定理保证若一个函数,fxy关于y满足Lipschitz条件，它所对应的微分方程的初值问题就有唯一解.在解的存在唯一性得到保证的前提下，自然要考虑方程的求..解问题．求解微分方程虽然有多种解析方法，但根据工程和科学实践问题所得到的微分方程往往很复杂，在很多情况下不能或很难给出解析解，有时即使能求出形式解，也往往因形式过于复杂或计算量太大而不实用，因此从实际问题中归结出来的微分方程主要依靠数值解法．定义2微分方程数值解：对初值问题（1）寻求数值解就是寻求解yx在一系列离散节点0121nnxxxxx上的近似解0121,,,,,,nnyyyyy，相邻两个节点的间距1nnnhxx称为步长.在一般情况下假定0,1,ihhi为常数，这时节点为0,0,1,2,nxxnhn．要求微分方程数值解，首先要建立数值算法，即对初值问题（1）中的方程离散化，建立求解数值解法的递推公式．一类是计算1ny时只用到前一点的值ny，称为单步法；另一类是用到1ny前面k点的值11,,,nnnkyyy称为k步法.对初值问题（1）式的单步法可用一般形式表示为11(,,,)nnnnnyyhxyyh，其中多元函数与,fxy有关，当含有1ny时，方法是隐式的；若中不含1ny，则为显式方法，所以显式单步法可表示为1(,,)nnnnyyhxyh.（2）设yx是初值问题（1）的准确解，称11(,,)nnnnnTyxyxhxyxh为显式单步法（2）的局部截断误差.若存在最大正整数p,使显式单步法（2）式的局部截断误差满足11,,pnTyxhyxhxyhOh，则称（2）式有p阶精度.1.2几种常用的数值解法及其分析、比较1.2.1欧拉法与后退欧拉法1）欧拉法：欧拉曾简单地用差分代替微分，即利用公式1()()()(,)nnnnnyxyxyxfxyh将初值问题（1）离散化，则问题（1）可化为..1(,),nnnnyyhfxy0nxxnh，（3）此方法称为欧拉法.欧拉方法的几何意义在数值计算公式中体现了出来.在xy平面上，一阶微分方程的解yyx称作它的积分曲线.积分曲线上一点,xy的切线斜率等于函数,fxy，按函数,fxy在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，从初始点000(,)Pxy出发，先依方向场在该点的方向上推进到1xx上一点1P，再从1P依方向场的方向推进到2xx上一点2P，循环前进便作出一条折线012PPP，因此欧拉方法又称为折线法.若初值0y已知，则由（3）式可逐步算出10002111(,),(,),yyhfxyyyhfxy为了分析计算公式的精确度，通常可用泰勒展开将1nyx在nx处展开，则有2''11,,.2nnnnnnnnhyxyxhyxyxhyxx在nnyyx的前提下，,,.nnnnnfxyfxyxyx可得欧拉法（3）的误差为2211.22nnnnhhyxyyyx容易看出，欧拉法（3）式具有一阶精度.2）向后欧拉方法：如果对微分方程（1）从nx到1nx积分，得11,nnxnnxyxyxftytdx，（4）如果（4）式右端积分用右矩形公式11,nnhfxyx近似，则得到另一个公式111,nnnnyyhfxy，（5）称为后退欧拉法...值得一提的是:后退欧拉法与欧拉公式有着本质的区别，后者是关于1ny的直接计算公式，它是显式的，而（5）式的右端含有关于1ny的表达式，它是隐式的.在利用后退欧拉法时，我们通常利用迭代法求解，实质就是逐步显示化.具体迭代过程如下：首先利用欧拉公式(0)1(,)nnnnyyhfxy给出迭代初值(0)1ny，把它代入（5）式的右端，使之转化为显式，直接计算得11(1)(0)1(,)nnnnyyhfxy.如此反复进行，得(1)()111(,)kknnnnyyhfxy0,1,k，则得到后退欧拉法的迭代公式(0)1(1)()111(,)(,)nnnnkknnnnyyhfxyyyhfxy，可以看出，后退欧拉法具有一阶精度，且计算比较麻烦.1.2.2梯形方法为得到比欧拉法精确度高的计算公式，在等式（4）式右端积分中若用梯形求积公式近似，并用ny代替nyx，1ny代替1nyx，则得111,,2nnnnnnhyyfxyfxy，（6）称其为梯形方法.梯形方法与后退欧拉法一样，都是隐式单步法，可用迭代法求解，其迭代公式为(0)1(1)()111(,),,2nnnnkknnnnnnyyhfxyhyyfxyfxy.（7）为了分析梯形公式的收敛性，将（6）与（7）式相减，得(1)()111111,,2kknnnnnnhyyfxyfxy，0,1,2,k因为,fxy满足Lipschitz条件，于是有(1)()11112kknnnnhLyyyy，其中L为,fxy关于y的Lipschitz常数.如果选取h充分小，使得12hL，则当k时有(1)11knnyy，这说明迭代过程（7）式是收敛的[4].容易推导得出梯形法（7）式是二阶方法...经分析，梯形方法虽然提高了精度，但是以增加计算量为代价的.从上述的迭代公式可以看出，每迭代一次都要重新计算,fxy的值，而且迭代又要进行若干次，计算相当的复杂.为此，有没有比较简便的计算方法呢？下面给出改进的欧拉方法.1.2.3改进的欧拉方法由前面的讨论可知，梯形法计算相对复杂，现对上面的梯形法进行简化，具体方法是只计算一两次就转入下一步的计算，先用欧拉公式（3）求得一个初步的近似解1ny，称为预测值，再利用公式（6）把它校正一次，这样建立的预测-校正系统通常称为改进的欧拉公式.具体公式如下11,,,2nnnnnnnnhyyfxyfxyhfxy（8）改进的欧拉法与梯形法一样，是二阶方法.1.2.4Runge-Kutta方法由前面讨论可知，从（4）式11,nnxnnxyxyxfxyxdx可以看出，若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必然要增加求积积点，为此将（4）式的右端用求积公式表示为11,,nnrxininixifxyxdxhcfxhyxh，（9）一般来说，点数r越多，精度越高，上式右端相当于增量函数,,xyh，为得到便于计算的显式方法，将公式（9）表示为：1,,,nnnnyyhxyh（10）其中1111,,,,,2,rnniiinniininijjjxyhcKKfxyKfxhyhKir（11）这里,,iiijc均为常数.ic为加权因子，iK为第i段斜率，共有r段.我们把（10）和..（11）称为r级显式Runge-Kutta法，简称为R-K方法.下面给出其中最经典最常用的一个公式：11234121324322,6,,,22,,22,.nnnnnnnnnnhyyKKKKKfxyhhKfxyKhhKfxyKKfxhyhK（12）Runge-Kutta方法作为一种重要的单步方法，具有很高的实用价值，它关于初值是稳定的，其解连续地依赖于初值，是一类便于应用的单步法，为了计算1ny，只用到前面一步的值ny即可，因此每步的步长可以独立取定.常用的Runge-Kutta方法精度较高，为了达到预定的精度，与欧拉方法与梯形法相比，步长h可取得大些，求解区间上的总步数可以少些.但Runge-Kutta方法也有些缺点，比如四阶Runge-Kutta方法每算一步需要四次计算,fxy的值，计算量较大（对于复杂的,fxy而言）．2数值方法的应用实例[5-9]例1对于初值问题1001yyy，分别用欧拉法、改进的欧拉法，梯形法求1y的近似值.解：易得该方程的解析解10xyxe，14.5400e-005y，为比较，将按不同数值计算方法所得结果列表如下：..表1三种不同方法的数值结果h欧拉法改进的欧拉法梯形法0.2-1110.109.7656E-0045.4994E-0050.012.6561E-0054.6223E-0054.5026E-0050.0014.3717E-0054.5408E-0054.5396E-0050.00014.5173E-0054.5400E-0054.5400E-00500.20.40.60.8100.0050.010.0150.02欧拉法误差改进欧拉法误差梯形法误差图1三种不同方法数值解与精确解的误差曲线从表1中可以看出：当0.2h时，三种方法均不稳定，计算结果严重偏离精确值；0.1h时，改进后的欧拉和梯形法均稳定，但欧拉法效果很差；当0.01h时，三种方法均稳定，但精确度有区别．可以看出，h越小，计算结果越好，要想计算结果充分接近于解析解还须取较小的h值．图1反映的步长0.01h时，三种数值方法的所得数值解与解析解在[0,1]区间的误差曲线，由图可知，在步长相同的情况下，梯形法的精确度略高于改进的欧拉法；改进的欧拉法和梯形法精确度都明显高于欧拉法．例2用欧拉法、改进的欧拉法和Runge-Kutta法求解初值问题2,0,101dyxyxdxyy并比较三种方法的结果．..解：方程为1n的伯努利方程，可求得解析解为2442xxxyeexe现用MATLAB软件编程，用题目要求的方法求解，可得如下图示结果：00.20.40.60.8111.21.41.61.8解析解Runge-kutta法图2（a）步长为0.2时R-K法和解析解比较00.20.40.60.8111.21.41.61.8解析解改进的Euler法图2（b）步长为0.2时改进的Euler法和解析解比