double pendulum animation in matlab

main.m
%% Define parameters
m1 = 1;   % mass of pendulum 1
m2 = 1;   % mass of pendulum 2
L1 = 1;   % length of pendulum 1
L2 = 1;   % length of pendulum 2
g = 9.81; % gravitational constant

%% Define initial conditions
q0 = [pi/2; pi/2]; % initial angles
p0 = [0; 0];       % initial angular velocities

%% Define simulation parameters
t_span = [0 20];   % time span of simulation
dt = 0.05;         % time step

%% Define simulation function
odefun = @(t, y) double_pendulum_odefun(t, y, m1, m2, L1, L2, g);

%% Run simulation
[t, q] = ode45(odefun, t_span, [q0; p0], odeset('RelTol', 1e-10, 'AbsTol', 1e-10));

%% Create animation
figure;
p1 = [L1*sin(q(:,1)), -L1*cos(q(:,1))];
p2 = [L2*sin(q(:,1)+q(:,2)), -L2*cos(q(:,1)+q(:,2))];

for i = 1:length(t)
    plot([0 p1(i,1) p2(i,1)], [0 p1(i,2) p2(i,2)], '-o', 'LineWidth', 2, 'MarkerSize', 10, 'MarkerFaceColor', 'r');
    xlim([-2 2]);
    ylim([-2 2]);
    title(sprintf('Double Pendulum Animation (t = %0.2f s)', t(i)));
    drawnow;
end
989 chars
34 lines

main.m
function [dydt] = double_pendulum_odefun(t, y, m1, m2, L1, L2, g)
theta1 = y(1);
theta2 = y(2);
p1 = y(3);
p2 = y(4);
dtheta1 = p1 / (m1 * L1^2) - (m1 + m2) * g / (m1 * L1) * sin(theta1) - m2 / (m1 * L1^2) * ...
    (p2 / p1 * cos(theta1 - theta2) * sin(theta1 - theta2) + (L2 / L1) * dtheta2^2 * sin(theta1 - theta2));
dtheta2 = p2 / (m2 * L2^2) - g / L2 * sin(theta2) / L1 * sin(theta1 - theta2) ...
    - p1 / (m2 * L2^2) * cos(theta1 - theta2) * sin(theta1 - theta2);
dp1 = -(m1 + m2) * g * L1 * sin(theta1) - m2 * L1 * L2 * (dtheta1^2 * sin(theta1 - theta2) + dtheta2^2 * sin(theta2)) * cos(theta1 - theta2);
dp2 = -m2 * g * L2 * sin(theta2) + m2 * L1 * L2 * (dtheta1^2 * sin(theta1 - theta2) + dtheta2^2 * sin(theta2)) * cos(theta1 - theta2);
dydt = [dtheta1; dtheta2; dp1; dp2];
end
790 chars
14 lines