rungekutta_4.m

以lorenz吸引子为例说明四阶定步长Runge-Kutta算法 详情见程序说明

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x=lorenz_series()
odefile='lorenz';
tspan=[0 100 0.01];
y=[0; 1; 0];
[tout,yout]=RungeKutta_4(odefile,tspan,y);
x=yout(:,1);
x=x';
x=x(5000:8000);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tout,yout]=RungeKutta_4(odefile,tspan,y)
%-------------------------------------------------------
%四阶定步长Runge-Kutta算法:
%  k1=h*f(tk,xk)
%  k2=h*f(tk+h/2,xk+k1/2)
%  k3=h*f(tk+h/2,xk+k2/2)
%  k4=h*f(tk+h,xk+k3)          其中h计算步长
%下一个步长的状态变量值为：
%   x(k+1)=xk +(k1+2*k2+2*k3+k4)/6
%-------------------------------------------------------
t0=tspan(1); th=tspan(2);
if length(tspan)<=3,h=tspan(3);
else, h=tspan(2)-tspan(1); th=tspan(2); end
tout=[t0:h:th]'; yout=[];
for t=tout'
    k1=h*eval([odefile '(t,y)']);
    k2=h*eval([odefile '(t+h/2,y+0.5*k1)']);
    k3=h*eval([odefile '(t+h/2,y+0.5*k2)']);
    k4=h*eval([odefile '(t+h,y+k3)']);
    y=y+(k1+2*k2+2*k3+k4)/6;
    yout=[yout;y'];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out1 = lorenz(tspan,y)
out1 = [-1*16*(y(1)-y(2)); -1*y(1)*y(3)+45.92*y(1)-y(2); y(1)*y(2)-4*y(3)];