第四讲-龙格-库塔方法

龙格-库塔方法法

3.2Runge-Kutta

3.

2.1显式Rimge・Kutta法的一般形式上节已给出与初值问题（

1.

2.1）等价的积分形式y（Xn+i）-y

（4）=（

3.

2.1）只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题（L

2.1）的数值方法，若用显式单步法Yn+1=居八，我），附=01，…（

3.

2.2）当数:…%/）=政2%），即数值求积用左矩形公式，它就是Euler法（

3.

1.2）,方法只有一阶精度，若取0（4，%，h）=；［/6,片）（，〃+4/（/，”））］⑶

2.3）+/4+I乙就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的,若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将⑶

2.1）右端积分表示为/（x J（x））dx=+a㈤）+（为］）注意，右端f中双Xn+%h）还不能直接得到，需要像改进Euler法（

3.

1.11）一样，用前面已算得的f值表示为（

3.

2.3）,一般情况可将（

3.

2.2）的表示为r工刖，”Xi%5=

3.

2.4其中无1=/（/，八）2-1舄砧）0=23…/）这里6吗，C=12…,r，J=12…,i-1）均为待定常数，公式⑶

2.2）,⑶

2.4）称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值（即%，…，匕）,而ki由前面（i-D个已算出的片也，…,如表示,故公式是显式的.例如当r=2时，公式可表示为共计8个参数待设定Taylor展开分析局部截断误差，使得精度达到三级，即误差为（、）0于是，r=3,p=3的8个参数所满足的关系式（C3不能为0）为+勺+公=1々2=玩1，%=玩1+运21122利C=70Kutta法function y=DELGKT3_kuttaf,h,a,b,yO,varvec formatlong;b-a/h;zeros N+l,1;yD二yO;x二a:h:b;var=findsymf;for i=2:N+l[xi-l yi-1];KI Funval f,varvec,[xi-l+h/2yi-l+Kl*h/2];Funval f,varvec,K2[xi-l+h yi-l-h*Kl+K2*2*h];Funval f,varvec,y K⑴3yi-l+h*Kl+4*K2+K3/6;endformat short;Suen法function y=DELGKT3_suenf,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;y1=yO;x=a:h:b;var=findsymf;for i=2:N+lKI=Funval f,varvec,[xi-l y i-l];K2二Funval f,varvec,[xi-l+h/3y i-l+Kl*h/3];K3二Funval f,varvec,[x i-1+2*h/3y i-l+K2*2*h/3];y i=yi-l+h*Kl+3*K3/4;endformat short;

3.四阶四级R-K方法举例经典四级Rung-Kutta法h八二占+仁/+止+i=K+2+246勺=人人八）k=八Xx+］况八+暂占）k3=八A+W力.+鼻£）k4=Z（x+为，乂+为x公）程序function y=DELGKT4_Rungkuttaf,h,a,b,yO,varvec formatlong;N=b-a/h;y二zeros N+l,1;yD=yO；x=a:h:b;var=findsymf;for i=2:N+lKI=Funval f,varvec,[xi-l y i-l];K2=Funval f,varvec,[xi-l+h/2y iT+Kl*h/2];K3=Funval f,varvec,[xi-l+h/2y iT+K2*h/2];K4=Funval f,varvec,[xi-l+h y iT+h*K3];yi=yi-l+h*Kl+2*K2+2*K3+K4/6;end formatshort;基尔jer、Gear法function y=DELGKT4_jer f,h,a,b,yO,varvec formatlong;N=b-a/h;C2=sqrt⑵；y=zeros N+l,1;y1=yO;x.a:h:b;var=findsymf;for i=2:N+lKI=Funval f,varvec,[x i-l y i-l];K2=Funvalf,varvec,[xi-l+h/2yi-l+Kl*h/2];K3=Funval f,varvec,[xi-l+h/2y i-l+Kl*h*C2T/2+K2*h*2-C2/2];K4=Funvalf,varvec,[xi-1+h y i-1-K2*h*C2/2+K3*h*2+C2/2];yi=yi-l+h*Kl+2-C2*K2+2+C2*K3+K4/6;end formatshort;变型Rung-Kutta法function y=DELGKT4_qt f,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;y1=yO;x=a:h:b;var=findsymf;for i=2:N+lKI二Funval f,varvec,[xi-l y i-l];K2二Funval f,varvec,[xi-l+h/3y iT+Kl*h/3];K3=Funval f,varvec,[x i-1+2*h/3yiT-Kl*h/3+K2*h];K4=Funval f,varvec,[xi-l+h yi-l+h*K1-K2+K3];yi=yi-l+h*Kl+3*K2+3*K3+K4/8;end formatshort;说明四阶四级R-K方法的步长自动选取问题用四阶R-K方法求解初值问题（

