Serious bug in matlab's solve - what can I do?

Question

This post updates my question in Solve finds wrong solution?

Given the "simple" function in x:

f = (x + 3/10)^(1/2) - (9*(x + 3/10)^5)/5 - (x + 3/10)^6 + (x - 1)/(2*(x + 3/10)^(1/2));

Find the zeros by the call

solve(f,x)

This yields 3 zeros:

 ans =

 0.42846617518653978966562924618638
 0.15249587894102346284238111155954
 0.12068186494007759990714181154349

A simple look at the plot shows that the third root is nonsense:

I have a serious problem because I need to retrieve the smallest zero from the above vector. Calling min(ans) returns a wrong zero. What can I do to workaround?

Answer 1

This is a non-polynomial equation, and it will probably fallback to a numeric solver (non-symbolic). So there might be numerical errors, or the numeric algorithm might get stuck and report false solutions, I'm not sure...

What you can do is substitute the solutions back into the equation, and reject ones that are above some specified threshold:

% define function
syms x real
syms f(x)
xx = x+3/10;
f(x) = sqrt(xx) - 9/5*xx^5 - xx^6 + (x - 1)/(2*sqrt(xx));
pretty(f)

% find roots
sol = solve(f==0, x, 'Real',true)

% filter bad solutions
err = subs(f, x, sol)
sol = sol(abs(err) < 1e-9);  % this test removes the 2nd solution

% plot
h = ezplot(f, [0.1 0.5]);
line(xlim(), [0 0], 'Color','r', 'LineStyle',':')
xlabel('x'), ylabel('f(x)')

% programmatically insert data tooltips
xd = get(h, 'XData'); yd = get(h ,'YData');
[~,idx] = min(abs(bsxfun(@minus, xd, double(sol))), [], 2);
dcm = datacursormode(gcf);
pos = [xd(idx) ; yd(idx)].';
for i=1:numel(idx)
    dtip = createDatatip(dcm, h);
    set(get(dtip,'DataCursor'), 'DataIndex',idx(i), 'TargetPoint',pos(i,:))
    set(dtip, 'Position',pos(i,:))
end

We are left only with the two desired solutions (one is rejected by our test):

                   /      3 \5
                 9 | x + -- |
    /      3 \     \     10 /    /      3 \6         x - 1
sqrt| x + -- | - ------------- - | x + -- |  + ----------------
    \     10 /         5         \     10 /          /      3 \
                                               2 sqrt| x + -- |
                                                     \     10 /
sol =
 0.42846617518653978966562924618638
 0.12068186494007759990714181154349       % <== this one is dropped
 0.15249587894102346284238111155954

err(x) =
 -9.1835496157991211560057541970488e-41
   -0.058517436737550288309001512815475   % <==
  1.8367099231598242312011508394098e-40

---

I also tried using MATLAB's numeric solvers, which were able to find the two solutions given reasonable starting points:

• fzero function (a derivative-free root-finding algorithm)
• fsolve function (from the Optimization Toolbox)

(See this related question)

fcn = matlabFunction(f);  % convert symbolic f to a regular function handle
opts = optimset('Display','off', 'TolFun',1e-9, 'TolX',1e-6);

% FZERO
sol2(1) = fzero(fcn, 0.1, opts);
sol2(2) = fzero(fcn, 0.5, opts);
disp(sol2)    % [0.1525, 0.4285]

% FSOLVE
sol3(1) = fsolve(fcn, 0.0, opts);
sol3(2) = fsolve(fcn, 1.0, opts);
disp(sol3)    % [0.1525, 0.4285]

---

For comparison, I tried solving the equation directly into MuPAD, as well as in Mathematica.

Mathematica 9.0

MuPAD (R2014a)

Of course we can always directly call MuPAD from inside MATLAB:

% f is the same symbolic function we've built above

>> sol = evalin(symengine, ['numeric::solve(' char(f) ' = 0, x, AllRealRoots)'])
sol =
[ 0.15249587894102346284238111155954, 0.42846617518653978966562924618638]

The above call is equivalent to searching the entire range x = -infinity .. infinity (which can be slow!). We should assist numeric::solve by providing a specific search ranges when possible:

>> sol = feval(symengine, 'numeric::solve', f==0, 'x = 0 .. 1', 'AllRealRoots')
sol =
[ 0.15249587894102346284238111155954, 0.42846617518653978966562924618638]