Please read to the end here, because there is a lot of ground that I cover in explaining the real issue with your problem.
1. Fsolve is an optimizer, specifically, a rootfinder, a numerical optimization tool.
2. Numerical solvers find only ONE solution, even if multiple solutions exist, and that is if they do converge at all. So at MOST one solution.
3. The solution that is found depends on where you start the solver. That is called the basin of attraction, the locus of start points that will converge to any given solution for a given solver.
So if you start in the wrong place, then you will end where you will end.
Finally, problems with trig functions, roots, powers, etc., in them tend to have multiple solutions. While that appears not to be the case here, it is often true. This problem is just one of poor starting values. Ok, they don't seem poor. But in fact, they are poor, but in a subtle way.
Since you don't actually tell us what solution you expected, or if indeed it is truly a solution, it is impossible to know why you think you got the "wrong" answer. But keep reading, as I finally get to an interesting conclusion. You do get the correct answer here, just not what you expected.
del1=0.0710; del2=0.0897;
F1 = @(x1,x2) x2-600./(2.*pi.*x1).*asin(sqrt(del1./(1+4.*(x1-1).^2./(x1+1).^2.*(1-del1)./(1-(x1-1).^2./(x1+1).^2).^2)));
F2 = @(x1,x2) x2-540./(2.*pi.*x1).*asin(sqrt(del2./(1+4.*(x1-1).^2./(x1+1).^2.*(1-del2)./(1-(x1-1).^2./(x1+1).^2).^2)));
ezsurf(F1,[0 20, 0 20])
figure
ezsurf(F2,[0 20, 0 20])
Now, you apparently expected to see a solution in the vicinity of x0.
x0=[1.33,18.6];
First of all, PLOT these functions. ALWAYS PLOT EVERYTHING.
ezsurf(F1,[0 20, 0 20])
figure
ezsurf(F2,[0 20, 0 20])
So, two very similar looking surfaces. That is to be expected, since the forms of those functions are nearly identical, except for some subtle changes in a few constants. Lets see what fsolve found.
F12 = @(x12) [F1(x12(1),x12(2)),F2(x12(1),x12(2))];
[x12final,ffinal] = fsolve(F12,x0)
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 200 (the default value).
x12final =
11.4875348892883 0.400563824853909
ffinal =
0.00377412920813197 -0.00484216991838443
So fsolve thinks it has gone in the correct direction for a solution. Certainly it has decreased the function value at the end result.
So where does the truth lie? Again, we need to look more deeply.
[X1,X2] = meshgrid(linspace(0.1,20,500));
F1surf = F1(X1,X2);
F2surf = F2(X1,X2);
A nice thing about a contour plot, is you can view it as a root finder, when used properly.
F1locus = contour(X1,X2,F1surf,[0 0],'r');
F2locus = contour(X1,X2,F2surf,[0 0],'r');
[F1locus(:,2:10);F2locus(:,2:10)]'
ans =
1.2565 20 1.2724 20
1.2565 19.999 1.2745 19.96
1.2586 19.96 1.2766 19.92
1.2607 19.92 1.2787 19.88
1.2628 19.88 1.2808 19.84
1.2649 19.84 1.2829 19.801
1.267 19.801 1.285 19.761
1.2691 19.761 1.287 19.721
1.2712 19.721 1.2891 19.681
So columns 1 and 2 are the set of points that contour found, that solves for where F1 is zero. Likewise, columns 3 and 4 are the same result, but for function F2.
See how close they are, as curves. In fact, if you plotted those two curves on top of each other, they would virtually overlap!
But in fact, they are very different sets of numbers as far as F1 and F2 are concerned.
contour(X1,X2,F1surf,[0 0],'r')
hold on
contour(X1,X2,F2surf,[0 0],'g')
As you can see, the red and the green curves are almost indistinguishable.
Lets see what happens when I push those numbers through F1 and F2.
F1(F1locus(1,2:10),F1locus(2,2:10))
ans =
2.0293e-05 0 0.00086698 0.0016545 0.0023419 0.0029296 0.0034176 0.0038063 0.0040958
F2(F1locus(1,2:10),F1locus(2,2:10))
ans =
-0.30571 -0.30571 -0.30433 -0.30301 -0.3018 -0.30069 -0.29968 -0.29877 -0.29796
So the points from F1locus give pretty small numbers for F1, but relatively large numbers for F2.
F1(F2locus(1,2:10),F2locus(2,2:10))
ans =
0.30604 0.30566 0.30519 0.30462 0.30396 0.3032 0.30234 0.30139 0.30035
F2(F2locus(1,2:10),F2locus(2,2:10))
ans =
0.0042814 0.0044167 0.0044544 0.0043946 0.0042377 0.0039839 0.0036332 0.0031861 0.0026426
Conversely, the F2locus points do the opposite. So even though these two solution loci are extremely close, in fact, they are very different in what they do to F1 and F2.
Where did you start the solver? Since hold is still on, just add the start and end points.
plot(x0(1),x0(2),'bo')
plot(x12final(1),x12final(2),'bs')
So the o is the start, the square the end point for fsolve.
Conclusion? Your starting point SEEMED ok to you. Hey it seems not too far off!
In fact however, your start point was pure crap. Fsolve had to walk all along those two curves to find what it deems abetter solution. Even the final solution is not terribly great, and fsolve recognizes that, based on the message where fsolve thinks it has not converged. I see that you extended the default option there. That never really helps. It just lets fsolve walk a little further along those curves.
Here is what fsolve told us:
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 200 (the default value).
So eventually, fsolve gives up, because those two curves are growing both infinitesimally close and virtually parallel, but always subtly apart! In that region, your functions have a very large slope, so tiny changes in x1 and x2 will make a large difference.
If I now try it again, but this time letting fsolve go as far as you did, we get this answer:
options = optimset('MaxFunEvals',2000);
[x12final,ffinal] = fsolve(F12,x0,options)
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
x12final =
215.57 0.0011482
ffinal =
1.2105e-05 -1.2821e-05
fsolve is now happier, because I let it chase all down that pair of curves, until finally both functions give zero to within the required tolerances. Had I set a tighter tolerance, and given fsolve more time to converge, it would have just walked much further along those curves.
So, no, you only thought you had good starting values.
And if I had to guess, there os only one true solution to this system of equations, at x1=inf, x2=0. That is easy enough to test out.
F12([1e6,0])
ans =
-5.2799e-11 -5.3957e-11
F12([1e7,0])
ans =
-5.2799e-13 -5.3957e-13
F12([1e8,0])
ans =
-5.2799e-15 -5.3957e-15
MATLAB will give up on me if I go all the way.
F12([inf,0])
ans =
NaN NaN
So I can make both functions as small as I wish, simply by choosing a point closer and closer to [inf,0].
