Does fsolve depend on the stability of the system?

Hello all,

I am trying to solve a system of non linear equations using fsolve; lets say

F(x;lambda) = 0, where lambda is a vector of parameters, and x the vector I want to solve for.

I have 2 values of the parameter lambda, that I want to solve the system for. For the one value of lambda I get a solution, which seems alright.

For the other value of lambda I get a solution again (matlab exits with a flag of 1. However I know this is not an actual solution For example I know that some of the dimensions of x have to be equal to each other, and this is not the case in the solution I get from fsolve.

I have tried both trust-region and the levenberg-marquardt algorithm, and I am not getting any better results. (explicitly enforcing those x's to be the same, still seems to give solutions that are not consistent with what I would be expecting from the properties of the system)

My question is: do the algorithms used by fsolve depend on any kind of stability of the system? Could it be that changing the parameter lambda in the second case I mention above, I make the system unstable, and could that make fsolve unable to solve it correctly?

Thank you, George

Answer 1

[Matt J]

Well, the first formulation in your Comment to Alan definitely looks like it could be unstable if the g_i(x)--->0 as x becomes large. That would also mean that f_i(x)--->0 with large x and thus would be insensitive to changes in x in that region.

You could try running fsolve with your second formulation, but that formulation is non-differentiable and undefined for g_i(x)<0, so you may need to move to a constrained solver.

[George A.]

Thank you, good point.

Well the problem is in this second formulation, that for small x, (g_i(x)).^(-1/alpha)--> 0, and so for very small x's f(x) becomes insensitive to changes of x again. Using lsqnonlin and setting a lower bound away from 0, say 0.5, does not help; it prevents the solver from getting any closer to the solution than a point where |fval|=17.

[Matt J]

Hard to know what to do without knowing more about g_i(x). What about

 f(x) = log(g(x))+alpha*log(x)

Answer 2

[Alan Weiss]

If fsolve exited with an exitflag of 1, then the reported solution is a solution. You can verify it by putting it back into your function. If your function is F and your solution is xsol, you will see

 F(xsol)

is a vector with very small values.

What it sounds like you want is a different solution than the one fsolve already found. Generally, this means starting fsolve from a variety of initial points. See this section for ideas on how to take different initial points.

You can also try to use lsqnonlin instead of fsolve, but for this cse I doubt there would be any advantage in doing so, unless you use MultiStart to search for minima of lsqnonneg.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[George A.]

Thank you Allan. Actually my question was more on whether there are (mathematical) conditions to the system, that cause fsolve to perform poorly, e.g. the eigenvalues of F satisfying (or violating) a condition, so that know if there is a way to avoid any such "pathological" formulation.

What I do is the following. My F(x)=0 system looks like the following:

f(1) = x(1)^(-alpha) - g1(x)

f(2) = x(2)^(-alpha) - g2(x) ....

alpha>1;

fsolve returns a solution to that with an exit flag of 1, satisfying a tolerance of say 10^-10.

Now if I give this solution back to an equivalent formulation of the system, like:

f(1) = x(1) - g1(x)^(-1/alpha)

f(2) = x(2) - g2(x)^(-1/alpha) ...

it is far from recognized as a solution; it gives a residual (|fval|) of 10^3 or so. I understand that this makes sense in mathematical terms, because it is as if I raise the residuals to (-1/alpha), alpha>1, this yielding a high residual for the equivalent formulation at the "solution", but the problem remains that I supposedly have a solution to my system which is clearly (I believe) not one.