4 views (last 30 days)

Show older comments

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

John D'Errico
on 25 Sep 2021

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

...

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

Start Hunting!