Numerically find cubic polynomial roots where coefficients widely vary in magnitude

Question

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.

Answer 1

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}$

Answer 2

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}$