Numerically find cubic polynomial roots where coefficients widely vary in magnitude Consider the following polynomial:
$$
p(x) = x^3 + (C_b+K_a)x^2 - (C_aK_a + K_w)x - K_aK_w
$$
Where:


*

*$x, C_a, C_b$ are concentrations, positive real numbers typically within $[10^{-7};1]$. The unknown $x$ is the proton concentration and $C_a, C_b$ are supposed to be constant;

*$K_a$ is acid/base dissociation constant, typically within $[10^{-12},10^{-2}]$;

*$K_w$ is the water ionic product constant which is about $10^{-14}$.


I expect this polynomial to have at least a positive real solution because proton concentration physically exist. 
My goal is to find the real positive solution to the equation $p(x)=0$ using some numerical method (Newton, etc. or even Cardano formulae), but I face some normalization problem because of different magnitudes between constants or monomials.
For a typical setup:
$$K_a = 10^{-4.75},\quad C_a = 0.1,\quad C_b = 0.2$$
Polynomial coefficients widely varies in magnitude:
$$
a_3=1.0000,\quad a_2=0.2000,\quad a_1=-1.7792\cdot 10^{-11},\quad a_0=-1.7783\cdot 10^{-24}
$$
The discriminant of this polynomial for this setup is about $\Delta = 5.6905\cdot 10^{-26}$ which is really small, it could be anything: zero, positive or negative, who knows?
I have used both float and fixed decimal (with 40 decimals) arithmetics to compute those quantities. They slightly differ with the arithmetic kind, but the scaling problem remains. I think my problem is about numerical stability and normalization.
Wolfram seems to be able to solve it when I use symbolic values (answer is correct) instead of coefficients decimal development (answer is too small). 
I am not asking for a complete solution, rather for some hints or references to a method for roots finding for this kind of ill-scaled polynomials. So my questions are:


*

*How can I accurately solve this kind of polynomial with a numerical method?

*What kind of normalization must I apply before solving it?

*Is there a specific numerical method for this class of problem?


Nota
The polynomial $p(x)$ comes from the definition of the acid/base equilibrium constant:
$$
K_a = x\frac{b}{a}, \quad K_w = xy
$$
Where matter amount and charge balances have been injected:
$$
C_a + C_b = a + b,\quad b + y = x + C_b
$$
Update
I have managed to find a credible root using numpy.roots and decimal which complies with Wolfram Alpha. The following Python 3 code: 
import numpy as np
import decimal

Ka = decimal.Decimal("10")**decimal.Decimal("-4.75")
Kw = decimal.Decimal("1e-14")
Ca = decimal.Decimal("0.1")
Cb = decimal.Decimal("0.2")
coefs = [decimal.Decimal(1), (Cb+Ka), -(Ca*Ka + Kw), -Ka*Kw]

r0 = np.roots(coefs)

Returns:
array([-2.00026673e-01,  8.89021156e-06, -9.99999983e-14])

I am still interested to know if there are other care to take in order to properly solve this kind of problems. 
 A: If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (\text{Cb}+\text{Ka})-x (\text{Ca} \,\text{Ka}+\text{Kw})-\text{Ka}
  \, \text{Kw}$$
$$f'(x)=3 x^2+2 x (\text{Cb}+\text{Ka})-(\text{Ca} \text{Ka}+\text{Kw})$$
$$f''(x)=6 x+2  (\text{Cb}+\text{Ka})$$ you must not start at $x=0$ since
$$f(0) \,f''(0)=-2\, \text{Ka}  \,\text{Kw}  \,(\text{Cb}+\text{Ka})<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-\frac{\text{Ka} \, \text{Kw}}{\text{Ca}  \,\text{Ka}+\text{Kw}} <0$$
The first derivative cancels at
$$x_*=\frac{1}{3} \left(\sqrt{(\text{Cb}+\text{Ka})^2+3 (\text{Ca}\,
   \text{Ka}+\text{Kw})}-(\text{Cb}+\text{Ka})\right) > 0$$
 and the second derivative test shows that this is a minimum.
In order to generate the starting point $x_0$, perform around $x=x_*$ (were $f'(x_*)=0$), the Taylor expansion
$$0=f(x_*)+\frac 12 f''(x_*)\,(x-x_*)^2+O((x-x_*)^3 ) \implies x_0=x_*+ \sqrt{-2\frac {f(x_*)}{f''(x_*)}}$$
Using your numbers, Newton iterates would be
$$\left(
\begin{array}{cc}
 n & x_n \\
 0 & \color{blue}{8.8902}609452502354578 \times 10^{-6}\\
 1 & \color{blue}{8.89021155}71665711819 \times 10^{-6}\\
 2 & \color{blue}{8.8902115568921917511} \times 10^{-6}
\end{array}
\right)$$ This is a very reliable procedure.
Edit
I am almost sure that we could even skip the Newton procedure using a $[1,2]$ Padé approximant built around $x_1$. This would give
$$x_1=x_0+\frac 12\frac{ f(x_0) \left(f(x_0)\, f''(x_0)-2\, f'(x_0)^2\right)}{f(x_0)^2 +f'(x_0)^3- f(x_0)\,f'(x_0)\, f''(x_0)}$$
For the worked example, this would give $x_2=\color{blue}{8.890211556892191751}2 \times 10^{-6}$
A: According to the computation on Wolfram Alpha linked in the question, the task is to find the roots of the cubic equation 
$x^{3} + (0.2 + 1 \times 10^{-4.75})x^{2}-(0.1 \times 10^{-4.75} + 1 \times 10^{-14})x - (1 \times 10^{-4.75} \times 1 \times 10^{-14})$
If this is representative of the cubic equations that need to be solved, there seems to be no problem in tackling this with IEEE-754 double precision computation and a standard solver, such as Kahan's QBC solver:
William Kahan, "To Solve a Real Cubic Equation". Lecture notes, UC Berkeley, Nov. 1986 (online)
I am using Kahan's code pretty much verbatim, but have enhanced it in various places by the use of the FMA (fused multiply-add) operation which is readily accessible from C++ as the standard math function fma(). I specified the coefficients as follows:
double a3 =  1.0;
double a2 =  0.2 + exp10 (-4.75);
double a1 = -(0.1 * exp10 (-4.75) + exp10 (-14.0));
double a0 = -(exp10 (-4.75) * exp10 (-14.0));

The results returned by the QBC solver are:
-2.0002667300555729e-001
-9.9999998312876059e-014 + i *  0.0000000000000000e+000
 8.8902115568921925e-006 + i *  0.0000000000000000e+000

These match the results computed by Wolfram Alpha: $-0.200027, -1.\times10^{-13}, 8.89021\times10^{-6}$
