Solving these two equations using Fsolve is giving the wrong answer don't understand why ?

if true
function F = myfun(x)
global Xpos
l1 = 0.20000000;
l2 = 0.10000000;
eq1 = l1*sin(x(1,1)) + l2*sin(x(1,2))-0.2;
eq2 = l1*cos(x(1,1)) + l2*cos(x(1,2)) - Ypos;
F=[eq1;eq2];
end
and then..........
options=optimset('Display', 'off', 'TolFun', 1.0e-54, 'TolX',1.0e-54);
       theta(i,1:2)= fsolve(@myfun,[1 2.8],options);

I need a mass to follow a specific path dependent using two links which are defined by theta1 and theta 2. My above code is trying to get a point to follow a specific path i.e. keep x set as 0.2m hence the -0.2 in eq1 and Ypos which is a global variable which changes its value depending on the change in the y direction, which changes with each time step, hence Y updates itself. However it is yielding some bad values which are quite a way off. I was wondering whether anyone could help?

[John D'Errico]

I'm surprised it works at all. I'll just copy the two relevant lines...

global Xpos
eq2 = l1*cos(x(1,1)) + l2*cos(x(1,2)) - Ypos;

So, you pass in Xpos, then use Ypos. I bet that code fails, telling you that Ypos is undefined. :)

Answer 1

[Star Strider]

I would define your function as:

function F = myfun(x, XPos, YPos)

with the fsolve call changed to:

theta(i,1:2)= fsolve(@(x) myfun(x, XPos, YPos), [1 2.8], options);

and delete the global call. (In general global variables are to be avoided. They make code more difficult to debug, and since anything can change them, can cause unpredictable behaviour.)

You do not use ‘XPos’ in your code, and since your function is a function file and not an anonymous function, it will not pick up ‘YPos’ from your workspace. You have to pass your parameters as variables. Fortunately, this is easy to do with the anonymous function I used in my revision of your fsolve call.

I cannot run your code, but these changes should work.

Answer 2

[John D'Errico]

Assuming that the confusion between Xpos and Ypos is resolved, AND that your problem is not due to poor use of Global variables (a terrible way to pass things around, because it makes for buggy code.) Then...

Since sin and cos are periodic functions, you need to recognize that there are infinitely many solutions. fsolve will return one of them, but perhaps not the solution you want. So my claim is that fsolve is not giving you the wrong answer, just not the answer that you want to see.

By the way, trying to set the tolerance as low as 1e-54 is an insane waste of time. You will NEVER get that accurate of a solution, so you will probably force fsolve to continue trying until it gives up due to to many iterations, or some similar reason.

Regardless, lets look at how I might resolve the "problem".

syms u1 u2
l1 = 0.20000000;
l2 = 0.10000000;
Ypos = 0.03; % chosen arbitrarily as an example.

eq1 = l1*sin(u1) + l2*sin(u2) - 0.2 == 0;
eq2 = l1*cos(u1) + l2*cos(u2) - Ypos == 0;
res = solve(eq1,eq2);

vpa(res.u1)
ans =
  1.9242743377708882328399117968457
 0.91953842059991050444967080312206

vpa(res.u2)
ans =
 0.12396842138987707464108063974351
  2.7198443369809216626485019602243

So, there are at least TWO solutions. We can verify that as fact.

simplify(subs(eq1,{u1,u2},{res.u1,res.u2}))
ans =
   TRUE
   TRUE

And, we can add any integer multiple of 2*pi to that either u1 or u2 or both, and get an equally valid pair of solutions.

simplify(subs(eq1,{u1,u2},{res.u1 + 2*pi,res.u2}))
ans =
   TRUE
   TRUE

If I had to guess at what you are doing, is you are getting a different (EQUALLY VALID!!!!!) solution from fsolve from what you expect. But since neither you nor fsolve can control which solution it returns, you see a "problem".

Perhaps using solve is a better choice here. At least it is capable of returning both primary solutions, and you can choose the one that you prefer. And you can always recognize that multiples of 2*pi can be added or subtracted.

Finally, I could PROBABLY come up with a scheme to find the OTHER primary solution from that which fsolve found. But is it really worth the effort? Just use solve.

Oh, what the heck. Let me do some pencil and paper work.

(sin(u1))^2 = (0.2 - l2*sin(u2))^2 / l1^2
(cos(u1))^2 = (Ypos - l2*cos(u2))^2 / l1^2

Add the two together, then recognize a very classical identity for trig functions.

1 = (0.2 - l2*sin(u2))^2 / l1^2 + (Ypos - l2*cos(u2))^2 / l1^2

So, this is now a problem in only ONE unknown, u2. Just a bit more playing around...

l1^2 = (0.2 - l2*sin(u2))^2 + (Ypos - l2*cos(u2))^2

Now, expand those squares, and remember that nice identity again...

sin(u2)^2 + cos(u2)^2 = 1

so...

l2^2 - sin(u2))*24/5 + Ypos^2 - l1^2 - 24*Ypos*cos(u2) + 0.04 == 0

In essence, we have a simple equation of the form

a*sin(u2) + b*cos(u2) + c == 0

This will generally have two solutions in the interval [0,2*pi], and then we can add multiples of 2*pi to those solutions.

Answer 3

[Alan Weiss]

Your tolerances don't make much sense. See the documentation, where you find that setting any value less than eps ~ 2e-16 is the same as setting it to zero, which might cause convergence issues.

I don't see a global Ypos statement in your function, but maybe that is just a typo.

It is possible that a solution does not exist locally, and the "bad values" you talk about might indicate that the function itself has to jump to reach a solution. But I don't know, partly because I am not sure what you mean by "bad values."

One minor point: you call x(1,1) and x(1,2) within your objective function. It is less error-prone to call x(1) and x(2); that way, if a vector has an unexpected orientation, you won't get an error.

Alan Weiss

MATLAB mathematical toolbox documentation