Listingul MATLAB al modelului pentru rezolvarea sistemului de ecuatii Lang

1- Rutina de calcul a arborelui cu bifurcatii pentru modelul microeconomic (MATLAB)

clear

Stabilizare=100;

N=1000;

a=0.2;b=0.8;c=100;alfa=0.46;beta=0.8; %1 > alfa, beta > 0

x_=0.2;

y_=0.12;

cc=0:1:110;

axis([cc(1) cc(size(cc,2)) 0 0.3]);title('','FontSize',18);hold on;

for k=1:size(cc,2)

c=cc(k);

for n=1:Stabilizare

x_=(1-alfa)*x_+a/(1+exp(-c*(x_-y_)));

y_=(1-beta)*y_+b/(1+exp(-c*(x_-y_)));

end

x(1)=x_;y(1)=y_;

for n=1:N

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

end

media(k)=mean(x(1:N));

imprastierea(k)=std(x(1:N));

plot(c,x,'ro','LineWidth',0.5,'MarkerEdgeColor','k','MarkerFaceColor','r','MarkerSize',3)

end

figure;plot(1:1:size(media,2),media)

figure;plot(1:1:size(imprastierea,2),imprastierea)

2-Rutina de calcul seriei temporale pentru modelul microeconomic (MATLAB)

clear

N=1000;

a=0.16;b=0.9;c=105;alfa=0.46;beta=0.7; %1 > alfa, beta > 0

m=1; K=0; %fara control K=-0.1...

x(1)=0.02;

y(1)=0.12;

for n=1:m

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

end

for n=m+1:N

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

%control

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_1_N1=mean(x(1:N))

plot(1:1:size(x,2),x,'*')

%ATENTIE perioadele orbitelor "O" depind foarte mult de m si k (extrem de mult!) de ex k=0.1 si m=1 are O=4

3-Rutina de calcul a punctelor fixe pentru modelul microeconomic (MATLAB)

%rezolvarea unui sistem de ecuatii neliniare pentru punctele fixe modelul BehrensFeichtinger

clear

options=optimset('Display','iter'); % Option to display output

x0 =[0; 0]; % prima aproximare

options.TolFun=1e-12;options.TolX=1e-7;

x0 = [0.01182215476882; 0.04369975066333]; %dupa prima ghicire

options.TolFun=1e-20; options.TolX=1e-20; %scadem mult si tolerantele!

[x,fval] = fsolve(@functieBehrensFeichtinger,x0,options) % Call optimizer

%verificare

N=100;

a=0.16;b=0.9;c=105;alfa=0.46;beta=0.7; %1 > alfa, beta > 0

X(1)=x(1);

Y(1)=x(2);

for n=1:N

X(n+1)=(1-alfa)*X(n)+a/(1+exp(-c*(X(n)-Y(n))));

Y(n+1)=(1-beta)*Y(n)+b/(1+exp(-c*(X(n)-Y(n))));

end

plot(1:1:size(X,2),X,'*',1:1:size(Y,2),Y,'.')

4-Rutina de calcul a coef. Lyapunov, model map 2D (MATLAB)

clear

a=0.16;b=0.9;c=105;alfa=0.46;beta=0.7; %1 > alfa, beta > 0

x=0.02;y=0.12;y_p=0; %valorile initiale

N=10^4;

NN=50;dx=1/NN;dy=1/NN;

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

for q=1:1:NN;

xx(q)=q*dx;

for qq=1:1:NN;

yy(qq)=qq*dy;

x=xx(q);y=yy(qq);

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

lam=0;det=0;

for n=2:N

y0=y;x0=x; %vechile puncte

%noile puncte, functie de vechile puncte

fdx = 1-alfa+a/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y)); %calculul derivatelor cu punctele vechi pentru a calcula noul y_p

fdy = -a/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

gdx = b/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

gdy = 1-beta-b/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

y_p=(gdx + gdy*y_p)/(fdx + fdy*y_p); %merege si in functie de x si y deiarece pana aici x0==x si y0==y

x=(1-alfa)*x0+a/(1+exp(-c*(x0-y0)));

y=(1-beta)*y0+b/(1+exp(-c*(x0-y0)));

fdx = 1-alfa+a/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y)); %calculul derivatelor cu punctele noi pentru a calcula noul lam

fdy = -a/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

gdx = b/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

gdy = 1-beta-b/(1+exp(-c*(x-y)))^2*c*exp(-c*(x-y));

lam=lam+log(((fdx + fdy*y_p)^2+(gdx + gdy*y_p)^2)/(1+y_p^2));

det=det+fdx*gdy-fdy*gdx; %determinantul Jacobianului...

end

lambda1(q,qq)=lam/2/N; %cel mai mare coeficient Lyapunov

det=det/2/N;

lambda2(q,qq)=log(abs(det))-lambda1(q,qq); %al doilea coficint Lyapunov

end

[X,Y]=meshgrid(xx,yy);

surf(X,Y,lambda1)

5-Rutina de control a manifestării haotice prin controlul investitiilor, model (MATLAB)

clear

N1=100;N2=200;N3=400;N4=600;N5=700;N6=900;

a=0.16;b=0.9;c=105;alfa=0.46;beta=0.7; %alfa>0 si beta<1

m=1;

x(1)=0.01;

y(1)=0.12;

for n=1:m+1

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

end

K=0;m=1; %nu mai conteaza m

for n=m+1:N1

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_1_N1=mean(x(1:N1))

K=0.5;m=1;

for n=N1+1:N2

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_N1_N2=mean(x(N1:N2))

K=0;m=1; %nu mai conteaza m

for n=N2+1:N3

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_N2_N3=mean(x(N2:N3))

K=-0.1;m=2;

for n=N3+1:N4

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_N3_N4=mean(x(N3:N4))

K=0;m=2; %nu mai conteaza m

for n=N4+1:N5

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

m_N4_N5=mean(x(N4:N5))

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

K=-0.5;m=1;

for n=N5+1:N5

x(n+1)=(1-alfa)*x(n)+a/(1+exp(-c*(x(n)-y(n))));

y(n+1)=(1-beta)*y(n)+b/(1+exp(-c*(x(n)-y(n))));

KK(n+1)=K;

F(n+1)=KK(n+1)*(y(n)-y(n-m));

y(n+1)=y(n+1)+F(n+1);

mm(n+1)=m;

end

%m_N5_N6=mean(x(N5:N6))

figure:plot(x(3:N5),'--','LineStyle','none','Marker','x');title('X')

figure:plot(F(3:N5),'--','LineStyle','none','Marker','x');title('X')

figure:plot(KK(3:N5),'--','LineStyle','none','Marker','x');title('X')

figure:plot(mm(3:N5),'--','LineStyle','n