function sticks % Chaotic Swinging Sticks % https://www.the-swinging-sticks.com/en % https://blogs.mathworks.com/cleve/2024/02/20/chaotic-swinging-sticks/ % 02/20/2024 function M = swingmass(~,u) % Mass matrix c = cos(u(1)-u(2)); M = [1 0 0 0; 0 1 0 0; 0 0 5 c; 0 0 c 1]; end % ------------------------ function f = swingrhs(~,u) % Driving force g = 1; s = sin(u(1)-u(2)); f = [u(3); u(4); -g*sin(u(1))-s*u(4)^2; -g*sin(u(2))+s*u(3)^2]; end % ------------------------ function theta = swinginit(x,y) r = norm([x,y]); if r > 2 alpha = 0; else alpha = acos(r/2); end beta = atan2(y,x) + pi/2; theta = [beta+alpha; beta-alpha]; end % ------------------------ function status = swingplot(t,u,task) if ~isequal(task,'done') theta = u(1:2); x = sin(theta); y = -cos(theta); switch task case 'init' % Initialize plot plt = plot( ... [-x(1) 3*x(1) NaN 3*x(1)-x(2) 3*x(1)+2*x(2)], ... [-y(1) 3*y(1) NaN 3*y(1)-y(2) 3*y(1)+2*y(2)], ... '-','markersize',14,'linewidth',3); plt(2) = line(plt.XData(2),plt.YData(2), ... 'marker','.','markersize',24); line(0,0,'marker','.','markersize',24) axis(6*[-1 1 -1 1]) axis square set(gca,'xtick','','ytick','') titl = title(sprintf('t = %8.1f',t)); set(gcf,'userdata',{plt,titl}) shg gif_frame('sticks_gif.gif',inf') case '' % Update plot fud = get(gcf,'userdata'); [plt,titl] = deal(fud{:}); set(plt(1), ... 'xdata',[-2*x(1) 3*x(1) NaN 3*x(1)-x(2) 3*x(1)+2*x(2)], ... 'ydata',[-2*y(1) 3*y(1) NaN 3*y(1)-y(2) 3*y(1)+2*y(2)]) set(plt(2),... 'xdata',plt(1).XData(2), ... 'ydata',plt(1).YData(2)) % Display time in title set(titl,'string',sprintf('t = %8.1f',t)); gif_frame % writeVideo(VidObj,getframe(gcf)) pause(0.025) end status = false; end end % ------------------------ % Initial position x0 = 0.1; y0 = pi/2; theta0 = swinginit(x0,y0); % Initialize figure and buttons. clf reset shg set(gcf,'name','swinging sticks','menu','none','numbertitle','off') % Start the ode solver % VidObj = VideoWriter('sticks.mp4','mpeg-4'); % open(VidObj) delt = 0.1; u0 = [theta0; 0; 0]; tspan = [0 10]; opts = odeset('mass',@swingmass,'outputfcn',@swingplot,'maxstep',delt); ode23(@swingrhs,tspan,u0,opts); gif_frame('wrap') % close(VidObj) end