(→lorenz.m) |
|||
(14 versioni intermedie di uno stesso utente non sono mostrate) | |||
Riga 1: | Riga 1: | ||
+ | == L'attrattore di Lorenz == | ||
+ | [[Immagine:Lorenz.png|right|thumb|400px|''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 | %Equazione differenziale Ordinaria dell'attrattore di lorenz | ||
%----------------------------------------------------------- | %----------------------------------------------------------- | ||
− | % dx/dt = -10 x | + | % dx/dt = -10 x + 10 y |
% dy/dt = x - y - (sqrt(27)+x)z | % dy/dt = x - y - (sqrt(27)+x)z | ||
% dz/dt = sqrt(27)*(x+y) - z + xy | % dz/dt = sqrt(27)*(x+y) - z + xy | ||
− | |||
function dy = lorenz(t,x) | function dy = lorenz(t,x) | ||
dy = zeros(3,1); | dy = zeros(3,1); | ||
Riga 11: | Riga 17: | ||
dy(3) = sqrt(27)*(x(1)+x(2))-x(3)+x(1)*x(2); | dy(3) = sqrt(27)*(x(1)+x(2))-x(3)+x(1)*x(2); | ||
end | 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 == | ||
+ | |||
+ | [http://mathworld.wolfram.com/LorenzAttractor.html Mathword: Lorenz Attractor] |
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).
%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
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
%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