sandbox/easystab/lorenz.m
Numerical integration of the Lorenz system
The Lorenz system is a classical dimension-3 dynamical system modelling thermal convection in a rectangular cell.
The system is defined as follows :
\displaystyle \frac{ dx}{dt} = P (y-x) \displaystyle \frac{ dy}{dt} = r x - y - x z; \displaystyle \frac{ dz}{dt} = x y - bz.
r is the control parameter and represents the temperature difference between upper and lower surfaces (Rayleigh number)
P reprents the Prandtl number (classical value 10)
b is a geometrical parameter related to the aspect ratio of the cell (classical value 8/3)
This program will draw two trajectories with slightly different initial conditions.
What you can modify when playing with the program is :
control parameter r (line 47)
initial conditions X0 and Xb0 : change definition in lines 55-56 or uncomment lines 61-63
function [] = lorenz()
Definition of the parameters and initialization of the figures
% control parameter
r= 24;
% Physical parameters
P = 10;
b = 8/3;
%% Initial conditions
%% X0 is initial condition of "blue" trajectory and Xb is initial condition of "red" trajectory
% First idea : two points close to [1,1,1]
X0 = [1;1;1];
Xb0 = [1+1e-7;1;1] ;
%% Second idea : two points close to the equilibrium point
%% the system has an equilibrium point [Xf,Xf,Xf^2/b] ; with Xf = sqrt((r-1)*b) ;
%% we will use initial condition close to this point
% Xf = sqrt((r-1)*b);
% X0 = [Xf+1e-2;Xf;Xf^2/b+1e-3]
% Xb0 = [Xf;Xf+1e-2;Xf^2/b-1e-3]
% Numerical parameters
T0 = 0; % initial time
Tstep = 1; % time-step
Nt = 100; % max-number of integrations
Tpause = .3; %time (in seconds between two plots)
% Initialisation of the figures
close all;
figure(1);
xlabel('x'), ylabel('y'), zlabel('z');
title('Trajectory in phase-space for two initial conditions');
grid;
view(20,10); % set angle for 3D view
figure(2);
subplot(2,1,1);
xlabel('t');ylabel('x (-) , y(--) , z (..)');
title('Time-history of the solutions');
subplot(2,1,2);
title('Distance between the two trajectories');
xlabel('t');ylabel('|X(blue)-X(red)|');
disp('Position the figures and press enter to launch');
pause;
Performing the temporal integration
for N = 1:Nt
Dt = [T0 : Tstep/100 : T0+Tstep] ; % bornes d'integration
[t, X] = ode45(@(t,X)(f_lorenz(X,r,P,b)), Dt, X0); % integration for "blue"
[tb, Xb] = ode45(@(t,X)(f_lorenz(X,r,P,b)), Dt, Xb0); % integration for "red"
figure(1); hold on;%
plot3(X(:,1), X(:,2), X( :,3),'b-');
plot3(Xb(:,1), Xb(:,2), Xb(:,3),'r-');
X0 = X(end,:); Xb0 = Xb(end,:); T0 = t(end);% initial condition for next run
figure(2);
subplot(2,1,1);hold on;
plot(t,X(:,1),'b-',t,X(:,2),'b--',t,X(:,3),'b:', tb,Xb(:,1),'r-',tb,Xb(:,2),'r--',tb,Xb(:,3),'r:');
Xbb = interp1(tb,Xb,t);
diff = sqrt( (Xbb(:,1)-X(:,1)).^2 + (Xbb(:,2)-X(:,2)).^2 + (Xbb(:,3)-X(:,3)).^2);
subplot(2,1,2);hold on;
semilogy(t,diff,'k');
pause(Tpause);
end
Results
figure(1);
set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','Lorenz1.png');
figure(2);
set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','Lorenz2.png');
end
Definition of the function governing the dynamical system
function f = f_lorenz(X,r,P,b)
% cette fonction definit le systeme de Lorenz
x = X(1); y = X(2); z = X(3);
f = zeros(3,1); % vecteur colonne
f(1) = P*(y-x);
f(2) = r*x - y - x*z;
f(3) = x*y - b*z;
end
