微分方程常用的两种数值解法：欧拉方法与龙格-库塔法

-PAGE.z.四川师范大学本科毕业论文微分方程常用的两种数值解法：欧拉方法与龙格—库塔法学生姓名***院系名称数学与软件科学学院专业名称信息与计算科学班级2006级4班学号20060640**指导教师***四川师范大学教务处二○一○年五月微分方程常用的两种数值解法：欧拉方法与龙格—库塔法学生姓名：***指导教师：**【内容摘要】微分方程是最有生命力的数学分支，在自然科学的许多领域中，都会遇到常微分方程的求解问题。当前计算机的开展为常微分方程的应用及理论研究提供了非常有力的工具，利用计算机解微分方程主要使用数值方法，欧拉方法和龙格——库塔方法是求解微分方程最典型常用的数值方法。本文详细研究了这两类数值计算方法的构造过程，分析了它们的优缺点，以及它们的收敛性，相容性，及稳定性。讨论了步长的变化对数值方法的影响和系数不同的同阶龙格—库塔方法的差异。通过编制C程序在计算机上实现这两类方法及对一些典型算例的结果分析比拟，能更深切体会它们的功能，优缺点及适用场合，从而在实际应用中能对不同类型和不同要求的常微分方程会选取适当的求解方法。关键词：显式单步法欧拉〔Euler〕方法龙格—库塔〔Runge—Kutta〕方法截断误差收敛性Twocommonlyusednumericalsolutionofdifferentialequations：EulermethodandRunge-KuttamethodStudentName:*iongShiyingTutor：ZhangLi【Abstract】Thedifferentialequationisthemostvitalitybranchinmathematics.Inmanydomainsofnaturalscience,wecanmeettheordinarydifferentialequationsolutionquestion.Currently,thedevelopmentofcomputerhasprovidedthee*tremelypowerfultoolfortheordinarydifferentialequationapplicationandthefundamentalresearch,thecomputersolvingdifferentialequationmainlyusesvaluemethod.TheEulermethodandtheRunge—Kuttamethodarethemosttypicalcommonlyvaluemethodtosolvethedifferentialequation.Thisarticledissectsthestructureprocessofthesetwokindsofvaluescommonlyvaluemethodtosolvetheanalysestheirgoodandbadpoints,totheirastringency,thecompatibility,andthestabilityhasmadetheproof.Atthesametime,thearticlediscussthelengthofstridetothenumericalmethodchanginginfluenceandthedifferenceofthecoefficientdifferentsamestepRunge—kuttamethod.ThroughestablishingCprogramonthecomputercanrealizethesetwokindofmethods,Anglicizingsomemodelsofcalculatee*ampleresultcansincerelyrealizetheirfunction,theadvantageanddisadvantagepointsandthesuitablesituation,thusthesuitablesolutionmethodcanbeselectedtosolvethedifferenttypeandthedifferentrequestordinarydifferentialequationinthepracticalapplication.Keywords:E*plicitsingle-stepprocessEulermethodRunge—Kuttamethodtruncationerrorconvergence目录-PAGE.z.微分方程常用的两种数值解法：欧拉方法与龙格—库塔法前言常微分方程的形成与开展是和力学、天文学、物理学以及其他科学技术的开展密切相关的。数学其他分支的新开展，如复变函数、群、组合拓扑学等，都对常微分方程的开展产生了深刻的影响，当前计算机的开展更是为常微分方程的应用及理论研究提供了非常有力的工具。牛顿在研究天体力学和机械力学的时候，利用了微分方程这个工具，从理论上得到了行星运动规律。后来，法国天文学家勒维烈和英国天文学家亚当斯使用微分方程各自计算出那时尚未发现的海王星的位置。这些都使数学家更加深信微分方程在认识自然、改造自然方面的巨大力量，微分方程也就成了最有生命力的数学分支。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值，利用计算机解微分方程主要使用数值方法。 本文主要讨论一阶常微分方程初值问题〔1.1〕在区间上的数值解法，其中为关于，的函数，为给定的初始值，将上述问题的准确解记为。该问题常用的数值解法有：欧拉〔Euler〕方法、龙格—库塔〔Runge—Kutta〕方法及一些常用的线性多步法。本文重点介绍欧拉〔Euler〕方法和龙格—库塔〔Runge—Kutta〕方法。并对这两种方法编制程序，体会它们的功能、优缺点及适用场合，对不同类型常微分方程会选取适当的求解方法。１根本概念和准备知识一阶常微分方程初值问题是： 其中是平面上*个区域上的连续函数，式〔〕的微分方程一般有无穷多个解，式〔1.1.2〕是确定解的初始条件，如果一元函数对一切满足〔1〕；〔2〕；〔3〕存在而且；则称是初值问题〔1.1〕在区间上的解。误差：假定在计算时，用到的前一步的值是准确的即，把用计算得到的近似值记为，估计误差：=y(*n+1)，这种误差称为局部截断误差。如果不作这一假定，在每一步计算时除局部截断误差以外，还有由于前一步不准确而引起的误差，称为总体截断误差。收敛性：对于解初值问题的数值方法，我们希望它产生的数值解收敛于初值问题的准确解，"收敛性〞的一般定义为：对于所有满足引理1.1条件的初值问题〔1.1〕，如果有一种显式单步法：产生的近似解，对于任意固定的，均有，则称该显式单步法是收敛的。相容性：显式单步法〔〕称为与原微分方程相容，如果〔〕成立。并称式〔〕为相容性条件。稳定性：在实际计算中，一方面初始值不一定准确，往往带有一定误差，同时由于计算机字长有限，在计算过程中有舍入误差，而且这种误差在式〔〕的递推过程中会传递下去，对以后的结果产生影响。因此要考虑舍入误差的积累是否会得到控制，也即要考虑数值方法的稳定性。当时，假设舍入误差引起的后果是有限的，则可以认为该方法是数值稳定的。2欧拉方法2.1欧拉方法简介对常微分方程初值问题〔1.1〕用数值方法求解时，我们总是认为〔1.1〕的解存在且唯一。 欧拉方法是解初值问题的最简单、最原始的数值方法，它是显式单步法。下面介绍几种导出欧拉法的途径，每个途径皆可以推导出更为有效的数值方法。〔1〕Taylor展开在点将作Taylor展开得：当充分小时略去误差项，并注意到，得，以近似代替，以近似代替，且用"=〞代替"〞得差分方程初值问题：，，()用式()求解初值问题〔1.1〕的方法称为欧拉方法。〔2〕数值微分由导数的定义知，对于充分小的整理得，对此作相似的处理也可以得到欧拉方法()。〔3〕数值积分在区间对积分得〔〕用数值积分的左矩形公式计算式〔〕右端的积分，得，于是同样可以得到欧拉方法(2.1.1)。〔4〕多项式插值利用解和其导数在点的值，作一次埃尔米特插值，得到关于的插值多项式：，用近似代替就得到欧拉方法。〔5〕待定系数法在第步，和，利用这两个值估计出下一步的，将的值与估计值作线性组合：，其中，为待定系数。为确定这两个参数，要求这个估计值对和〔为常数〕准确成立。如果，则，得到方程：，得。如果，则，这样有：，说明，这样估计值为：，即为欧拉方法。欧拉方法的几何意义：由点斜式得切线方程等步长为，则，可由切线方程算出：，逐步计算出在点的值：，，用分段的折线逼近函数为"折线法〞而非"切线法〞，除第一个点是曲线上的切线，其它都不是。00图1欧拉方法的几何意义2.2欧拉方法的截断误差，收敛性，相容性，稳定性设，把在处展开成泰勒级数，即再由欧拉方法：两式相减得欧拉方法的局部截断误差为：，假设在上充分光滑，且令，则，故欧拉方法是一阶方法或具有一阶精度。欧拉方法的增量函数就是，由引理1.3、引理1.4知当满足Lipschitz条件时欧拉方法是收敛的而且是相容的。用欧拉方法求解典型方程〔〕的计算公式为：，有。要让，必须有，因此欧拉方法的绝对稳定域为：，当为实数时，绝对稳定区间为。在复平面上，是以1为半径、以为圆心的内部。3龙格—库塔法3.1龙格—库塔法的根本思想为了导出龙格—库塔法的一般公式，我们取如下的线性组合形式： 〔〕其中 〔〕即,；；a21,a31除c1=0外均为待定系数。显然用公式〔〕每计算一个新值要计算函数的值s次，又因每个都能以一种明显的方式由,计算出来，故将公式〔3.2.3〕称为s级显式龙格—库塔法。s级显式龙格—库塔法又可以写成下面既简洁又直观的阵列形式：03.2二、三、四级龙格—库塔法常见的二级二阶龙格—库塔法有：中点方法〔取〕二阶Heun方法〔取〕取s=3，完全仿照上述方法推导出三阶龙格—库塔公式。这时参数满足以下条件这个方程组解也不唯一，从而可以得到不同的三级三阶龙格—库塔法。两个较为常见的三级三阶龙格—库塔法是：Kutta方法〔取=，=1，a21=，a31=，a32=2，=，=，=〕将它代入〔〕得()三阶Heun方法〔取=，=，a21=，a31=0，a32=，=，=，=〕将它代入〔〕得通常人们所说的龙格—库塔法是指四阶而言的。取s=4，我们同样可以仿照二阶的情形推导出此公式。这样就得到如下含12个未知数但仅有10个方程的方程组此方程组的解同样不是唯一的，从而有许多不同的四级四阶龙格—库塔法，最常用的两个四阶公式是：经典龙格—库塔法〔〕RKG(Runge-Kutta-Gill)方法()经典龙格—库塔方法的几何解释：将〔〕式中的改记为，从而得到经典龙格—库塔方法的计算公式为只考虑区间上的解曲线。由()式可知，曲线在处的斜率，在处的斜率，而和则是在中点处的斜率的两个近似值，如图2：图2图2另一方面有：，用Simpson积分公式逼近此式右端的积分，得上式右端要计算的三个值。由前面的结果，可取，，，见图3：图3图3代入后得到〔〕另外，从()式容易看出，当与无关时，公式(3.3.5)实际上就是用标准的Simpson积分公式计算积分得到的，因此，可以将经典龙格——库塔方法视为是用Simpson积分公式的推广形式计算积分而得到的数值方法。经典龙格—库塔法的稳定性：增量函数：，其中代入后得于是有=从而得出龙格—库塔方法的绝对稳定区间为1经典龙格—库塔方法的算法框图：开始开始读入输出=结束图4经典龙格—库塔法算法框图4应用举例建模及其分析4.1欧拉方法解题及其数学模型问题提出例假定*公司的净资产因资产本身产生了利息而以4%的年利率增长，同时，该公司以每年100万的数额支付职工工资。净资产的微分方程为〔t以年为单位〕分别以初始值问题：用Euler公式预测公司24后的净资产趋势。模型建立分析：这是求微分方程的数值积分，为的是预测公司24年后的净资产趋势。确定变量：设净资产是时间的微分函数，不妨设变量t为时间〔以年为单位〕。设为第n年的净资产，为第n+1年的净资产，利息以每年4%的速度增长，且公司每年支付职工工资为100万，则第n+1年的净资产增长数额为，由于已得第n年的净资产为，则第n年的净资产加上第n+1年的净资产增长数额就得到第n+1年的净资产。归纳公式：〔4-1〕确定其微分方程的标准形式，这就是例的微分方程模型。解决问题分别以代入，*表示利息赢利低于工资支出的数额，y表示利息赢利与工资支出平衡的数额，z表示利息赢利高于工资支出的数额，计算结果见表1.表1nn11460.25003540..13834.92625004165.0721418.425003581.614768.32425004231.6831375.1425003624.8615699.05625004300.9441330.1425003669.8616627.01925004372.9851283.3525003716.6517552.125004447.961234.6825003765.32184748271184.0725003815.93193938581131.4325003868.5720308.87725004691.1291076.6925003923.3121221.23225004778.77101019.7625003980.2422130.08125004869.9211960.54625004039.452335.284525004964.7212898.96825004101.0324-63.304225005063.3从表1可以看到当利息赢利低于工资的支出，公司的净资产逐年减少，以致净资产为负值；当利息赢利与工资的支出平衡时，公司的净资产每年保持不变；当利息赢利超过工资的支出，公司的净资产稳步增长。再用欧拉方法求解，每计算一步，依次需要计算1次、2次和4次函数的值，为了比拟在计算量一样的条件下近似解的精度，步长分别取0.025、0.05、0.1。我对欧拉的C程序做了一些优化，参加了计算误差，计算结果如下〔表示近似解、表示准确解、表示误差〕：ThisisEulermethod:pleaseinputabandy0：011pleaseinputN：40*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.025000y[1]=1.025000y1[1]=1.024695e[1]=-0.000305*[2]=0.050000y[2]=1.049405y1[2]=1.048809e[2]=-0.000597*[3]=0.075000y[3]=1.073258y1[3]=1.072381e[3]=-0.000878*[4]=0.100000y[4]=1.096596y1[4]=1.095445e[4]=-0.001151*[5]=0.125000y[5]=1.119451y1[5]=1.118034e[5]=-0.001417*[6]=0.150000y[6]=1.141854y1[6]=1.140175e[6]=-0.001679*[7]=0.175000y[7]=1.163832y1[7]=1.161895e[7]=-0.001937*[8]=0.200000y[8]=1.185410y1[8]=1.183216e[8]=-0.002194*[9]=0.225000y[9]=1.206609y1[9]=1.204159e[9]=-0.002450*[10]=0.250000y[10]=1.227451y1[10]=1.224745e[10]=-0.002706*[11]=0.275000y[11]=1.247954y1[11]=1.244990e[11]=-0.002964*[12]=0.300000y[12]=1.268134y1[12]=1.264911e[12]=-0.003223*[13]=0.325000y[13]=1.288009y1[13]=1.284523e[13]=-0.003486*[14]=0.350000y[14]=1.307593y1[14]=1.303841e[14]=-0.003753*[15]=0.375000y[15]=1.326900y1[15]=1.322876e[15]=-0.004024*[16]=0.400000y[16]=1.345941y1[16]=1.341641e[16]=-0.004300*[17]=0.425000y[17]=1.364730y1[17]=1.360147e[17]=-0.004583*[18]=0.450000y[18]=1.383278y1[18]=1.378405e[18]=-0.004873*[19]=0.475000y[19]=1.401594y1[19]=1.396424e[19]=-0.005170*[20]=0.500000y[20]=1.419689y1[20]=1.414214e[20]=-0.005475*[21]=0.525000y[21]=1.437572y1[21]=1.431782e[21]=-0.005790*[22]=0.550000y[22]=1.455251y1[22]=1.449138e[22]=-0.006113*[23]=0.575000y[23]=1.472735y1[23]=1.466288e[23]=-0.006447*[24]=0.600000y[24]=1.490032y1[24]=1.483240e[24]=-0.006792*[25]=0.625000y[25]=1.507149y1[25]=1.500000e[25]=-0.007149*[26]=0.650000y[26]=1.524093y1[26]=1.516575e[26]=-0.007518*[27]=0.675000y[27]=1.540872y1[27]=1.532971e[27]=-0.007900*[28]=0.700000y[28]=1.557490y1[28]=1.549193e[28]=-0.008297*[29]=0.725000y[29]=1.573955y1[29]=1.565248e[29]=-0.008708*[30]=0.750000y[30]=1.590273y1[30]=1.581139e[30]=-0.009134*[31]=0.775000y[31]=1.606449y1[31]=1.596872e[31]=-0.009577*[32]=0.800000y[32]=1.622489y1[32]=1.612452e[32]=-0.010037*[33]=0.825000y[33]=1.638397y1[33]=1.627882e[33]=-0.010515*[34]=0.850000y[34]=1.654180y1[34]=1.643168e[34]=-0.011013*[35]=0.875000y[35]=1.669842y1[35]=1.658312e[35]=-0.011530*[36]=0.900000y[36]=1.685388y1[36]=1.673320e[36]=-0.012068*[37]=0.925000y[37]=1.700823y1[37]=1.688194e[37]=-0.012629*[38]=0.950000y[38]=1.716151y1[38]=1.702939e[38]=-0.013212*[39]=0.975000y[39]=1.731376y1[39]=1.717556e[39]=-0.013820*[40]=1.000000y[40]=1.746504y1[40]=1.732051e[40]=-0.014453thisisimproveEulermethodpleaseinputabandy0:：011pleaseinputN：20*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.050000y[1]=1.048869y1[1]=1.048809e[1]=-0.000060*[2]=0.100000y[2]=1.095561y1[2]=1.095445e[2]=-0.000116*[3]=0.150000y[3]=1.140345y1[3]=1.140175e[3]=-0.000169*[4]=0.200000y[4]=1.183437y1[4]=1.183216e[4]=-0.000221*[5]=0.250000y[5]=1.225017y1[5]=1.224745e[5]=-0.000272*[6]=0.300000y[6]=1.265236y1[6]=1.264911e[6]=-0.000325*[7]=0.350000y[7]=1.304219y1[7]=1.303841e[7]=-0.000378*[8]=0.400000y[8]=1.342074y1[8]=1.341641e[8]=-0.000433*[9]=0.450000y[9]=1.378896y1[9]=1.378405e[9]=-0.000492*[10]=0.500000y[10]=1.414766y1[10]=1.414214e[10]=-0.000553*[11]=0.550000y[11]=1.449756y1[11]=1.449138e[11]=-0.000618*[12]=0.600000y[12]=1.483927y1[12]=1.483240e[12]=-0.000687*[13]=0.650000y[13]=1.517337y1[13]=1.516575e[13]=-0.000762*[14]=0.700000y[14]=1.550035y1[14]=1.549193e[14]=-0.000841*[15]=0.750000y[15]=1.582066y1[15]=1.581139e[15]=-0.000927*[16]=0.800000y[16]=1.613472y1[16]=1.612452e[16]=-0.001020*[17]=0.850000y[17]=1.644289y1[17]=1.643168e[17]=-0.001121*[18]=0.900000y[18]=1.674551y1[18]=1.673320e[18]=-0.001230*[19]=0.950000y[19]=1.704288y1[19]=1.702939e[19]=-0.001349*[20]=1.000000y[20]=1.733529y1[20]=1.732051e[20]=-0.001478结论这道题用普通的微分方程也能列式求解，关键在于如何预测假设干年后的净资产趋势，用普通的微分方程就无法进展预测，且人工计算量相当大，这里我使用欧拉方法可以计算出精度以及误差，通过电脑运行程序就可以预测出假设干年后的净资产情况，欧拉方法的使用使得解题更加方便且准确。4.2龙格—库塔解题及其数学模型问题提出例两种果树寄生虫，其数量分别是其中一种寄生虫以吃另一种寄生虫为生，两种寄生虫的增长函数如以下常微分方程组所示：问题：预测3年后这一对寄生虫的数量。模型建立分析：这是一个典型的常微分方程例题，要求预测出3年后这对寄生虫的数量。确定变量：假定时间是寄生虫数量的积分函数，不妨设变量t为时间〔以年为单位〕。由题知有两种寄生虫u和v，u寄生虫年增长函数为，v寄生虫年增长函数为，初始值：u寄生虫为1.6，v寄生虫为1.2，由于其中一种寄生虫以吃另一种寄生虫为生，我们可建立u和v的关联函数f(u,v),g(u,v)。归纳后得到公式：〔4-2〕〔4-2〕即为例所述问题的微分方程模型。解决问题在此题中用Euler预估—校正公式取h=1，计算结果如表2.表2t/年u(t)v(t)11.61.221.024571.2683430.6409121.336640.3912111.41077我把龙格—库塔的C程序进展编辑后得到结果:〔e[i]表示误差〕pleaseinputabandy0:：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.095446y1[1]=1.095445e[1]=-0.000000*[2]=0.200000y[2]=1.183217y1[2]=1.183216e[2]=-0.000001*[3]=0.300000y[3]=1.264912y1[3]=1.264911e[3]=-0.000001*[4]=0.400000y[4]=1.341642y1[4]=1.341641e[4]=-0.000001*[5]=0.500000y[5]=1.414215y1[5]=1.414214e[5]=-0.000002*[6]=0.600000y[6]=1.483242y1[6]=1.483240e[6]=-0.000002*[7]=0.700000y[7]=1.549196y1[7]=1.549193e[7]=-0.000003*[8]=0.800000y[8]=1.612455y1[8]=1.612452e[8]=-0.000004*[9]=0.900000y[9]=1.673324y1[9]=1.673320e[9]=-0.000004*[10]=1.000000y[10]=1.732056y1[10]=1.732051e[10]=-0.000005结论从上面的结果可以分析，用每一种方法计算节点=0.1、0.20.9、1.0，上的值都需要计算4次的值，即它们的计算量根本一样，其结果是经典的龙格—库塔方法的精度最好，龙格—库塔方法的推导是基于Taylor级数方法的，因而在使用高阶龙格—库塔方法计算时可以很准确的推算出寄生虫每一年的数量。4.3分别使用二阶、三阶龙格—库塔法解初值问题对一些特殊的微分方程，使用欧拉方法和低阶的龙格—库塔方法也能到达很高的精度，例如：微分方程的解析解是一次函数，则用欧拉方法求得的数值解与准确解相符，微分方程的解析解是二次函数，则用二阶龙格—库塔方法求得的数值解与准确解相符。微分方程的解析解是三次多项式，用三阶龙格—库塔方法求得的数值解与准确解相符。建立模型建立初值问题模型：(4.1)用不同的龙格—库塔法解〔4.1〕初值问题模型〔1〕用二阶龙格库塔方法求解初值问题〔4.1〕结果如下：thisisthesecond-orderRunge-Kutta(Heun)methodpleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.110000y1[1]=1.110000e[1]=0.000000*[2]=0.200000y[2]=1.240000y1[2]=1.240000e[2]=0.000000*[3]=0.300000y[3]=1.390000y1[3]=1.390000e[3]=0.000000*[4]=0.400000y[4]=1.560000y1[4]=1.560000e[4]=0.000000*[5]=0.500000y[5]=1.750000y1[5]=1.750000e[5]=0.000000*[6]=0.600000y[6]=1.960000y1[6]=1.960000e[6]=0.000000*[7]=0.700000y[7]=2.190000y1[7]=2.190000e[7]=0.000000*[8]=0.800000y[8]=2.440000y1[8]=2.440000e[8]=0.000000*[9]=0.900000y[9]=2.710000y1[9]=2.710000e[9]=0.000000*[10]=1.000000y[10]=3.000000y1[10]=3.000000e[10]=0.000000thisisimproveEulermethodpleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.110000y1[1]=1.110000e[1]=0.000000*[2]=0.200000y[2]=1.240000y1[2]=1.240000e[2]=0.000000*[3]=0.300000y[3]=1.390000y1[3]=1.390000e[3]=0.000000*[4]=0.400000y[4]=1.560000y1[4]=1.560000e[4]=0.000000*[5]=0.500000y[5]=1.750000y1[5]=1.750000e[5]=0.000000*[6]=0.600000y[6]=1.960000y1[6]=1.960000e[6]=0.000000*[7]=0.700000y[7]=2.190000y1[7]=2.190000e[7]=0.000000*[8]=0.800000y[8]=2.440000y1[8]=2.440000e[8]=0.000000*[9]=0.900000y[9]=2.710000y1[9]=2.710000e[9]=0.000000*[10]=1.000000y[10]=3.000000y1[10]=3.000000e[10]=0.000000thisisthesecond-orderRunge-Kutta(middle)methodpleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.110000y1[1]=1.110000e[1]=0.000000*[2]=0.200000y[2]=1.240000y1[2]=1.240000e[2]=0.000000*[3]=0.300000y[3]=1.390000y1[3]=1.390000e[3]=0.000000*[4]=0.400000y[4]=1.560000y1[4]=1.560000e[4]=0.000000*[5]=0.500000y[5]=1.750000y1[5]=1.750000e[5]=0.000000*[6]=0.600000y[6]=1.960000y1[6]=1.960000e[6]=0.000000*[7]=0.700000y[7]=2.190000y1[7]=2.190000e[7]=0.000000*[8]=0.800000y[8]=2.440000y1[8]=2.440000e[8]=0.000000*[9]=0.900000y[9]=2.710000y1[9]=2.710000e[9]=0.000000*[10]=1.000000y[10]=3.000000y1[10]=3.000000e[10]=0.000000〔2〕用三阶龙格—库塔方法求解初值问题模型：〔4.2〕结果如下：thisisthethird-orderHeunmethod:pleaseinputabandy0：021pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.248000y1[1]=1.248000e[1]=0.000000*[2]=0.400000y[2]=1.624000y1[2]=1.624000e[2]=0.000000*[3]=0.600000y[3]=2.176000y1[3]=2.176000e[3]=0.000000*[4]=0.800000y[4]=2.952000y1[4]=2.952000e[4]=0.000000*[5]=1.000000y[5]=4.000000y1[5]=4.000000e[5]=-0.000000*[6]=1.200000y[6]=5.368000y1[6]=5.368000e[6]=0.000000*[7]=1.400000y[7]=7.104000y1[7]=7.104000e[7]=-0.000000*[8]=1.600000y[8]=9.256000y1[8]=9.256000e[8]=-0.000000*[9]=1.800000y[9]=11.872000y1[9]=11.872000e[9]=-0.000000*[10]=2.000000y[10]=15.000000y1[10]=15.000000e[10]=-0.000000thisiskuttamethodpleaseinputabandy0：021pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.248000y1[1]=1.248000e[1]=0.000000*[2]=0.400000y[2]=1.624000y1[2]=1.624000e[2]=0.000000*[3]=0.600000y[3]=2.176000y1[3]=2.176000e[3]=0.000000*[4]=0.800000y[4]=2.952000y1[4]=2.952000e[4]=0.000000*[5]=1.000000y[5]=4.000000y1[5]=4.000000e[5]=0.000000*[6]=1.200000y[6]=5.368000y1[6]=5.368001e[6]=0.000000*[7]=1.400000y[7]=7.104000y1[7]=7.104000e[7]=-0.000000*[8]=1.600000y[8]=9.256001y1[8]=9.256001e[8]=0.000000*[9]=1.800000y[9]=11.872001y1[9]=11.872001e[9]=0.000000*[10]=2.000000y[10]=15.000001y1[10]=15.000000e[10]=-0.000001结论如果常微分方程的解析解是任意的次多项式，总存在阶龙格—库塔方法，用它求解该微分方程的所得的数值解与准确解相符。证明：龙格—库塔方法的推导是基于Taylor级数方法的，阶龙格——库塔方法的局部截断误差为：，其中：，为常数，如果常微分方程的解析解是任意的次多项式，则有，因此局部截断误差等于零。从而上述结论成立。5步长算法模型及其分析下面我们分析步长对不同数值方法的影响。用欧拉方法解初值问题建立初值模型：〔4.3〕当h=0.2时thisisEulermethodpleaseinputabandy0：011pleaseinputN：5*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.200000y1[1]=1.183216e[1]=-0.016784*[2]=0.400000y[2]=1.373333y1[2]=1.341641e[2]=-0.031693*[3]=0.600000y[3]=1.531495y1[3]=1.483240e[3]=-0.048256*[4]=0.800000y[4]=1.681085y1[4]=1.612452e[4]=-0.068633*[5]=1.000000y[5]=1.826948y1[5]=1.732051e[5]=-0.094898当h=0.1时thisisEulermethod:pleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.100000y1[1]=1.095445e[1]=-0.004555*[2]=0.200000y[2]=1.191818y1[2]=1.183216e[2]=-0.008602*[3]=0.300000y[3]=1.277438y1[3]=1.264911e[3]=-0.012527*[4]=0.400000y[4]=1.358213y1[4]=1.341641e[4]=-0.016572*[5]=0.500000y[5]=1.435133y1[5]=1.414214e[5]=-0.020920*[6]=0.600000y[6]=1.508966y1[6]=1.483240e[6]=-0.025727*[7]=0.700000y[7]=1.580338y1[7]=1.549193e[7]=-0.031145*[8]=0.800000y[8]=1.649784y1[8]=1.612452e[8]=-0.037332*[9]=0.900000y[9]=1.717780y1[9]=1.673320e[9]=-0.044460*[10]=1.000000y[10]=1.784771y1[10]=1.732051e[10]=-0.052720用改良的欧拉方法解初值问题模型〔4.3〕。当h=0.2时thisisimproveEulermethod:pleaseinputabandy0：011pleaseinputN：5*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.186667y1[1]=1.183216e[1]=-0.003451*[2]=0.400000y[2]=1.348312y1[2]=1.341641e[2]=-0.006672*[3]=0.600000y[3]=1.493704y1[3]=1.483240e[3]=-0.010464*[4]=0.800000y[4]=1.627861y1[4]=1.612452e[4]=-0.015410*[5]=1.000000y[5]=1.754205y1[5]=1.732051e[5]=-0.022154当h=0.1时thisisimproveEulermethod:pleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.095909y1[1]=1.095445e[1]=-0.000464*[2]=0.200000y[2]=1.184097y1[2]=1.183216e[2]=-0.000881*[3]=0.300000y[3]=1.266201y1[3]=1.264911e[3]=-0.001290*[4]=0.400000y[4]=1.343360y1[4]=1.341641e[4]=-0.001719*[5]=0.500000y[5]=1.416402y1[5]=1.414214e[5]=-0.002188*[6]=0.600000y[6]=1.485955y1[6]=1.483240e[6]=-0.002715*[7]=0.700000y[7]=1.552514y1[7]=1.549193e[7]=-0.003320*[8]=0.800000y[8]=1.616474y1[8]=1.612452e[8]=-0.004023*[9]=0.900000y[9]=1.678166y1[9]=1.673320e[9]=-0.004846*[10]=1.000000y[10]=1.737867y1[10]=1.732051e[10]=-0.005816用二级二阶龙格—库塔方法〔二阶Heun方法〕解初值问题模型〔4.3〕。当h=0.2时thisisthesecond-orderRunge-Kuttamethod:pleaseinputabandy0：011pleaseinputN：5*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.184706y1[1]=1.183216e[1]=-0.001490*[2]=0.400000y[2]=1.344644y1[2]=1.341641e[2]=-0.003003*[3]=0.600000y[3]=1.488062y1[3]=1.483240e[3]=-0.004822*[4]=0.800000y[4]=1.619653y1[4]=1.612452e[4]=-0.007201*[5]=1.000000y[5]=1.742496y1[5]=1.732051e[5]=-0.010445当h=0.1时thisisthesecond-orderRunge-Kuttamethod:pleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.095625y1[1]=1.095445e[1]=-0.000180*[2]=0.200000y[2]=1.183572y1[2]=1.183216e[2]=-0.000356*[3]=0.300000y[3]=1.265449y1[3]=1.264911e[3]=-0.000538*[4]=0.400000y[4]=1.342374y1[4]=1.341641e[4]=-0.000733*[5]=0.500000y[5]=1.415162y1[5]=1.414214e[5]=-0.000948*[6]=0.600000y[6]=1.484431y1[6]=1.483240e[6]=-0.001191*[7]=0.700000y[7]=1.550663y1[7]=1.549193e[7]=-0.001470*[8]=0.800000y[8]=1.614246y1[8]=1.612452e[8]=-0.001794*[9]=0.900000y[9]=1.675493y1[9]=1.673320e[9]=-0.002173*[10]=1.000000y[10]=1.734671y1[10]=1.732051e[10]=-0.002620用三级三阶龙格—库塔方法〔Kutta方法〕解初值问题模型〔4.3〕。当h=0.2时thisiskuttamethod:pleaseinputabandy0：011pleaseinputN：5*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.183244y1[1]=1.183216e[1]=-0.000028*[2]=0.400000y[2]=1.341729y1[2]=1.341641e[2]=-0.000088*[3]=0.600000y[3]=1.483408y1[3]=1.483240e[3]=-0.000168*[4]=0.800000y[4]=1.612727y1[4]=1.612452e[4]=-0.000275*[5]=1.000000y[5]=1.732472y1[5]=1.732051e[5]=-0.000421当h=0.1时thisiskuttamethod:pleaseinputabandy0：011pleaseinputN：10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.095445y1[1]=1.095445e[1]=0.000001*[2]=0.200000y[2]=1.183217y1[2]=1.183216e[2]=-0.000001*[3]=0.300000y[3]=1.264915y1[3]=1.264911e[3]=-0.000004*[4]=0.400000y[4]=1.341648y1[4]=1.341641e[4]=-0.000007*[5]=0.500000y[5]=1.414225y1[5]=1.414214e[5]=-0.000011*[6]=0.600000y[6]=1.483256y1[6]=1.483240e[6]=-0.000016*[7]=0.700000y[7]=1.549214y1[7]=1.549193e[7]=-0.000021*[8]=0.800000y[8]=1.612479y1[8]=1.612452e[8]=-0.000027*[9]=0.900000y[9]=1.673355y1[9]=1.673320e[9]=-0.000034*[10]=1.000000y[10]=1.732094y1[10]=1.732051e[10]=-0.000043用四级四阶龙格—库塔方法〔经典Rung--Kutta方法〕解初值问题模型〔4.3〕。当h=0.2时thisisRunge-Kuttamethod:pleaseinputabandy0：011pleaseinputN：5*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.200000y[1]=1.183229y1[1]=1.183216e[1]=-0.000013*[2]=0.400000y[2]=1.341667y1[2]=1.341641e[2]=-0.000026*[3]=0.600000y[3]=1.483281y1[3]=1.483240e[3]=-0.000042*[4]=0.800000y[4]=1.612514y1[4]=1.612452e[4]=-0.000062*[5]=1.000000y[5]=1.732142y1[5]=1.732051e[5]=-0.000091当h=01时thisisRunge-Kuttamethod:pleaseinputabandy0:011pleaseinputN:10*[0]=0.000000y[0]=1.000000y1[0]=1.000000e[0]=0.000000*[1]=0.100000y[1]=1.095446y1[1]=1.095445e[1]=-0.000000*[2]=0.200000y[2]=1.183217y1[2]=1.183216e[2]=-0.000001*[3]=0.300000y[3]=1.264912y1[3]=1.264911e[3]=-0.000001*[4]=0.400000y[4]=1.341642y1[4]=1.341641e[4]=-0.000001*[5]=0.500000y[5]=1.414215y1[5]=1.414214e[5]=-0.000002*[6]=0.600000y[6]=1.483242y1[6]=1.483240e[6]=-0.000002*[7]=0.700000y[7]=1.549196y1[7]=1.549193e[7]=-0.000003*[8]=0.800000y[8]=1.612455y1[8]=1.612452e[8]=-0.000004*[9]=0.900000y[9]=1.673324y1[9]=1.673320e[9]=-0.000004*[10]=1.000000y[10]=1.732056y1[10]=1.732051e[10]=-0.000005从上面的结果可以看出，对不同的数值方法步长越少准确度越高，但步长的减少将导致计算量的增加，因此为了提高计算的效率，除了选择适当的计算方法外，还需考虑步长的恰中选取，使得数值解能到达预定的精度，计算效率也尽可能高。在单步法中每一步积分步长实际上是独立的，步长的选择具有相当的灵活性，根据具体情况合理地选择每一步步长，既可以保证精度要求，又可以减少计算量，从而也就能减少舍入误差。下面我们以四级四阶龙格—库塔方法为例，说明如何选取适宜的步长。分析：假设从点出发，以为步长经过一步计算