How to store all of the converged roots found by Matlab's fsolve algorithm?
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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,
0 commentaires
Réponse acceptée
Walter Roberson
le 3 Oct 2020
Initialize
converged_roots = [];
and replace
converged_roots = my_root;
with
converged_roots{end+1} = my_root;
17 commentaires
Walter Roberson
le 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
le 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.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Solver Outputs and Iterative Display dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!