# 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.

• You could try perturbation theory, using the smallest coefficient $a_0$ as the small parameter. Substitute $x = x_0 + a_0 y$ and solve for $y$, where $x_0$ is the solution to the quadratic obtained when $a_0 = 0$. – John Barber Jul 28 '18 at 16:15
• For the coefficients given, $(1.0, 0.2, -1.7792\cdot10^{-11}, -1.7783\cdot10^{-24})$, Wolfram Alpha returns solutions $(-0.2, -9.98318\cdot10^{-14}, 8.90598\cdot10^{-11})$. I get very similar solutions out of Kahan's QBC solver, using a C++ code that uses double arithmetic enhanced by the use of FMA (fused multiply-add): -2.0000000008896002e-1, -9.9837370721755525e-14, 8.9059837331107879e-11 – njuffa Jul 28 '18 at 19:40
• When the coefficients vary rapidly and monotonically it can be improved by scaling your variable. Here the ratio between the cubic term and the linear term is about $10^{-24}$, so if you replace $x$ by $10^{-8}y$ all the terms will be about $10^{-24}$ and you can divide it out. Now your quadratic term will be about $10^8$ and the linear one about $10^4$ while the constant and cubic are about $1$. – Ross Millikan Jul 28 '18 at 20:39
• Another similar approach is to note that the value of $x$ will be quite small, so the value of $x^3$ will be much smaller than $0.2x^2$, We can just ignore the cubic term and find that $x$ is of order $10^{-12}$, getting the correct value from the quadratic formula. The cube of this is $10^{-36}$, which is negligible. If you really want to incorporate it, write your equation as x=5\sqrt{stuff}$where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic. – Ross Millikan Jul 28 '18 at 20:45 • I did not want to edit once more my answer but, please, notice that my last formula simplifies since$f'''(x_1)=6$. So,$6$can factor in denominator and$3$in numerator makes the formula quite nice. – Claude Leibovici Jul 29 '18 at 6:55 ## 2 Answers 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}$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}\$

• I don't think that we need to solve the cubic since only the positive root is required. – Claude Leibovici Jul 29 '18 at 6:56
• @ClaudeLeibovici I am not disagreeing, but I think our perspectives may differ: As a software engineer, I tend to look for an existing proven approach that can deliver the requested result (call it the cookbook approach, if you will), even if it delivers a bit more than what is strictly needed. Only if I cannot find that do I embark on crafting a custom solution. It seemed to me that in this case, the former approach would work, but it is hard to tell for sure with only a single test case provided. – njuffa Jul 29 '18 at 7:19