How to code double pendulum by using rk4
10 次查看(过去 30 天)
显示 更早的评论
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 个)
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).
- Implement the RK4 method for solving the ODEs.
- Simulate the double pendulum for the given initial conditions.
- Vary the dt to find a reliable time step.
- 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
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!