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

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)

Answer 1

[David Goodmanson]

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

Answer 2

[David Goodmanson]

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.

[Emily Friedman]

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.

Answer 3

[Torsten]

% 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