阶RungeKutta法求解一阶常微分方程.docx

《MATLAB语言及应用》
大作业

姓名：学号：学院：班级：题目编号：2013年10月13

4阶Runge-Kutta法求解一阶常微分方程。
Runge-Kutta法的数学理论
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。(Yn+1-Yn)/h=f(Xn,Yn)推出（近似等于，极限为Yn'）Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理，存在0<t<1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th))这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式：K1＝f(Xn,Yn);K2=f(Xn+h/2,Yn+(h/2)*K1);K3=f(Xn+h/2,Yn+(h/2)*K2);K4=f(Xn+h,Yn+h*K3);Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6)
Runge-Kutta的算法和流程图
在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为0(h5),被广泛应用于解微分方程的初值问题。其算法公式为:

流程图：
(1)、四阶龙格-库塔方法流程图:

（2）、实例求解流程图：

Runge-Kutta的Matlab实现
function[x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数forii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解end
Runge-Kutta的算例实现
例:求解常微分方程d(y)/d(x)=-2*y+2*x*x+2*x，0≤x≤0.5，y(0)=1.
编辑m文件Untitled3：
fun=inline('-2*y+2*x*x+2*x');
[x,y]=ode45(fun,[0,0.5],1)
>>Untitled3

x=

0
0.0125
0.0250
0.0375
0.0500
0.0625
0.0750
0.0875
0.1000
0.1125
0.1250
0.1375
0.1500
0.1625
0.1750
0.1875
0.2000
0.2125
0.2250
0.2375
0.2500
0.2625
0.2750
0.2875
0.3000
0.3125
0.3250
0.3375
0.3500
0.3625
0.3750
0.3875
0.4000
0.4125
0.4250
0.4375
0.4500
0.4625
0.4750
0.4875
0.5000

y=

1.0000
0.9755
0.9519
0.9291
0.9073
0.8864
0.8663
0.8471
0.8287
0.8112
0.7944
0.7785
0.7633
0.7489
0.7353
0.7224
0.7103
0.6989
0.6883
0.6783
0.6690
0.6605
0.6526
0.6454
0.6388
0.6329
0.6277
0.6231
0.6191
0.6157
0.6130
0.6109
0.6093
0.6084
0.6080
0.6083
0.6091
0.6104
0.6124
0.6148
0.6179