Puede usar esta función para ser reescrito para mejorar la estabilidad numérica?

Question

Estoy escribiendo un programa que necesita para evaluar la función $$f(x) = \frac{1 - e^{-ux}}{u}$$ often with small values of $u$ (i.e. $u \ll x$). In the limit $u \a 0$ we have $f(x) = x$ using L'Hôpital's rule, so the function is well-behaved and non-singular in this limit, but evaluating this in the straightforward way using floating-point arithmetic causes problems when $u$ is small: then the exponential is quite close to 1, and when we subtract it from 1, the difference doesn't have great precision. Is there a way to rewrite this function so it can be computed more accurately both in the small-$u$ limit and for larger values of $u$ (i.e. $u \aprox$x)?

Por supuesto, yo podría añadir una prueba a mi programa que cambia a la informática sólo $f(x) = x$, o quizás $f(x) = x - ux^2/2$ (de segundo orden de la aproximación) al $u$ es menor que un cierto umbral. Pero tengo curiosidad por ver si hay una manera de hacer esto sin la introducción de una prueba.

Answer 1

Una vez tenido esto (o al menos algo muy similar) en mi primer análisis numérico de la clase:

Allí fue propuesto para calcular el $f$ en términos de $y = \exp(-ux)$ lugar. I. e. primer set $y$ y, a continuación, calcular $$f(x) = \frac{1-e^{-ux}}{u} = -x\frac{1-y}{\log(y)}$$

Nunca he entendido por qué esto debería tener ningún efecto sobre el cálculo, sin embargo... (Que sólo puede tener un efecto, si hemos de hacer algunas suposiciones sobre cómo los valores se calculan dentro del ordenador).

Edición: J. M. dio un muy buen enlace en los comentarios, en donde se explica por qué el anterior funciona!

Si puedes hacer algunos tests numéricos, que me haga saber lo que se obtiene (me interesaría si esto realmente hace lo que afirmó que lo haría! ;)

Answer 2

La razón por la que el error es debido a la precisión de la máquina de aritmética de punto flotante.

Lo que está sucediendo no es tanto que $1-e^{ux}$ tiene problemas con la precisión de la máquina, sino que se están dividiendo, por algo que es bastante pequeño, así.

Esencialmente, si $e^{ux}$ es de menos de aproximadamente $2.2 \times 10^{-16}$ (la máquina de $\varepsilon$) de doble precisión de la aritmética, entonces la cantidad total $1-e^{ux}$ regresará exactamente 1. Sin embargo, esto no es tanto un problema... 1 veces algo es casi la misma cosa como 0.99999999999999... los tiempos de algo.

El problema es que luego dividiendo por un número muy pequeño: $1/u$. Si su $u < 2.2\times 10^{-16}$, $1/u$ es muy grande, y el numerador es incapaz de "cancelar" para este efecto. Por otra parte, y no recuerdo lo que la OFL límite para el dobles están fuera de la parte superior de mi cabeza, si $1/u$ es mayor que algo como $2^32$, estás de suerte, de todos modos.

Una cosa que puedes hacer es calcular $\log f(x)$ primera: $\log f(x) = \log(1-e^{ux})-\log(u)$. Usted todavía tiene la máquina epsilon problema (es decir, los valores pequeños de a $u$ esencialmente hacer el registro de devolución 0, pero esto va a ser casi cierto para los valores en un barrio de todos modos), pero su denominador será más manejable.

Edit: por supuesto, otro método es utilizar la fiel ol' expansión de Taylor para $e^{ut}$. Entonces, usted tiene $f(x) = \frac{1}{u}+\sum_i \frac{(ux)^i}{i!u}$. El numerador es pequeño, y el denominador es pequeño, por lo que en cierto grado se cancelan uno al otro. Calcular cada término en la suma de forma independiente, almacenarlos, y luego en añadir.

Answer 3

Como J. M., señaló en un comentario a otra respuesta, algunas plataformas y lenguajes de programación tienen una primitiva que calcula el $2^x-1$ o $e^x-1$ en un solo paso con sin pérdida de precisión para las pequeñas $x$. Por ejemplo, en C99 <math.h> no es la expm1() función, que debe asignar a la F2XM1 de la instrucción en la plataforma x86 (a través de una escala de $x$$\frac{1}{\log2}$).

Answer 4

Su problema es que no hay una función de biblioteca que conserva la precisión de que el argumento para la función exponencial en sus resultados (en relación a uno) cuando el argumento se encuentra en el barrio de cero. Pero se puede escribir como una función de sí mismo, aunque de forma rudimentaria:

public static double ExpMinusOne(double x)
{
    return Math.Abs(x) < 1e5 ? x + x * x / 2 + x * x * x / 6 : Math.Exp(x) - 1;
}

que debe dar alrededor de diez dígitos de precisión en $f(x)$ con IEEE de punto flotante, independientemente de lo pequeño que $u$ es. Funciona mediante el uso de la normal función exponencial cuando el argumento es lo suficientemente grande y el uso de los cuatro primeros términos de la expansión en series de Taylor de otra manera. Debido a que la serie converge rápidamente cuando el argumento es pequeño, usted todavía obtener una razonable precisión.

Answer 5

Todavía otra posibilidad: si su entorno de computación de elección apoya el seno hiperbólico, entonces también puede utilizar la fórmula

$$f(x)=\frac{2\exp\left(-\frac{u x}{2}\right)\sinh\left(\frac{u x}{2}\right)}{u}$$

Suponiendo que el seno hiperbólico de rutina en su entorno de computación está bien escrito para muy pequeños argumentos, los errores en el cálculo del numerador y el denominador deben cancelar a cabo muy bien.