Cuando se quiere calcular derivadas numéricas, el método presentado por Bengt Fornberg aquí (y reportado aquí ) es muy conveniente (preciso y simple de implementar). Como el documento original data de 1988, me gustaría saber si hay una mejor alternativa hoy (como (o casi como) simple y más precisa).

Visión general

Buena pregunta. Hay un documento titulado "Mejora de la precisión del método de diferenciación de matriz para puntos de colocación arbitrarios" por R. Baltensperger. No es gran cosa en mi opinión, pero tiene un punto (que ya era conocido antes de la aparición en 2000): destaca la importancia de una representación precisa del hecho de que la derivada de la función constante $f(x)=1$ debería ser cero (esto se cumple exactamente en el sentido matemático, pero no necesariamente en la representación numérica).

Es simple ver que esto requiere que las sumas de filas de las n-ésimas matrices derivadas $D^{(n)}$ sean cero. Es común hacer cumplir esta restricción ajustando la entrada diagonal, es decir, estableciendo

$$\tag 1 D^{(n)}_{jj} := -\sum_{\substack{i=1\\i \neq j }}^N D_{ij}\,.$$ Está claro que esta característica no se cumple exactamente cuando se trabaja en una computadora debido a errores de redondeo en los cálculos de coma flotante. Lo que es más sorprendente es que estos errores son aún más graves cuando se utilizan las fórmulas analíticas para la matriz derivada (que están disponibles para muchos puntos de colocación clásicos, por ejemplo, Gauss-Lobatto).

Ahora, el documento (y las referencias en él) establece que el error de la derivada está en el orden de la desviación de las sumas de la fila desde cero. Por lo tanto, el objetivo es hacer que estos sean numéricamente lo más pequeños posible.

Pruebas numéricas

Lo bueno es que el procedimiento de Fornberg parece ser bastante bueno a este respecto. En la siguiente imagen, he comparado el comportamiento de la matriz de primera derivada exacta, es decir, analítica, y la derivada por el algoritmo de Fornberg, para un número variable de puntos de colocación de Chebyshev-Lobatto.

Nuevamente, creyendo en la declaración en el documento citado, esto implica que el algoritmo de Fornberg producirá resultados más precisos para la derivada.

$$\tag 2 f(x) = \frac{1}{1+x^2}\,.$$ $$\tag 3 E_n = \max_{i\in \{0,\ldots,n\}} \bigg| f^\prime(x_i) - \sum_{j=1}^n D_{ij} f(x_j) \bigg|\,.$$ $$\tag 4 \tilde D_{jj} = D_{jj} - \left(\sum_{i=1}^n D_{ji} \right), \qquad \text{for all}\ j\,.$$

Conclusión

En conclusión, el método de Fornberg parece ser bastante preciso, en el caso de incluso en aproximadamente 3 órdenes de magnitud más precisos que los resultados de las fórmulas analíticas. Esto debería ser lo suficientemente preciso para la mayoría de las aplicaciones. Además, esto es notable porque Fornberg no parece incluir explícitamente este hecho en su método (al menos no se menciona en los dos documentos de Fornberg).$N=512$

Se puede obtener otro orden de magnitud para este ejemplo a través de una inclusión directa de la ecuación (4). Como este es un enfoque bastante simple y se aplica solo una vez para cada derivado, no veo ninguna razón para no usarlo.

El método del artículo de Baltensperger, que utiliza un enfoque más sofisticado para evaluar la suma en la ecuación (1) para reducir los errores de redondeo, produce aproximadamente el mismo orden de magnitud para el error. Entonces, al menos para este ejemplo, es más o menos equivalente al método "Fornberg ajustado" anterior.

Suponiendo que está tratando de diferenciar una implementación numérica de una función continua, hay una gran cantidad de métodos:

1) Diferenciación automática. El método más preciso y general. Doloroso de codificar, que requiere sobrecarga del operador y búsqueda dependiente de argumentos. Pone una carga en el usuario para comprender estos conceptos. También lucha con singularidades removibles, como diferenciar sinc en .$x= 0$

2) Una transformación de Chebyshev. Proyecte su función en un lapso de polinomios de Chebyshev y diferencie la recurrencia de tres términos. Súper rápido, muy preciso. Pero requiere que tenga un dominio de interés compatible de forma compacta; fuera del dominio seleccionado , la recurrencia de tres términos es inestable.$[a,b]$

3) Diferenciación finita. Subestimado en 1D; ver los consejos y trucos de Nick Higham en computación numérica . La idea es que si equilibra el error de truncamiento y el error de redondeo, no necesita seleccionar un tamaño de paso; Se puede elegir automáticamente. En Boost, esta idea se usa para recuperar (por defecto) 6/7 de los dígitos correctos para el tipo. (Higham solo muestra la idea para el caso más simple de la mitad de los dígitos correctos, pero la idea se extiende fácilmente.) Los coeficientes son de la tabla equiespaciada de Fornberg, pero el tamaño de pasos se elige bajo el supuesto de que la función se puede evaluar a 1ULP exactitud. La desventaja es que requiere 2 evaluaciones de funciones para recuperar la mitad de los dígitos del tipo, 4 para recuperar 3/4 de los dígitos, etc. En 1D, no es un mal negocio. En dimensiones superiores, es catastrófico.

