刚性微分方程组隐式龙格库塔方法.doc

《刚性微分方程组隐式龙格库塔方法.doc》由会员分享，可在线阅读，更多相关《刚性微分方程组隐式龙格库塔方法.doc（36页珍藏版）》请在三一办公上搜索。

1、毕 业 设 计题 目：刚性系统的隐式RK方法摘要本文主要介绍单步隐式方法,简要的介绍了型隐式方法、型隐式方法和型隐式方法。并利用这些根本的隐式方法来对刚性方程组进展数值求解，并将隐式方法与显式经典方法求解的结果进展比照，说明两种数值解法的优缺点。关键词：刚性系统 隐式方法 单步方法 迭代法AbstractThis paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-K

2、utta method and Lobatto implicit Runge-Kutta method. These basic Implicit Runge-Kutta methods are used to solve the stiff equations. These implicit Runge-Kutta methods iare pared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of nu

3、merical methods.Keywords:Stiff systemImplicitRunge-Kutta method One step method Newton iterative method目录1.绪论1刚性方程1隐式RK方法的研究意义2 RK方法的研究现状32.单步RK方法的收敛性和稳定性5单步RK方法的收敛性5单步RK方法的稳定性6RK方法8引言8 Gauss型隐式RK方法93.3 Radau型隐式RK方法103.4 Lobatto型隐式RK方法114隐式RK方法的实现13非线性系统的改良13简化的Newton迭代法135数值实验与结果分析15参考文献18附录211. 绪论

4、1.1刚性方程对于一般的线性常系数系统为的矩阵，特征值为。定义123假如一个系统满足12其中为刚性比，如此这个系统称为刚性系统。定义227 假如线性系统或非线性系统的矩阵或矩阵的特征值满足如此其是刚性的。定理1解的存在性与唯一性1对于所有，函数是连续的；2对于任何，存在常数，是函数满足如此初值问题有唯一解。其中，。其中被称为常数定义3如果一个常微分系统的常数很大(大于20)，如此它是刚性的。1.2隐式RK方法的研究意义在常微分方程与常微分方程组的数值解法中，方法是目前应用最为广泛的数值解法之一，同时又具有误差小，精度高的特点。尽管显式方法能够非常准确、快速的给出大局部常微分方程组的数值解。但是

5、在化学、自动控制电力系统等领域中，会出现一些病态的常微分方程组，也就是刚性方程组。刚性方程组对于数值解法的稳定性要求苛刻，比如方程组将其表示为矩阵形式：令发现特征值为：，刚性比。方程组的解为：快瞬态 慢瞬态解由快瞬态和慢瞬态两局部构成。由于慢瞬态的局部，衰减变得十分缓慢。当自变量变到时，函数值还未下降到初值的，求解区间至少取为。另一方面，由于快瞬态的局部，衰减的非常快，因此步长要取得非常小。从绝对稳定性的方面来看，如果用四阶显式经典方法求解，绝对稳定区间要求，如此要求。这样，在上就要计算步，计算量巨大，因此计算区间内的解时，舍入误差积累会特别严重。例如取求解区间为，用不同步长来计算和的值。利用

6、四阶显式经典方法求解如下：2.9802322e+172.9802322e+17真值很显然，保存八位有效数字的情况下，要保持良好的精度，步长要取得非常小，这就增加了计算量。而随着步长的加大，误差也会越来越大。当步长加大到绝对稳定与之外时即，计算结果就完全不可信了。对于刚性方程组，显式方法已经远远不能胜任了，一般采用绝对稳定性更好的方法如隐式方法进展求解，本文采用单步隐式方法，而对于隐式方法中的级值得求解，本文采用迭代法。1.3RK方法的研究现状研究基于标准模型方程的方法的常见形式为：显式 隐式因为方法是比拟成熟的常微分方程数值解法。所以如今主要是对于经典的方法进展完善和扩大，在一定的条件下，提高

7、级数以提高精度。或者是将方法与某些领域结合使用。在年，费景高1给出了一种显式并行方法，从而实现方法在多处理机上的应用。年，Enenkel和Jackson2实现了方法的对角隐式并行改良。年，廖文远和李庆扬5给出了一类求解刚性常微分方程的半隐式多步方法。年，X诚坚和余洪兵3针对非线性延迟系统构造了一类并行预校算法，给出其算法的局部误差估计，数值实验明确该算法是有效的，且具一定的可比性。年，李爱雄4通过对传统单支方法的计算格式进展改造，得到了解的两类单支并行算法单支并行预校算法和单支并行算法，并对方法的收敛性和稳定性做出了分析。年，庞丽君和朱永忠6给出了一类随机微分方程方法的指数稳定性。RK方法的收

