How to code double pendulum by using rk4

Please help me with this problem.

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

[James Tursa]

What have you done so far? What specific problems are you having with your code?

[numpy]

The simulation of double pendulum was successful, but RK4 was not applied

[numpy]

%% 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

Answer 1

[Milan Bansal]

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