Matlab中龙格-库塔(Runge-Kutta)方法原理及实现.doc

函数功能ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(x)。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解. www.iLoveM使用方法T,Y=ode45(odefun,tspan,y0) odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 www.iLoveM tspan 是区间 t0tf或者一系列散点t0,t1,.,tf book.iLoveM y0 是初始值向量 www.iLoveMT 返回列向量的时间点 www.iLoveM Y 返回对应T的求解列向量 www.iLoveMT,Y=ode45(odefun,tspan,y0,options) Simulink与信号处理 options 是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 Simulink与信号处理 T,Y,TE,YE,IE=ode45(odefun,tspan,y0,options) Matlab中文论坛在设置了事件参数后的对应输出 www.iLoveM TE 事件发生时间 book.iLoveM YE 事件解决时间 Matlab中文论坛IE 事件消失时间 www.iLoveMsol=ode45(odefun,t0tf,y0.) book.iLoveM sol 结构体输出结果 www.iLoveM应用举例1 求解一阶常微分方程程序:一阶常微分方程odefun=(t,y)(y+3*t)/t2;%定义函数 Simulink与信号处理tspan=14;%求解区间y0=-2;%初值 Simulink与信号处理 t,y=ode45(odefun,tspan,y0);plot(t,y)%作图title(t2y=y+3t,y(1)=-2,1t4)legend(t2y=y+3t) xlabel(t) Simulink与信号处理 ylabel(y)%精确解%dsolve(t2*Dy=y+3*t,y(1)=-2)%ans= 一阶求解结果图%(3*Ei(1)-2*exp(1)/exp(1/t)-(3*Ei(1/t)/exp(1/t) book.iLoveM2 求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y,y.y(n-1),t)=0用变量替换,y1=y,y2=y.注意odefun方程定义为列向量dxdy=y(1),y(2).Simulink与信号处理程序:functionTestode45tspan=3.94.0;%求解区间y0=28;%初值t,x=ode45(odefun,tspan,y0);plot(t,x(:,1),-o,t,x(:,2),-*)legend(y1,y2)title(y=-t*y+et*y+3sin2t) book.iLoveMxlabel(t)ylabel(y)functiony=odefun(t,x)y=zeros(2,1);%列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend高阶求解结果图相关函数ode23,ode45,ode113,ode15s,ode23s,ode23t,ode23tbMatlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的“理论基础”来源于“泰勒公式”和“使用斜率近似表达微分”。它在积分区间内多计算几个点的斜率，然后进行加权平均，用作下一点的依据，从而构造出了精度更高的数值积分计算方法。如预先求两个点的斜率就是二阶龙格库塔法，如预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y=f(x,y),使用差分概念。(Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn）Yn+1=Yn+h*f(Xn,Yn)另外根据微分中值定理，存在0t1,使得Yn+1=Yn+h*f(Xn+th,Y(Xn+th)这里Kf(Xn+th,Y(Xn+th)称为平均斜率，龙格库塔方法就是求得K的一种算法。利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h5)的四阶龙格库塔公式：K1f(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);所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数：function x,y=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor(b-a)/h);%求步数x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数for ii=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调用的子函数以及其调用语句：function dy=test_fun(x,y)dy = zeros(3,1);%初始化列向量dy(1) = y(2) * y(3);dy(2) = -y(1) + y(3);dy(3) = -0.51 * y(1) * y(2);对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下：T,F = ode45(test_fun,0 15,1 1 3);subplot(121)plot(T,F)%Matlab自带的