What şs wrong ?

% Mechanism parameters

AB = [0.09, 0.15, 0.22, 0.16]; % Link lengths AB (m)

BC = [0.40, 0.67, 1.00, 0.70]; % Link lengths BC (m)

CE = [0.25, 0.45, 0.65, 0.50]; % Link lengths CE (m)

CD = [0.12, 0.22, 0.35, 0.25]; % Link lengths CD (m)

EF = [0.21, 0.32, 0.60, 0.48]; % Link lengths EF (m)

a = [0.22, 0.33, 0.55, 0.90]; % Link lengths a (m)

b = [0.35, 0.60, 0.40, 0.60]; % Link lengths b (m)

C = [0.40, 0.65, 1.20, 0.70]; % Link lengths C (m)

n = [60, 135, 240]; % Angular speeds (rpm)

rpm_to_rad_per_sec = 2*pi/60; % Conversion factor from RPM to rad/s

% Input link's angular position values

input_angles = linspace(0, 2*pi, 100); % Assuming 100 points for a full rotation

% Initialize arrays to store results

angular_displacement = zeros(length(input_angles), length(AB));

angular_velocity = zeros(length(input_angles), length(AB));

angular_acceleration = zeros(length(input_angles), length(AB));

linear_displacement = zeros(length(input_angles), length(AB));

linear_velocity = zeros(length(input_angles), length(AB));

linear_acceleration = zeros(length(input_angles), length(AB));

% Numerical solution

for i = 1:length(input_angles)

% Calculate kinematics for the given input angle

for j = 1:length(AB)

% Angular displacement

angular_displacement(i, j) = acos((AB(j)^2 + BC(j)^2 - CE(j)^2) / (2 * AB(j) * BC(j)));

% Angular velocity (derivative of angular displacement)

angular_velocity(i, j) = -AB(j) * sin(angular_displacement(i, j)) * n(j) * rpm_to_rad_per_sec;

% Angular acceleration (derivative of angular velocity)

angular_acceleration(i, j) = -AB(j) * cos(angular_displacement(i, j)) * n(j)^2 * rpm_to_rad_per_sec^2;

% Linear displacement

linear_displacement(i, j) = AB(j) * sin(angular_displacement(i, j));

% Linear velocity (derivative of linear displacement)

linear_velocity(i, j) = AB(j) * cos(angular_displacement(i, j)) * angular_velocity(i, j);

% Linear acceleration (derivative of linear velocity)

linear_acceleration(i, j) = -AB(j) * sin(angular_displacement(i, j)) * angular_velocity(i, j)^2 ...

+ AB(j) * cos(angular_displacement(i, j)) * angular_acceleration(i, j);

end

end

% Plotting

figure;

% Plot Angular Displacement

subplot(3, 2, 1);

plot(input_angles, rad2deg(angular_displacement));

title('Angular Displacement (degrees)');

% Plot Angular Velocity

subplot(3, 2, 2);

plot(input_angles, rad2deg(angular_velocity));

title('Angular Velocity (degrees/s)');

% Plot Angular Acceleration

subplot(3, 2, 3);

plot(input_angles, rad2deg(angular_acceleration));

title('Angular Acceleration (degrees/s^2)');

% Plot Linear Displacement

subplot(3, 2, 4);

plot(input_angles, linear_displacement);

title('Linear Displacement (m)');

% Plot Linear Velocity

subplot(3, 2, 5);

plot(input_angles, linear_velocity);

title('Linear Velocity (m/s)');

% Plot Linear Acceleration

subplot(3, 2, 6);

plot(input_angles, linear_acceleration);

title('Linear Acceleration (m/s^2)');

Answer 1

[Rishi]

Hi Furkan,

I understand from your query that you want to know the mistake in the provided code.

The array 'n' is defined as [60, 135, 240] in the following line:

n = [60, 135, 240];              % Angular speeds (rpm)

The error occurs in the following line, when the value of 'j' becomes 4. During this iteration, the code tries to acces the value of 'n(4)', which is wrong and hence, gives the following error:

% Angular velocity (derivative of angular displacement)
angular_velocity(i, j) = -AB(j) * sin(angular_displacement(i, j)) * n(j) * rpm_to_rad_per_sec;

Index exceeds the number of array elements. Index must not exceed 3.

To get rid of the error, you can define the array 'n' with 4 elements. For e.g, you can redefine 'n' as follows:

n = [60, 135, 240, 350];              % Angular speeds (rpm) (350 is a random value, replace it with the one that suits your case)

Hope this helps!

[Furkan]

thank you

[Furkan]

% Define physical parameters of the mechanism

AB = 0.09; % m, length of arm 1

BC = 0.4; % m, length of arm 2

CD = 0.7; % m, length of arm 3

La = 0.75; % m

Lb = 0.35; % m

% Define simulation parameters

Nmax = 100; % Maximum number of iterations

