why does this trig iteration generates complex numbers?

7 次查看（过去 30 天）
Mick
Mick 2014-9-9
评论： Roger Stafford 2014-9-10
im trying to write a program that analyses angles in a four bar mechanism given by three members d1&d2 (a kinked bar- kink at angle beta), d3 and d4. d1 and d4 are hinged at the origin and a point (h,v) respectively, where v is below the origin. The driving angle between the horiz and d1 is alpha, the angles between the horiz and d3, and the horiz and d4 are phi and theta respectively. I want to output theta, but i need to iterate to find theta and phi for a given set of d's. currently, my program is outputting complex values of theta and i cant figure out why... my suspicion is that the iteration doesnt work. THis is the code:
d1 = 2;
d2 = 0.5;
v = 2;
h = 5;
beta = pi()/4;
%set minimum length with flap and hinge to hinge lengths
h2h_len = sqrt(v*v + h*h)
flap_len = sqrt(d1*d1 + d2*d2)
min_len = h2h_len + flap_len
d3 = 1;
%loop for later calculations...
% while (d3 < min_len + 3)
d4 = min_len + 2 - d3;
alpha = pi()/4;
%while (alpha < pi()/2)
Ah = -d1*cos(alpha) + d2*cos(alpha - beta);
Av = d1*sin(alpha) - d2*sin(alpha - beta);
theta = pi()/4;
%ITERATION SECTION STARTS HERE
count = 1;
while (count <= 50)
theta = asin((sqrt(d3*d3 - ((h-Ah)-d4*cos(theta))^2) + v + Av)/d4);
count = count + 1;
end
theta
phi = acos(((h - Ah) - d4*cos(theta))/d3)
%alpha = alpha + 0.0873;
%end
% d3 = d3 + 1;
%end
end
  1 个评论
José-Luis
José-Luis 2014-9-9
For real elements of x outside the range [-1,1], asin(X) is complex.
Therefore your function might evaluation to one or more such values.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Roger Stafford
Roger Stafford 2014-9-10
You are apparently trying to solve the trigonometric equation
sin(theta) = (sqrt(d3*d3 - ((h-Ah)-d4*cos(theta))^2) + v + Av)/d4
but the iterative method you are using is definitely not the way to accomplish this. You have no reason to believe that your iteration would ever converge to a solution even if you started with a close approximation. Other techniques are certainly called for here. By the appropriate algebraic manipulation that equation can be put in this form:
d4*sin(theta)-v-Av = sqrt(d3^2-(h-Ah-d4*cos(theta))^2)
By squaring both sides and further manipulation it becomes
2*(v+Av)*d4*sin(theta)+2*((h-Ah)*d4*cos(theta)
= d4^2-d3^2+(v+Av)^2+(h-Ah)^2
This equation is of the form
A*sin(theta)+B*cos(theta) = C
where A = 2*(v+Av)*d4, B = 2*((h-Ah)*d4, and C = d4^2-d3^2+(v+Av)^2+(h-Ah)^2. Its solution can be found using matlab's 'atan2' function. If we evaluate
eta = atan2(A,B)
it will satisfy sin(eta) = A/sqrt(A^2+B^2) and cos(eta) = B/sqrt(A^2+B^2) so that
cos(theta-eta) = cos(theta)*cos(eta)+sin(theta)*sin(eta)
= (B*cos(theta)+A*sin(theta))/sqrt(A^2+B^2)
= C/sqrt(A^2+B^2)
This finally allows a solution
theta = eta + acos(C/sqrt(A^2+B^2)) = atan2(A,B) + acos(C/sqrt(A^2+B^2))
Of course you can also use iteration with the 'fzero' function, but that is iteration of a very different kind than you have attempted here.
  1 个评论
Roger Stafford
Roger Stafford 2014-9-10
I should have added that
theta = atan2(A,B) - acos(C/sqrt(A^2+B^2))
is also a solution, as well as with the addition or subtraction of any multiple of 2*pi.
Also it should be noted that while all of these are certainly solutions to the equation
2*(v+Av)*d4*sin(theta)+2*((h-Ah)*d4*cos(theta)
= d4^2-d3^2+(v+Av)^2+(h-Ah)^2
they may not all be solutions to
sin(theta) = (sqrt(d3*d3-((h-Ah)-d4*cos(theta))^2)+v+Av)/d4
but some would instead satisfy
sin(theta) = (-sqrt(d3*d3-((h-Ah)-d4*cos(theta))^2)+v+Av)/d4

请先登录，再进行评论。

the cyclist
the cyclist 2014-9-9
You could set a condition breakpoint that checks for the first time theta goes complex, and trace back through the value that go into the calculation, and what went awry. Maybe you have a tiny bit of roundoff error that takes you out of the expected input range.
See this page for details on debugging, setting breakpoints, etc, if you don't know how to do that.

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by