Back to Gallery
  • Remix
  • Share
  • New Entry
on 15 Nov 2023
  • 26
  • 74
  • 0
  • 4
  • 1556
drawframe(1);
Write your drawframe function below
function drawframe(f)
persistent phis
persistent l
persistent tpts
numP = 50;
if isempty(tpts)
[phis, l, tpts] = computeDyn(numP);
end
figure(1)
col='rbgcmyk';
txt="50 chaotic pendula";
title(txt)
rmax=1.1*sum(l);
xlim(rmax*[-1,1])
ylim(rmax*[-1,1])
axis equal
ti = f;
for jj=1:numP
y=squeeze(phis(jj,ti,:));
end
plot(0,0,'sq')
hold on
xlim(rmax*[-1,1])
ylim(rmax*[-1,1])
for jj=1:numP
y=squeeze(phis(jj,ti,1:2));
[x1,y1,x2,y2]=state_2_cartesian(y,l);
line([0 x1],[0 y1],'Color',col(rem(jj-1,5)+1))
plot(x1,y1,['o',col(rem(jj-1,5)+1)],'MarkerSize',4,"MarkerFaceColor",col(rem(jj-1,5)+1))
line([x1 x2],[y1 y2],'Color',col(rem(jj-1,5)+1))
plot(x2,y2,['o',col(rem(jj-1,5)+1)],'MarkerSize',7,"MarkerFaceColor",col(rem(jj-1,5)+1))
end
hold off
drawnow
title(txt)
xlabel("x_{lateral}")
ylabel("height")
end
function dydt=dp(~,y,m,l)
dphi=y(1)-y(2);
a=(m(1)+m(2))*l(1);
b=m(2)*l(2)*cos(dphi);
c=m(2)*l(1)*cos(dphi);
d=m(2)*l(2);
e=-m(2)*l(2)*y(4)^2*sin(dphi)-10*(m(1)+m(2))*sin(y(1));
f=m(2)*l(1)*y(3)^2*sin(dphi)-m(2)*10*sin(y(2));
dydt=[y(3);y(4);(e*d-b*f)/(a*d-c*b);(a*f-c*e)/(a*d-c*b)];
end
function [x1,y1,x2,y2]=state_2_cartesian(y,l)
x1=l(1)*sin(y(1));
y1=-l(1)*cos(y(1));
x2=x1+l(2)*sin(y(2));
y2=y1-l(2)*cos(y(2));
end
function [phis, l, tpts] = computeDyn(numP)
m=[1,5];
l=[1.25,1];
dydt=@(t,y)dp(t,y,m,l);
phi0=pi;
state_0=[phi0*(1+0.01*(rand(numP,2)-0.5)),zeros(numP,2)];
tSp=[0,2.5];
tsMax=0.05;
tSol=tSp(1):tsMax:tSp(2);
tpts=numel(tSol);
for jj=1:numP
sol(jj)=ode45(dydt,tSp,state_0(jj,:));
phis(jj,:,:)=deval(sol(jj),tSol)';
r1(jj,:,:)=l*[sin(phis(jj,:,1));-cos(phis(jj,:,1))];
dr2(jj,:,:)=l(2)*[sin(phis(jj,:,2));-cos(phis(jj,:,2))];
r2(jj,:,:)=r1(jj,:,:)+dr2(jj,:,:);
end
end
Animation
Remix Tree