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

1、 创新实验论文 题 目： 刚性常微分问题的数值解法 课程名称： 创新实验 学 院： 理学院 专 业： 数学与应用数学 年 级： 应数131 学 号： 1307010239 234 236 姓 名： 袁蕊 张蕾 刘霖 指导教师： 罗贤兵 2015年07月14日目录第一章 绪论3 选题背景3 刚性问题的算法3 引言4第二章 刚性问题5第三章 预备知识8第四章 计算实验15附页20第1章 绪论 自然界和工程技术的很多现象,其数学模型是常微分方程（组）的初值问题，普通的常微分方程的数值解法已经比较成熟，理论比较完整，也有许多方法可供选择。但，有一类常微分方程组，求解值时遇到相当大的困难，这类常微分方程

2、组解的分量有的变化很快，有的变化缓慢，常常出现这种现象：变化快的分量很快趋于它的稳定值，而变化慢的分量缓慢趋于它的稳定值。从数值解的观点看来，当变化快时应该用小步长积分，当变化快的分量已趋于稳定，或者说已经没有变化快的分量时就应该用较大的步长积分，但是理论和实际都说明，很多方法特别是显示方法的步长任不能放大，否者便出现数值不稳定现象，即误差急剧增加，已致掩盖了真解，使求解过程无法继续进行。常微分方程组的这种性质叫做刚性， 我们考虑一阶常微分方程初值问题的数值解法。 (1.1)常微分方程的解能用初等函数、特殊函数或它们的级数与积分形式表达的非常之少，用解析办法只能求解线性常系数等特征类型的常微分

3、方程。在实际问题中归结出来的求解微分方程的方法只要依靠数值解法。所谓数值解法，就是通过某种离散化办法，将微分方程转化为差分方程来求解。求方程（2-1）的数值解，即对一系列离散节点建立求的近似值的递推格式，由此求得解y(x)在各节点的近似值，n=0，1,2，。相邻两个节点的间距称为步长，这时节点。因此，这样得到的数值解法也称为差分方法。初值问题（1.1）的数值解法，可区分为两大类： （1）单步法：此类方法在计算上的近似值时只用到了前一点上的信息。如Euler法、Runge-Kutta法和Taylor级数法这就是这类方法的典型代表。 （2）多步法：此类方法在计算时，除了需要点的信息外，还需要前面若

4、干个点上的信息。线性多步法是这类方法的典型代表。 本文讨论的是几种隐式方法（向后欧拉，梯形公式，改进欧拉）和隐式Runge-Kutta第二章 刚性问题 选取的步长h必须很小，满足h1/100,才能保证绝对稳定性要求。对于非线性常微分方程初值问题 若初值问题是稳定的，即dy/dx 0.用欧拉法进行数值求解时，h应满足 。若 ，h应满足h=2/M。在方程组的情况下，例如一阶常系数线性方程组 这里的A=，y=.记A的特征值为，对稳定的初值问题应满足Re.用欧拉法数值求解时，为了保证计算的稳定性，h的选取应满足 由下面的例子可以知道，当比值很大时，h很小，这时计算不熟很多，耗时很长，给实际计算带来了很

5、大的困难。例，某一物理现象可归结为一个线性方程组 （1.2）其中x为时间变量，而 A的特征值分别为。 方程（2-11）的解为 （1.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(01-11.5);holdon;h2=ezplot(1/2)*exp(-2*x)-(1/2)*(cos(40*x)+sin(40*x

6、)*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)holdon legend(y1,y2,y3) 图像：图一 由于在开始的一段时间量x，解曲线变化激烈，对方程进行数值求解时，自然要求数值有较高的精度，而对较大的时间量x,解曲线变化缓慢，因此，对数值方法的精度不必苛刻的要求，但就数值方法的稳定性而言，它并不随时间量x的大小而改变。例如对初值问题（1.2）用欧拉折现法，步长必须满足h1称式（1.2）为刚性方程

7、，比值R为刚性比。 第三章 预备知识 隐式方法11欧拉公式（显示） 考虑初值问题为了求得它在等距离散点上的数值解，首先将（1.1）离散化。设，将式（1.4）离散化的办法有Taylor展开法、数值微分法及数值积分法如果在点将做Taylor展开，得那么当h充分小时，略去误差项，用近似值代替近似代替，并注意到，便得 （1.7）其中 用（1.4）求解（1.1）的方法称为Euler方法。 如果利用差商代似替代微商，那么可得 （1.8）在（1.8）若用近似替代近似替代，同样得到递推公式（1.7）. 如果在上对积分，得 （1.9）那么对上式右端积分用左矩形求积公式，并用近似替代近似替代，也可得递推公式（1.

8、7）. 我们知道，在xOy平面上，微分方程的解称为积分曲线，积分曲线上一点（x,y）的切线斜率等于函数f（x,y）D的值。如果D中每一点，都画上一条以f（x,y）在这点的值为斜率并指向x增加方向的有向线段，即在D上作出了一个由方程确定的方向场，那么方程的解f=y（x）.从几何上看，就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向一致。 从初始点出发，过这点的积分曲线为y=y(x)，斜率为.设在附近y(x)可用过点的切线近似表示，切线方程为.当时，的近似值为，并记为，这就是得到时计算的近似公式.当时，的近似值为，并记为.于是就得到当的近似公式 .重复上面方法，一般可得当的计算的

