Implementación de la interpolación cúbica monótona

Question

Estoy en la necesidad de implementar Interpolación cúbica monótona para interpolar una secuencia de puntos.
La información que tengo sobre los puntos es x,y y timestamp.

Soy mucho más un informático que un matemático, así que busco un ejemplo de implementación. Lo que necesito hacer con las funciones resultantes es almacenarlas para futuros análisis.

Mi sintaxis favorita es PHP, Python, Java, Delphi, VB o un lenguaje algorítmico genérico de su elección.

Answer 1

No es muy difícil de hacer; de hecho, incluso si te limitas a la interpolación Hermite, hay varios métodos. Voy a hablar de los tres con los que tengo más experiencia.

Recordemos que dados los puntos $(x_i,y_i),\quad i=1\dots n$ y asumiendo que no hay dos $x_i$ son iguales, se puede ajustar a los datos una interpolante de Hermite cúbica a trozos. (He dado la forma de la cúbica de Hermite en esta respuesta anterior .)

Para utilizar la notación de esa respuesta, ya tienes $x_i$ y $y_i$ y requieren una estimación del $y_i^\prime$ a partir de los datos dados. Hay al menos tres esquemas para hacer esto: Fritsch-Carlson, Steffen y Stineman.

(En las fórmulas siguientes, utilizo la notación $h_i=x_{i+1}-x_i$ y $d_i=\frac{y_{i+1}-y_i}{h_i}$ .)

El método de Fritsch-Carlson calcula una media armónica ponderada de las pendientes:

$$y_i^\prime = \begin{cases}3(h_{i-1}+h_i)\left(\frac{2h_i+h_{i-1}}{d_{i-1}}+\frac{h_i+2h_{i-1}}{d_i}\right)^{-1} &\text{ if }\mathrm{sign}(d_{i-1})=\mathrm{sign}(d_i)\\ 0&\text{ if }\mathrm{sign}(d_{i-1})\neq\mathrm{sign}(d_i)\end{cases}$$

el método de Steffen se basa en una media ponderada (alternativamente, un ajuste parabólico dentro del intervalo):

$$y_i^\prime = (\mathrm{sign}(d_{i-1})+\mathrm{sign}(d_i))\min\left(|d_{i-1}|,|d_i|,\frac12 \frac{h_i d_{i-1}+h_{i-1}d_i}{h_{i-1}+h_i}\right)$$

y el método de Stineman se ajusta a los círculos:

$$y_i^\prime = \frac{h_{i-1} d_{i-1}h_i^2(1+d_i^2)+h_i d_i h_{i-1}^2(1+d_{i-1}^2)}{h_{i-1} h_i^2(1+d_i^2)+h_i h_{i-1}^2(1+d_{i-1}^2)}$$

Las fórmulas que he dado sólo son aplicables a los puntos "internos"; tendrás que consultar esos documentos para las fórmulas de pendiente para manejar los puntos finales.

---

Como demostración de estos tres métodos, consideremos estos dos conjuntos de datos debidos a Akima :

$$\begin{array}{|l|l|} \hline x&y\\ \hline 1&10\\2&10\\3&10\\5&10\\6&10\\8&10\\9&10.5\\11&15\\12&50\\14&60\\15&95\\ \hline \end{array}$$

y Fritsch y Carlson :

$$\begin{array}{|l|l|} \hline x&y\\ \hline 7.99&0\\8.09&2.7629\times 10^{-5}\\8.19&4.37498\times 10^{-3}\\8.7&0.169183\\9.2&0.469428\\10&0.94374\\12&0.998636\\15&0.999919\\20&0.999994\\ \hline \end{array}$$

Aquí están los gráficos de estos dos conjuntos de datos:

Aquí están los gráficos de los ajustes de splines cúbicos a estos dos conjuntos:

Obsérvese la inestabilidad que no estaba presente en los datos originales; este es el precio que se paga por la continuidad de la segunda derivada de la que goza el spline cúbico.

A continuación se muestran los gráficos de los interpolantes utilizando los tres métodos mencionados anteriormente.

Este es el resultado de Fritsch-Carlson:

Este es el resultado de Steffen:

Este es el resultado de Stineman:

(El resultado no demasiado bueno para el conjunto de datos de Fritsch-Carlson podría deberse al uso de una interpolante cúbica de Hermite en lugar de la interpolante racional que Stineman recomendaba utilizar con su prescripción de derivadas).

---

Como he dicho, he tenido buena experiencia con estos tres; sin embargo, tendrás que experimentar en tu entorno sobre cuál de ellos se adapta mejor a tus necesidades.

Answer 2

NB: Si hubiera podido hacer este comentario, probablemente lo habría hecho.

Hay un Implementación de Haskell de interpolación cúbica natural, aunque no parece especificar la monotonicidad. Puede que merezca la pena mirarlo si tienes algo de experiencia al respecto. Si no es así, puede servir como punto de partida. También hay alguna discusión relativa a un motor de renderizado 3D (cuya discusión toca la variación de los métodos de interpolación, muy ligeramente) aquí . Ponerse en contacto con algunos de los implicados puede conducirle a mejores recursos.