Mixing
Numerically find cubic polynomial roots where coefficients widely vary in magnitude
Clash Royale CLAN TAG#URR8PPP
up vote
4
down vote
favorite
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.7792cdot 10^-11,quad a_0=-1.7783cdot 10^-24
$$
The discriminant of this polynomial for this setup is about $Delta = 5.6905cdot 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 = xfracba, 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.
polynomials numerical-methods roots machine-precision
asked Jul 26 at 8:11
jlandercy
210212
 |Â
show 1 more comment
up vote
4
down vote
favorite
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.7792cdot 10^-11,quad a_0=-1.7783cdot 10^-24
$$
The discriminant of this polynomial for this setup is about $Delta = 5.6905cdot 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 = xfracba, 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.
polynomials numerical-methods roots machine-precision
asked Jul 26 at 8:11
jlandercy
210212
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 at 16:15
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-11)$. I get very similar solutions out of Kahan'sQBCsolver, using a C++ code that usesdoublearithmetic enhanced by the use of FMA (fused multiply-add):-2.0000000008896002e-1, -9.9837370721755525e-14, 8.9059837331107879e-11
– njuffa
Jul 28 at 19:40
1
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^-8y$ 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 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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 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 at 6:55
 |Â
show 1 more comment
up vote
4
down vote
favorite
up vote
4
down vote
favorite
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.7792cdot 10^-11,quad a_0=-1.7783cdot 10^-24
$$
The discriminant of this polynomial for this setup is about $Delta = 5.6905cdot 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 = xfracba, 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.
polynomials numerical-methods roots machine-precision
asked Jul 26 at 8:11
jlandercy
210212
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.7792cdot 10^-11,quad a_0=-1.7783cdot 10^-24
$$
The discriminant of this polynomial for this setup is about $Delta = 5.6905cdot 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 = xfracba, 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.
polynomials numerical-methods roots machine-precision
asked Jul 26 at 8:11
jlandercy
210212
edited Jul 26 at 9:01
asked Jul 26 at 8:11
jlandercy
210212
asked Jul 26 at 8:11
jlandercy
210212
asked Jul 26 at 8:11
jlandercy
210212
210212
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 at 16:15
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-11)$. I get very similar solutions out of Kahan'sQBCsolver, using a C++ code that usesdoublearithmetic enhanced by the use of FMA (fused multiply-add):-2.0000000008896002e-1, -9.9837370721755525e-14, 8.9059837331107879e-11
– njuffa
Jul 28 at 19:40
1
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^-8y$ 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 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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 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 at 6:55
 |Â
show 1 more comment
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 at 16:15
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-11)$. I get very similar solutions out of Kahan'sQBCsolver, using a C++ code that usesdoublearithmetic enhanced by the use of FMA (fused multiply-add):-2.0000000008896002e-1, -9.9837370721755525e-14, 8.9059837331107879e-11
– njuffa
Jul 28 at 19:40
1
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^-8y$ 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 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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 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 at 6:55
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 at 16:15
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 at 16:15
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-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 at 19:40
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-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 at 19:40
1
1
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^-8y$ 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 at 20:39
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^-8y$ 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 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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 at 20:45
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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 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 at 6:55
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 at 6:55
 |Â
