# How to code double pendulum by using rk4

73 views (last 30 days)
numpy on 27 May 2024 at 4:51
Answered: Milan Bansal on 6 Jun 2024 at 20:27
It has to satisfy these conditions below.
Simulate the motion of the double pendulum for the following two initial conditions and observe the difference in motion (butterfly effect)
Initial conditions 1: Initial angles theta=pi/2, omg=pi/2 (initial speeds are all 0)
Initial conditions 2: Initial angles theta=pi/2, omg=pi/2+0.001 (initial speeds are all 0)
Precautions 1: Use RK4
Precautions 2: fps should be 30 frames per second
Precautions 3: Calculate by changing dt=1/30/50, 1/30/100, 1/30/200, 1/30/400, etc. and find a reliable size of dt
Precautions 4: Simulate 25 seconds of exercise
##### 3 CommentsShow 1 older commentHide 1 older comment
numpy on 28 May 2024 at 3:18
The simulation of double pendulum was successful, but RK4 was not applied
numpy on 28 May 2024 at 3:19
Edited: Sam Chak on 28 May 2024 at 4:59
%% clear everything section
clear;
close all;
clc;
%% parameter section
g = 9.81;
L1 = 1.0;
L2 = 1.0;
m1 = 1.0;
m2 = 1.0;
theta1_0 = pi / 2;
theta2_0 = pi / 2;
omega1_0 = 0;
omega2_0 = 0;
%% solve ode section
tspan = [0 25];
dt = 0.001;
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t, y] = ode45(@(t, y) doublePendulumODE(t, y, g, L1, L2, m1, m2), tspan, [theta1_0, omega1_0, theta2_0, omega2_0], options);
%% make animation section
theta1 = y(:, 1);
theta2 = y(:, 3);
x1 = L1 * sin(theta1);
y1 = -L1 * cos(theta1);
x2 = x1 + L2 * sin(theta2);
y2 = y1 - L2 * cos(theta2);
frames = [];
animationFig = figure;
animationLimits = [-3 3 -3 1];
axis(animationLimits);
for i = 1:length(t)
plot([0, x1(i)], [0, y1(i)], 'b', 'LineWidth', 2);
hold on;
plot(x1(i), y1(i), 'bo', 'MarkerSize', 20, 'MarkerFaceColor', 'b');
plot([x1(i), x2(i)], [y1(i), y2(i)], 'r', 'LineWidth', 2);
plot(x2(i), y2(i), 'ro', 'MarkerSize', 20, 'MarkerFaceColor', 'r');
xlabel('X Position (m)');
ylabel('Y Position (m)');
title(['Double Pendulum Animation - Time: ', num2str(t(i), '%.2f'), 's']);
ylim(animationLimits(3:4));
xlim(animationLimits(1:2));
frame = getframe(animationFig);
frames = [frames, frame];
if i < length(t)
cla(animationFig);
end
end
close(animationFig);
figure;
movie(frames, 1, 30);
%% dynamics of double pendulum section
function dydt = doublePendulumODE(t, y, g, L1, L2, m1, m2)
theta1 = y(1);
omega1 = y(2);
theta2 = y(3);
omega2 = y(4);
dydt = zeros(4, 1);
dydt(1) = omega1;
dydt(2) = (-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)));
dydt(3) = omega2;
dydt(4) = (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)));
end

