double pendulum in matlab

main.m
%% Parameters
g = 9.81; % acceleration due to gravity (m/s^2)
L1 = 1; % length of first pendulum (m)
L2 = 1; % length of second pendulum (m)
m1 = 1; % mass of first pendulum (kg)
m2 = 1; % mass of second pendulum (kg)

%% Initial conditions
theta1_0 = pi/2; % initial angle of first pendulum (rad)
theta2_0 = pi/2; % initial angle of second pendulum (rad)
omega1_0 = 0; % initial angular velocity of first pendulum (rad/s)
omega2_0 = 0; % initial angular velocity of second pendulum (rad/s)
y0 = [theta1_0; theta2_0; omega1_0; omega2_0];

%% Simulation
tspan = [0 10]; % simulation time interval (s)
opts = odeset('RelTol',1e-6,'AbsTol',1e-6); % ODE solver options
[t, y] = ode45(@(t,y) double_pendulum_odefun(t, y, g, L1, L2, m1, m2), tspan, y0, opts); % solve ODEs

%% Plotting
theta1 = y(:,1);
theta2 = y(:,2);

x1 = L1*sin(theta1);
y1 = -L1*cos(theta1);

x2 = x1 + L2*sin(theta2);
y2 = y1 - L2*cos(theta2);

figure;
plot(x1, y1, 'r', x2, y2, 'b');
axis equal;
xlabel('x (m)');
ylabel('y (m)');
legend('Pendulum 1', 'Pendulum 2');

function dy = double_pendulum_odefun(t, y, g, L1, L2, m1, m2)
theta1 = y(1);
theta2 = y(2);
omega1 = y(3);
omega2 = y(4);

dtheta1 = omega1;
dtheta2 = omega2;

domega1 = (-g*(2*m1 + m2)*sin(theta1) - m2*g*sin(theta1 - 2*theta2) - 2*sin(theta1 - theta2)*m2*(omega2^2*L2 + omega1^2*L1*cos(theta1 - theta2))) / (L1*(2*m1 + m2 - m2*cos(2*theta1 - 2*theta2)));
domega2 = (2*sin(theta1 - theta2) * (omega1^2*L1*(m1+m2) + g*(m1+m2)*cos(theta1) + omega2^2*L2*m2*cos(theta1-theta2))) / (L2*(2*m1+m2-m2*cos(2*theta1-2*theta2)));

dy = [dtheta1; dtheta2; domega1; domega2];
end
1602 chars
51 lines