How to store all of the converged roots found by Matlab's fsolve algorithm?
3 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Walter Roberson
2020-10-3
Initialize
converged_roots = [];
and replace
converged_roots = my_root;
with
converged_roots{end+1} = my_root;
17 个评论
Walter Roberson
2020-10-4
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
2020-10-4
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!