fsolve yields wrong answer

Hello

It is an example of many equations in my lengthy code. I want to use fsolve to solve it as it is more efficient; however, it yields wrong answer:

F = @(x)((x-20.6667)^2+57.89);
opts = optimoptions('fsolve', 'Display', 'off');
    x(1) = fsolve(@(x) F(x), 0, opts);

in this case it yields x(1)=20.6667 ; howver, it should not yield any solution in real domain!!!!I want to get real domain roots by fsolve.

What should i do? I dont want to use solve as it results in longer run times.

Bests

Pooneh

Answer 1

[Star Strider]

Give fsolve a complex initial estimate to get a complex root —

F = @(x)((x-20.6667).^2+57.89);
opts = optimoptions('fsolve', 'Display', 'off');
x0 = 1 + 1i;
x(1) = fsolve(@(x) F(x), x0, opts)

x = 20.6667 + 7.6085i

.

[Pooneh Shah Malekpoor]

Is there any other way? I really dont want complex numbers in my code as I have to change many things...is there anyway to get only real roots from fsolve?

[Walter Roberson]

For real-valued x, then (x-20.6667) is real-valued, and (x-20.6667).^2 is always non-negative (could be 0).

57.89 is always positive.

The sum of a non-negative number (x-20.6667).^2 and a positive number is always positive.

fsolve() tries to find values for which the given expression is 0. But we just demonstrated that for any real value x, the expression is always positive. Therefore no real roots exist for the function.

If you look at Star Strider's answer, you might notice that the real part is the same as the 20.6667 being subtracted, so you get a 0 real contribution from x-20.6667 . And the imaginary part + 7.6085i is sqrt(-57.89) so 7.6085i ^2 is -57.89 and so when you add +57.89 to that you get a total of 0. So 20.6667 + 7.6085i is a correct solution, and no real-only solutions exist.

[Pooneh Shah Malekpoor]

Yes, you are right and it is clear...let me ask how is it possible that fsolve understands that this equation does not have a real root? as you already mentioned this equation sums two positive values!

[Walter Roberson]

format long g
F = @(x)((x-20.6667).^2+57.89);
opts = optimoptions('fsolve', 'Display', 'off');
[x, fval, errorflag, output] = fsolve(@(x) F(x), 0, opts);
x

x = 
          20.6666998416185

fval

fval = 
                     57.89

errorflag

errorflag = 
    -2

output

output = struct with fields:
       iterations: 29
        funcCount: 48
        algorithm: 'trust-region-dogleg'
    firstorderopt: 1.33568026341535e-06
          message: '↵No solution found.↵↵fsolve stopped because the last step was ineffective. However, the vector of function↵values is not near zero, as measured by the value of the function tolerance. ↵↵<stopping criteria details>↵↵fsolve stopped because the sum of squared function values, r, changed by 8.814899e-15↵relative to its initial value; this is less than max(options.FunctionTolerance^2,eps) = 1.000000e-12.↵However, r = 3.351252e+03, exceeds sqrt(options.FunctionTolerance) = 1.000000e-03.↵↵'

Notice the errorflag is negative, indicating a failure. Notice the output message says that the last steps it tried did not do any good, but the function value was not near 0.

If the question is how to find out whether fsolve() solved the problem, then look at errorflag

If the question is how fsolve() figures out that it is not going to be able to solve the problem, then the answer is approximately that it did numeric gradient estimation to get to a place where it could see that the gradient was 0 but that the function was not 0.

[Star Strider]

The fsolve function will find a complex root if given a complex initial estimate, or if a purely imaginary root is the only root.

One option to eliminate complex numbers could be to take their absolute values —

x = 20.6667 + 7.6085i;
xa = abs(x)

xa = 22.0228

and the other option of course is just to use the real part.

However changing the code to accommodate complex numbers may be the best option, because it will likely provide the most accurate results.

.

[Matt J]

Give fsolve a complex initial estimate to get a complex root —

But note that this will only lead to a proper complex-domain solution under certain conditions:

https://www.mathworks.com/help/optim/ug/complex-numbers-in-optimization-toolbox-solvers.html

Answer 2

[Matt J]

I want to use fsolve to solve it as it is more efficient; however, it yields wrong answer:

In your example, the function is a polynomial. Is this always the case for you? If so, it is more efficient to use roots(). Moreover, roots() will find all solutions, both real and complex, whereas fsolve can find only one solution.

[Pooneh Shah Malekpoor]

Thank you so much Matt J! your comment made a huge positive impact on my code and its run time has decreased significantly. This is highly appreciated!

Answer 3

[Matt J]

parts=@(z) [real(z);imag(z)];
complx=@(x) complex(x(1),x(2));

F = @(x)parts((complx(x)-20.6667)^2+57.89);

opts = optimoptions('fsolve', 'Display', 'off');
   
x = complx( fsolve(F, [0,0], opts))

x = 20.6667 + 7.6085i