show 1 more comment
2 Answers
2
active
oldest
votes
up vote
3
down vote
accepted
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (textCb+textKa)-x (textCa ,textKa+textKw)-textKa
, textKw$$
$$f'(x)=3 x^2+2 x (textCb+textKa)-(textCa textKa+textKw)$$
$$f''(x)=6 x+2 (textCb+textKa)$$ you must not start at $x=0$ since
$$f(0) ,f''(0)=-2, textKa ,textKw ,(textCb+textKa)<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-fractextKa , textKwtextCa ,textKa+textKw <0$$
The first derivative cancels at
$$x_*=frac13 left(sqrt(textCb+textKa)^2+3 (textCa,
textKa+textKw)-(textCb+textKa)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-2frac f(x_*)f''(x_*)$$
Using your numbers, Newton iterates would be
$$left(
beginarraycc
n & x_n \
0 & colorblue8.8902609452502354578 times 10^-6\
1 & colorblue8.8902115571665711819 times 10^-6\
2 & colorblue8.8902115568921917511 times 10^-6
endarray
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 12frac f(x_0) left(f(x_0), f''(x_0)-2, f'(x_0)^2right)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=colorblue8.8902115568921917512 times 10^-6$
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
add a comment |Â
up vote
1
down vote
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.89021times10^-6$
answered Jul 28 at 20:29
njuffa
6761813
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
add a comment |Â
2 Answers
2
active
oldest
votes
2 Answers
2
active
oldest
votes
active
oldest
votes
active
oldest
votes
up vote
3
down vote
accepted
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (textCb+textKa)-x (textCa ,textKa+textKw)-textKa
, textKw$$
$$f'(x)=3 x^2+2 x (textCb+textKa)-(textCa textKa+textKw)$$
$$f''(x)=6 x+2 (textCb+textKa)$$ you must not start at $x=0$ since
$$f(0) ,f''(0)=-2, textKa ,textKw ,(textCb+textKa)<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-fractextKa , textKwtextCa ,textKa+textKw <0$$
The first derivative cancels at
$$x_*=frac13 left(sqrt(textCb+textKa)^2+3 (textCa,
textKa+textKw)-(textCb+textKa)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-2frac f(x_*)f''(x_*)$$
Using your numbers, Newton iterates would be
$$left(
beginarraycc
n & x_n \
0 & colorblue8.8902609452502354578 times 10^-6\
1 & colorblue8.8902115571665711819 times 10^-6\
2 & colorblue8.8902115568921917511 times 10^-6
endarray
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 12frac f(x_0) left(f(x_0), f''(x_0)-2, f'(x_0)^2right)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=colorblue8.8902115568921917512 times 10^-6$
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
add a comment |Â
up vote
3
down vote
accepted
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (textCb+textKa)-x (textCa ,textKa+textKw)-textKa
, textKw$$
$$f'(x)=3 x^2+2 x (textCb+textKa)-(textCa textKa+textKw)$$
$$f''(x)=6 x+2 (textCb+textKa)$$ you must not start at $x=0$ since
$$f(0) ,f''(0)=-2, textKa ,textKw ,(textCb+textKa)<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-fractextKa , textKwtextCa ,textKa+textKw <0$$
The first derivative cancels at
$$x_*=frac13 left(sqrt(textCb+textKa)^2+3 (textCa,
textKa+textKw)-(textCb+textKa)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-2frac f(x_*)f''(x_*)$$
Using your numbers, Newton iterates would be
$$left(
beginarraycc
n & x_n \
0 & colorblue8.8902609452502354578 times 10^-6\
1 & colorblue8.8902115571665711819 times 10^-6\
2 & colorblue8.8902115568921917511 times 10^-6
endarray
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 12frac f(x_0) left(f(x_0), f''(x_0)-2, f'(x_0)^2right)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=colorblue8.8902115568921917512 times 10^-6$
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
add a comment |Â
up vote
3
down vote
accepted
up vote
3
down vote
accepted
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (textCb+textKa)-x (textCa ,textKa+textKw)-textKa
, textKw$$
$$f'(x)=3 x^2+2 x (textCb+textKa)-(textCa textKa+textKw)$$
$$f''(x)=6 x+2 (textCb+textKa)$$ you must not start at $x=0$ since
$$f(0) ,f''(0)=-2, textKa ,textKw ,(textCb+textKa)<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-fractextKa , textKwtextCa ,textKa+textKw <0$$
The first derivative cancels at
$$x_*=frac13 left(sqrt(textCb+textKa)^2+3 (textCa,
textKa+textKw)-(textCb+textKa)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-2frac f(x_*)f''(x_*)$$
Using your numbers, Newton iterates would be
$$left(
beginarraycc
n & x_n \
0 & colorblue8.8902609452502354578 times 10^-6\
1 & colorblue8.8902115571665711819 times 10^-6\
2 & colorblue8.8902115568921917511 times 10^-6
endarray
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 12frac f(x_0) left(f(x_0), f''(x_0)-2, f'(x_0)^2right)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=colorblue8.8902115568921917512 times 10^-6$
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
If you consider using Newton method to find the zero of
$$f(x)=x^3+x^2 (textCb+textKa)-x (textCa ,textKa+textKw)-textKa
, textKw$$
$$f'(x)=3 x^2+2 x (textCb+textKa)-(textCa textKa+textKw)$$
$$f''(x)=6 x+2 (textCb+textKa)$$ you must not start at $x=0$ since
$$f(0) ,f''(0)=-2, textKa ,textKw ,(textCb+textKa)<0$$ (think about Darboux-Fourier theorem). If you use $x_0=0$, then the first iterate would be
$$x_1=-fractextKa , textKwtextCa ,textKa+textKw <0$$
The first derivative cancels at
$$x_*=frac13 left(sqrt(textCb+textKa)^2+3 (textCa,
textKa+textKw)-(textCb+textKa)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-2frac f(x_*)f''(x_*)$$
Using your numbers, Newton iterates would be
$$left(
beginarraycc
n & x_n \
0 & colorblue8.8902609452502354578 times 10^-6\
1 & colorblue8.8902115571665711819 times 10^-6\
2 & colorblue8.8902115568921917511 times 10^-6
endarray
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 12frac f(x_0) left(f(x_0), f''(x_0)-2, f'(x_0)^2right)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=colorblue8.8902115568921917512 times 10^-6$
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
edited Jul 30 at 3:48
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
answered Jul 28 at 16:04
Claude Leibovici
111k1055126
111k1055126
add a comment |Â
add a comment |Â
up vote
1
down vote
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.89021times10^-6$
answered Jul 28 at 20:29
njuffa
6761813
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
add a comment |Â
up vote
1
down vote
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.89021times10^-6$
answered Jul 28 at 20:29
njuffa
6761813
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
add a comment |Â
up vote
1
down vote
up vote
1
down vote
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.89021times10^-6$
answered Jul 28 at 20:29
njuffa
6761813
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.89021times10^-6$
answered Jul 28 at 20:29
njuffa
6761813
answered Jul 28 at 20:29
njuffa
6761813
answered Jul 28 at 20:29
njuffa
6761813
answered Jul 28 at 20:29
njuffa
6761813
6761813
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
add a comment |Â
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 at 6:56
I don't think that we need to solve the cubic since only the positive root is required.
– Claude Leibovici
Jul 29 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 at 7:19
@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 at 7:19
add a comment |Â
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
StackExchange.ready(
function ()
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2863186%2fnumerically-find-cubic-polynomial-roots-where-coefficients-widely-vary-in-magnit%23new-answer', 'question_page');
);
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Sign up or log in
StackExchange.ready(function ()
StackExchange.helpers.onClickDraftSave('#login-link');
);
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
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 at 16:15
For the coefficients given, $(1.0, 0.2, -1.7792cdot10^-11, -1.7783cdot10^-24)$, Wolfram Alpha returns solutions $(-0.2, -9.98318cdot10^-14, 8.90598cdot10^-11)$. I get very similar solutions out of Kahan's
QBCsolver, using a C++ code that usesdoublearithmetic enhanced by the use of FMA (fused multiply-add):-2.0000000008896002e-1, -9.9837370721755525e-14, 8.9059837331107879e-11– njuffa
Jul 28 at 19:40
1
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^-8y$ 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 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=5sqrtstuff$ where stuff comes from the non-quadratic terms and do fixed point iteration starting from the root of the quadratic.
– Ross Millikan
Jul 28 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 at 6:55