How to code double pendulum by using rk4

29 次查看（过去 30 天）
numpy
numpy 2024-5-27
回答： Milan Bansal 2024-6-6
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
  3 个评论
显示 1更早的评论
numpy
numpy 2024-5-28
The simulation of double pendulum was successful, but RK4 was not applied
numpy
numpy 2024-5-28
编辑：Sam Chak 2024-5-28
%% 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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Milan Bansal
Milan Bansal 2024-6-6
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

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by