In nested for-loops, how and where should I set the counting index (indices?) correctly, for doing numerical root-finding using fsolve?

Hi,

Let's say I have a function, f, that maps R^3 to R^3, and I want to find its multivariable roots numerically, using fsolve.

What's a good way to use the counting index, i, in the loops?

For instance, if f(x,y,z) = ( f_1(x,y,z), f_2(x,y,z), f_3(x,y,z) ), and I want to solve f = (0,0,0), then I could write nested for-loops that do something like this:

for x_guess = linspace(-10,10,10)   % first, fix a guess for x
for y_guess = linspace(-10,10,10)   % now for a fixed x, also fix a guess for y 
for z_guess = linspace(-10,10,10)   % now for a fixed x and y, check all values of z, before fixing a new x and y
    
    i = ... 
        
end

end

end

[ multivariable_roots, FVAL, EXITFLAG, OUTPUT, JACOB ] = fsolve( f, [x_guess, y_guess, z_guess] )

Where's a good place to put the counting index? At the innermost loop? Could I do so at the outermost loop?

My mentor showed me a code that uses the counting index in the innermost loop (the code structure that I provided above), but I currently don't understand it, so I'd like to write my own nested for-loops in order to understand what it does.

Also, would using more than one counting index be helpful? I'm thinking of using counting indices i, j, k, if it makes the nested for-loops easier to understand. For instance, in a double summation, we would typically fix i and sum through j -- and then fix another i and sum through j again, etc. Could I do something analogous in nested for-loops? I think I might actually prefer using a different counting index for each loop, since things might look more explicit that way.

Thanks,

[Matt J]

Could I do something analogous in nested for-loops?

Every for-loop has a loop index (otherwise, what action is it performing), so if you are going to use 3 loops, you will inevitably have 3 loop indices.

[Noob]

Hi Matt,

Yes, I understand what you're getting at, which is the reason why I want to write my own nested for-loops. However, the code I have from my mentor only has one loop index -- at the innermost loop. The code runs fine too, and the solutions seem valid. I've edited my question to show you the code structure of what I mean.

Thanks,

Answer 1

[Matt J]

Well, first of all, you would not use fsolve when doing a gridded root search. In general, the search could be done like this:

[x_guess,y_guess,z_guess]  = ndgrid(linspace(-10,10,10));

fvals=nan(size(x_guess));

for i=1:numel(x_guess)
    
    fvals(i) = f( [x_guess(i), y_guess(i), z_guess(i)] );
        
end

iroot=abs(fvals)<=something_small;

solutions=[x_guess(iroot),y_guess(iroot),z_guess(iroot)];

Note however that you would only ever resort to a for-loop if your objective function f() is absolutely not vectorizable. For a vectorizable function like f(x,y,z)=x+y+z, the above could be simplified to something that calls f() only once (and hence is much more efficient):

[x_guess,y_guess,z_guess]  = ndgrid(linspace(-10,10,10));

fvals=f(x_guess, y_guess, z_guess);

iroot=abs(fvals)<=something_small;

solutions=[x_guess(iroot),y_guess(iroot),z_guess(iroot)];

Once you have done one of the above search schemes, and assuming the solution space contains only a finite set of isolated points, you could then loop through solutions(i), and try using fsolve to calculate more exact roots. Note, however, that there is generally no systematic way of choosing initial guesses for fsolve that guarantees it will converge to a root you haven't already found previously. There is no numeric method that will guarantee that you find all roots to a general function.

[Noob]

Hi Matt,

Thanks so much for your answer. Why would we then loop through solutions(i) and use fsolve to calculuate more exact roots? For instance, couldn't I simply set the value, something_small, in the variable, iroot, to be as small as the function tolerance in fsolve? Or is something gained by using fsolve?

Is the geometric aspects of fsolve the reason why it may find roots that are different from the gridded root search? That would be interesting. I wrote two simple solvers for gradient descent and Newton's second-order method; and, fixing the same initial guess for each solver, they may converge to different roots -- with gradient descent sometimes zig-zagging a great deal towards the minimizer and requiring many more iterations to converge.

Thanks,

[Matt J]

If you're happy with the solutions you get from the grid search, you may elect not to refine them. However, the precision of the solutions you get from the grid search are limited by the fineness of the grid. You could make the grid super-fine, but then you would have an exponentially increasing number of points to search. You also don't know in advance how fine a grid you would need to achieve a certain tolerance on f(), if that's what matters to you.

fsolve, by contrast, uses local gradients and Jacobians to narrow in on solutions a lot more efficiently than a grid search. fsolve is, on the other hand, more vulnerable to getting stuck at local solutions than is a grid search and needs to satisfy regularity conditions like differentiability.

[Noob]

Ah, interesting. I'll think more about this over dinner now.

Thanks again!

[Noob]

Hi Matt,

What if the function has a parameter c (a constant value) that I would also like to vary -- and time T?

Could I expand the n-dimensional grid?

The range of guesses would be much smaller for c though, say, .5 to 1.5. I could also use linspace for c, but that seems very inefficient to do.

And would it make sense to loop through time T, or manually change it but keep it fixed, like, say, T = 20 seconds?

Thanks,

[Matt J]

Yes, ndgrid() will let you make grid sample arrays in any dimension, assuming you have adequate memory to hold them. However, if your functions are vectorizable, it may be more memory-efficient to use ndgridVecs(),

https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs