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

조회 수: 4(최근 30일)
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)

채택된 답변

David Goodmanson
David Goodmanson 2022년 10월 9일
편집: David Goodmanson 2022년 10월 9일
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

추가 답변(2개)

David Goodmanson
David Goodmanson 2022년 10월 1일
편집: David Goodmanson 2022년 10월 6일
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
Emily Friedman
Emily Friedman 2022년 10월 1일
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.

댓글을 달려면 로그인하십시오.


Torsten
Torsten 2022년 10월 6일
편집: Torsten 2022년 10월 6일
% 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

범주

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!

Translated by