0
$\begingroup$

I know that in a cubic equation $ax^3+bx^2+cx+d=0$ holds $a,b>0$, and $d<0$, so there exists exactly one positive real root. I can solve for all of the roots, no problem, but I want to numerically find only the positive root fast and reliably (millions of times in parallel, on each iteration of an algorithm). Any naïve tests (take a root with a positive real part, and "not too large" imaginary part), turn out not to be numerically reliable. Just taking the second root in the code below appears to work in my overall algorithm better that such tests, but is there any reason why this should be so? Or is there a better way to get just the positive real root in this special case that I have? 

function res=cubicsolve(a,b,c,d)
    delta0=b.^2-3.*a.*c;
    delta1=2*b.^3-9*a.*b.*c+27*a.^2.*d;
    xx=delta1.^2-4*delta0.^3;
    c=((delta1+sqrt(xx))/2).^(1/3);
    zeta=-1/2+1i/2*sqrt(3);

    x(:,1)=-(b+c+delta0./c)./(3*a);
    x(:,2)=-(b+zeta*c+delta0./(zeta*c))./(3*a);
    x(:,3)=-(b+zeta^2*c+delta0./(zeta^2*c))./(3*a);

    % This "appears" to work, but I don't know why it should be so
    res=real(x(:, 2));
end
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ What's the problem: numerical stability, or selecting the correct solution among the three? $\endgroup$
    – Federico Poloni
    Commented Mar 22, 2017 at 17:50
  • $\begingroup$ Selecting the correct solution: any naive selection rules are not stable themselves (because of the inaccuracy of the calculations). $\endgroup$
    – user68635
    Commented Mar 22, 2017 at 17:55
  • $\begingroup$ Another approach that appears to be stable in the overall algorithm is to take the largest orthogonal projection of the roots to the positive real axis. But I don't know why this should be the correct root (it may depend on my specific numbers a, b, c, and d). $\endgroup$
    – user68635
    Commented Mar 22, 2017 at 18:08

1 Answer 1

Reset to default
3
$\begingroup$

One possibility is numerical rootfinding. A few iterations of Newton's method should get you home, and I imagine that for a third-degree polynomial one can construct (by sketching a plot of the function) a starting point for which the sufficient condition for convergence is guaranteed to hold (the one with $f f'' > 0$).

Even better if you already have a good solution guess available (for instance, if your millions of problems are slowly varying).

Improve this answer
$\endgroup$