جاذب لورنز lorenz attractor

در این بخش برنامه جاذب لورنز را با استفاده از دو روش رانگ کوتای مرتبه4 و ode45 با توجه به توضیحات داده شده مینویسم.

جاذب لورنز یک سیستم آشوبناک است که برای نختسن بار توسط ادوارد لورنز از انستیتو تکنولوژی ماساچوست معرفی شد.

http://mathworld.wolfram.com/LorenzAttractor.html

http://en.wikipedia.org/wiki/Lorenz_attractor

http://ozviz.wasp.uwa.edu.au/~pbourke/fractals/lorenz

سیستم معادلات لورنز عبارتند از

نمی خواهم در مورد نحوه استخراج این معادلات بحث کنم و تنها به این بسنده می­کنم که اساس آنها مباحث فیزیکی در ارتباط با جابجایی اتمسفراست. متغیرهای x، y و z نشان دهنده کمیت ها فیزیکی همچون دما و سرعت سیال هستند و پارامترها معرفی شده مربوط به ویژگیهای جو می باشند که با داده های زیر معین می شوند،

در ابتدا این معادلات را با استفاده از روش رانگ کوتا حل می کنیم. تابع Lorenz_rk اینکار را انجام می دهد. حاصل کار بصورت یک انیمیشن نمایش داده می­شود.

function lorenz_rk

%-----------------

N=5000;

sigma = 10;

rho = 25;

beta = 8/3.;

t=linspace(0,50,N); % time between 0-5s

h=t(2)-t(1);   % time step

f=@(x,y) sigma*(y-x);

g=@(x,y,z) x*(rho-z)-y;

u=@(x,y,z) x*y-beta*z;

x(1)=1;

y(1)=0;

z(1)=0;

for i=1:length(t)-1

kx1=f(x(i),y(i));

ky1=g(x(i),y(i),z(i));

kz1=u(x(i),y(i),z(i));

%------------------------

kx2=f(x(i)+h/2*kx1,y(i)+h/2*ky1);

ky2=g(x(i)+h/2*kx1,y(i)+h/2*ky1,z(i)+h/2*kz1);

kz2=u(x(i)+h/2*kx1,y(i)+h/2*ky1,z(i)+h/2*kz1);

%------------------------

kx3=f(x(i)+h/2*kx2,y(i)+h/2*ky2);

ky3=g(x(i)+h/2*kx2,y(i)+h/2*ky2,z(i)+h/2*kz2);

kz3=u(x(i)+h/2*kx2,y(i)+h/2*ky2,z(i)+h/2*kz2);

%------------------------

kx4=f(x(i)+h*kx3,y(i)+h*ky3);

ky4=g(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);

kz4=u(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);

%------------------------

x(i+1)=x(i)+h/6*(kx1+2*kx2+2*kx3+kx4);

y(i+1)=y(i)+h/6*(ky1+2*ky2+2*ky3+ky4);

z(i+1)=z(i)+h/6*(kz1+2*kz2+2*kz3+kz4);

end

body = line( ...

'color','b', ...

'LineStyle','none', ...

'Marker','.', ...

'erase','none', ...

'xdata',[],'ydata',[],'zdata',[]);

axis([ min(y)-eps max(y)+eps...

min(z)-eps max(z)+eps])

for j=2:N

set(body,'XData',y(j),'YData',z(j))

drawnow

end

plot(y,z,'.')

xlabel('y')

ylabel('z')

 

 

 

با استفاده از تابع ode45 این سیستم را میتوانیم بشکل زیر حل کنیم که با نام Lorenz_ode مشخص شده است.

function Lorenz_ode

close all

x0=1;

y0=0;

z0=0;

N=10000;

[t,Y] = ode45(@lorenz_at,linspace(0,50,N),[x0 y0 z0]);

body = line( ...

'color','b', ...

'LineStyle','none', ...

'Marker','.', ...

'erase','none', ...

'xdata',[],'ydata',[],'zdata',[]);

axis([ min(Y(:,2))-eps max(Y(:,2))+eps...

min(Y(:,3))-eps max(Y(:,3))+eps])

for j=2:N

set(body,'XData',Y(j,2),'YData',Y(j,3))

drawnow

end

plot(Y(:,2),Y(:,3),'.')

xlabel('y')

ylabel('z')

function dy=lorenz_at(t,y)

dy=zeros(3,1);

SIGMA = 10;

RHO = 25;

BETA = 8/3.;

dy(1)=SIGMA*(y(2)-y(1));

dy(2)=y(1)*(RHO-y(3))-y(2);

dy(3)=y(1)*y(2)-BETA*y(3);

نوشتن دیدگاه


تصویر امنیتی
تصویر امنیتی جدید