8、敛性和稳定性2.1单步RK方法的收敛性对于常微分初值问题的单步显式公式局部截断误差可以表示为定理216：设为的解，为的解。如果：1存在常数，使得2存在，使得其中记，如此当时，有证:由得将与相减由知道，在时，成立。现在假设时也是成立的。在时有：进而其中是和之间的数。于是定理结合条件与式，可以得到从而因为与，得到所以当时也是成立的。证毕对于单步显式格式而言，当和在内连续时，那么定理1的条件2总是满足的。所以满足上述条件的单步显式公式，只要也满足相容性条件那么它一定是收敛的。2.2单步RK方法的稳定性定义4 对于初值问题,假如是由得到，是下面扰动问题的解：如果存在正常数，与，使所有的，当时，有就称是

9、稳定的或者零稳定的。定理316在定理1的条件下，单步显式公式是稳定的。RK方法3.1引言对于初值问题，隐式方法的一般形式为其中，为数轴上的离散点列；为步长，为解的近似值；，称为隐式方法的步长；，为权系数。令，。称为系数矩阵，因此上式可以简化表示为使用阵的记号，上表可以表示为可以看到，如果是一个主对角元素均为零的下三角矩阵，那么上表就可以表示一个显式方法。如果是一个主对角线非零的下三角矩阵，相应的方法就是半隐式方法，式等号左右就含有级值，计算时要解一个包含r个未知量的非线性方程组。如果是满足不全为零，如此相应的方法就是隐式方法。假如 系统是m维的，对于隐式Runge-Kutta方法实现的每一个递

10、推步，均需要求解一个维的非线性方程组，一般用牛顿迭代法求解。例如这是Kutta得到的级阶显式公式。而与分别是级阶的半隐式公式和隐式公式。要求出具体的公式，一般有两种方法。一种是将在点处进展泰勒展开，并且代入到中，再与在处的泰勒展开式进展比照，从而确定参数，。另一种方法就是将微分方程转化为等价的积分方程，用数值积分得到表达式，再进展比照得到参数。基于后一种方法，可以构造得到以下三种不同的隐式Runge-Kutta方法。3.2 Gauss型隐式RK方法17设，为的根，这里是上的次多项式，。求级的阶的Gauss型公式的参数，需要以下步骤:1求出次的多项式的个零点，；2由线性方程组确定系数，；3计算系

11、数矩阵在这个根底上，就给出了，的一系列型隐式方法。定理4910 如果一个数值方法应用于方程,其中为复常数，如此对于所有，和对于固定步长，当时，得到的数值解趋于零，如此称这种方法是稳定的。定理5910级型隐式方法是阶的，其稳定函数是近似且是稳定的。以如下出，的型隐式方法之一：3.3 Radau型隐式RK方法17这种方法的参数，需要下面几个步骤求得：1求多项式的零点，并指定，；2计算系数，3计算系数矩阵，例如，的型方法之一为：定理631级型方法和级型方法是阶的，稳定函数是次对角近似。这两种都是稳定的。3.4 Lobatto型隐式RK方法17这种方法的参数，需要下面几个步骤求得：1求多项式的零点，并

12、指定，；2计算系数，3计算系数矩阵，列出4级6阶型隐式方法定理731级，型隐式方法是阶的，和型方法的方法的稳定函数是对角近似，型隐式方法是次对角近似。所以这些方法都是稳定的。4隐式RK方法的实现4.1非线性系统的改良对于级隐式隐式方法令于是变为当得到的解，时，由显式可以很快得到解。假如隐式方法的系数矩阵式非奇异的，那么可以写为所以可以看成与等价的，其中比如，的型方法中，。4.2简化的Newton迭代法对于一般的非线性微分方程，系统必须用迭代的方法求解。本文采用迭代法。迭代法应用于系统，每次迭代都需要一个系数矩阵的线性方程组。定义5 假如阶矩阵，阶矩阵，且有如此称运算是与的积。为了简化的计算，我

13、们用近似值代替矩阵的值，于是方程的简化迭代法变为这里是解的第k次近似，是增量，是的缩写。每一次的迭代需要次由函数的求值和一个维线性方程组的求解。因为矩阵对于所有的迭代是一样的，所以其分解在一个积分步内只要进展一次，进而减少了计算时间。5数值实验与结果分析问题：写成矩阵形式就是：很明显该方程组的刚性比,因此是个刚性方程组。取步长为，在区间内分别用四级四阶经典显式方法和二级四阶隐式方法进展求解。得到结果如下表一：真值结果 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e-01 1.0049958e+00 1.0999384e+00 2.000

