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

2 views (last 30 days)
Hi,
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;
else
disp('not a converged root')
end
end
end
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 ...
Thanks,

Accepted Answer

Walter Roberson
Walter Roberson on 3 Oct 2020
Initialize
converged_roots = [];
and replace
converged_roots = my_root;
with
converged_roots{end+1} = my_root;
  17 Comments
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)

Products


Release

R2017a

Community Treasure Hunt

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

Start Hunting!