1.

2.1）精度已经足够高了，但要从理论上给出误差限-y（当）1的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h通常是通过不同步长在计算机上的计算结果近似估计.设贝＞“）在％”处的值此=（%“），当兀+1=尤+/7时，y（%”+i）的近似为约+/），于是由四阶R-K方法有%=收+-斓仁好）1h若以5为步长，计算两步到醺小则有1+1=M/+1）一月处2分（红乙于是得丁（，稼+知一朗）3+1y i-淄

3.

2.13—\y2-y{k\l e它给出了误差的近似估计.如果15川（e为给定精度，如-中2_|『勺_再|1E-3）,则认为以2为步长的计算结果八十】满足精度要求，若1512口则还可放大步长,因此⑶

2.13）提供了自动选择步长的方法.二隐式半隐式R-K方法举例一种半隐式举例，该格式可以通过除法化为显式单步法格式，受到除法复杂度限制，容易造成耗时，不稳，相对精度降低罗赛布诺克半隐式法券+i=K+卬内+吗，k、=h[l-ha,f x,y]-1fx,y n n nk=h[l-hafx+,y+ck}~}fx+dh,y+dk2y n n21n xn2xfunction y=DELSBRKf,h,a,b,yO,varvecformat long;al=

1.40824829;a2=

0.59175171;cl=

0.17378667;c2=cl;wl=-

0.41315432;w2=

1.41315432;N=b-a/h;y=zeros N+l,1;yD=y0;x=a:h:b;var=findsymf;dy=diff f,varvec2;for i=2:N+lfl=Funval f,varvec,[xi-l y i-l];dyl=Funvaldy,varvec,[xi-1yi-l];kl=h*f1/1-h*al*dyl;dy2=Funvaldy,varvec,[xi-l+cl*h y i-l+c2*kl];f2=Funvalf,varvec,[xi-l+cl*h yi-l+c2*kl];k2=h*f2/l-h*a2*dy2;yi=yi-l+wl*kl+w2*k2;endformat short;单步法的收敛性与绝对稳定性

3.

33.

3.1单步法的收敛性定义

3.1设y（x）是初值问题（

1.

2.1）的精确解，八是单步法⑶

2.2）在X+.处产生的近似解，若瘟〃二y（Q则称方法⑶

2.2）产生的数值解收敛于）㈢）.实际上,定义中/是一固定点，当h-o时n-8,n不是固定的.因“=1二血n显然方法收敛，则在固定点x=/处的整体误差/=尸（％）一入=°（川），当p、i时40%=U.下面定理给出方法⑶

2.2）收敛的条件.定理

3.1设初值问题（L

2.1）的单步法⑶

2.2）是p阶方法（p》l）,且函数步对y满足Lipschitz条件，即存在常数L0,使对均有I”（x,y,h）-^（x,z,h）|L|y-z|则方法⑶

2.2）收敛，且％=尸（/）-片=伊）.定理证明略.绝对稳定性用单步法⑶

2.2）求数值解Ji，…，居，…，由于原始数据及计算过程舍入误差影响，实际得到的不是必而是其=居+外，其中5是误差，再计算下一步得到区+1=5+”0（巧，见，如以Euler法为例，若令a=%一匕,则Px+1=A+/（/,忆）-/（乙,〉*）］=［1+切（/4）］（

3.

3.1）♦如果11+lif|1,则从外计算到八+i误差不增长，它是稳定的.但如果条件不y满足就不稳定.例

3.4y=T00y,y

（0）=l,精确解为y（x）=°*用Euler法求解得居+i=入+蛆醺/=-100加入”=0,1，…若取h=

0.025,则为+1=一1・5/，当胃—8,胃一人|=8,而lim\_lim八一lOOx_AXT9灵-°,显然计算是不稳定的.X\T9,如果用后退Euler法

