www.pudn.com > zishiyingjiami.rar > RKK.m, change:2011-12-27,size:1615b


% Runge_Kutta_4 
% 方程表达式  
% dx/dt = a*(y-x)  
% dy/dt = (c-a)*x - x*z + c*y  
% dz/dt = x*y - b*z  
function Z=RKK(x_0,y_0,z_0) 
%参数值  
a = 35;  
b = 3;  
c = 28;  
 
h = 0.01;             % 积分时间步长  
step1 = 10000;        % 前面的迭代点数  
step2 = 5000;         % 后面的迭代点数  
  
X = [];  
Y = [];  
Z = [];  
  
for(i = 1:1:(step1 + step2))  
    %e1  
    x_e1 = -a*x_0 + a*y_0;  
    y_e1 = c*y_0 + (c - a)*x_0 - x_0*z_0;  
    z_e1 = -b*z_0 + x_0*y_0;  
    %e2  
    y_h = y_0 + 0.5*h*y_e1;  
    x_e2 = -a*(x_0 + 0.5*h*x_e1) + a*y_h;  
      
    x_h = x_0 + 0.5*h*x_e1;  
    z_h = z_0 + 0.5*h*z_e1;  
    y_e2 = c*(y_0 + 0.5*h*y_e1) + (c - a)*x_h - x_h*z_h;  
      
    z_e2 = -b*(z_0 + 0.5*h*z_e1) + x_h*y_h;  
    %e3  
    y_h = y_0 + 0.5*h*y_e2;  
    x_e3 = -a*(x_0 + 0.5*h*x_e2) + a*y_h;  
      
    x_h = x_0 + 0.5*h*x_e2;  
    z_h = z_0 + 0.5*h*z_e2;  
    y_e3 = c*(y_0 + 0.5*h*y_e2) + (c - a)*x_h - x_h*z_h;  
      
    z_e3 = -b*(z_0 + 0.5*h*z_e2) + x_h*y_h;  
    %e4  
    y_h = y_0 + h*y_e3;  
    x_e4 = -a*(x_0 + h*x_e3) + a*y_h;  
      
    x_h = x_0 + h*x_e3;  
    z_h = z_0 + h*z_e3;  
    y_e4 = c*(y_0 + h*y_e2) + (c - a)*x_h - x_h*z_h;  
      
    z_e4 = -b*(z_0 + h*z_e2) + x_h*y_h;  
    %叠代  
    x_1 = x_0 + 1/6*h*(x_e1 + 2*x_e2 +2*x_e3 + x_e4);  
    y_1 = y_0 + 1/6*h*(y_e1 + 2*y_e2 +2*y_e3 + y_e4);  
    z_1 = z_0 + 1/6*h*(z_e1 + 2*z_e2 +2*z_e3 + z_e4);  
      
    X = [X,x_1];  
    Y = [Y,y_1];  
    Z = [Z,z_1];  
      
    x_0 = x_1;  
    y_0 = y_1;  
    z_0 = z_1;  
end  
X = X(step1+1:end);  
Y = Y(step1+1:end);  
Z = Z(step1+1:end);