My code is outputting very illogical answers for the second block (my second vector loop equation).

1 view (last 30 days)
The first block of code (E1) runs and comes out with a correct answer (validated in excel). It has only one input and two outputs. The second block (E2) of code runs, but comes out with an illogical answer. This block has two inputs that are dependent on the previous block's outputs. I'm expecting the result to be some sort of sine or cosine wave as this code is the position analysis of an ornithoptor.
"fsolve" in the in E2 also says "no solution found", but provides an answer. I do not understand why.
Why won't my code work properly for the E2 block?
% given values
R1 = 59.7;
R2 = 18.8;
R3 = 41.0;
R4 = 40.1;
R5 = 136.4;
R6 = 11.9;
R7 = 15.6;
R8 = 11.7;
Rgi = 13.1;
Rec = 12.7;
Rfc = 137.2;
theta1 = 0;
theta2_all = linspace(0, 2*pi, 500); % all possible input values
% E1
for i = 1:500
VLE1 = @(x0) [(R2.*cos(theta2_all(i)) + R3.*cos(x0(1)) - R4.*cos(x0(2)) - R1.*cos(theta1));
(R2.*sin(theta2_all(i)) + R3.*sin(x0(1)) - R4.*sin(x0(2)) - R1.*sin(theta1))];
x0 = [pi/2; pi/2];
x = fsolve(VLE1, x0);
theta3(i) = x(1);
theta4(i) = x(2);
end
theta3
theta4
% E2
theta_ec_off = 16.5*(pi/180);
theta_fc_off = 8.5*(pi/180);
theta_ec = theta3 + theta_ec_off;
theta_fc = theta4 + theta_fc_off;
for i = 1:500
VLE2 = @(x0) [(Rec*cos(theta_ec(i)) + Rfc*cos(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)));
(Rec*sin(theta_ec(i)) + Rfc*sin(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)))];
x0 = [pi/2; pi/2];
x = fsolve(VLE2, x0);
theta5(i) = x(1);
theta6(i) = x(2); % storing values in a vector
end
theta5
theta6
% validating
plot(theta2_all, theta5)
hold on
plot(theta2_all, theta6)

Accepted Answer

David Goodmanson
David Goodmanson on 9 Oct 2022
Edited: David Goodmanson on 9 Oct 2022
Hi Emily,
No one seems to have mentioned that there are two independent solutions for the angle pair angle3 and angle4, since if a side of S of the triangle S R3 R4 lies along, say, the x axis, then the triangle can lie either above or below the x axis, as determined by the angle between S and R3 being either positive or negative.
For each of those solutions there are two independent solutions for the angle pair angle5 and angle6, so there are four independent solutions in all for the set of angles 3 through 6. Of course only one of those sets will correspond to the real world geometry that you have.
Not having access to fsolve, the code below consists of algebraic vectorized solutions. The angles are stndard, positive angle being ccw from the x axis.
% VLE1 = @(x0) [(R2.*cos(theta2(i)) + R3.*cos(x0(1)) - R4.*cos(x0(2)) - R1.*cos(theta1));
% (R2.*sin(theta2(i)) + R3.*sin(x0(1)) - R4.*sin(x0(2)) - R1.*sin(theta1))];
% VLE2 = @(x0) [(Rec*cos(theta_ec(i)) + Rfc*cos(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)));
% (Rec*sin(theta_ec(i)) + Rfc*sin(theta_fc(i)) - R6*sin(x0(2)) - R5*sin(x0(1)))];
R1 = 59.7;
R2 = 18.8;
R3 = 41.0;
R4 = 40.1;
R5 = 136.4;
R6 = 11.9;
R7 = 15.6;
R8 = 11.7;
Rgi = 13.1;
Rec = 12.7;
Rfc = 137.2;
theta1 = 0;
theta2 = linspace(0, 2*pi, 500); % all possible input values
eps1 = [1 1 -1 -1];
eps2 = [1 -1 1 -1];
for n = 1:4 % loop over solutions
% E1
R7 = abs(R2*iexp(theta2)-R1*iexp(theta1));
theta7 = angle(R2*iexp(theta2)-R1*iexp(theta1));
u1 = eps1(n)*acos((R4^2-R7.^2-R3^2)./(2*R7*R3)); % pos or neg
theta3 = u1 + theta7;
theta4 = angle(R7 +R3*iexp(u1)) + theta7;
theta3 = unwrap(theta3);
theta4 = unwrap(theta4);
theta7 = unwrap(theta7);
% check complex version of VLE1, should be small
max(abs(R7.*iexp(theta7) + R3*iexp(theta3)-R4*iexp(theta4)))
% E2
theta_ec_off = 16.5*(pi/180);
theta_fc_off = 8.5*(pi/180);
theta_ec = theta3 + theta_ec_off;
theta_fc = theta4 + theta_fc_off;
R8 = abs(Rec*iexp(theta_ec) + Rfc*iexp(theta_fc));
theta8 = angle(Rec*iexp(theta_ec) + Rfc*iexp(theta_fc));
u2 = eps2(n)*acos((R8.^2+R6^2-R5^2)./(2*R8*R6)); % pos or neg
theta6 = u2 + theta8;
theta5 = angle(R8 -R6*iexp(u2)) + theta8;
theta5 = unwrap(theta5);
theta6 = unwrap(theta6);
theta8 = unwrap(theta8);
% check complex version of VLE2, should be small
max(abs(R8.*iexp(theta8) - R5*iexp(theta5) - R6*iexp(theta6)))
figure(1)
subplot(2,2,n)
plot(theta2,[theta3;theta4;theta5;theta6])
grid on
end