x = [20*pi/180, 320*pi/180]; % Adjusted initial guess for th2, th3 (in radians)

xe = 0.0001 * abs(x); % Slightly relaxed error tolerance

dth = 5 * pi / 180; % Step size for thq

thQ = 0:dth:2*pi; % Range of thq (in radians)

w1 = -0.2944 * ones(size(thQ)); % Constant angular velocity w1

acc1 = zeros(size(thQ)); % Constant angular acceleration acc1

% Preallocate arrays for results

th2 = zeros(size(thQ));

th3 = zeros(size(thQ));

w2 = zeros(size(thQ));

w3 = zeros(size(thQ));

acc2 = zeros(size(thQ));

acc3 = zeros(size(thQ));

% Main loop for each th2 value

for k = 1:length(thQ)

x_temp = x; % Temporary variable for the current iteration

kerr = 1; % Flag for convergence

% Iterative solver for position analysis

for n = 1:Nmax

th2_temp = x_temp(1);

th3_temp = x_temp(2);

J = [-BC*sin(th2_temp), -CD*sin(th3_temp); BC*cos(th2_temp), CD*cos(th3_temp)];

% Check for singularity in Jacobian

if cond(J) > 1e10

disp(['Warning: Jacobian is near singular at th2 index: ', num2str(k)]);

break;

end

f = [-(AB*cos(thQ(k)) + BC*cos(th2_temp) + CD*cos(th3_temp) - La); ...

-(AB*sin(thQ(k)) + BC*sin(th2_temp) + CD*sin(th3_temp) + Lb)];

eps = J\f; % More efficient than inv(J)*f

x_temp = x_temp + eps';

if all(abs(eps) < xe)

kerr = 0;

break;

end

end

if kerr == 1

disp(['Error: Solution did not converge at th2 index: ', num2str(k), ' (th2 = ', num2str(thQ(k) * 180/pi), ' degrees)']);

else

th2(k) = x_temp(1);

th3(k) = x_temp(2);

% Velocity analysis

fv = [-AB*sin(thQ(k))*w1(k); AB*cos(thQ(k))*w1(k)];

vel = J\fv;

w2(k) = vel(1);

w3(k) = vel(2);

% Acceleration analysis

fa = [-AB*cos(thQ(k))((w1(k))^2) - BC*cos(th2(k))((w2(k))^2) - CD*cos(th3(k))*((w3(k))^2) - AB*sin(thQ(k))*acc1(k); ...

AB*cos(thQ(k))acc1(k) - BC*sin(th2(k))((w2(k))^2) - CD*sin(th3(k))((w3(k))^2) - AB*sin(thQ(k))((w1(k))^2)];

acc = J\fa;

acc2(k) = acc(1);

acc3(k) = acc(2);

end

end

% Convert angles to degrees for plotting

th2d = th2 * 180/pi;

th3d = th3 * 180/pi;

% Plotting results

figure(1);

subplot(3,3,1), plot(thQ, th2d, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\theta_2 (deg)'), grid on;

subplot(3,3,2), plot(thQ, w2, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\omega_2 (rad/s)'), grid on;

subplot(3,3,3), plot(thQ, acc2, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\alpha_2 (rad/s^2)'), grid on;

subplot(3,3,4), plot(thQ, th3d, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\theta_3 (deg)'), grid on;

subplot(3,3,5), plot(thQ, w3, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\omega_3 (rad/s)'), grid on;

subplot(3,3,6), plot(thQ, acc3, 'r', 'LineWidth', 2), xlabel('\theta_Q (rad)'), ylabel('\alpha_3 (rad/s^2)'), grid on;

Can I ask you to look at this?

[Rishi]

The error occurs in the following line of code:

 % Acceleration analysis
        fa = [-AB*cos(thQ(k))((w1(k))^2) - BC*cos(th2(k))((w2(k))^2) - CD*cos(th3(k))*((w3(k))^2) - AB*sin(thQ(k))*acc1(k); ...
              AB*cos(thQ(k))acc1(k) - BC*sin(th2(k))((w2(k))^2) - CD*sin(th3(k))((w3(k))^2) - AB*sin(thQ(k))((w1(k))^2)];       

I am assuming you want to multiply the terms within the parantheses. To do so, you need to include the asterisk ('*') sign between the terms.

Here is the line after correction:

% Acceleration analysis
fa =[-AB*cos(thQ(k))*((w1(k))^2) - BC*cos(th2(k))*((w2(k))^2) - CD*cos(th3(k))*((w3(k))^2) - AB*sin(thQ(k))*acc1(k); AB*cos(thQ(k))*acc1(k) - BC*sin(th2(k))*((w2(k))^2) - CD*sin(th3(k))*((w3(k))^2) - AB*sin(thQ(k))*((w1(k))^2)];

Hope this helps!

[Furkan]

I'm sorry to bother you, but I'm trying to find a solution for this mechanism. Please help me. If you want, I would like to reach you via e-mail.