3.

1.5解此例,仍取h=

0.025,则丁/1=+%i=^«-2・5j\+i gn=2+，川外+1g,八显然当I入=O,计算是稳定的.由此看到稳定性与方法有关，也与力有关，在此例中《=70°.在研究方法的稳定性时，通常不必对一般的fx,y进行讨论，而只针对模型方程y1=2y,Re

203.

3.2这里兀可能为复数.规定Re“°是因为ReU0时微分方程⑶

3.2本身是不稳定的，而讨论数值方法

3.

2.2的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题L

2.1,若将fx,y在仆,Yn处线性展开，可得Rx,y处f醺，〃+£；“x一4+£x yy-yXf xx于是方程L

2.1可近似表示为y=移+g3n，yn以=£；Xn，yn它表明用模型方程⑶

3.2是合理的，至于模型方程

3.

3.2中所以用复数人是因为初值问题

1.

2.1如果是方程组，即yeRmjeKj则〃是mXm阶矩阵，其特征值可能是复数.当然对单个方程，入就是实数，此时只要规定丸0即可.用单步法⑶

2.2解模型方程⑶

3.2可得到%+i=

3.

3.3其中下㈤依赖所选方法，如用Euler法，则以+i=〃+〃必=1+hZy,以吟=1+U

3.

3.4n此时由⑶

3.1看到误差方程也为=矶研,与

3.

3.4是一样的.因此对一般单步法⑶

2.2误差方程也与⑶

