fsolve/fzero: No solution found, appears regular

Question

I am trying to do the following algorithm using fsolve or fzero:

K5=8.37e-2
P=1

Choose an A
    S2=(4*K5/A)^(2/3)
    S6=3*S2
    S8=4*S2

    SO2 = (5*P)/149 - (101*S2)/149 - (293*S6)/149 - (389*S8)/149
    H2O = (40*P)/149 + (556*S2)/447 + (636*S6)/149 + (2584*S8)/447
    H2S = 2*SO2

    newA = (H2O)^2/(SO2)^3

Repeat until newA=oldA

The main thing to solve is K5=1/4 * A * S2^3/2. It is from this that S2 is calculated in the first place.

So here is what I did in Matlab:

function MultipleNLEexample
clear, clc, format short g, format compact

Aguess = 300000; % initial guess

options = optimoptions('fsolve','Display','iter','TolFun',[1e-9],'TolX',[1e-9]); % Option to display output
xsolv=fsolve(@MNLEfun,Aguess,options); 

[~,ans]=MNLEfun(xsolv)

%- - - - - - - - - - - - - - - - - - - - - -
function varargout = MNLEfun(A); 
    K5 = 8.37e-2;

    S2 = (4*K5/A)^(2/3);
    S6 = 3*S2;
    S8 = 4*S2;

    P=1; %atm

    SO2 = (5*P)/149 - (101*S2)/149 - (293*S6)/149 - (389*S8)/149;
    H2O = (40*P)/149 + (556*S2)/447 + (636*S6)/149 + (2584*S8)/447;
    newA=H2O^2/SO2^3;
    fx=1/4*newA*S2^(3/2)-K5;

    varargout{1} = fx;
    if nargout>1
        H2S = 2*SO2;
        varargout{2} = ((2*S2+6*S6+8*S8)/(2*S2+6*S6+8*S8+H2S+SO2)*100);
    end

I cannot get my code to run, I get the following error: No solution found.

fsolve stopped because the problem appears regular as measured by the gradient, but the vector of function values is not near zero as measured by the selected value of the function tolerance.

I have tried setting the tolerances as low as 1e-20 but that did not change anything.

Answer 1

The way your system is set up, it is actually convenient to plot it and observe its behavior. I vectorized your function and plottedf(x) = MLNEfun(x)-x, where the output of MLNE(x) is newA. Effectively, you are interested in a fixed point of your system.

What I observe is this:

There is a singularity, and a root crossing, at A ~ 3800. We can use fzero since it is a bracketed root solver, and give it very tight bounds on the solution fzero(@(x)MLNEfun(x)-x, [3824,3825]) which produces 3.8243e+03. This is a couple order of magnitudes from your starting guess. There is no solution to your system near ~3e5.

Update In my haste, I failed to zoom in on the plot, which shows another (well-behaved) root at 1.3294e+04. It is up to you to decide which is the physically meaningful one(s). Everything I say below still applies. Just start your guess near the solution you're interested in.

In response to comment Since you want to perform this for varying values of K, then your best bet is to stick with fzero so long as you are solving for one variable, instead of fsolve. The reason behind this is that fsolve uses variants of Newton's method, which are not bracketed and will struggle to find solutions at singular points like this. fzero on the other hand, uses Brent's method which is guaranteed to find a root (if it exists) within a bracket. It is also much more well behaved near singular points.

MATLAB's implementation of fzero also searches for a bracketing interval before carrying out Brent's method. So, if you provide a single starting guess sufficiently close to the root, it should find it for you. Sample output from fzero below:

fzero(@(x)MLNEfun(x)-x, 3000, optimset('display', 'iter'))

Search for an interval around 3000 containing a sign change:
 Func-count    a          f(a)             b          f(b)        Procedure
    1            3000       -616789          3000       -616789   initial interval
    3         2915.15       -433170       3084.85       -905801   search
    5            2880       -377057          3120  -1.07362e+06   search
    7         2830.29       -311972       3169.71  -1.38274e+06   search
    9            2760       -241524          3240  -2.03722e+06   search
   11         2660.59       -171701       3339.41  -3.80346e+06   search
   13            2520       -109658          3480  -1.16164e+07   search
   15         2321.18      -61340.4       3678.82   -1.7387e+08   search
   17            2040      -29142.6          3960   2.52373e+08   search

Search for a zero in the interval [2040, 3960]:
 Func-count    x          f(x)             Procedure
   17            2040      -29142.6        initial
   18         2040.22      -29158.9        interpolation
   19         3000.11       -617085        bisection
   20         3480.06  -1.16224e+07        bisection
   21            3960   2.52373e+08        bisection
   22         3720.03  -4.83826e+08        interpolation
....
   87         3824.32  -5.46204e+48        bisection
   88         3824.32   1.03576e+50        bisection
   89         3824.32   1.03576e+50        interpolation

Current point x may be near a singular point. The interval [2040, 3960] 
reduced to the requested tolerance and the function changes sign in the interval,
but f(x) increased in magnitude as the interval reduced.

ans =

   3.8243e+03