4) La derivada del paso complejo. Use . Tome como redondeo de la unidad y esto se recuperará casi por completo. Sin embargo, es un poco engañoso, porque generalmente es más difícil implementar una función en el plano complejo que codificar manualmente su derivada real. Sigue siendo una idea genial y útil en ciertas circunstancias.$f'(x) \approx \Im(f(x+ih))$$h$

No estoy al tanto de que nadie haya mejorado el algoritmo de Fornberg (Véase también su poco más reciente de papel ). Por otro lado, me parece que mirar su algoritmo como una forma de calcular derivadas numéricas no es correcto. Todo lo que ha hecho es derivar un algoritmo eficiente para calcular los pesos de los métodos de diferencias finitas. La ventaja de su método es que le brinda los pesos de todos los derivados hasta el derivado deseado de una sola vez.

Un esquema mas simple

Además de mi otra respuesta, que trata más sobre una extensión del método Fornberg, abordaré aquí la pregunta para obtener alternativas más simples.

Para esto bosquejo un esquema alternativo que produce los coeficientes derivados de la interpolación lagrangiana más directamente. Su implementación requiere solo unas pocas líneas de código, funciona para cuadrículas arbitrarias y, según mis primeros experimentos, es tan precisa como la de Fornberg.

La base de la implementación es la derivada de paso imaginario donde es una variable en el orden de precisión de máquina. Se sabe que la derivada de paso imaginario produce valores derivados de forma estable y no sufre la inestabilidad numérica de una implementación de diferencia finita con .

$$f'(x) \ = \ \frac{1}{h}\,\text{Im}\big(f(x+ih)\big)\,,$$ $h$$h\rightarrow 0$

El segundo ingrediente es el polinomio de interpolación de Lagrange en la cuadrícula evaluado con una de las formas , por ejemplo, donde Para usar la derivada de pasos complejos, uno debe asegurarse de que estas fórmulas también funcionen para complejos argumentos . Además, para una función dada f (x) y un vector de coeficientes , denotamos el polinomio de interpolación a través de por $L_i(x)$$\{x_1,\ldots, x_N\}$

$$L_i(z) \ = \ \begin{cases}1 & \text{if }z = x_i\\\ \frac{\frac {\mu_i}{z-x_k}}{ \sum_k\frac {\mu_k}{z-x_k} } & \text{otherwise}\,. \end{cases}$$ $\mu_i = \frac{1}{\prod_{k\neq i} (x_i - x_k)}$$z$$f_i = f(x_i)$$(x_1,f_i)$ $$P(x;f) = \sum_{i=1}^N f_i L_i(x)\,.$$

---

Algoritmo

El algoritmo se bosqueja a continuación. Tiene los mismos parámetros de entrada y salida que los de Fornberg, pero es mucho más obvio.

Entrada:

• $x$ : una cuadrícula con N puntos de cuadrícula distintos
• $ord$ : orden derivado
• z: un punto en el que se debe evaluar la derivada
• Posiblemente: una función o sus valores de función en los puntos de la cuadrícula (solo se requiere para la variante de salida 2.)$f(x)$$f_i$

Inicialización

• Inicialice los polinomios de Lagrange mediante interpolación baricéntrica.$L$
• Inicialice una matriz de -matrices , para$N\times N$$D^{(o)}$$o=0, \ldots, ord$
• Establezca , es decir, la matriz unitaria de dimensión$D^{(0)} := I_N$$N\times N$
• Establecer$o := 0$

Algoritmo

Mientras :$o < ord$

• Calcule por la derivada de paso complejo para todo y . Aquí, denota la -ésima fila de .$D^{(o+1)}_{ik} = CSD \bigg( P\big(x_k; D^{(o)}_{i}\big) \bigg)$$CSD$$i$$k$
  $D^{(o)}_{i}$$i$$D^{(o)}$

• Establecer o = o + 1;

Decidir qué producir :

1. El vector de los coeficientes de diferencia finita en el punto , donde . Esto es lo que hace Fornberg.$d^{(ord)}$$z$$d_i = P(z;D^{(ord)}_i)$

2. La función de interpolación a la derivada de orden . Para esto, debe ingresar una función resp. la función valora en para el algoritmo.$p^{(ord)}(x) = \sum_{i=1}^N f(x_i) P(x;D^{(ord)}_i)$$f^{(ord)}(x)$$ord$$f_i$$x_i$

3. Una metafunción que devuelve la función de interpolación de la variante 2., pero para una función arbitraria que se debe interpolar en los puntos de la cuadrícula.$\text{diff}(\text{int } ord, \text{function } g)$$g$

Personalmente, me gusta más la variante 3.

---

Análisis del algoritmo.

Como el de Fornberg, este algoritmo es . Publicaré resultados más empíricos con respecto a la precisión, estabilidad, etc. si encuentro el tiempo.$\mathcal O(ord \cdot N^2)$

Para aumentar la precisión de la diferenciación numérica, haga lo siguiente:

1) Elija su método "estándar" favorito de alta precisión basado en algún tamaño de paso h .

2) Calcule el valor de la derivada con el método elegido en 1) muchas veces con tamaños de paso diferentes pero razonables h . Cada vez que puede elegir h como un número aleatorio del intervalo (0.5 * H / 10, 1.5 * H / 10) donde H es un tamaño de paso apropiado para el método que utiliza.

3) Promedio de los resultados.

Su resultado puede ganar 2-3 órdenes de magnitud en el error absoluto wrt. El resultado no promediado.

https://arxiv.org/abs/1706.10219