I have a problem in solving and animation

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]

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?

Answer 1

[Sulaymon Eshkabilov]

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
...

[Mehrdad Nasirshoaibi]

Thanks man. It works.

[Sulaymon Eshkabilov]

Most welcome!