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

毕业设计（论文）毕 业 设 计题 目：刚性系统的隐式RK方法学 院： 数理学院 专业名称： 信息与计算科学 学 号： 201241210127 学生姓名： 丁楠 指导教师： 汪玉霞 2016年05月15日摘要本文主要介绍单步隐式Runge Kutta方法,简要的介绍了Gauss型隐式Runge Kutta方法、Radau型隐式Runge Kutta方法和Lobatto型隐式Runge Kutta方法。并利用这些基本的隐式Runge Kutta方法来对刚性方程组进行数值求解，并将隐式Runge Kutta方法与显式经典Runge Kutta方法求解的结果进行对比，说明两种数值解法的优缺点。关键词：刚性系统 隐式Runge Kutta方法 单步方法 Newton迭代法AbstractThis paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-Kutta 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 compared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of numerical methods.Keywords: Stiff system Implicit Runge-Kutta method One step method Newton iterative method目录1.绪论11.1刚性方程11.2隐式RK方法的研究意义21.3 RK方法的研究现状32.单步RK方法的收敛性和稳定性52.1单步RK方法的收敛性52.2单步RK方法的稳定性63.三类隐式RK方法83.1引言83.2 Gauss型隐式RK方法93.3 Radau型隐式RK方法103.4 Lobatto型隐式RK方法114隐式RK方法的实现134.1非线性系统的改进134.2简化的Newton迭代法135数值实验与结果分析15参考文献18附录211. 绪论1.1刚性方程对于一般的线性常系数系统y=Ay+(t)A为mm的矩阵，特征值为i(i=1,2,m)。定义123 若一个系统满足（1）Rei0, i=1,2,m（2）maxiReiminiRei=R1其中R为刚性比，则这个系统称为刚性系统。定义227 若线性系统y=Ay x0,T或非线性系统y=f(x,y) x0,T的矩阵A或Jacobi矩阵fy的特征值i满足max1imRei1则其是刚性的。定理1（解的存在性与唯一性）（1）对于所有x,yD，函数fx,y是连续的；（2）对于任何x,y,(x,y*)D，存在常数L，是函数满足fx,y-f(x,y*)Ly-y*则初值问题y=fx,y axbya= 有唯一解。其中y=(y1,y2,ym)T，D=x,y|axb,-ab;-yi,i=1,2,m 。其中L被称为Lipschitz常数定义3 如果一个常微分系统的Lipschitz常数L很大(大于20)，则它是刚性的。1.2隐式RK方法的研究意义在常微分方程及常微分方程组的数值解法中，Runge Kutta方法是目前应用最为广泛的数值解法之一，同时又具有误差小，精度高的特点。尽管显式Runge Kutta方法能够非常准确、快速的给出大部分常微分方程组的数值解。但是在化学、自动控制电力系统等领域中，会出现一些病态的常微分方程组，也就是刚性方程组。刚性方程组对于数值解法的稳定性要求苛刻，比如方程组y1=-0.01y1-99.99y2, y10=2y2=-100y2, y20=1将其表示为矩阵形式：y1y2=-0.01-99.990-100y1y2令A=-0.01-99.990-100发现A特征值为：1=-100，2=-0.01，刚性比s=12=100001。方程组的解为：y1x=e-100x+e-0.01xy2x=e-100x 快瞬态 慢瞬态解由快瞬态和慢瞬态两部分构成。由于慢瞬态的部分，y1x衰减变得十分缓慢。当自变量变到x=391时，函数值还未下降到初值的1%，求解区间至少取为(0,391)。另一方面，由于快瞬态的部分，y2x衰减的非常快，因此步长要取得非常小。从绝对稳定性的方面来看，如果用四阶显式经典RK方法求解，绝对稳定区间要求h(-2.78,0)，则要求h0.0278），计算结果就完全不可信了。对于刚性方程组，显式方法已经远远不能胜任了，一般采用绝对稳定性更好的方法（如隐式Runge Kutta方法）进行求解，本文采用单步隐式Runge Kutta方法，而对于隐式方法中的级值得求解，本文采用Newton迭代法。1.3 RK方法的研究现状研究基于标准模型方程的Runge Kutta方法的常见形式为：yi+1=yi+hj=1rjkj k1=fxi,yi kj=f(xi+ih,yi+hl=1j-1jlkl) (j=2,3,r) （显式）yi+1=yi+hj=1rjkj kj=f(xi+ih,yi+hl=1rjlkl) (j=1,2,3,r) （隐式）因为Runge Kutta方法是比较成熟的常微分方程数值解法。所以如今主要是对于经典的Runge Kutta方法进行完善和扩充，在一定的条件下，提高级数以提高精度。或者是将Runge Kutta方法与某些领域结合使用。在1994年，费景高1给出了一种显式Runge Kutta并行方法，从而实现Runge Kutta方法在多处理机上的应用。1997年，Enenkel和Jackson2实现了Runge Kutta方法的对角隐式并行改进。1999年，廖文远和李庆扬5给出了一类求解刚性常微分方程的半隐式多步Runge Kutta方法。2000年，张诚坚和余洪兵3针对非线性延迟系统构造了一类并行预校算法，给出其算法的局部误差估计，数值实验表明该算法是有效的，且具一定的可比性。2003年，李爱雄4通过对传统单支方法的计算格式进行改造，得到了解DDEs的两类单支并行算法单支并行预校算法和单支并行算法，并对方法的收敛性和稳定性做出了分析。2008年，庞丽君和朱永忠6给出了一类随机微分方程Runge Kutta方法的指数稳定性。2.单步RK方法的收敛性和稳定性2.1单步RK方法的收敛性对于常微分初值问题 y=fx,y axbya= (1)的单步显式公式 yi+1=yi+h(xi,yi,h) i=0,1,n-1y0= (2)局部截断误差可以表示为 Ri+1=yxi+1-yxi+hxi,yi,h i=0,1,n-1 3定理216：设y(x)为(1)的解，yii=0n为(2)的解。如果：（1）存在常数c，使得Ri+1chp+1 (i=0,1,2,n-1)（2）存在a0，使得maxx,yD0ah(x,y,h)yL其中D=x,y|axb,yx-yyx+记c0=cLeLb-a-1，则当hmina,pc0 时，有E(h)c0hp证:由(3)得yxi+1=yxi+hxi,y(xi),h+Ri+1 (i=0,1,n-1) (4)将(4)与(2)相减yxi+1-yxi=yxi-yi+hxi,yxi,h-xi,yi,h+ Ri+1 i=0,1,n-1 由yx0-y0=0知道，在i=0时，yxi-yic0hp成立。现在假设0ik-1时也是成立的。在hpc0时有：yxi-yic0pc0p= (i=0,1,n-1)进而 xi,yxi,h-xi,yi,h =xi,i,hyyxi-yiLyxi-yi 0ik-1其中i是yxi和yi之间的数。于是定理结合条件与(4)式，可以得到 yxi+1-yi+1yxi-yi+h xi,yxi,h-xi,yi,h+ Ri+1yxi-yi+Lhyxi-yi+chp+1=1+Lhyxi-yi+chp+1 0ik-1 从而y(xk-yk)1+Lhkyx0-y0+1+Lhk-11+Lh-1chp+1因为yx0-y0=0及1+LhkeLkheL(b-a)，得到yxk-ykcLeLb-a-1hp=c0hp所以当i=k时yxi-yic0hp也是成立的。（证毕）对于单步显式格式而言，当f(x,y)和f(x,y)y在D内连续时，那么定理1的条件（2）总是满足的。所以满足上述条件的单步显式公式，只要也满足相容性条件x,yx,0=f(x,y(x)那么它一定是收敛的。2.2单步RK方法的稳定性定义4 对于初值问题(1),若yii=0n是由(2)得到，zii=0n是下面扰动问题的解：zi+1=zi+hxi,zi,h+i+1 (i=0,1,2,n-1)z0=+0 如果存在正常数C，0及h0，使所有的(0,0，h(0,h0，当max0ini时，有max0inyi-ziC就称(2)是稳定的或者零稳定的。定理316 在定理1的条件下，单步显式公式(2)是稳定的。3.三类隐式RK方法3.1引言对于初值问题(1)，隐式Runge Kutta方法的一般形式为yi+1=yi+hj=1rbjkj (5) kj=fxi+cih,yi+hl=1rajlkl j=1,2,3,r (6)其中，xi=x0+ih，n=0,1,2,，为数轴上的离散点列；h为步长，yi为解y(xi)的近似值；c1，c2，cr称为隐式Runge Kutta方法的步长；b1，b2，br为权系数。令A=ajl，j=1,2,3,r。称A为系数矩阵，因此上式可以简化表示为c1c2a11a1rar1arrb1br使用Butcher阵的记号，上表可以表示为cAbT可以看到，如果A是一个主对角元素均为零的下三角矩阵，那么上表就可以表示一个显式Runge Kutta方法。如果A是一个主对角线非零的下三角矩阵，相应的Runge Kutta方法就是半隐式方法，(5)式等号左右就含有级值k1，k2，kr，计算ki时要解一个包含r个未知量ki的非线性方程组。如果A是满足Aij(ij)不全为零，则相应的方法就是隐式方法。若 系统是m维的，对于隐式Runge-Kutta方法实现的每一个递推步，均需要求解一个mr维的非线性方程组，一般用牛顿迭代法求解。例如01200001200-120162316这是Kutta得到的3级3阶显式Runge Kutta公式。而012100014140010162316与012100052413-124162316162316分别是3级4阶的半隐式Runge Kutta公式和隐式Runge Kutta公式。要求出具体的Runge Kutta公式，一般有两种办法。一种是将(6)在点xi,yi处进行泰勒展开，并且代入到(5)中，再与y(xi+h)在xi处的泰勒展开式进行对比，从而确定参数c，b，A。另一种方法就是将微分方程转化为等价的积分方程，用数值积分得到表达式，再进行对比得到参数。基于后一种方法，可以构造得到以下三种不同的隐式Runge-Kutta方法。3.2 Gauss型隐式RK方法17设c1，c2，cs为Ps2c-1=0的根，这里Ps(t)是0,1上的s次Legendre多项式，0ci1，i=1,s。求s级的2s阶的Gauss型Runge Kutta公式的参数c，b，A需要以下步骤:（1）求出s次的Legendre多项式Ps2c-1=0的s个零点c1，c2，cs；（2）由线性方程组j=1sbjcjk-1=1k k=1,s确定系数bj，j=1,s；（3）计算系数矩阵A=CW-1在这个基础上，就给出了s=1,2,5，p=2s的一系列Gauss型隐式Runge Kutta方法。定理4910 如果一个数值方法应用于方程y=y,其中为复常数，则对于所有，Re0，当n时，得到的数值解趋于零，则称这种方法是A稳定的。定理5910 s级Gauss型隐式Runge Kutta方法是2s阶的，其稳定函数是s,sPade近似且是A稳定的。以下列出s=3，p=6的Gauss型隐式Runge Kutta方法之一：12-15101212+151053629-1515536-1530536+152429536-1524536+153029+1515536518495183.3 Radau型隐式RK方法17这种方法的参数c，b，A需要下面几个步骤求得：（1）求多项式Pst-Ps-1(t)的零点c1，c2，cs，并指定c10，cs=1；（2）计算系数bk，bk=01Pst-Ps-1(t)t-ckPsck-Ps-1(ck)dt k=1,s（3）计算系数矩阵A，A=CW-1例如s=3，p=5的RadauA型Runge Kutta方法之一为：4-6104+610188-76360296-16961800-2+36225296+1696180088+76360-2-3622516-63616+6361916-63616+63619定理631 s级RadauA型Runge Kutta方法和s级RadauA型Runge Kutta方法是2s-1阶的，稳定函数是(s-1,s)次对角Pade近似。这两种都是A稳定的。3.4 Lobatto型隐式RK方法17这种方法的参数c，b，A需要下面几个步骤求得：（1）求多项式Pst-Ps-2(t)的零点c1，c2，cs，并指定c1=0，cs=1；（2）计算系数bk，bk=01Pst-Ps-2(t)t-ckPsck-Ps-2(ck)dt k=1,s（3）计算系数矩阵A，A=D-1(WT)-1N-CTD列出4级6阶RadauA型隐式Runge Kutta方法05-5105+5101011+5120025-512011-512011225+135120512025-1351200-1+512025+5120512-1-5120112112512512112定理731 s级LobattoA，B，C型隐式Runge Kutta方法是2s-2阶的，LobattoA和B型Runge Kutta方法的方法的稳定函数是(s-1,s-1)对角Pade近似，LobattoC型隐式Runge Kutta方法是(s-2,s)次对角Pade近似。所以这些方法都是A稳定的。4隐式RK方法的实现4.1非线性系统的改进对于s级隐式隐式Runge Kutta方法Yi=yn+hj=1saijfxn+cjh,Yj i=1,s (7)yn+1=yn+hj=1sbjfxn+cjh,Yj (8)令 zi=Yi-y (9)于是(7)变为 zi=j=1saijfxn+cjh,yn+zj (10)当得到的(10)解z1，z1时，由显式(8)可以很快得到解yn+1。若隐式Runge Kutta方法的系数矩阵式非奇异的，那么(10)可以写为z1zs=Ahf(xn+c1h,yn+z1)hf(xn+csh,yn+zs) (11)所以(8)可以看成与yn+1=yn+i=1sdizi等价的，其中d1,ds=b1,bsA-1 (12)比如s=3，p=5的RadauA型Runge Kutta方法中，d=0,0,1。4.2简化的Newton迭代法对于一般的非线性微分方程，系统(10)必须用迭代的方法求解。本文采用Newton迭代法。Newton迭代法应用于系统(10)，每次迭代都需要一个系数矩阵I-ha11fyxn+c1h,yn+z1ha1sfyxn+csh,yn+zs-has1fyxn+c1h,yn+z1hassfyxn+csh,yn+zs的线性方程组。定义5 若s阶矩阵B=bij，m阶矩阵A=aij，且有BA=b11Ab1sAbs1AbssA则称运算是B与A的Kronecker积。为了简化(10)的计算，我们用近似值Jfyxn,yn代替Jacobi矩阵fyxn+cih,yn+zi的值，于是方程(10)的简化Newton迭代法变为I-hAJZk=-Zk+h(AI)F(Z(k) Z(k+1)=Z(k)+Z(k) (13)这里Z(k+1)=(z1k,zs(k)T是解的第k次近似，Zk=Z1k,Zs(k)T是增量，F(Z(k)是FZk=fxn+c1h,yn+z1(k),fxn+csh,yn+zs(k)T的缩写。每一次的迭代需要s次由函数f的求值和一个NS维线性方程组的求解。因为矩阵I-hAJ对于所有的迭代是相同的，所以其LU分解在一个积分步内只要进行一次，进而减少了计算时间。5数值实验与结果分析问题：y1=10-7y1+sinx y10=1y2=10-3y2+cosx y20=1 写成矩阵形式就是：y1y2=10-70010-3y1y2+sinxcosx很明显该方程组的刚性比R=10-310-7=1041,因此是个刚性方程组。取步长为0.1，在区间0,1内分别用四级四阶经典显式Runge Kutta方法yi+1=yi+16k1+2k2+2k3+k4k1=fxi,yi k2=fxi+12h,yi+12hk1 k3=fxi+23h,yi+23hk2 k4=fxi+h,yi+hk3 和二级四阶隐式Runge Kutta方法yi+1=yi+12k1+k2 k1=f(xi+3-36h,yi+14k1h+3-2312k2h)k2=f(xi+3+36h,yi+14k2h+3+2312k1h)进行求解。得到结果如下表一：真值结果xi y1 y2 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-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-01 1.3783901e+00 1.7846058e+00 1.0000000e+00 1.4596978e+00 1.8429313e+00表二：隐式结果xi y1 y2 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-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表三：显式结果xi y1 y2 0.0000000e+00 1.0000000e+00 1.0000000e+00 1.0000000e-01 1.0049958e+00 1.0999336e+00 2.0000000e-01 1.0199334e+00 1.1988802e+00 3.0000000e-01 1.0446635e+00 1.2958521e+00 4.0000000e-01 1.0789390e+00 1.3898814e+00 5.0000000e-01 1.1224175e+00 1.4800297e+00 6.0000000e-01 1.1746645e+00 1.5653972e+00 7.0000000e-01 1.2351579e+00 1.6451319e+00 8.0000000e-01 1.3032934e+00 1.7184382e+00 9.0000000e-01 1.3783901e+00 1.7845846e+00 1.0000000e+00 1.4596978e+00 1.8429111e+00很明显的看到，在保留八位小数的情况下，二级四阶隐式Runge Kutta方法与真值无异，能够精确到小数点后七位以上。而经典Runge Kutta只能精确到小数点后四位。因此对于刚性方程组，相同步长下，隐式Runge Kutta方法的精度比显式Runge Kutta方法高得多。并且由于步长小，在实际的求解区间里面，涉及的递推次数少，从而函数的计算量就小。参考文献1 费景高. 并行显式 RungeKutta 公式的实现J. 计算机工程与设计, 1994 (5): 34-40.2 Enenkel R F, Jackson K R. DIMSEMs-diagonally implicit single-eigenvalue methods for the numerical solution of stiff ODEs on parallel computersJ. Advances in Computational Mathematics, 1997, 7(1-2): 97-133.3 张诚坚, 余红兵. 非线性延迟系统的一类并行预校算法J. 华中理工大学学报, 2000, 28(11): 104-106.4 李爱雄, 刘伟丰, 张诚坚, 等. 求解 DDEs 的多导龙格库塔方法的渐近稳定性 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).7 Curtiss C F, Hirschfelder J O. Integration of stiff equationsJ. Proc. Nat. Acad. Sci, 1952, 38(235): 1.8 Gear C W. 常微分方程初值问题的数值解法J. 1978.9 Butcher J C. Implicit runge-kutta processesJ. Mathematics of Computation, 1964, 18(85): 50-64.10 Ehle B L. High order A-stable methods for the numerical solution of systems of DEsJ. BIT Numerical 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 Statistical Computing, 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-implicit mid-point rule for stiff systems of ordinary differential equationsJ. Numerische Mathematik, 1983, 41(3): 373-398.16 孙志忠, 袁慰平, 闻震初. 数值分析M. 东南大学出版社, 2002.17 刘德贵, , 费景高. 刚性大系统数字仿真方法M. 河南科学技术出版社, 1996.18 Cockburn B, Shu C W. The RungeKutta discontinuous Galerkin method for conservation laws V: multidimensional systemsJ. Journal of Computational 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 Discretizations 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, 7(2): 433-437.22 Zhong X, Shu C W. A simple weighted essentially nonoscillatory limiter for RungeKutta discontinuous Galerkin methodsJ. Journal of Computational Physics, 2013, 232(1): 397-415.23 Lambert J D, Lambert J D. Computational methods in ordinary differential equationsM. London: Wiley, 1973.24 Zhu J, Zhong X, Shu C W, et al. RungeKutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshesJ. Journal of Computational Physics, 2013, 248: 200-220.25 Jayakumar T, Maheshkumar D, Kanagarajan K. Numerical solution of fuzzy differential equations by Runge-Kutta method of 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. Numerical 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 (ICTAP) 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 Computation, 2014, 247: 589-606.30 Jimenez J C, Sotolongo A, Sanchez-Bornot J M. Locally linearized runge kutta method of dormand and princeJ. Applied Mathematics and Computation, 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=zeros(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*k21,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第五章程序：syms 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+c