《刚性常微分问题的数值解法及编程解读》由会员分享，可在线阅读，更多相关《刚性常微分问题的数值解法及编程解读（17页珍藏版）》请在金锄头文库上搜索。
1、创新实验论文题 目： 刚性常微分问题的数值解法课程名称:创新实验学院：理学院专 业：数学与应用数学年级：应数131学 号： 1307010239 234 236姓 名：袁蕊张蕾刘霖指导教师：罗贤兵2015年07月14日18目录第一章绪论3选题背景3刚性问题的算法3弓I言4第二章刚性问题5第三章预备知识8第四章计算实验15附页20第一章绪论自然界和工程技术的很多现象，其数学模型是常微分方程(组)的初值问题， 普通的常微分方程的数值解法已经比较成熟，理论比较完整，也有许多方法可供 选择。但，有一类常微分方程组，求解值时遇到相当大的困难，这类常微分方程 组解的分量有的变化很快，有的变化缓慢， 常常出
2、现这种现象：变化快的分量很 快趋于它的稳定值，而变化慢的分量缓慢趋于它的稳定值。从数值解的观点看来， 当变化快时应该用小步长积分，当变化快的分量已趋于稳定，或者说已经没有变 化快的分量时就应该用较大的步长积分， 但是理论和实际都说明，很多方法特别 是显示方法的步长任不能放大，否者便出现数值不稳定现象，即误差急剧增加， 已致掩盖了真解，使求解过程无法继续进行。常微分方程组的这种性质叫做刚性，我们考虑一阶常微分方程初值问题的数值解法。；x=f(x，y)外 y(a) = v。(1.1)常微分方程的解能用初等函数、特殊函数或它们的级数与积分形式表达的非常之 少，用解析办法只能求解线性常系数等特征类型的
3、常微分方程。在实际问题中归结出来的求解微分方程的方法只要依靠数值解法。所谓数值解法，就是通过某种离散化办法，将微分方程转化为差分方程来求解。求方程(2-1)的数值解，即对一系列离散节点a =x。X1 Xn Xn 1 建立求y(Xn)的近似值yn的递推格式，由此求得解y(x)在各节点的近似值，n=0, 1,2 ,。相邻两个节点的间距h =xn书-xn称为步长，这时节点xn=x0+nh。因此，这样得到的数值解法也称为差分方法。初值问题(1.1 )的数值解法，可区分为两大类：(1)单步法：此类方法在计算人书上的近似值外书时只用到了前一点小上的信 息。如Euler法、Runge-Kutta法和Tayl
4、or级数法这就是这类方法的典型代表。(2)多步法：此类方法在计算外书时，除了需要 点的信息外，还需要xnxn前面若干个点上的信息。线性多步法是这类方法的典型代表。本文讨论的是几种隐式方法(向后欧拉，梯形公式，改进欧拉)和隐式Runge-Kutta选取的步长h必须很小, 线性常微分方程初值问题若初值问题是稳定的，即第二章刚性问题满足h 1/100,才能保证绝对稳定性要求。对于非V= f(x,y(x),a xb,=、y(x。)= yo,yo 亡 Rm, 5f./1 + h 1 o M 例在方程组的情况下,=maxcf,h应满足h=2/M。例如一阶常系数线性方程组型二Ay dxy(a) = yo这里
5、的A=faj屋,y=(yi, y2ys T.记a的特征值为%,，h ，对稳定的初值问题应满足ReA 0.用欧拉法数值求解时，为了保证计算的稳定性，h的选取、一2应满足h -2- max 九 i 1s由下面的例子可以知道，当比值maxRe i半艮大时，h很小，这时计算不熟min Re i1is很多，耗时很长，给实际计算带来了很大的困难。例，某一物理现象可归结为一个线性方程组dy 二 Aydxy 0 = 1,0,-1T(1.2)-2119-20”其中x为时间变量，而A= 19-212040-40 40A的特征值分别为九 1 =-1, % =-40(1 *i )% =40(1 i )。方程(2-11
6、 )的解为dy/dx 0.用欧拉法进行数值求解时，h应满足一 1, C、，1 ,、一、y1 (x )= -exp( 一2x) +- (cos40x +sin 40x)exp( -40 x)(1.3)，、1，八、1，八，一、，、y2 (x) = exp(-2x) - (cos40x +sin 40x)exp(-40x) 22y3(x) =-(cos40xsin 40x)exp(-40x)进入稳态，对应于我们对解式（1.3）编程作图，可以看出这组解在开始时刻变化激烈，随后逐渐任,，的分量在解中的作用随时间x的推移越来越显得无足轻重。解式(1.3)程序和图像曲线(图1)如下面所示解程序：h1=ezp
7、lot(1/2)*exp(-2*x)+(1/2)*(cos(40*x)+sin(40*x)*exp(-40*x);axis(0 1 -11.5);hold on;h2=ezplot(1/2)*exp(-2*x)-(1/2)*(cos(40*x)+sin(40*x)*exp(-40*x);axis(0 1 -11.5);set(h2,Color,r)hold on;h3=ezplot(-(cos(40*x)-sin(40*x)*exp(-40*x);axis(0 1 -11.5);set(h3,Color,k)hold on legend(yT,y2,y3)图像：图一由于在开始的一段时间量X,解
8、曲线变化激烈，对方程进行数值求解时，自 然要求数值有较高的精度，而对较大的时间量X,解曲线变化缓慢，因此，对数值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量 的大小而改变。例如对初值问题（1.2）用欧拉折现法，步长必须满足 h1min Re i称式(1.2)为刚性方程，比值R为刚性比。第三章预备知识隐式方法1.1 欧拉公式(显示) 考虑初值问题字=f(x,y)(1.4)dxy(xo) = y。(1.5)为了求得它在等距离散点x1x2，xn一上的数值解，首先将(1.1)离散化。设卜=为x-(i=1,2)，将式(1.4)离散化的办法有Taylor展开法、数值微分法及数值积分法
9、如果在点xn将y(4书)=y(xn+h)做Taylor展开，得(1.6)yn书近似代替. h2 .y(xn 1) = y(xn) hy (xn)2! y ( n), n 函凡.1)h2那么当h充分小时，略去误差项一y”(二)，用yn近似值代替y(xn)、2y(xn 书)，并注意到 y(xn) = f (xn, y(xn),使得(1.7)。= y(xo),=Jn+ =yn +hf (xn,yn),0 = 0,1；b - a其中xn =x。+nh,h =.用(1.4 )求解(1.1 )的方法称为Euler方法。N如果利用差商代似替代微商，那么可得y(xn 1) - y(xn):y(xn)= f (
10、xn, y(xn).(1.8)在(1.8)若用yn近似替代y(xn)、yn由近似替代y(xn书)，同样得到递推公式(1.7).如果在xn，xn书上对y= f (x, y(x)积分，得y(xn41)-y(xn) = f f (x, y(x)dx.(1.9)那么对上式右端积分用左矩形求积公式，并用yn近似替代y(xn)、yn由近似替代y(xn+),也可得递推公式(1.7).我们知道，在xOy平面上，微分方程dy = f(x,y)的解称为积分曲线，积分dx曲线上一点(x,y)的切线斜率等于函数f (x,y) D的值。如果D中每一点，者B 画上一条以f (x,y)在这点的值为斜率并指向x增加方向的有向
11、线段，即在 D上作出了一个由方程dy = f(x, y)确定的方向场，那么方程的解f=y (x).从几何 dx上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方 向一致。从初始点B(xo，y。)出发，过这点的积分曲线为y=y(x),斜率为f(xo,y0).设在x = %附近y(x)可用过B点的切线近似表示，切线方程为 y = y0+f(x,yo)(xx0).当 x=x1时,y(x1)的近似值为 y0 + f (%, yo)(x1 -xO),并 记为yi,这就是得到x=x时计算y(xj的近似公式y = y f(xi,yi)(x-xi).当x = x2时，y(x2)的近似值为y
12、I + f (x，y1)(x2 -x1)，并记为y2.于是就得到当x = x2的近似公式N2 f f (x1，yi)(x2 -xi).重复上面方法，一般可得当x = %书的计算丫(书)的近似公式yn 1 =Vn f (xn，yn)(xn i - ) .如果h=X -x-(i =1,2)，则上面公式就是(1.7).将PoF,，Pn连续起来， 就得到一条折线，所以Euler方法又称为折线法。由公式(1.4)看出，已知y0便可算出yi.已知yi ,便可算出 刈,如此继续下去，这种只用前一步的值yk便可计算出yk书的递推公式称为单步法。1.2 向后欧拉公式(隐式)yn替代若在(1.6)中，其右边的积分
13、由数值积分的右矩形公式近似，并用y(xn)替代y(xn书)，则可得到(1.10)“。=y(xo) yn 1- ” hf(xn 1, yn 1)并称(1.10)为后退Euler公式。1.3 梯形公式(隐式)yn替代若在公式(1.6)中，其右边积分用数值梯形积分公式近似，并用y(xn) , yn由替代f (xn书),则可得到梯形方法公式(1.11)h.yn1 y n 2f(xn，yn) f (xn 1，yn 1)通常采用方法的迭梯形方法同头退Euler方法一样都是隐式单步法。对于隐式方法, 迭代法。对后退Euler方法，先Euler方法计算y。并将它作为初值yn) = y0+hf (xn，yn)，再把它代入(1.7)的右端，便得到后退Euler代公式为(0), I、(1.12)yn 44 = yn hf ( xn , yn)(k +)，-,(k)、，yn m=yn+hf(Xn 4,ynH1)k=0,1同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为yn0)1=yn hf (Xn,yn)* h(1.13)yn；) =yn +万f (Xn, yn) + f (Xn 书,y) k = 0,1：1.4 改进欧拉(隐式)yn 1 =yn hf (Xn,yn) J/一 、一，一、yn+ = yn +2(f (Xn，yn)+ f (Xn+fn 书)1.5Runge