14、0000e-01 1.0199334e+00 1.1988893e+00 3.0000000e-01 1.0446635e+00 1.2958649e+00 4.0000000e-01 1.0789390e+00 1.3898974e+00 5.0000000e-01 1.1224175e+00 1.4800481e+00 6.0000000e-01 1.1746644e+00 1.5654174e+00 7.0000000e-01 1.2351579e+00 1.6451531e+00 8.0000000e-01 1.3032934e+00 1.7184598e+00 9.0000000e-

15、01 1.3783901e+00 1.7846058e+00 1.0000000e+00 1.4596978e+00 1.8429313e+00表二：隐式结果 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e-01 1.0049958e+00 1.0999384e+00 2.0000000e-01 1.0199334e+00 1.1988893e+00 3.0000000e-01 1.0446635e+00 1.2958649e+00 4.0000000e-01 1.0789390e+00 1.3898974e+00 5.0000000e

16、-01 1.1224175e+00 1.4800481e+00 6.0000000e-01 1.1746644e+00 1.5654173e+00 7.0000000e-01 1.2351579e+00 1.6451531e+00 8.0000000e-01 1.3032934e+00 1.7184598e+00 9.0000000e-01 1.3783901e+00 1.7846058e+00 1.0000000e+00 1.4596978e+00 1.8429313e+00表三：显式结果 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000

17、e-01 1.0049958e+001.0999336e+00 2.0000000e-01 1.0199334e+001.1988802e+00 3.0000000e-01 1.0446635e+001.2958521e+00 4.0000000e-01 1.0789390e+001.3898814e+00 5.0000000e-01 1.1224175e+001.4800297e+00 6.0000000e-01 1.1746645e+001.5653972e+00 7.0000000e-01 1.2351579e+001.6451319e+00 8.0000000e-01 1.303293

18、4e+001.7184382e+00 9.0000000e-01 1.3783901e+001.7845846e+00 1.0000000e+00 1.4596978e+001.8429111e+00很明显的看到，在保存八位小数的情况下，二级四阶隐式方法与真值无异，能够准确到小数点后七位以上。而经典只能准确到小数点后四位。因此对于刚性方程组，一样步长下，隐式方法的精度比显式方法高得多。并且由于步长小，在实际的求解区间里面，涉与的递推次数少，从而函数的计算量就小。参考文献1费景高. 并行显式 RungeKutta 公式的实现J. 计算机工程与设计, 1994 (5): 34-40.2 Enenk

19、el R F, Jackson K R. DIMSEMs-diagonally implicit single-eigenvalue methods for the numerical solution of stiff ODEs on parallel putersJ. Advances in putational Mathematics, 1997, 7(1-2): 97-133.3 X诚坚, 余红兵. 非线性延迟系统的一类并行预校算法J. 华中理工大学学报, 2000, 28(11): 104-106.4 李爱雄, X伟丰, X诚坚, 等. 求解 DDEs 的多导龙格库塔方法的渐近稳定性

20、 Asymptotic stability of multi-derivative Runge-Kutta method for delay differential equationsJ. 华中科技大学学报 (自然科学版), 2002, 6: 108-110.5 廖文远, 李庆扬. 一类求解刚性常微分方程的半隐式多步 RK 方法J. 清华大学学报: 自然科学版, 1999, 39(6): 1-4.6庞立君, 朱永忠. 类随机微分方程 Runge. Kutta 方法的指数稳定性J. 河海大学学报 (自然 科学版), 2008, 36(3).7Curtiss C F, Hirschfelder

21、J O. Integration of stiff equationsJ. Proc. Nat. Acad. Sci, 1952, 38(235): 1.8Gear C W. 常微分方程初值问题的数值解法J. 1978.9 Butcher J C. Implicit runge-kutta processesJ. Mathematics of putation, 1964, 18(85): 50-64.10 Ehle B L. High order A-stable methods for the numerical solution of systems of DEsJ. BIT Numer

22、ical Mathematics, 1968, 8(4): 276-278.11 Burrage K, Butcher J C, Chipman F H. An implementation of singly-implicit Runge-Kutta methodsJ. BIT Numerical Mathematics, 1980, 20(3): 326-340.12 Shampine L F. Implementation of implicit formulas for the solution of ODEsJ. SIAM Journal on Scientific and Stat

