Numerically find cubic polynomial roots where coefficients widely vary in magnitude

machine precisionnumerical methodspolynomialsroots

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.

Best Answer

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

Related Question