clear all;
global U V k1 k2 k3 k4 k21 k22 k23 k24 a t0 m1 m2 m3 m4 mu d0 d2 m;
t0=0.1;U=1.0;m=0.0;n = 2^4;k1 = 1.1;k2 =0.1;k3 = 0.1;k4 =0.1;k21=0.0;k22=0.0;k23=0.0;k24=0.0;
a=0.01;N0=2^4;mu=0.0;
for l =1:5
U=0.2.*(l-1)
for N0=1:2^4
clear y u w1 w wt w2 wt2 w12 sum1 w3 b1 b2 b3 sum1;
for k = 1:n
r2=0.001.*rand(1);
r3=0.001.*exp(i.*2.*pi.*rand(1));
r4=0.001.*exp(i.*2.*pi.*rand(1));
d=1./3;
%g=[d,d,d,d,0.0+r2,0.0+r2,0.0+r2,0.0+r2];
%[t,u]=ode45('loreq30',[N0],g');
g=[d,d,d,0.0+r3,0.0+r3,0.0+r3]';
[t,u]=ode45('loreq49',[N0],g);
xi = 0.0:0.1:N0;
y(:,1,k) = interp1(t,u(:,1),xi);
y(:,2,k) = interp1(t,u(:,2),xi);
y(:,3,k) = interp1(t,u(:,3),xi);
y(:,4,k) = interp1(t,u(:,4),xi);
y(:,5,k) = interp1(t,u(:,5),xi);
y(:,6,k) = interp1(t,u(:,6),xi);
%y(:,7,k) = interp1(t,u(:,7),xi);
%y(:,8,k) = interp1(t,u(:,8),xi);
end
N=N0./0.1;
for k=1:n
w(:,k)=fft(y(:,1,k)+y(:,2,k)+y(:,3,k));
wt(:,k)=conj(fft(y(:,1,k)+y(:,2,k)+y(:,3,k)));
w2(:,k)=fft(y(:,4,k)+y(:,5,k)+y(:,6,k));
wt2(:,k)=conj(fft(y(:,4,k)+y(:,5,k)+y(:,6,k)));
end
for i=1:N
for j=1:N
% clear w w1 wt wt2;
for k=1:n
w1(i,j,k)=w(i,k).*wt(j,k)./n;
w12(i,j,k)=w2(i,k).*wt2(j,k)./n;
end
end
end
for k=1:n
w3(:,k)=fft(y(:,1,k)+y(:,2,k)+y(:,3,k)).*conj(fft(y(:,1,k)+y(:,2,k)+y(:,3,k)));
end
for i=1:N
for j=1:N
b1(i,j)=sum(w1(i,j,:));
b2(i,j)=sum(w12(i,j,:));
b3(i)=sum(w3(i,:));
end
end
Entropy(N0,l)=trace(b1)./(N.^2);
Entropyph(N0,l)=trace(b2)./(N.^2);
sum1=abs(sum(b3));
rate(N0,l)=N./(sum1-pi);
end
end
x= 1:N0;
rate1=rate(:,1);
rate2=rate(:,2);
rate3=rate(:,3);
rate4=rate(:,4);
rate5=rate(:,5);
Entropy1=Entropy(:,1);
Entropy2=Entropy(:,2);
Entropy3=Entropy(:,3);
Entropy4=Entropy(:,4);
Entropy5=Entropy(:,5);
Entropyph1=Entropyph(:,1);
Entropyph2=Entropyph(:,2);
Entropyph3=Entropyph(:,3);
Entropyph4=Entropyph(:,4);
Entropyph5=Entropyph(:,5);
plot(rate1,log(2.*pi.*abs(Entropy1)),rate2,log(2.*pi.*abs(Entropy2)),rate3,log(2.*pi.*abs(Entropy3)),rate4,log(2.*pi.*abs(Entropy4)),rate5,log(2.*pi.*abs(Entropy5)))
%plot(x,rate(:,1),x,rate(:,2),x,rate(:,3),x,rate(:,4),x,rate(:,5))
%plot(x,log(2.*pi.*abs(Entropyph(:,1)))./2,x,log(2.*pi.*abs(Entropyph(:,2)))./2,x,log(2.*pi.*abs(Entropyph(:,3)))./2,x,log(2.*pi.*abs(Entropyph(:,4)))./2,x,log(2.*pi.*abs(Entropyph(:,5)))./2)
%plot(real(C))
