MATLAB Answers

I have a problem in solving and animation

4 views (last 30 days)
I have written a code for this problem to solve the equations and, at the end, get the equations, but I get an error.
Can anybody help me?
% Simulation of coupled pendulum by Lagrangian mechanics
close all; clear; clc
% generalized coordinates
syms t dum_
theta = str2sym('theta(t)');
phi = str2sym('phi(t)');
% constants, length, mass, g, geometry
L_1 = 7.5;
L_2 = 7.5;
L_3 = 7.5;
m_1 = 4;
m_2 = 4;
m_3 = 4;
g = 9.81;
d_0 = 15; % rest length spring
% positions and velocities as function of the generalized coordinates=
x1 = L_1 * cos(theta);
y1 = L_1 * sin(theta);
x2 = 2 * L_2 * cos(theta);
y2 = 0;
x3 = 2*L_1 * cos(theta)+L_3*sin(phi);
y3 = L_3 * cos(phi);
x1_dot = diff(x1, t);
x2_dot = diff(x2, t);
y1_dot = diff(y1, t);
y2_dot = diff(y2, t);
x3_dot = diff(x3, t);
y3_dot = diff(y3, t);
% kinetic and potential energy
T = m_1/2 * (x1_dot^2 + y1_dot^2) + m_2/2 * (x2_dot^2 + y2_dot^2)+ m_3/2 * (x3_dot^2 + y3_dot^2);
k = 0.5;
V = -m_1 * g * y1 - m_3 * g * y3 +...
1/2 * k * (d_0-2*L_1*cos(theta))^2;
% determine for which theta = alpha and phi = beta the system is at rest
% alpha = sym('alpha');
% beta = sym('beta');
% v = subs(V, {theta, phi}, {alpha, beta});
% [alpha, beta] = vpasolve([diff(v, alpha), diff(v, beta)], [alpha, beta]);
% Lagrangian
L = T - V;
% dL/d(qdot)
dL_dthetadot = subs(diff(subs(L, diff(theta, t), dum_), dum_), dum_, diff(theta, t));
dL_dphidot = subs(diff(subs(L, diff(phi, t), dum_), dum_), dum_, diff(phi, t));
% dL/dq
dL_dtheta = subs(diff(subs(L, theta, dum_), dum_), dum_, theta);
dL_dphi = subs(diff(subs(L, phi, dum_), dum_), dum_, phi);
% dFdq
k = 0.25; % dissipation constant
% generalized equations of motion
deq_1 = diff(dL_dthetadot, t) - dL_dtheta;
deq_2 = diff(dL_dphidot, t) - dL_dphi;
% abbreviation of variables
variables = {theta, phi, diff(theta, t), diff(phi, t), diff(theta, t, 2), diff(phi, t, 2)};
variables_short = arrayfun(@str2sym, {'x(1)', 'x(2)', 'x(3)', 'x(4)', 'thetaddot', 'phiddot'});
deq_1 = subs(deq_1, variables, variables_short);
deq_2 = subs(deq_2, variables, variables_short);
% solve for thetaddot, phiddot
solution = solve(deq_1, deq_2, str2sym('thetaddot'), str2sym('phiddot'));
THETADDOT = solution.thetaddot;
PHIDDOT = solution.phiddot;
% solve non linear ode system
time = linspace(0, 60, 2000);
% initial conditions [theta, phi, thetadot, phidot]
x_0 = [-pi/4 pi/6 0 0];
str = ['x_dot = @(t, x)[x(3); x(4);', char(THETADDOT), ';', char(PHIDDOT), '];'];
eval(str);
[t, q] = ode45(x_dot, time, x_0);
% Calculute positions as function of generalized coordinates
X1 = L_1 * cos(q(:, 1));
Y1 = L_1 * sin(q(:, 1));
X2 = 2* L_2 * cos(q(:, 2));
Y2 = 0;
X3= 2*L_1 * cos(q(:, 3))+L_1*sin(q(:, 3));
Y3 = -L_3 * cos(q(:, 3));
% plot solution
set(gcf, 'color', 'w')
set(gcf, 'position', [10, 100, 750, 750])
h = plot([]);
hold on
box on
axis equal
for i = 1 : numel(time)
if ~ishghandle(h)
break
end
cla
plot([0, X1(i)], [0, Y1(i)], 'k', 'Linewidth', 2);
plot(X1(i), Y1(i), 'o', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k', 'MarkerSize', 4 * m_1);
plot([X1(i), X2(i)], [Y1(i), Y2(i)], 'k', 'Linewidth', 2);
plot(X2(i), Y2(i), 'o', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k', 'MarkerSize', 4 * m_2);
plot([X2(i), X3(i)], [Y2(i), Y3(i)], 'k', 'Linewidth', 2);
plot(X3(i), Y3(i), 'o', 'MarkerFaceColor', 'k', 'MarkerEdgeColor', 'k', 'MarkerSize', 4 * m_3);
axis([-12, 60, -10, 20]);
h = draw_spring_2D([0; 0], [X2(i); 0], 12, 0.5);
drawnow
end
function h = draw_spring_2D(A, B, number_of_coils, y_amplitude)
persistent t
normalvector_AB = (B - A) / norm(B - A);
offset_A = A + 1.25 * normalvector_AB;
offset_B = B - 1.25 * normalvector_AB;
distance_between_offsets = norm(offset_B - offset_A);
t = linspace(-pi, number_of_coils * 2 * pi, 500);
x_coordinate_between_offsets = distance_between_offsets * linspace(0, 1, numel(t));
% ratio between x amplitude and y
ratio_X_div_Y = 0.5;
x = x_coordinate_between_offsets + ratio_X_div_Y * y_amplitude * cos(t);
y = y_amplitude * sin(t);
coil_positions = [x; y];
rotation_matrix = [normalvector_AB, null(normalvector_AB')];
rotated_coil_positions = rotation_matrix * coil_positions;
h = plot([A(1), offset_A(1) + rotated_coil_positions(1,:), B(1)], ...
[A(2), offset_A(2) + rotated_coil_positions(2,:), B(2)], 'k');
end
  1 Comment
John D'Errico
John D'Errico on 25 Sep 2021
Never just say you got an error. Show the COMPLETE error message, thus everything in red. If all you do is say you got an error, how can we know what you did wrong? Perhaps you don't even know how to execute the script you wrote?

Sign in to comment.

Accepted Answer

Sulaymon Eshkabilov
Sulaymon Eshkabilov on 25 Sep 2021
There was one potential size related err in Y2 that is fixed. Plug in this into your code and then your simulation works OK.
...
X2 = 2* L_2 * cos(q(:, 2));
Y2 = zeros(size(X2)); % Size of Y2 must match with the one of X2
...

More Answers (0)

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!