MATLAB Answers

How to store all of the converged roots found by Matlab's fsolve algorithm?

1 view (last 30 days)
I'm practicing with Matlab's fsolve function, and am calling it with a few nested for loops that provide fsolve with various initial guesses.
How can I store all of the converged roots? Should I first store all of the roots, and then use another variable to store only the converged roots that satisfy the function tolerance?
For instance, I can write
% converged_roots = zeros(2, nguesses^2) % Pre-allocate a 2 x nguesses^2 matrix of zeros to store the converged roots
nguesses = 20
for x_guess = linspace(-10, 10, nguesses)
for y_guess = linspace(-5, 5, nguesses)
initial_guess = [ x_guess, y_guess ];
[ my_root, fval, exit_flag, output, Jacobian ] = fsolve( myfunction, initial_guess )
if norm( fval ) < 1e-5
converged_roots = my_root;
disp('not a converged root')
But, if I have, say, 5 converged roots found, the matrix 'converged_roots' only stores the last converged root -- so it's being overwritten each time in the loop.
How could I fix this?
If I pre-allocate a matrix of zeros (first line of my code, commented out), nothing happens either -- 'converged_roots' will still only return the last converged root. I think it might be because I'm not using a loop index i, j, anywhere in the code ...


Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 3 Oct 2020
converged_roots = [];
and replace
converged_roots = my_root;
converged_roots{end+1} = my_root;


Noob on 4 Oct 2020
Hi Walter,
I just saw your last comment, after I submitted an additional comment up there ^.
So if x^2 + 1e-50 is indistinguishable from x^2, then, numerically, x = 0 is considered a root of the function.
Is that correct to say?
And are these misleading "roots" still valuable to numerical analysts, say, from an optimization point of view?
Although they aren't true roots analytically, they would still give function evaluations that might satisfy the tolerances that are set.
Walter Roberson
Walter Roberson on 4 Oct 2020
x^2 +1e-50 is indistinguishable from x^2 for abs(x) > 1e-17 . But there are a lot of values between 0 and 1e-17 in double precision, so if fsolve happens to look in the area at all, it could potentially have some meaning.
The question is whether it will look in the area at all. For example at f = eps then f(x) is about 1e-31, so fsolve might have already decided that it had a good enough root. If the function tolerance were 1e-13 then the range +/- 0.00000007 satisfies that tolerance.
Thus, we should not be surprised if fsolve() decides that there is a root anywhere in +/- 0.00000007 unless we adjust the tolerances.
But no matter what absolute tolerance we use, fsolve() is going to have problems.
x^2 + sqrt(eps(realmin))
does not become distinguishable from x^2 until about 1e-72.
Walter Roberson
Walter Roberson on 4 Oct 2020
We could hypothesize that there might be different algorithms for looking at the tolerance of a root determined by jacobian (root has an even multiplicity so function does not cross 0 there), versus a root determined by looking for a zero crossing. But then we just substitute the function (x^2 - 1e-50) * (x-20) which does have 0 crossings (true roots) at +/- 1e-25 and ask ourselves whether we can realistically expect a solver that works by looking for zero crossings to be able to find such a root. Sure, a solver might deliberately try exactly 0, but we can obviously shift the root slightly ((x-epsilon)^2 - 1e-50)*(x-20) and have the same problem.

Sign in to comment.

More Answers (0)





Community Treasure Hunt

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

Start Hunting!