Considere el siguiente polinomio:
$$ 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$ se supone que ser constante;
- $K_a$ es ácido/base de la constante de disociación, normalmente dentro de $[10^{-12},10^{-2}]$;
- $K_w$ es el agua iónica producto constante que se acerca $10^{-14}$.
Espero que este polinomio para tener al menos una solución real positiva debido a que la concentración de protones existen físicamente.
Mi objetivo es encontrar la verdadera solución positiva de la ecuación de $p(x)=0$ mediante el uso de algún método numérico (Newton, etc. o incluso fórmulas de Cardano), pero me enfrento a algunos de normalización problema debido a la diferencia de magnitudes entre las constantes o monomials.
Para una instalación típica:
$$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}$, que es muy pequeño, puede ser cualquier cosa: cero, positivo o negativo, ¿quién sabe?
He utilizado tanto float y fixed decimal (con 40 decimales) aritmética para calcular estas cantidades. Ellos difieren ligeramente con la aritmética, sino de la ampliación problema sigue. Creo que mi problema es sobre numérico de la estabilidad y la normalización.
Wolfram parece ser capaz de resolverlo cuando yo uso simbólico de valores (la respuesta es correcta) en lugar de los coeficientes de desarrollo decimal (la respuesta es demasiado pequeño).
No estoy pidiendo una solución completa, en lugar de pistas o referencias a un método para encontrar raíces de este tipo de enfermos en escala de polinomios. Así que mis preguntas son:
- ¿Cómo puedo resolver con precisión este tipo de polinomio con un método numérico?
- ¿Qué tipo de normalización debe aplicar antes de resolverlo?
- Hay un específico método numérico para esta clase de problema?
Nota
El polinomio $p(x)$ proviene de la definición del ácido/base de la constante de equilibrio:
$$ 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 $$
Actualización
He logrado encontrar una creíble root usando numpy.roots y decimal que cumple con Wolfram Alpha. El siguiente Python 3 código:
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)
Devuelve:
array([-2.00026673e-01, 8.89021156e-06, -9.99999983e-14])
Todavía estoy interesado en saber si hay otros cuidados a tomar para solucionar este tipo de problemas.
