刚性常微分问题的数值解法及编程

1题目：刚性常微分问题的数值解法课程名称：创新实验学院：理学院专业：数学与应用数学年级：应数131学号：1307010239234236姓名：袁蕊张蕾刘霖指导教师：罗贤兵2015年07月14日创新实验论文2目录第一章绪论·············································3选题背景·········································3刚性问题的算法···································3引言·············································4第二章刚性问题·········································5第三章预备知识·········································8第四章计算实验·········································15附页····················································203第一章绪论自然界和工程技术的很多现象,其数学模型是常微分方程（组）的初值问题，普通的常微分方程的数值解法已经比较成熟，理论比较完整，也有许多方法可供选择。但，有一类常微分方程组，求解值时遇到相当大的困难，这类常微分方程组解的分量有的变化很快，有的变化缓慢，常常出现这种现象：变化快的分量很快趋于它的稳定值，而变化慢的分量缓慢趋于它的稳定值。从数值解的观点看来，当变化快时应该用小步长积分，当变化快的分量已趋于稳定，或者说已经没有变化快的分量时就应该用较大的步长积分，但是理论和实际都说明，很多方法特别是显示方法的步长任不能放大，否者便出现数值不稳定现象，即误差急剧增加，已致掩盖了真解，使求解过程无法继续进行。常微分方程组的这种性质叫做刚性，我们考虑一阶常微分方程初值问题的数值解法。0)(),(yaybxayxfdxdy(1.1)常微分方程的解能用初等函数、特殊函数或它们的级数与积分形式表达的非常之少，用解析办法只能求解线性常系数等特征类型的常微分方程。在实际问题中归结出来的求解微分方程的方法只要依靠数值解法。所谓数值解法，就是通过某种离散化办法，将微分方程转化为差分方程来求解。求方程（2-1）的数值解，即对一系列离散节点1n10nxxxxa建立求)(nxy的近似值ny的递推格式，由此求得解y(x)在各节点的近似值，n=0，1,2，…。相邻两个节点的间距nnxxh1称为步长，这时节点nhxxn0。因此，这样得到的数值解法也称为差分方法。初值问题（1.1）的数值解法，可区分为两大类：（1）单步法：此类方法在计算1nx上的近似值1ny时只用到了前一点nx上的信息。如Euler法、Runge-Kutta法和Taylor级数法这就是这类方法的典型代表。（2）多步法：此类方法在计算1ny时，除了需要nx点的信息外，还需要,,21nnxx前面若干个点上的信息。线性多步法是这类方法的典型代表。本文讨论的是几种隐式方法（向后欧拉，梯形公式，改进欧拉）和隐式Runge-Kutta4第二章刚性问题选取的步长h必须很小，满足h＜1/100,才能保证绝对稳定性要求。对于非线性常微分方程初值问题,,)(,)),(,(000mRyyxybxaxyxfy若初值问题是稳定的，即dy/dx0.用欧拉法进行数值求解时，h应满足11hyf。若yfMmax，h应满足h=2/M。在方程组的情况下，例如一阶常系数线性方程组0)(yyayAdxdy这里的A=ijass，y=Tsyyy21,.记A的特征值为s,,21，对稳定的初值问题应满足Re0i.用欧拉法数值求解时，为了保证计算的稳定性，h的选取应满足isi1max2h由下面的例子可以知道，当比值isiisiReminRemax11很大时，h很小，这时计算不熟很多，耗时很长，给实际计算带来了很大的困难。例，某一物理现象可归结为一个线性方程组TyAydxdy1,0,10（1.2）其中x为时间变量，而404040202119201921AA的特征值分别为i140,i140-1-321，。方程（2-11）的解为5)40exp()40sin40(cosx)40exp()40sin40(cos21)2exp(21)()40exp()40sin40(cos21)2exp(21y321xxxyxxxxxyxxxxx）（（1.3）我们对解式（1.3）编程作图，可以看出这组解在开始时刻变化激烈，随后逐渐进入稳态，对应于3,2，的分量在解中的作用随时间x的推移越来越显得无足轻重。解式（1.3）程序和图像曲线（图1）如下面所示解程序：h1=ezplot('(1/2)*exp(-2*x)+(1/2)*(cos(40*x)+sin(40*x))*exp(-40*x)');axis([01-11.5]);holdon;h2=ezplot('(1/2)*exp(-2*x)-(1/2)*(cos(40*x)+sin(40*x))*exp(-40*x)');axis([01-11.5]);set(h2,'Color','r')holdon;h3=ezplot('-(cos(40*x)-sin(40*x))*exp(-40*x)');axis([01-11.5]);set(h3,'Color','k')holdonlegend('y1','y2','y3')图像：6图一由于在开始的一段时间量x，解曲线变化激烈，对方程进行数值求解时，自然要求数值有较高的精度，而对较大的时间量x,解曲线变化缓慢，因此，对数值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量x的大小而改变。例如对初值问题（1.2）用欧拉折现法，步长必须满足h035.0402max2i1is,这样小的步长对于较大的求解区间是难以接受的。我们看到，步长主要受特征值24032的限制，如前所述，正是这两个特征值，在微分方程中随时间量x的增大而显得作用减小，这种矛盾完全是比值isiisiReminRemax11过大造成的。定义1.若线性系统式（2-11）中A的特征值i满足条件（1）mii,,1,0,0Re（2）R=isiisiReminRemax111称式（1.2）为刚性方程，比值R为刚性比。7第三章预备知识隐式方法1．1欧拉公式（显示）考虑初值问题)5.1()()4.1(),(ooyxyyxfdxdy为了求得它在等距离散点21nxxx上的数值解，首先将（1.1）离散化。设)2,1(1ixxhii，将式（1.4）离散化的办法有Taylor展开法、数值微分法及数值积分法如果在点nx将)()(1hxyxynn做Taylor展开，得)6.1(),(),(!2)()()(121nnnnnnnxxyhxyhxyxy那么当h充分小时，略去误差项)(22nyh，用ny近似值代替1)(nnyxy、近似代替)(1nxy，并注意到))(,()(nnnxyxfxy，便得,1,0),,(),(10nyxhfyyxyynnnno（1.7）其中.,0Nabhnhxxn用（1.4）求解（1.1）的方法称为Euler方法。如果利用差商代似替代微商，那么可得)).(,()()()(1nnnnnxyxfxyhxyxy（1.8）在（1.8）若用ny近似替代1)(nnyxy、近似替代)(1nxy，同样得到递推公式（1.7）.如果在],[1nnxx上对))(,(xyxfy积分，得.))(,()()(1dxxyxfxyxynn（1.9）那么对上式右端积分用左矩形求积公式，并用ny近似替代1)(nnyxy、近似替代)(1nxy，也可得递推公式（1.7）.我们知道，在xOy平面上，微分方程),(yxfdxdy的解称为积分曲线，积分曲线上一点（x,y）的切线斜率等于函数f（x,y）D的值。如果D中每一点，都画上一条以f（x,y）在这点的值为斜率并指向x增加方向的有向线段，即在D8上作出了一个由方程),(yxfdxdy确定的方向场，那么方程的解f=y（x）.从几何上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向一致。从初始点),(000yxP出发，过这点的积分曲线为y=y(x)，斜率为),(00yxf.设在0xx附近y(x)可用过0P点的切线近似表示，切线方程为))(,(0000xxyxfyy.当1xx时，)(1xy的近似值为))(,(01000xxyxfy，并记为1y，这就是得到1xx时计算)(1xy的近似公式))(,(1111xxyxfyy.当2xx时，)(2xy的近似值为))(,(12111xxyxfy，并记为2y.于是就得到当2xx的近似公式))(,(121112xxyxfyy.重复上面方法，一般可得当1nxx的计算)(1nxy的近似公式))(,(11nnnnnnxxyxfyy.如果)2,1(1ixxhii，则上面公式就是（1.7）.将NPPP,,,10连续起来，就得到一条折线，所以Euler方法又称为折线法。由公式（1.4）看出，已知0y便可算出1y.已知1y，便可算出2y，如此继续下去，这种只用前一步的值ky便可计算出1ky的递推公式称为单步法。1.2向后欧拉公式（隐式）若在（1.6）中，其右边的积分由数值积分的右矩形公式近似，并用ny替代)(nxy替代)(1nxy，则可得到),()(1110nnnnoyxhfyyxyy（1.10）并称（1.10）为后退Euler公式。1.3梯形公式（隐式）若在公式（1.6）中，其右边积分用数值梯形积分公式近似，并用ny替代)(nxy，1ny替代)(1nxf，则可得到梯形方法公式9)],(),([2111nnnnnnyxfyxfhyy(1.11)梯形方法同头退Euler方法一样都是隐式单步法。对于隐式方法，通常采用迭代法。对后退Euler方法，先Euler方法计算1ny，并将它作为初值)0(1ny，即),(0)0(1nnfnyxhyy，再把它代入（1.7）的右端，便得到后退Euler方法的迭代公式为.1,0),(),()(11)1(1)0(1kyxhfyyyxhfyyknnnknnnnn(1.12)同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为,1,0)],(),([2),()(11)1(1)0(1kyxfyxfhyyyxhfyyknnnnnknnnnn(1.13)1.4改进欧拉（隐式）),(),((2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy1.5Runge-Kutta法龙格-库塔方法的基本思路是想办法计算f(x,y)在某点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式和解的泰勒展开式相比较，便得前面的若干项相吻合，从而达到较高的精度，一般的显示R-K方法的形式如下：)2(),,(),(11111rikyhxhfkyxhfkkyyijjijnninnriiinn(1.14)当式（1.14）的布局截断误差达到)(1phO时成公式（1.14）为p阶r段R-K方法。类似于显示方法我