Evaluación numéricamente estable de $\frac{x^a - 1}{x^a (x - 1)}$

Question

Quiero evaluar la siguiente función como parte de un programa de ordenador: $$f(x) = \frac{x^a - 1}{x^a (x - 1)}, x \in \mathbb{R}^+, a \in \mathbb{R}_0^+$$

Matemáticamente, esta función se comporta bastante bien, con la singularidad en $x = 1$ siendo una singularidad removible. Sin embargo, aunque la singularidad puede tratarse por separado, la evaluación es numéricamente inestable en la vecindad de $x = 1$ debido a que tanto el numerador como el denominador son muy pequeños.

Aunque $x = 1$ es una raíz del numerador, he sido incapaz de expresar algebraicamente el numerador como un producto $(x - 1)r$ para anular el factor lineal $x - 1$ para eliminar la singularidad. Sé que es posible que $a \in \mathbb{N}$ mediante las fórmulas de la progresión geométrica, pero el polinomio resultante es costoso de evaluar y en mi caso, $a$ puede ser un número real no negativo arbitrario, por lo que esto no se aplica de todos modos.

¿Hay alguna forma de expresar algebraicamente una función que es $f$ ¿pero sin la singularidad? Si no es así, ¿hay alguna forma de evaluar $f$ evitando la inestabilidad numérica en torno a $x = 1$ ?

Answer 1

Quiero destacar el comentario de Lutz Lehmann sobre la expm1 que es omnipresente. Su propósito es calcular eficientemente $e^x -1$ para $x \approx 0$ .

Lutz da la fórmula

expm1(a*log(x))/pow(x,a)/(x-1)

que me parece correcto, y de hecho parece ser estable para $x \approx 1$ . Debe comprobar $x = 1$ para evitar la división por cero.

Sospecho que todas las respuestas que utilizan series de Taylor son esencialmente reinventar la rueda aquí.

Answer 2

Para mayor comodidad $x=1+s$ entonces cerca de $s=0$ se puede utilizar la serie de Taylor

$$ \eqalign{f(1+s) &= \sum_{n=0}^\infty \frac{(-1)^n \Gamma(n+a+1)}{\Gamma(a) (n+1)!} s^n \cr &= a - \dfrac{a(a+1)}{2} s + \dfrac{a(a+1)(a+2)}{6} s^2 - \dfrac{a(a+1)(a+2)(a+3)}{24} s^3 + \ldots}$$

Answer 3

Presumiblemente, los problemas no se deben a que los valores cercanos a $1$ son pequeños o grandes pre se, sino por cancelación cuando $x\approx 1$ . $\def\e{\varepsilon}$

Así que cuando $x\approx 1$ puede intentar aproximar $f$ por una función más simple como la expansión de Taylor en torno a $1$ por ejemplo

$$f(x)\approx \frac ax-\frac12 a(a-1)(x-1)$$

En $x\approx 1$ escribámoslo como $x=1+\e$ :

$$\begin{align} f(x)=\frac{x^a-1}{x^a(x-1)} &= \frac1{x^a}\underbrace{\frac{(1+\e)^a-1^a}{\e}}_{\textstyle\approx g'(1)+\dfrac\e2 g''(1)} \\ \end{align}$$ donde $g(x) = x^a$ y así $$g'(x) = ax^{a-1}$$ $$g''(x) = a(a-1)x^{a-2}$$

He aquí una Parcela Desmos acercado a $x\approx 1$ . Obsérvese que las funciones del gráfico tienen en realidad $a$ restados para que todos alcancen $(1,0)$ independientemente del valor de $a$ . Puede jugar con $a$ arrastrando el control deslizante.

La función roja es $f(x)$ que se vuelve loco cerca de 1. La función violeta es la aproximación cerca de 1.

Por lo que sé, Desmos utiliza IEEE-754 double, es decir, coma flotante binaria "ordinaria" de 64 bits, por lo que se obtienen algunas estimaciones de calidad visual.

Answer 4

Utilizando una aproximante simple de Padé

$$\frac{x^a - 1}{x^a (x - 1)}=a \frac{1+\frac{7-a}{10} (x-1)+\frac{(a-1)(a-2)}{60} (x-1)^2 } {1+\frac{2(a+3)}{5} (x-1)+\frac{(a+2)(a+3)}{20} (x-1)^2 }$$

Intentando $a=\pi$ y $x=1+10^{-k}$ algunos resultados

$$\left( \begin{array}{ccc} k & \text{approximation} & \text{exact value} \\ 1 & \color{red}{2.58756}3328300637403099800 & 2.587562502050549708509754 \\ 2 & \color{red}{3.0776347615}76885361735203 & 3.077634761563679134659739 \\ 3 & \color{red}{3.13509818768015}6052084880 & 3.135098187680155913340110 \\ 4 & \color{red}{3.14094220521706779548}2107 & 3.140942205217067795480713 \\ 5 &\color{red}{ 3.141527598719473976878852} & 3.141527598719473976878852 \end{array} \right)$$ y podríamos hacerlo mucho mejor.