Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@nolte
nolte authored Jan 29, 2021
1 parent aa6d589 commit 57a1f9f
Show file tree
Hide file tree
Showing 23 changed files with 2,194 additions and 0 deletions.
162 changes: 162 additions & 0 deletions DoublePend.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
function DoublePend
format compact
c(1) = 'r';
c(2) = 'b';
c(3) = 'g';
c(4) = 'k';
c(5) = 'c';
c(6) = 'm';
c(7) = 'y';

h = colormap(lines);

E = 0.8; % 1.0 1.2; % Kinetic Energy

h2 = figure(2);
h2.Position = [11 -31 560 420];
dum = set(h2);

h3 = figure(3);
h3.Position = [10 381 560 420];
dum = set(h3);

h5 = figure(5);
h5.Position = [229 1 1008 804];
dum = set(h5);
clf
hold on


repnum = 128; % 32
mulnum = 64/repnum;
for reploop = 1:repnum

clear xvar yvar vy vx px py

vx0 = 2*(0.5 - rand)*sqrt(E);
vy0 = -vx0 + sign(randn)*sqrt(2*E - vx0.^2);

v = (vx0 + vy0);
KE = 0.5*(vx0.^2 + v.^2);

x0 = 0.0; % 0.1
y0 = -0.0; % -0.2

y0 = [x0 y0 vx0 vy0];

tspan = [1 3000]; % 500

options = odeset('RelTol',1e-9);

[t,y] = ode45(@f1,tspan,y0,options);


siz = size(t);
% y1 = mod(y(1:siz,1) + pi,2*pi)-pi; % angle-1
% y2 = mod(y(1:siz,2) + pi,2*pi)-pi; % angle-2
% y3 = mod(y(1:siz,3) + pi,2*pi)-pi; % angular speed-1
% y4 = mod(y(1:siz,4) + pi,2*pi)-pi; % antular speed-2
y1 = y(1:siz,1); % angle-1
y2 = y(1:siz,2); % angle-2
y3 = y(1:siz,3);
y4 = y(1:siz,4);

figure(1)
plot(t(1:siz),y1,t(1:siz),y2)

figure(2)
plot(t(1:siz),y3,t(1:siz),y4)
xlabel('Time')
legend('speed')

figure(3)
plot(y1,y2)
xlabel('x')
ylabel('y')

KinEnergy = 0.5*y3.^2 + 0.5*(y3.^2 + y4.^2 + 2*y3.*y4.*cos(y2-y1));
PotEnergy = - 2*cos(y1) - cos(y2);
Energy = KinEnergy + PotEnergy;
figure(4)
plot(t(1:siz),KinEnergy,t(1:siz),PotEnergy,t(1:siz),Energy)
title('Energy')
legend('Kin','Pot','E')

% Power Spectrum
Pow = 0;
if Pow == 1
yf = y(:,1);
F = fft(yf);
Pow = F.*conj(F);

figure(4)
plot(Pow)
end

trig = y2; % Trigger for crossing Poincare section: try y1 and switch xvar/yvar and vx/vy

%First Return Map
Fst = 1;
if Fst == 1
cnt = 0;
last = trig(1);
for loop = 2:siz
if (last < 0)&&(trig(loop) > 0)
cnt = cnt+1;
del1 = -trig(loop-1)./(trig(loop) - trig(loop-1));

vx(cnt) = y3(loop-1) + del1*(y3(loop)-y3(loop-1));
vy(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1));
xvar(cnt) = y1(loop-1) + del1*(y1(loop)-y1(loop-1));
yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1));

px(cnt) = 2*vx(cnt) + vy(cnt)*cos(yvar(cnt)-xvar(cnt));
py(cnt) = vy(cnt) + vx(cnt)*cos(yvar(cnt)-xvar(cnt));

last = trig(loop);
else
last = trig(loop);
end
end

figure(5)
plot(xvar,px,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:))
xlabel('y')
ylabel('vy')
end

pause(0.001)

end % end reploop

figure(5)
hold off
print -dtiff DoublePendPrint


% Model function
function dy = f1(t,y)

A = y(4).^2.*sin(y(2)-y(1));
B = -2*sin(y(1));
C = y(3).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1));
D = sin(y(2)).*cos(y(2)-y(1));
EE = 2 - (cos(y(2)-y(1))).^2;