3.3一致.下面再考虑二阶R-K方法有h=[1++-税2M,EhZ=1+税+士/ay22对四阶R-K方法，可得Eh2=1+LI+1hl2+-hl3+2/a«定义

3.2将单步法⑶

2.2用于解模型方程

3.

3.2,若得到⑶

3.3中的I11则称方法是绝对稳定的.在复平面上复变量以满足I EhN的区域，称为方法⑶

2.2的绝对稳定域，它与实轴的交点称为绝对稳定区间.例如对Euler法,|Eh；l|=|l+h；ql在复平面h兄上是以-1,0为圆心，以1为半径的单位圆域内部，当兄为实数时，则得绝对稳定区间为-2%元0,2因兄0,故有〈在例

3.4中“三=

0.02时方法稳定，而例中取一儿.100h=

0.025故不稳定.对后退Euler法

3.

1.5,片+1=以+加15rlm吟二因兄0,故口一川》1,其绝对稳定域是以1,0为圆心的单位圆外部，绝对稳定区间为-8^0,即对任何h0方法都是绝对稳定的.二阶R-K方法的绝对稳定区间为-2^

20.三阶R-K方法的绝对稳定区间为-

2.51粗

0.四阶R-K方法的绝对稳定区间为-2785K

0.例

3.5用经典四阶R-K方法计算初值问题y=-2oxo x i,xo=i步长取h=o.1及o.2,给出计算误差并分析其稳定性.解本题直接按R-K方法⑶

2.12的公式计算.因精确解为x=e-20x,其计算误差上7—|如表所示

0.

20.

40.

60.

81.0h=

0.

10.

0927950.

0120100.

0013660.

0001520.000017h=

0.

24.

9825.

0125.

0625.

03125.0从计算结果看到，h=

0.2时误差很大，这是由于在入=-20,h=

0.2时X h=-4,而四阶R-K方法的绝对稳定区间为［-

2.785,0］,故h=

0.2时计算不稳定，误差很大.而h=

0.1时此二-2,其值在绝对稳定区间［-

2.785,0］内，计算稳定，故结果是可靠的.讲解由于微分方程初值问题数值解公式求出的解尢，为，…，儿…是一个逐次递推的过程，因此原始数据误差及计算过程舍入误差对解的影响就是数值方法绝对稳定性研究的问题，如果由〃计算〃+1误差不增长，方法就是绝对稳定的.为使问题得到简化通常就是将方法用于解模型方程（

3.

3.2）,对于单步法得到的差分方程为片+1=双”㈤〃，由于模型方程的/阮川二由，代入Euler法分+1二%+也=居+板力=0+%㈤然，得颐瓦1）=1+成，对二阶方法，例R—K如，用改进Euler法=A[加+双耳+W«]=入+1=%+小//,”+勺+兄居+/,〃]♦于是1alzEQtZ=1+力兄+士瓦I2+上瓦13对三阶R-K方法有23!EhZ=1+以兄+1的兄2+1h6+-烟对四阶R-K方法有

2、/小只要欧刈＜1方法，就是绝对稳定的，这时的值当n增大式是减少的,故计算稳定.这时舍入误差影响可忽略不计，而当忸（“㈤1＞1,则卜」增大，方法不稳定，计算结果是不可靠的.因此用显式单步法必须使I矶“刈＜1,也就是步长选择要满足这一要求.对于隐式的梯形公式将模型方程矽，即/=由代入得+1=+51电」加11+6_2于是〃+1=—i一入注意兀〈°,于是有1」函2这就表明对任意步长h,梯形法都是绝对稳定的.对VA e0,8成立.练习题

3.4dy兀，，3兀—=sm y+cosx,—x——

2.对于一阶微分方程初值问题产22,用Euler法,吟=，、乙改进Euler法，二阶R-K方法求解，并作图比较dy

3.对于一阶微分方程初值问题区二一k+1°L5,用Eg法，改、y0=1进Euler法，二阶R-K法，三阶R-K法求解，并与真解作图比较，列出误差对比表格dy

4.对于一阶微分方程初值问题区=lnx+l，°*,、用三阶、四阶y0=1R-K法与真解作图对比

5.对于一阶微分方程初值问题2孙’°-x-1,R-K法的半隐式、四阶y0=l与真解作图对比——Q v2V Y

46.对一阶微分方程初值问题公一进行稳定性分析丁⑵二1a.得出二阶R-K法中步长的稳定区间域；b.自选两个步长（hl为区间内的数（如

0.1）,h2为区间外的数（如

0.3））得到两个不同数值解和精确解作图比较，并列出误差表格=〃+城+右后

3.

2.5其中公=〃/月=/4+夕用儿+如城

1.改进Euler法

3.

1.11就是一个二级显式R-K方法.参数5勺，与，与1取不同的值，可得到不同公式.

二、三级显式R・K方法对片2的显式『1方法方.

2.5,要求选择参数1，0，叼，与I,使公式的精度阶P尽量高，由局部截断误差定义T=yx-yx-h[cjx,yx+cf x+ah,y+b hk]

3.n+}n+i n nn2n2n2]{

2.6令％=双,对

3.

2.6式在%，%处按Taylor公式展开,由于/必4一-+1=式/+坊=%+沙」百3「百丁二+〃.=/%〃=右Z=2，片=/区,入+£4,〃以+a抱〃+41妁=///+///必+％i峭+*将上述结果代入⑶

2.6得h

2.Tn+i=g3■[]4，〃+力/必]+*-Mq工+与一+小以M乙+%区区,人力+薪]=1-q-与乜+8-匕22//区,其+；-匕血1我乐外力+薪要使公式

3.

2.5具有的阶p=2,即Tx+i=◎*,必须11l-q-c=0,--c^=0,--c Z=032722221即c】+G=1，”=呆次1=尚由此三式求心勺，与也I的解不唯「因r=2,由

3.

2.5式可知勺H0,于是有解1q=1-勺色=砥=⑶

2.8它表明使⑶

2.5具有二阶的方法很多，只要勺h0都可得到二阶精度R-K方法.若取0=L则q=；,町=1，则得改进Euler法

3.

1.11,若取巧722则得9=0，“2=瓯=5,此时

3.

2.5为乙yn+1=入++；氏乂占

3.

2.9乙乙其中％=//,居称为中点公式.改进的Euler法

3.

1.11及中点公式

3.

2.9是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为厂2只有4个参数，要达到p=3则在

3.

2.6的展开式中要增加3项，即增加三个方程，加上

3.

2.7的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法

3.

2.5当勺HO取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对厂3的情形，要计算三个k值，即^x,y,h=c/i+q后+q网nn其中占=fg以劣=fg+与阮儿+瓦/Q心=/+电比乙+%也占+与，月将匕出按二元函数在/，％处按Taylor公式展开，然后代入局部截断误差表达式，可得1+1+4-丁4-加4,八M=必可得三阶方法，其系数共有8个，所应满足的方程为+0+个=1%=瓦1,出=打1+421122eg的+Q%=-

