% Integration par la methode de Runge-Kutta
% du systeme de Lorenz.
% On s'interesse ici aux trajectoires :
% trace de (x,y,z) en fonction de t, 
% trace dans des plans de phase

clear all
close all

% on fixe les parametres du systeme
r=35;
s=16;
b=4;

% on fixe le pas de temps, le temps jusqu'auquel 
% on calcule, la condition initiale
dt=0.01
Tmax=50;
x(1)=1;
y(1)=1;
z(1)=1;
tps(1)=0;
k=1;
% iteration
  for t=0:dt:Tmax,
        k=k+1;
        dx1=s*(y(k-1)-x(k-1));
        dy1=r*x(k-1)-y(k-1)-x(k-1)*z(k-1);
        dz1=x(k-1)*y(k-1)-b*z(k-1);
        
        x2=x(k-1)+dt*dx1/2;
        y2=y(k-1)+dt*dy1/2;
        z2=z(k-1)+dt*dz1/2;
        
        dx2=s*(y2-x2);
        dy2=r*x2-y2-x2*z2;
        dz2=x2*y2-b*z2;
        
        x3=x(k-1)+dt*dx2/2;
        y3=y(k-1)+dt*dy2/2;
        z3=z(k-1)+dt*dz2/2;
        
        dx3=s*(y3-x3);
        dy3=r*x3-y3-x3*z3;
        dz3=x3*y3-b*z3;
        
        x4=x(k-1)+dt*dx3;
        y4=y(k-1)+dt*dy3;
        z4=z(k-1)+dt*dz3;
        
        
        dx4=s*(y4-x4);
        dy4=r*x4-y4-x4*z4;
        dz4=x4*y4-b*z4;
        
        x(k)=x(k-1)+dt*(dx1/6+dx2/3+dx3/3+dx4/6);
        y(k)=y(k-1)+dt*(dy1/6+dy2/3+dy3/3+dy4/6);
        z(k)=z(k-1)+dt*(dz1/6+dz2/3+dz3/3+dz4/6);
        
        tps(k)=t;
        
  end

       
plot(tps,x,tps,y,tps,z)
legend('x','y','z')
disp('pressez une touche pour continuer')
pause 
clf;
plot3(x,y,z)
xlabel('x(t)')
ylabel('y(t)')
zlabel('z(t)') 
disp('pressez une touche pour continuer')
pause 
clf;
plot3(z,x,y)
xlabel('z(t)')
ylabel('x(t)')
zlabel('y(t)') 
disp('pressez une touche pour continuer')
pause 
clf;
plot(x,y)
xlabel('x(t)')
ylabel('y(t)')
disp('pressez une touche pour continuer')
pause 
clf;
plot(y,z)
xlabel('y(t)')
ylabel('z(t)')
disp('pressez une touche pour continuer')
pause 
clf;
plot(x,z)
xlabel('x(t)')
ylabel('z(t)')
disp('pressez une touche pour continuer')