FF = y(4).^2.*sin(y(2)-y(1)).*cos(y(2)-y(1));
G = -2*sin(y(1)).*cos(y(2)-y(1));
H = 2*y(3).^2.*sin(y(2)-y(1));
I = 2*sin(y(2));
JJ = (cos(y(2)-y(1))).^2 - 2;

dy = zeros(4,1);
dy(1) = y(3);
dy(2) = y(4);
dy(3) = (A + B + C + D)./EE;
dy(4) = (FF + G + H + I)./JJ;

end % end f5


end % end ltest

161 changes: 161 additions & 0 deletions Heiles.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
function Heiles

c(1) = 'r';
c(2) = 'b';
c(3) = 'g';
c(4) = 'k';
c(5) = 'c';
c(6) = 'm';
c(7) = 'y';

h = colormap(lines);


E = 0.1; % 1 0.1

epsE = 1; % 0.3425 1

prms = sqrt(E);
pmax = sqrt(2*E);

% Potential Function
for xloop = 1:100
x = -2 + 4*xloop/100;
for yloop = 1:100
y = -2 + 4*yloop/100;

V(yloop,xloop) = 0.5*x^2 + 0.5*y^2 + epsE*(x^2*y - 0.33333*y^3);

end
end
figure(6)
colormap(gray)
imagesc(V)
hold on
contour(V,30,'LineColor','k')
hold off
caxis([-0.5 1])
colorbar

%keyboard

h2 = figure(2);
h2.Position = [11 -31 560 420];
dum = set(h2);

h3 = figure(3);
h3.Position = [10 381 560 420];
dum = set(h3);


h5 = figure(5);
h5.Position = [229 1 1008 804];
dum = set(h5);
clf
hold on
repnum = 256; % 32
mulnum = 64/repnum;
for reploop = 1:repnum

clear yvar py

% px = 2*(rand-0.499)*pmax;
% yp1 = 2*(rand-0.499)*sqrt(2*(E-px^2/2));
% py = 2*(rand-0.499)*sqrt(2*(E-px^2/2-yp1^2/2));
% xp1 = sign(rand-0.5)*real(sqrt(2*(E-px^2/2-yp1^2/2-py^2/2)));

%Etest = px^2/2 + py^2/2 + xp1^2/2 + yp1^2/2

px = 2*(rand-0.499)*pmax;
py = sign(rand-0.5)*real(sqrt(2*(E-px^2/2)));
xp1 = 0;
yp1 = 0;

y0 = [xp1 yp1 px py];



tspan = [1 2000]; % 500

% figure(1)
%options = odeset('OutputFcn',@odephas3);
options = odeset('RelTol',1e-9);

[t,y] = ode45(@f1,tspan,y0,options);

siz = size(t);
y1 = y(:,1);
y2 = y(:,2);
y3 = y(:,3);
y4 = y(:,4);

figure(2)
plot(t(1:siz),y(1:siz,1))
xlabel('Time')
legend('speed')

figure(3)
plot(y(1:siz,1),y(1:siz,2))
xlabel('x')
ylabel('y')


% Power Spectrum
Pow = 0;
if Pow == 1
yf = y(:,1);
F = fft(yf);
Pow = F.*conj(F);

figure(4)
plot(Pow)
end


%First Return Map
Fst = 1;
if Fst == 1
cnt = 0;
last = y1(1);
for loop = 2:siz
if (last < 0)&(y1(loop) > 0)
cnt = cnt+1;
del1 = -y1(loop-1)./(y1(loop) - y1(loop-1));
py(cnt) = y4(loop-1) + del1*(y4(loop)-y4(loop-1));
yvar(cnt) = y2(loop-1) + del1*(y2(loop)-y2(loop-1));
last = y1(loop);
else
last = y1(loop);
end
end

figure(5)
plot(yvar,py,'o','MarkerSize',2,'Color',h(floor(0.9*mulnum*reploop)+1,:),'MarkerFaceColor',h(floor(0.9*mulnum*reploop)+1,:))
xlabel('y')
ylabel('py')
end

pause(0.001)

end % end reploop

figure(5)
hold off
print -dtiff HeilesPrint


% Model function
function dy = f1(t,y)

B = 0.00; %0.002
dy = zeros(4,1);
dy(1) = y(3);
dy(2) = y(4);
dy(3) = -y(1) - epsE*(2*y(1)*y(2) + B*y(3));
dy(4) = -y(2) - epsE*(y(1)^2 - y(2)^2 + B*y(4));

end % end f5


end % end ltest

Loading

0 comments on commit 57a1f9f

Please sign in to comment.