Consider the following polynomial:

$$ p(x) = x^3 + (C_b+K_a)x^2 - (C_aK_a + K_w)x - K_aK_w $$


  • $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?


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


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)


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.

  • $\begingroup$ 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$. $\endgroup$ – John Barber Jul 28 '18 at 16:15
  • $\begingroup$ 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 $\endgroup$ – njuffa Jul 28 '18 at 19:40
  • 1
    $\begingroup$ 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$. $\endgroup$ – Ross Millikan Jul 28 '18 at 20:39
  • $\begingroup$ 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. $\endgroup$ – Ross Millikan Jul 28 '18 at 20:45
  • $\begingroup$ 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. $\endgroup$ – Claude Leibovici Jul 29 '18 at 6:55

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.


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:

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

  • $\begingroup$ I don't think that we need to solve the cubic since only the positive root is required. $\endgroup$ – Claude Leibovici Jul 29 '18 at 6:56
  • $\begingroup$ @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. $\endgroup$ – njuffa Jul 29 '18 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.