3.

2.10J乙内%用2=76这是8个未知数6个方程的方程组，解也是不唯一的，通常c3Ho.一种常见的三级三阶R-K方法是下面的三级Kutta方法h z、=〃+三占+/+月

3.

2.1161h十井八后=//+§片乙乙=/4+瓦〃-嫡+2g附R-K的三级Kutta方法程序如下function y=DELGKT3_kutaf,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;y1=yO;x=a:h:b;var=findsymf;for i=2:N+lKI=Funval f,varvec,[x i-1yi-1];K2=Funval f,varvec,[xiT+h/2yi-l+Kl*h/2];K3=Funval f,varvec,[xi-l+h y iT-h*Kl+K2*2*h];yi=yi-l+h*Kl+4*K2+K3/6;%满足cl+c2+c3=l,1/64/61/6endformat short;323四阶R・K方法及步长的自动选择利用二元函数Taylor展开式可以确定⑶

2.4中r=4,p=4的R-K方法，其迭代公式为hc、k\+c2k2+c3k3+c4kJ%+I=%+其中％]=/,％,k=fx+ah,y+b hfx,%,而2n2n2[n+bhk%3=/%〃+〃3〃,n2k=fx+ah,y+b hk+b hk+b hk4n4n4{x422433共计13个参数待定，Taylor展开分析局部截断误差，使得精度达到四阶，即误差为心于是，r=4,p=4的13个参数c4不能为0引出了多种方案和挑战，如参数优化使阶数增加到5阶，得到四阶五阶R-K方法，mat lab中有程序ode45；四级四阶R-K方法的步长自动选取；结合新算法的应用算法构造；适应于新的领域实现求解；经典的四阶R-K方法是卜〃+1=入+—肉+2上/2与+/）（

3.

2.12）6%=/Ed）M=/（^+；力，Vx+2用）.=/（4+鼻，X.+=/⑸+4乂+肽3）上4其中也需满足仿+2+3+4=1，这里为（1/62/62/61/6）.它的局部截断误差1+1=°的‘），故p=4,这是最常用的四阶R-K方法，数学库中都有用此方法求解初值问题的软件.这种方法的优点是精度较高，缺点是每步要算4个右端函数值,计算量较大.例

3.3用经典四阶R-K方法解例

1.1的初值问题=一y+x+1）

（0）=1,仍取h=

0.1,计算到、5=

0.5,并与改进Euler法、梯形法在=.5处比较其误差大小.解用四阶R-K方法公式

3.

2.12,此处/x，y=-y+x+l,x0=0,九=1收=

0.1,于是当口加时=两Jo=_%+而+1=0,“1,h..%=/两+5瓦y+5/=-1+9戊+1-如+1=

0.05乙乙.”1,h、、玲=/5+5，为+5舄乙乙=7+人+勺》0+1-g+1=

0.0475t乙乙S=」凝+瓦Vo+5与乙以2人3以2人3=7+4-+»0+1-4+3-7两+1=

0.09525乙乙于是必=0+缄月+2抬+2后+%=1+^x

0.29025=

1.00483750，按公式

663.

2.12可算出为=1,01873090,^=1,04081842%=

1.07032029,^=1,10653091/一双々）1=

2.8x1此方法误差:|y-yx|=

5.5xio-4改进Euler法误差:55梯形法误差|y-y（x）|=

2.5x10-55可见四阶R-K方法的精度比二阶方法高得多.用四阶R-K方法求解初值问题（

1.

2.1）精度较高，但要从理论上给出误差网-y（瑞的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长h通常是通过不同步长在计算机上的计算结果近似估计.设双A）在演处的值口=双居），当演+1=醺时，y（Xn+l）的近似为乂T，于是由四阶R-K方法有h若以5为步长，计算两步到醺小则有乙1+1=M/+1）一月处2分©5乙

3.

2.13于是得1它给出了误差的近似估计,如果行|y焉-八|«£（£为给定精度），则认为以《为步长的计算结果y中满足精度要求，若』力;一1|廉则还可放大步2八十115长.因此⑶

2.13）提供了自动选择步长的方法.附经典的4级也是4阶R-K方法的程序function y=DELGKT4_Rungkutaf,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;y1=y0;x=a:h:b;var=findsymf;KI=Funval f,varvec,[xi-l yi-1];Funval f,varvec,K2=[xi-l+h/2yi-l+Kl*h/21;Funval f,varvec,K3=[xi-l+h/2yi-l+K2*h/2];Funval f,varvec,K4=[x i-l+h yi-l+h*K3];=yi-l+h*Kl+2*K2+2*K3+K4/6;yifor i=2:N+lendformat short;讲解求初值问题（

1.

2.1）的单步法主要是指Runge-Kutta法，本节主要讨论显式R—K方法，建立具体的计算公式使用的是Taylor展开，形如（

3.

2.4）的显式R—K方法，当r=l时就是Euler法，因此只要讨论‘r之2的计算公式，在「确定后如何推导公式都是一样的，只是r越大计算越复杂，为了掌握了解公式来源，只要以r=2为例推导计算公式即可.因此本节重点就是用Taylor展开求出r=2的显式R-K方法的计算公式，由于方法的局部截断误差为（

3.

2.6）,4+i的右端有〃演+町瓦入+砥枕）的项，要对它做Taylor展开，就要用到二元函数/（xj）的Taylor展开，按照二元函数Taylor级数=锯（吗+呜）””）+叱）//\r/XJ A/xj」

1、2A「/A的嗤3安+・・・

3.

2.14将它用到（

3.

2.6）的Q1的展开式中，即可得到按治升幕整理出的结果，对r=2的公式只能得到P=2阶的公式，即＜+1=°（加），于是2级R-K方法（

3.

2.5）的系数巧色，”如必须满足（

3.

2.7）给出的方程,它的解由（

3.

2.8）给出,只要右w0,求出的公式都是r=2的2阶R-K方法.而常用的就是白改进Euler法（

3.

1.11）和J=1得到的中点公式（

3.

2.9）.我们知道，显式单步法的数值公式为yi=准+卜我％八漓门=°\,・，n+如果取定右边的r沙区,九上）内JI则称该数值公式为显式Runge-Kutta方法.

（一）显式Runge-Kutta方法

1.二阶二级R-K方法%+]=/+力左2）其中4=/区，笫），k=f（x+a/,yb2M（4,y口.四参数满足nG=1,G=0，4==—举例中点公式四个参数.一2,即h h券+i=%+比七+不，先+5/（%，笫）乙乙程序如下function y=DELGKT2_mid f,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;yl=y0;x=a:h:b;var=findsymf;for i=2:N+1vl=Funval f,varvec,[x i-l yi-l];t=yi-1+h*vl/2;v2=Funvalf,varvec,[xi-l+h/2t];yi=yiT+h*v2;endformat short;休恩suen法cl=l/4,c2=3/4,a2=b21=2/3function y=DELGKT2_suenf,h,a,b,yO,varvec formatlong;N=uintl6b~a/h;y=zeros N+l,1;y1=yO;x=a:h:b;var=findsymf;for i=2:N+1KI=Funvalf,varvec,[xi-l yi-l];K2二Funval f,varvec,[xi-l+2*h/3yiT+Kl*2*h/3];yi=yi-l+h*K1+3*K2/4;end formatshort;改进的Euler法cl=l/2,c2=l/2,a2=b21=lfunction y=DEModifEuler f,h,a,b,yO,varvec formatlong;N=b-a/h;y=zeros N+l,1;y1=yO;x=a:h:b;var=findsymf;for i=2:N+lvl=Funval f,varvec,[xi-l yi-l];t=yi-l+h*vl;v2=Funvalf,varvec,[xi t];yi=yi-l+h*vl+v2/2;end formatshort;实验对比参数取法不同得到不同的二阶R-K方法精确程度各异这个已经可以做（同学们已做？！）注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数,要达到p=3则在（

3.

2.6）的展开式中要增加3项，即增加三个方程,加上（

3.

2.7）的三个方程，共计六个方程求4个待定参数，验证得出是无解的.

2.三阶三级R-K方法^（x,y,h）=c/i+q后+q网nn其中K=/（Z，%）,%2=/（Z+%/〃+（Z，K））,而心=/（/+电凡乂也占+4心小）。