How to find all solutions to a system of two nonlinear equations?

I have a system of two nonlinear equations (f(x,y)=0 and g(x,y)=0) to which I want to find all roots over a region (say x from -5 to 5 and y from -5 to 5).

The problem with using fsolve is that I need to supply the initial guesses which may not be known easily. If the region where I am looking for solutions is small, I guess I could take a lot of initial guesses and get the solutions. However, I want to know if a more robust method exists.

One approach is tried was to find the values of (x,y) where f(x,y) and g(x,y) simultaneously change their sign and feeding those as the initial conditions. This approach works reasonably well until you supply equations which don't change sign after a zero crossing. (For example f(x,y) = x^2)

Any help/suggestion is appreciated!

Thanks,

Mohit

Answer 1

[John D'Errico]

Here is the rub - you are in trouble IF you really want to find ALL solutions of a system of totally general nonlinear equations, regardless of the domain.

First, any problem can have infinitely many solutions, even in a finite interval. Consider sin(1/x), for example, with infinitely many roots in any finite interval that contains zero. And while you can claim those solutions are describable analytically, it is easy enough to create a problem with roots that are not so easily describable. So finding all roots of any problem is therefore impossible.

Do I really need to cook up an example that has no simple analytical description, yet has infinitely many roots? Quickly, this should work:

fplot(@(x) x/10 + besselj(1,1./x),[0,1])

Pick any interval that contains zero, even an interval that does not explicitly contain zero, if the interval is something like (0,eps]. Thus infinitely many solutions, and since x appears both inside and outside the bessel function, I'll conjecture it won't have an analytical solution.

Next, you talk about the use of fsolve as being a problem, because tools like that provide only a single solution as a function of the start point. Change the start, and potentially get a different result. For this, you may want to do some reading about global optimization methods, and how they differ from more classical ones.

But even there, a global method cannot find all roots. As we just said before, since there can trivially be infinitely many roots, it must give up eventually. And then which roots did it miss? The issue is if you can bound the derivatives of your function. If no bounds are available, then you can provably never find all roots using a numerical method and finite precision arithmetic, since it is trivially easy to create a function that is constant everywhere, except for some immeasurably small interval where it bounces across zero.

In the end, that leaves you at best with probabilistic measures of whether you may have found any given solution, or all of the solutions, since you generally will not know exactly how many solutions truly exist. A good idea may be to use multi-start methods, then use clustering schemes to decide how many solutions you have actually found. Remember that after each solve, two solutions that may in theory be the same will typically only be as close as the convergence tolerance. And that will be impacted by the local shape of the objectives.

Again, do some reading. But don't expect to find a solver that will magically yield all of infinitely many solutions, or even a solver that will robustly always find all of a finite number of solutions on any general problem.

[Mohit Kumar]

Thanks for the answer! It was a nice read.

"But don't expect to find a solver that will magically yield all of infinitely many solutions, or even a solver that will robustly always find all of a finite number of solutions on any general problem."

This is definitely something I wanted to know. I didn't want to code a rough solution and end up finding that a shiny solution already exists, which has (unsurprisingly) happened more than a few times.

For all my needs and purposes, I can use the method mentioned in Ameer Hamza's answer or simply take several initial guesses or use a combination of both followed by checking it graphically if I obtain unexpected results. And for learning's sake, I'll read up on the topics you mentioned.

Now I'm in a fix which answer to accept ;)

[Ameer Hamza]

I would accept John's answer. Anyone searching documentation can find the thing I mentioned, but John's answer provides a more in-depth insight.

[John D'Errico]

Honestly, I don't care which one you accept, as I never chase site rep. We both made a serious effort to help you. If we did help you, then this is a good thing, and is all that matters to me. :)

So perhaps you might accept one based on a coin flip, and give an upvote to the other, thus a tip of the hat to both answers.

By the way, years ago, when I was actively reading about global solvers and the state of the art at the time, I did find there was a lot to be read. That was 30+ years ago. I'm sure there has been more written since then.

[Mohit Kumar]

Thanks guys!

The MATLAB community is awesome due to people like you :)

Answer 2

[Ameer Hamza]

If you have the symbolic toolbox, then the preferable approach will be to use solve(): https://www.mathworks.com/help/symbolic/solve.html. But that will only work if an analytical solution exists.

The next best approach is to use vpasolve(): https://www.mathworks.com/help/symbolic/vpasolve.html. This will probably work even if an analytical solution does not exists.

If both of these approaches are not feasible, then you can create a mesh and find the points where the value of function close to zero. Use those points as an initial guess to get the actual solutions.

[Mohit Kumar]

Thanks for your answer. About using solve - not feasible as the equations I'm looking to solve won't have analytical solutions in general. About using vpasolve - I had tried this a week or so back. If I remember correctly, I tried solving sin(x)=0.1x which clearly has more than one solution. However, I only got one solution depending on the initial guess. About using the mesh for initial guesses - In this method, how would I determine what is "small"? My initial idea would be something like say 1/10th of the average absolute values. Maybe I'll try this and update a couple of days later (currently I'm busy :( )

[Ameer Hamza]

Yes, choosing the threshold for small value can be difficult. What about this strategy. If you already know that there are k roots, then select the k smallest values from the mesh and use the corresponding points as the initial guess. The strategy is not perfect, so maybe try increasing the number of points, say 2*k minimum values.