23、istical puting, 1980, 1(1): 103-118.13 Lindberg B. On smoothing and extrapolation for the trapezoidal ruleJ. BIT Numerical Mathematics, 1971, 11(1): 29-52.14 Lindberg B. IMPEX 2-a procedure for solution of systems of stiff differential equationsR. CM-P00069411, 1973.15 Bader G, Deuflhard P. A semi-i

24、mplicit mid-point rule for stiff systems of ordinary differential equationsJ. Numerische Mathematik, 1983, 41(3): 373-398.16 孙志忠, 袁慰平, 闻震初. 数值分析M. 东南大学, 2002.17 X德贵, , 费景高. 刚性大系统数字仿真方法M. 某某科学技术, 1996.18 Cockburn B, Shu C W. The RungeKutta discontinuous Galerkin method for conservation laws V: multid

25、imensional systemsJ. Journal of putational Physics, 1998, 141(2): 199-224.19 Monovasilis T, Kalogiratou Z, Simos T E. Exponentially fitted symplectic rungeKuttaNystrm methodsJ. Appl. Math. Inf. Sci, 2013, 7(1): 81-85.20 Hochbruck M, Pazur T. Implicit Runge-Kutta Methods and Discontinuous Galerkin Di

26、scretizations for Linear Maxwells EquationsJ. SIAM Journal on Numerical Analysis, 2015, 53(1): 485-507.21 Papadopoulos D F, Simos T E. A modified Runge-Kutta-Nystrm method by using phase lag properties for the numerical solution of orbital problemsJ. Applied Mathematics & Information Sciences, 2013,

27、 7(2): 433-437.22 Zhong X, Shu C W. A simple weighted essentially nonoscillatory limiter for RungeKutta discontinuous Galerkin methodsJ. Journal of putational Physics, 2013, 232(1): 397-415.23 Lambert J D, Lambert J D. putational methods in ordinary differential equationsM. London: Wiley, 1973.24 Zh

28、u J, Zhong X, Shu C W, et al. RungeKutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshesJ. Journal of putational Physics, 2013, 248: 200-220.25 Jayakumar T, Maheshkumar D, Kanagarajan K. Numerical solution of fuzzy differential equations by Runge-Kutta method o

29、f order fiveJ. International Journal of Applied Mathematical Science, 2012, 6: 2989-3002.26 Dimarco G, Pareschi L. Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equationsJ. SIAM Journal on Numerical Analysis, 2013, 51(2): 1064-1087.27 Miranker W L, Miranker A. Num

30、erical Methods for Stiff Equations: And Singular Perturbation ProblemsM. Springer Science & Business Media, 2001.28 Hadi M, Anderson M, Husein A. Using 4th order Runge-Kutta method for solving a twisted Skyrme string equationC/THE 4TH INTERNATIONAL CONFERENCE ON THEORETICAL AND APPLIED PHYSICS (ICTA

31、P) 2014. AIP Publishing, 2016, 1719: 030005.29 Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta method of dormand and princeJ. Applied Mathematics and putation, 2014, 247: 589-606.30 Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta method of do

32、rmand and princeJ. Applied Mathematics and putation, 2014, 247: 589-606.31 Wanner G, Hairer E. Solving ordinary differential equations IIM. Springer-Verlag, Berlin, 1991.附录绪论程序：a=0;b=1; %求解区间f1=(t,x,y)(-0.01*x-99.99*y);f2=(t,x,y)(-100*y);h=0.0001; %步长T=zeros(1,(b-a)/h)+1);X=zeros(1,(b-a)/h)+1);Y=zer

33、os(1,(b-a)/h)+1);T=a:h:b;X(1)=2; %初值Y(1)=1;for j=1:(b-a)/h k11=feval(f1,T(j),X(j),Y(j); k12=feval(f2,T(j),X(j),Y(j); k21=feval(f1,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12); k22=feval(f2,T(j)+h/2,X(j)+h/2*k11,Y(j)+h/2*k12); k31=feval(f1,T(j)+h/2,X(j)+h/2*k21,Y(j)+h/2*k22); k32=feval(f2,T(j)+h/2,X(j)+h/2*k2

34、1,Y(j)+h/2*k22); k41=feval(f1,T(j)+h,X(j)+h*k31,Y(j)+h*k32); k42=feval(f2,T(j)+h,X(j)+h*k31,Y(j)+h*k32); X(j+1)=X(j)+h*(k11+2*k21+2*k31+k41)/6; Y(j+1)=Y(j)+h*(k12+2*k22+2*k32+k42)/6; R=T X Y; endsave 0.0001结果.txt R ascii真值程序：y1=exp(-100)+exp(-0.01);y2=exp(-100);R=1 y1 y2;save 真值结果.txt R ascii第五章程序：s

35、yms x;a=0;b=1; h=0.1; f1=dsolve(Dy=0.0000001*y+sin(x),y(0)=1,x);f2=dsolve(Dy=0.001*y+cos(x),y(0)=1,x);X=a:h:b;Y1=double(subs(f1,x,X);Y2=double(subs(f2,x,X);R1=X Y1 Y2;xlswrite(结果,R1,真值)format longf1=(x,y1)(0.0000001*y1+sin(x);f2=(x,y2)(0.001*y2+cos(x);a=0;b=1; h=0.1; X=zeros(1,(b-a)/h+1);Y1=zeros(1,

36、(b-a)/h+1);Y2=zeros(1,(b-a)/h+1);X=a:h:b;Y1(1)=1;Y2(1)=1;syms k1 k2;for j=1:(b-a)/hx10=1 1;x20=1 1; tol=1e-10; x1=x10(1);x2=x20(1);y1=x10(2);y2=x20(2);K1=feval(k1,k2)(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x1,y1);K2=feval(k1,k2)(k2-feval(f1,X(j)+h*(3+sqrt(3)/6),Y1(j

37、)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),x1,y1); % 迭代函数G1=feval(k1,k2)(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x2,y2);G2=feval(k1,k2)(k2-feval(f2,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),x2,y2); % 迭代函数F1=K1 K2;F2=G1 G2;J1=double(subs(jacobian(k1-feval(

38、f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f1,X(j)+h*(3+sqrt(3)/6),Y1(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),k1 k2),k1,k2,x1,y1); J2=double(subs(jacobian(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f2,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2

39、*sqrt(3)/12)*k1+h*(1/4)*k2),k1 k2),k1,k2,x2,y2);x11=x10-F1/J1;x21=x20-F2/J2;while (norm(x11-x10)tol)x10=x11;x1=x10(1);y1=x10(2);K1=feval(k1,k2)(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x1,y1);K2=feval(k1,k2)(k2-feval(f1,X(j)+h*(3+sqrt(3)/6),Y1(j)+h*(3+2*sqrt(3)/12)*k1

40、+h*(1/4)*k2),x1,y1); F1=K1 K2;J1=double(subs(jacobian(k1-feval(f1,X(j)+h*(3-sqrt(3)/6),Y1(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f1,X(j)+h*(3+sqrt(3)/6),Y1(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),k1 k2),k1,k2,x1,y1); x11=x10-F1/J1;end while (norm(x21-x20)tol)x20=x21;x2=x20(1);y2=x20(2);G1=feval(

41、k1,k2)(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2),x2,y2);G2=feval(k1,k2)(k2-feval(f2,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),x2,y2);F2=G1 G2;J2=double(subs(jacobian(k1-feval(f2,X(j)+h*(3-sqrt(3)/6),Y2(j)+h*(1/4)*k1+h*(3-2*sqrt(3)/12)*k2);k2-feval(f2

42、,X(j)+h*(3+sqrt(3)/6),Y2(j)+h*(3+2*sqrt(3)/12)*k1+h*(1/4)*k2),k1 k2),k1,k2,x2,y2);x21=x20-F2/J2;end Y1(j+1)=Y1(j)+h*(x11(1)+x11(2)/2; Y2(j+1)=Y2(j)+h*(x21(1)+x21(2)/2;endR2=X Y1 Y2; xlswrite(结果,R2,隐式结果)f1=(x,y1)(0.0000001*y1+sin(x);f2=(x,y2)(0.001*y2+cos(x);a=0;b=1; h=0.1; X=zeros(1,(b-a)/h+1);Y1=zeros(1,(b-a)/h+1);Y2=zeros(1,(b-a)/h+1);X=a:h:b;Y1(1)=1; Y2(1)=1;for j=1:(b-a)/h k11=feval(f1,X(j),Y1(j); k12=feval(f2,X(j),Y2(j); k21=feval(f1,X(j)+h/2,Y1(j)+h/2*k11); k22=feval(f2,X(j)+h/2,Y2(j)+h/2*k11); k31=feval(f1,X(j)+h/2,Y1(j)+h/2*k21); k32=feval(f2,X(j)+h/2,Y2(j)+h/2*