How to graph coupled differential equations

6 次查看（过去 30 天）
muhammad
muhammad 2025-4-10
编辑： Torsten 2025-4-10
Ran in:
I have the following set of equations of motions for 2 masses and I am trying to plot the x displacement for each of the masses.
The MATLAB code i have for it is:
% Parameters
m1 = 4.0;
m2 = 6.0;
L = 1.5;
k = 100.0; % Added missing parameter
g = 9.81;
F0 = 100;
tF = 1;
% Applied force F(t)
F = @(t) (t <= tF) * F0 .* sin(2*pi*t / tF);
% ODE function
odefun = @(t, y) [
y(2); % dx/dt
compute_ddx(t, y, m1, m2, k, L, g, F); % d²x/dt²
y(4); % dtheta/dt
compute_ddtheta(t, y, m1, m2, k, L, g, F); % d²theta/dt²
y(6); % dr/dt
compute_ddr(t, y, m1, m2, k, L, g, F) % d²r/dt²
];
% Initial state: [x, dx, theta, dtheta, r, dr]
y0 = [0; 0; 0; 0; L; 0];
% Duration
tspan = [0 3];
% Solve with tighter tolerances
options = odeset('RelTol',1e-6,'AbsTol',1e-8);
[t, Y] = ode45(odefun, tspan, y0, options);
% Extract states
x = Y(:,1); % x-displacement of m1
theta = Y(:,3);
r = Y(:,5); % Use actual r(t) instead of constant L
% x displacement of m2 (x_m1 + r*sinθ)
x_m2 = x + r .* sin(theta);
% Create two separate figures
figure(1);
plot(x, t, 'b', 'LineWidth', 2)
xlabel('x m1 [m]');
ylabel('t [s]');
title('m1 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');
grid on;
figure(2);
plot(x_m2, t, 'r', 'LineWidth', 2)
xlabel('x m2 [m]');
ylabel('t [s]');
title('m2 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');
grid on;
% Second derivative functions
function ddx = compute_ddx(t, y, m1, m2, k, L, g, F)
Ft = F(t);
theta = y(3); dtheta = y(4);
r = y(5); dr = y(6);
sinT = sin(theta); cosT = cos(theta);
% % Corrected numerator with proper operators
% numx = Ft - m2*(r*dtheta^2 + g*cosT + k*(L-r))*sinT - 2*m2*dr*dtheta*cosT + ...
% m2*r*(dtheta^2*sinT - (-2*dr*dtheta - g*sinT)*cosT/r);
%
% % Corrected denominator
% denx = m1 + m2 - m2*sinT^2 - m2*cosT^2;
% ddx = numx/denx;
ddx = (Ft - m2*k*(L-r)*sinT - 2*m2*r*dtheta^2*sinT)/( m1 + m2);
end
function ddtheta = compute_ddtheta(t, y, m1, m2, k, L, g, F)
theta = y(3);
r = y(5); dr = y(6);
ddx = compute_ddx(t, y, m1, m2, k, L, g, F);
dtheta = y(4);
% Corrected ddtheta calculation
ddtheta = (-ddx*cos(theta) - 2*dr*dtheta - g*sin(theta));
ddtheta = ddtheta/r;
end
function ddr = compute_ddr(t, y, m1, m2, k, L, g, F)
theta = y(3); dtheta = y(4);
r = y(5);
ddx = compute_ddx(t, y, m1, m2, k, L, g, F);
% Corrected ddr calculation
ddr = -ddx*sin(theta) + r*dtheta^2 + g*cos(theta) - k*(L-r)/m2;
end
Expectation
The graphs are supposed to come out as the green curves on these two graphs yet my code does not do that. I have checked that ddx, ddr and ddtheta have been computed correctly multiple times so the only thing i can think of it the way i am coding it.
Can someone pls assist?

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

回答（1 个）

Torsten
Torsten 2025-4-10
编辑：Torsten 2025-4-10
Ran in:
% Parameters
m1 = 4.0;
m2 = 6.0;
L = 1.5;
k = 100.0; % Added missing parameter
g = 9.81;
F0 = 100;
tF = 1;
% Initial state: [x, dx, theta, dtheta, r, dr]
y0 = [0; 0; 0; 0; L; 0 ];
% Duration
tspan = [0 3];
% Solve with tighter tolerances
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
[t, Y] = ode45(@(t,y)odefun(t,y,m1, m2, k, L, g, F0, tF), tspan, y0, options);
% Extract states
x = Y(:,1);
theta = Y(:,3);
r = Y(:,5);
x_m2 = x + r .* sin(theta);
% Create two separate figures
figure(1);
plot(x, t, 'b', 'LineWidth', 2)
xlabel('x m1 [m]');
ylabel('t [s]');
title('m1 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');
grid on;
figure(2);
plot(x_m2, t, 'r', 'LineWidth', 2)
xlabel('x m2 [m]');
ylabel('t [s]');
title('m2 Displacement vs Time');
xlim([-0.5, 3.5]);
ylim([0, 3]);
yticks(0:1:3)
set(gca, 'YDir', 'reverse');
grid on
function dydt = odefun(t,y, m1, m2, k, L, g, F0, tF)
x = y(1);
xdot = y(2);
theta = y(3);
thetadot = y(4);
r = y(5);
rdot = y(6);
F = (t <= tF) * F0 .* sin(2*pi*t / tF);
A = zeros(3);
b = zeros(3,1);
A(1,1) = m1+m2;
A(1,2) = m2*r*cos(theta);
A(1,3) = m2*sin(theta);
A(2,1) = m2*r*cos(theta);
A(2,2) = m2*r^2;
A(2,3) = 0;
A(3,1) = m2*sin(theta);
A(3,2) = 0;
A(3,3) = m2;
b(1) = F - 2*m2*rdot*thetadot*cos(theta) + m2*r*thetadot^2*sin(theta);
b(2) = -2*m2*r*rdot*thetadot - m2*g*r*sin(theta);
b(3) = m2*r*thetadot^2 + m2*g*cos(theta) + m2*k*(L-r);
ddot = A\b;
xddot = ddot(1);
thetaddot = ddot(2);
rddot = ddot(3);
dydt = [xdot;xddot;thetadot;thetaddot;rdot;rddot];
end

类别

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

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by