Milan Bansal on 6 Jun 2024 at 20:27
Hi numpy,
I understand that you want to simulate the motion of double pendulum in MATLAB using RK4 method but are facing issues with it. You also wish to run the simulation for different values of dt"
To simulate the motion of the double pendulum with the specified initial conditions using the RK4 (Runge-Kutta 4th order) method, it is required to modify the provided code. The code provided uses ode45, which is a built-in MATLAB ODE solver using a variable-step Runge-Kutta method, but it is not RK4 and does not allow explicit control over the time step as specified in the problem.
Implement the RK4 method manually and run the simulation for different time steps (dt).
1. Implement the RK4 method for solving the ODEs.
2. Simulate the double pendulum for the given initial conditions.
3. Vary the dt to find a reliable time step.
4. Generate and observe the animation for the specified initial conditions.
Here is how you can modify your code to implement the RK4 method.
%% clear everything section
clear;
close all;
clc;
%% parameter section
g = 9.81;
L1 = 1.0;
L2 = 1.0;
m1 = 1.0;
m2 = 1.0;
theta1_0 = pi / 2;
theta2_0 = pi / 2;
omega1_0 = 0;
omega2_0 = 0;
dt_values = [1/30/50, 1/30/100, 1/30/200, 1/30/400]; % Different dt values to test
t_total = 25; % Total simulation time
fps = 30; % Frames per second
num_frames = t_total * fps;
%% Main section to test different dt values and initial conditions
for dt = dt_values
fprintf('Simulating with dt = %.6f\n', dt);
simulate_pendulum(dt, g, L1, L2, m1, m2, theta1_0, theta2_0, omega1_0, omega2_0, t_total, num_frames, fps, 0);
simulate_pendulum(dt, g, L1, L2, m1, m2, theta1_0, theta2_0, omega1_0, omega2_0, t_total, num_frames, fps, 0.001);
end
Function to implement RK4 Method
%% RK4 method implementation
function y_next = rk4_step(f, t, y, dt)
k1 = f(t, y);
k2 = f(t + dt / 2, y + dt / 2 * k1);
k3 = f(t + dt / 2, y + dt / 2 * k2);
k4 = f(t + dt, y + dt * k3);
y_next = y + dt / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
end
Function to implement the dynamics (same as given in your code)
%% dynamics of double pendulum section
function dydt = doublePendulumODE(t, y, g, L1, L2, m1, m2)
theta1 = y(1);
omega1 = y(2);
theta2 = y(3);
omega2 = y(4);
dydt = zeros(4, 1);
dydt(1) = omega1;
dydt(2) = (-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)));
dydt(3) = omega2;
dydt(4) = (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)));
end
Function to simulate and animate the double pendulum with different initial conditions.
function simulate_pendulum(dt, g, L1, L2, m1, m2, theta1_0, theta2_0, omega1_0, omega2_0, t_total, num_frames, fps, initial_condition_variation)
%% Modification
t = 0:dt:t_total;
y = zeros(length(t), 4);
y(1, :) = [theta1_0, omega1_0, theta2_0 + initial_condition_variation, omega2_0];
for i = 1:length(t) - 1
y(i + 1, :) = rk4_step(@(t, y) doublePendulumODE(t, y, g, L1, L2, m1, m2), t(i), y(i, :)', dt);
end
theta1 = y(:, 1);
theta2 = y(:, 3);
x1 = L1 * sin(theta1);
y1 = -L1 * cos(theta1);
x2 = x1 + L2 * sin(theta2);
y2 = y1 - L2 * cos(theta2);
frames = [];
animationFig = figure;
animationLimits = [-3 3 -3 1];
axis(animationLimits);
for i = 1:round(length(t) / num_frames):length(t)
plot([0, x1(i)], [0, y1(i)], 'b', 'LineWidth', 2);
hold on;
plot(x1(i), y1(i), 'bo', 'MarkerSize', 20, 'MarkerFaceColor', 'b');
plot([x1(i), x2(i)], [y1(i), y2(i)], 'r', 'LineWidth', 2);
plot(x2(i), y2(i), 'ro', 'MarkerSize', 20, 'MarkerFaceColor', 'r');
xlabel('X Position (m)');
ylabel('Y Position (m)');
title(['Double Pendulum Animation - Time: ', num2str(t(i), '%.2f'), 's']);
ylim(animationLimits(3:4));
xlim(animationLimits(1:2));
frame = getframe(animationFig);
frames = [frames, frame];
if i < length(t)
cla(animationFig);
end
end
close(animationFig);
figure;
movie(frames, 1, fps);
end
Hope this helps