9、近似公式. 如果，则上面公式就是（1.7）.将连续起来，就得到一条折线，所以Euler方法又称为折线法。 由公式（1.4）看出，已知便可算出.已知，便可算出，如此继续下去，这种只用前一步的值便可计算出的递推公式称为单步法。1.2向后欧拉公式（隐式） 若在（1.6）中，其右边的积分由数值积分的右矩形公式近似，并用替代替代，则可得到 （1.10）并称（1.10）为后退Euler公式。1.3梯形公式（隐式） 若在公式（1.6）中，其右边积分用数值梯形积分公式近似，并用替代，替代，则可得到梯形方法公式 (1.11) 梯形方法同头退Euler方法一样都是隐式单步法。对于隐式方法，通常采用迭代法。 对后退

10、Euler方法，先Euler方法计算，并将它作为初值，即，再把它代入（1.7）的右端，便得到后退Euler方法的迭代公式为 (1.12) 同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为 (1.13)1.4改进欧拉（隐式）1.5Runge-Kutta法 龙格-库塔方法的基本思路是想办法计算f(x,y)在某点上的函数值，然后对这些函数值做数值线性组合，构造出一个近似的计算公式；再把近似的计算公式和解的泰勒展开式相比较，便得前面的若干项相吻合，从而达到较高的精度，一般的显示R-K方法的形式如下： (1.14)当式（1.14）的布局截断误差达到时成公式（1.14）为p阶r段R-K方法。 类

11、似于显示方法我们可以得到Runge-Kattu隐式方法 （i=） 常用隐式R-K方法：1级2阶中点公式：2级2阶梯形公式：2级4阶R-K公式： 在Matlab中，专门提供勒求解微分方程的ode函数，如ode45，ode23，ode113，ode15s，ode23s等等。Ode求解器分为变步长和固定布场两种求解模式，其中，ode45是采用R-K法的变步长求解器，和它一样的还有ode23。目前，ode45是使用最多的求解函数，主要用于求解非刚性常微分方程。（如果求解时长时间没有相应，则方程是刚性的，可以换用ode23求解），ode函数的调用方式大都类似，以ode45为例，常用的语法格式有：T,Y=

12、ode45(odefun,tspan,y0);函数问题类型精确度说明ode45非刚性中等采用算法为4-5阶Runge-Kutta法，大多数情况下首选的函数ode23非刚性低基于Bogacki-Shampine2-3阶Runge-Kutta公式，在精确要求不高的场合，以及对于轻度刚性方程，ode23的效率可能好于ode45ode113非刚性低到高基于变阶次Adams-Bashforth-MoutlonPECE算法。在对误差要求严格的场合或者输入参数doefun代表的函数本身计算量很大情况比ode45效率高。ode113可以看成一个多步解算器，因为它会利用前几次时间节点上的解计器当前时间节点的解。

13、因此它不适应与非连续系统。ode15s刚性低到中基于数值差分公式（后向差分公式，BDFs也叫Gear方法），因此效率不是很高。同ode113一样，ode15s也是一个多步计算器。当ode45求解失败，或者非常慢，并且怀疑问题是刚性的，或者求解微分代数问题时可以考虑用ode15sode23s刚性低基于修正的二阶Rosenbrock公式。由于是单步解算器，当精度要求不高时，它效率可能会高于ode15s。他可以解决一些ode15s求解起来效率不太高的刚性问题。ode23t适度刚性低ode23可以用来求解微分代数方程。ode23tb刚性低当方程是刚性的，并且求解要求精度不高时可以使用。第4章 实验计算

14、 例如:刚性常微分方程这个初值问题的解析解为，那么它在0.1处的解为？解:理论解：1.向后欧拉(隐式) 源程序functiony=Euler_2(fun,x0,y0,xN,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;forn=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);fork=1:3z1=y(n)+h*feval(fun,x(n+1),z0);ifabs(z1-z0)1e-3break;endz0=z1;endy(n+1)=z1;endT=x,yfunctionz=f

15、(x,y)z=-15*y;迭代100次迭代10000次源程序作图比较：（红色线是迭代的结果，蓝色线是准确值）functiony=Euler_2(fun,x0,y0,xN,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;forn=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);fork=1:3z1=y(n)+h*feval(fun,x(n+1),z0);ifabs(z1-z0)1e-3break;endz0=z1;endy(n+1)=z1;endT=x,y;plot(x,y,r)

16、由图可看出当取端点较小迭代次数较大时，误差较小。2.梯形公式 （隐式） 源代码functiony=Euler_3(fun,x0,y0,xN,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;forn=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);fork=1:3z1=y(n)+h/2*(feval(fun,x(n),y(n)+feval(fun,x(n+1),z0);ifabs(z1-z0)rk4(f,0,0.1,1000) 总结：经多次试验发现向后欧拉，梯形公式，改进欧拉公式都是取端