Contents

L'attrattore di Lorenz

(thumbnail)
L'attrattore di Lorenz

Utilizzando i due file per Matlab (testati con la versione 7.0) seguenti è possibile visualizzare una animazione (per t < tf) della posizione nel tempo di N soluzioni all'equazione differenziale di Lorenz. In praticolare le condizioni iniziali a partire dalle quali vengono calcolate le N soluzioni sono contenute in una sfera di raggio 0.1.

Dall'animazione è possibile notare come con l'avanzare del tempo le N soluzioni divergano anche a fronte di una piccolissima variazione delle condizioni iniziali, ma tuttavia mantenendosi confinate sull'attrattore di Lorenz (disegnato in celeste).

lorenz.m

 %Equazione differenziale Ordinaria dell'attrattore di lorenz
 %-----------------------------------------------------------
 % dx/dt = -10 x + 10 y
 % dy/dt = x - y - (sqrt(27)+x)z
 % dz/dt = sqrt(27)*(x+y) - z + xy
 function dy = lorenz(t,x)
    dy = zeros(3,1);
    dy(1) = -10*x(1) + 10 * x(2);
    dy(2) = x(1) - x(2) -(sqrt(27)+x(1))*x(3);
    dy(3) = sqrt(27)*(x(1)+x(2))-x(3)+x(1)*x(2);
 end

draw_lorenz.m

 function [t,x,y,out] = solve_lorenz(x_0,t_f,n)
   y = [ ] ;
   x = [ ] ;
   plot3(0,0,0);
   out = cell(n,1);
   hold;
   for i=1:n
       x_i = rand(1,3)*0.1 + x_0;
       x = [x x_i'];
       [t,temp]=ode45(@lorenz,[0 t_f], x_i);
       out{i} = temp;
       y_len = length(temp);
       y_f = temp(y_len,:);
       y = [ y y_f'];
       plot3(x_i(1),x_i(2),x_i(3),'ko');
       plot3(y_f(1),y_f(2),y_f(3),'kx');
   end
   min_size = 1e10;
   min_x = 1e10 ;
   min_y = 1e10 ;
   min_z = 1e10 ;
   max_x = 0;
   max_y = 0;
   max_z = 0;
   for i=1:n
     s = size(out{i});
     x_max = max(out{i}(:,1));
     x_min = min(out{i}(:,1));
     y_max = max(out{i}(:,2));
     y_min = min(out{i}(:,2));
     z_max = max(out{i}(:,3));
     z_min = min(out{i}(:,3));
     if s(1) < min_size
        min_size = s(1);
     end
     if x_max > max_x
        max_x = x_max;
     end
     if y_max > max_y
        max_y = y_max;
     end
     if z_max > max_z
        max_z = z_max;
     end
     if x_min < min_x
        min_x = x_min;
     end
     if y_min < min_y
        min_y = y_min;
     end
     if z_min < min_z
        min_z = z_min;
     end
   end
   %Plotta per ogni istante la posizione delle n soluzioni
   for j = 1:min_size
     clf;
     axis([min_x max_x min_y max_y min_z max_z]);
     plot3(out{n}(:,1),out{n}(:,2),out{n}(:,3),'c');
     hold on;
     for i=1:n
       v = out{i}(j,:);
       plot3(v(1),v(2),v(3),'ko');
     end
     pause(0.1);
     drawnow;
     hold off;
   end
 end

Riferimenti

Mathword: Lorenz Attractor