《刚性常微分问题的数值解法及编程》由会员分享，可在线阅读，更多相关《刚性常微分问题的数值解法及编程（16页珍藏版）》请在金锄头文库上搜索。
1、创新实验论文题 目：刚性常微分问题的数值解法课程名称学 院： 理学院专 业：数学与应用数学年 级： 应数131学 号：1307010239 234236姓 名：袁蕊蕾霖指导教师： 罗贤兵2015年07月14日目录第一章绪论选题背景3刚性问题的算法3引言4第二章刚性问题5第三章预备知识8第四章计算实验15附页20第一章绪论自然界和工程技术的很多现象，其数学模型是常微分方程（组）的初值问题，普 通的常微分方程的数值解法已经比较成熟，理论比较完整，也有许多方法可供选择。但，有一类常微分方程组，求解值时遇到相当大的困难，这类常微分方程组 解的分量有的变化很快，有的变化缓慢，常常出现这种现象：变化快的分
2、量很快 趋于它的稳定值，而变化慢的分量缓慢趋于它的稳定值。从数值解的观点看来， 当变化快时应该用小步长积分，当变化快的分量已趋于稳定，或者说已经没有变 化快的分量时就应该用较大的步长积分，但是理论和实际都说明，很多方法特别 是显示方法的步长任不能放大，否者便出现数值不稳定现象，即误差急剧增加， 已致掩盖了真解，使求解过程无法继续进行。常微分方程组的这种性质叫做刚性，我们考虑一阶常微分方程初值问题的数值解法。dy 、（1.1）=/（尤,y） dxy （。） = y”、0即对一系列离散节点常微分方程的解能用初等函数、特殊函数或它们的级数与积分形式表达的非常之 少，用解析办法只能求解线性常系数等特征
3、类型的常微分方程。在实际问题中归 结出来的求解微分方程的方法只要依靠数值解法。所谓数值解法，就是通过某种 离散化办法，将微分方程转化为差分方程来求解。求方程（2-1）的数值解，a = x x x x 01n n+1建立求y（x ）的近似值y的递推格式，由此求得解y（x）在各节点的近似值，n=0，1,2，。相邻两个节点的间距九=x -x称为步长，这时节点x = x + nh。 n+1 nn 0因此，这样得到的数值解法也称为差分方法。初值问题（1.1）的数值解法，可区分为两大类：（1）单步法：此类方法在计算x上的近似值y时只用到了前一点x上的信n+1n+1n息。如Euler法、Runge-Kutt
4、a法和Taylor级数法这就是这类方法的典型代表。（2）多步法：此类方法在计算yn+1时，除了需要xn点的信息外，还需要xn 1，xn 2,前面若干个点上的信息。线性多步法是这类方法的典型代表。本文讨论的是几种隐式方法（向后欧拉，梯形公式，改进欧拉）和隐式 Runge-Kutta第二章刚性问题选取的步长h必须很小，满足hl/100,才能保证绝对稳定性要求。对于非 线性常微分方程初值问题y = f(x,y(x),axb,y(x ) = y , j C Rm,I 000若初值问题是稳定的，即 dy/dx 0.用欧拉法进行数值求解时，h应满足 1 +堂九1。若M=max堂，h应满足h二2/M。dyd
5、y在方程组的情况下，例如一阶常系数线性方程组dy=Ay.记A的特征值为X ,X，人，对稳定的初值ij sxs12 s12 s问题应满足Re人0.用欧拉法数值求解时，为了保证计算的稳定性，h的选取 I2 应满足h maxX/lismaxlReX I由下面的例子可以知道，当比值心孔 ，很大时,h很小，这时计算不熟 min|ReXz| li其中x为时间变量，而人=仁2119、4019-21-40-2020-40A的特征值分别为方程A =-1, X =-4C)G + i)A = 4oG i) o123y G) = ； exp(-2x) + ： (cos40x + sin 40x) exp(-40x)1
6、 1(1. 3)y (x) = exp(-2x) - (cos40x + sin 40x) exp(-40x)2j (x) = -(cos40x一 sin40x)exp(-40x)3我们对解式(1.3)编程作图，可以看出这组解在开始时刻变化激烈，随后逐渐进入稳态，对应于气，气的分量在解中的作用随时间x的推移越来越显得无足轻重。解式(1.3)程序和图像曲线(图1)如下面所示解程序：h1二ezplot(1/2)*exp(-2*x) + (1/2)*(cos(40*x)+sin(40*x)*exp(-40*x);axis(0 1 -1 1.5);hold on;h2二ezplot(1/2)*exp(
7、-2*x)-(1/2)*(cos(40*x)+sin(40*x)*exp(-40*x);axis(0 1 -1 1.5);set(h2,Color,r)hold on;h3二ezplot(-(cos(40*x)-sin(40*x)*exp(-40*x);axis(0 1 -1 1.5);set(h3,Color,k)hold onlegend(y1,y2,y3)图像：图一由于在开始的一段时间量x，解曲线变化激烈，对方程进行数值求解时，自 然要求数值有较高的精度，而对较大的时间量x,解曲线变化缓慢，因此，对数 值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量x 的大小而改变。
8、例如对初值问题（1.2）用欧拉折现法，步长必须满足 看到，步长主要受特征值|气| = |气| = 4应的限制，如前所述，正是这两个特征 值，在微分方程中随时间量x的增大而显得作用减小，这种矛盾完全是比值 maxlRe Xh一-max人 ii 1= . 0.035，这样小的步长对于较大的求解区间是难以接受的。我们1企号 4过大造成的。min Re Xi1i s定义1.若线性系统式（2-11）中A的特征值Xi满足条件（1）ReX. 1minRe Xi|1i s称式（1.2）为刚性方程，比值R为刚性比。第三章预备知识隐式方法1. 1欧拉公式（显示）考虑初值问题(1.4)(1.5)dy、=/(x，y)
9、 dxy (x。) = y。为了求得它在等距离散点 x1 x2 x 上的数值解，首先将（1.1）离散化。设h = xt -xi1（i = 1,2.），将式（1.4）离散化的办法有Taylor展开法、数值(1.6)微分法与数值积分法如果在点x将y（x= y（x + h）做Taylor展开，得y(x = y(x) + hy，*) + 万 y”(&)，勺 (x,x那么当h充分小时，略去误差项h2y（&2），用y近似值代替y（x ）、y近似代替，nn n+1便得yJ，并注意到 y，(x) = f (x, y(x)，(1.7) y 近似替代y(x，同样得到递推公式(1.7).如果在x ,x上对y = f
10、 (x, y(x)积分，得y (x ) 一 y (x ) = j f (x, y (x)dx.(1.9)那么对上式右端积分用左矩形求积公式，并用yn近似替代y(xn)、yn+1近似替代y(x 1)，也可得递推公式(1.7).我们知道，在xOy平面上，微分方程空=f (x,y)的解称为积分曲线，积分 dx曲线上一点(x,y)的切线斜率等于函数f(x,y)D的值。如果D中每一点，都 画上一条以f(x,y)在这点的值为斜率并指向x增加方向的有向线段，即在D 上作出了一个由方程空=f (x, y)确定的方向场，那么方程的解f二y(x).从几何 dx上看，就是位于此方向场中的曲线，它在所经过的每一点都与
11、方向场的该点的方 向一致。从初始点P (x ,y )出发，过这点的积分曲线为y=y(x)，斜率为f (x ,y ).设 00000在x = x0附近y(x)可用过P0点的切线近似表示，切线方程为y = y + f (x ,y )(x-x ).当 x = x 时，y(x )的近似值为 y + f (x ,y )(x -x )，并 00001100010记为V，这就是得到x =气时计算y(气)的近似公式y = y.+f (x, y)(x -气).y2.于是就得到当时，y(x )的近似值为y +f (x ,y )(x -x )，并记为211121x = x2的近似公式y = y +f (x, y )
12、(x x).211121重复上面方法，一般可得当x =七1的计算y(x+1)的近似公式y = y +f (x , y )(x x ). n+1nn n n+1n如果h = x -x (i = 1,2)，则上面公式就是(1.7).将P,P,P连续起来，ii-10 1 N就得到一条折线，所以Euler方法又称为折线法。由公式(1.4)看出，已知y便可算出y .已知y，便可算出y，如此继续 0112下去，这种只用前一步的值七便可计算出七+1的递推公式称为单步法。1.2向后欧拉公式(隐式)若在(1.6)中，其右边的积分由数值积分的右矩形公式近似，并用y替代y(x )替代j(x 1)，则可得到(1.10
13、) y0 = y (xo)i y+1 = y* 询、jJ并称(1.10)为后退Euler公式。1.3梯形公式(隐式)若在公式(1.6)中，其右边积分用数值梯形积分公式近似，并用y替代y(x )，y 1替代f (x 1)，则可得到梯形方法公式h,、y+1 = y + 2 f (xn, y) + f (x +1, y(1.11)梯形方法同头退Euler方法一样都是隐式单步法。对于隐式方法，通常采用 迭代法。对后退Euler方法，先Euler方法计算y ，并将它作为初值y(。)，即 n+1n+1y(0) = y + h (x ,y )，再把它代入(1.7)的右端，便得到后退Euler方法的迭 n+10 f n n代公式为(1.12) yn。 = yn + hf (x, y)i yw) = y + hf (x , y( k) k = 0,1 .同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为y(1i = yn + hf (x, y)n+1 nn+1n+1y(5 = y + h f (x，