More Answers (2)

David Goodmanson
David Goodmanson on 1 Oct 2022
Edited: David Goodmanson on 6 Oct 2022
Hi Emily,
If VLE1 and VLE2 are lined up typographically, you have
VLE1 = @(x0) [(R2.*cos(theta2_all(i)) + R3.*cos(x0(1)) - R4.*cos(x0(2)) - R1.*cos(theta1));
(R2.*sin(theta2_all(i)) + R3.*sin(x0(1)) - R4.*sin(x0(2)) - R1.*sin(theta1))];
VLE2 = @(x0) [(Rec*cos(theta_ec(i)) + Rfc*cos(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)));
(Rec*sin(theta_ec(i)) + Rfc*sin(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)))];
Looking at just terms involving R3 and R4, VLE1 has cosines on the first line, sines on the second line. But for R5 and R6, VLE2 has cosines on both lines, meaning that the first line and the second line have identical terms. This does not seem right.
  1 Comment
Emily Friedman
Emily Friedman on 1 Oct 2022
Whoops! That was a typo. Those cos on the second line for VLE2 are supposed to be sin. The result from VLE2 with the correct equation is still illogical though. It makes a very, very weird plot.

Sign in to comment.


Torsten
Torsten on 6 Oct 2022
Edited: Torsten on 6 Oct 2022
% given values
R1 = 59.7;
R2 = 18.8;
R3 = 41.0;
R4 = 40.1;
R5 = 136.4;
R6 = 11.9;
R7 = 15.6;
R8 = 11.7;
Rgi = 13.1;
Rec = 12.7;
Rfc = 137.2;
theta1 = 0;
theta2_all = linspace(0, 2*pi, 500); % all possible input values
options = optimset('Display','off') ;
x0 = [pi/2 pi/2];
% E1
for i = 1:500
VLE1 = @(x0) [(R2.*cos(theta2_all(i)) + R3.*cos(x0(1)) - R4.*cos(x0(2)) - R1.*cos(theta1));
(R2.*sin(theta2_all(i)) + R3.*sin(x0(1)) - R4.*sin(x0(2)) - R1.*sin(theta1))];
x = fsolve(VLE1, x0, options);
x0 = x;
theta3(i) = x(1);
theta4(i) = x(2);
end
figure(1)
hold on
plot(theta2_all, theta3)
plot(theta2_all, theta4)
hold off
% E2
theta_ec_off = 16.5*(pi/180);
theta_fc_off = 8.5*(pi/180);
theta_ec = theta3 + theta_ec_off;
theta_fc = theta4 + theta_fc_off;
x0 = [pi/2 pi/2];
for i = 1:500
VLE2 = @(x0) [(Rec*cos(theta_ec(i)) + Rfc*cos(theta_fc(i)) - R6*cos(x0(2)) - R5*cos(x0(1)));
(Rec*sin(theta_ec(i)) + Rfc*sin(theta_fc(i)) - R6*sin(x0(2)) - R5*sin(x0(1)))];
x = fsolve(VLE2, x0, options);
x0 = x;
theta5(i) = x(1);
theta6(i) = x(2); % storing values in a vector
end
figure(2)
hold on
plot(theta2_all, theta5)
plot(theta2_all, theta6)
hold off

Categories

Find more on Simulink Functions in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!