Interpolación trigonométrica

Interpolación con polinomios trigonométricos

En matemáticas , la interpolación trigonométrica es la interpolación con polinomios trigonométricos . La interpolación es el proceso de encontrar una función que pasa por algunos puntos de datos dados . Para la interpolación trigonométrica, esta función tiene que ser un polinomio trigonométrico, es decir, una suma de senos y cosenos de períodos dados. Esta forma es especialmente adecuada para la interpolación de funciones periódicas .

Un caso especial importante es cuando los puntos de datos dados están igualmente espaciados, en cuyo caso la solución viene dada por la transformada de Fourier discreta .

Formulación del problema de interpolación

Un polinomio trigonométrico de grado K tiene la forma

Esta expresión contiene 2 K + 1 coeficientes, a ₀ , a ₁ , … a _K , b ₁ , …, b _K , y deseamos calcular esos coeficientes para que la función pase por N puntos:

${\displaystyle p(x_{n})=y_{n},\quad n=0,\ldots ,N-1.\,}$

Como el polinomio trigonométrico es periódico con período 2π, los N puntos se pueden distribuir y ordenar en un período como

${\displaystyle 0\leq x_{0}<x_{1}<x_{2}<\ldots <x_{N-1}<2\pi .\,}$

(Tenga en cuenta que, en general, no requerimos que estos puntos estén igualmente espaciados). El problema de interpolación es ahora encontrar coeficientes tales que el polinomio trigonométrico p satisfaga las condiciones de interpolación.

Formulación en el plano complejo

El problema se vuelve más natural si lo formulamos en el plano complejo . Podemos reescribir la fórmula para un polinomio trigonométrico como donde i es la unidad imaginaria . Si establecemos z = e ^ix , entonces esto se convierte en ${\displaystyle p(x)=\sum _{k=-K}^{K}c_{k}e^{ikx},\,}$

${\displaystyle q(z)=\sum _{k=-K}^{K}c_{k}z^{k},\,}$

con

${\displaystyle q(e^{ix})\triangleq p(x).\,}$

Esto reduce el problema de la interpolación trigonométrica al de la interpolación polinómica en el círculo unitario . La existencia y unicidad de la interpolación trigonométrica se desprende ahora inmediatamente de los resultados correspondientes de la interpolación polinómica.

Para obtener más información sobre la formulación de polinomios de interpolación trigonométrica en el plano complejo, consulte la página 156 de Interpolación utilizando polinomios de Fourier.

Solución del problema

En las condiciones anteriores, existe una solución al problema para cualquier conjunto dado de puntos de datos { x _k , y _k } siempre que N , el número de puntos de datos, no sea mayor que el número de coeficientes en el polinomio, es decir, N  ≤ 2 K +1 (una solución puede existir o no si N >2 K +1 dependiendo del conjunto particular de puntos de datos). Además, el polinomio de interpolación es único si y solo si el número de coeficientes ajustables es igual al número de puntos de datos, es decir, N  = 2 K  + 1. En el resto de este artículo, asumiremos que esta condición es verdadera.

Número impar de puntos

Si el número de puntos N es impar, digamos N=2K+1 , al aplicar la fórmula de Lagrange para interpolación polinómica a la formulación polinómica en el plano complejo se obtiene que la solución se puede escribir en la forma

dónde

${\displaystyle t_{k}(x)=e^{-iKx+iKx_{k}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K}{\frac {e^{ix}-e^{ix_{m}}}{e^{ix_{k}}-e^{ix_{m}}}}.}$

El factor en esta fórmula compensa el hecho de que la formulación del plano complejo también contiene potencias negativas de y, por lo tanto, no es una expresión polinómica en . La exactitud de esta expresión se puede verificar fácilmente observando que y que es una combinación lineal de las potencias correctas de . Al usar la identidad ${\displaystyle e^{-iKx+iKx_{k}}}$ ${\displaystyle e^{ix}}$ ${\displaystyle e^{ix}}$ ${\displaystyle t_{k}(x_{k})=1}$ ${\displaystyle t_{k}(x)}$ ${\displaystyle e^{ix}}$

El coeficiente se puede escribir en la forma ${\displaystyle t_{k}(x)}$

Número par de puntos

Si el número de puntos N es par, digamos N=2K , al aplicar la fórmula de Lagrange para interpolación polinómica a la formulación polinómica en el plano complejo se obtiene que la solución se puede escribir en la forma

dónde

Aquí, las constantes se pueden elegir libremente. Esto se debe al hecho de que la función de interpolación ( 1 ) contiene un número impar de constantes desconocidas. Una elección común es requerir que la frecuencia más alta tenga la forma una constante por , es decir, el término se anula, pero en general la fase de la frecuencia más alta se puede elegir como . Para obtener una expresión para , obtenemos mediante ( 2 ) que ( 3 ) se puede escribir en la forma ${\displaystyle \alpha _{k}}$ ${\displaystyle \cos(Kx)}$ ${\displaystyle \sin(Kx)}$ ${\displaystyle \varphi _{K}}$ ${\displaystyle \alpha _{k}}$

${\displaystyle t_{k}(x)={\frac {\cos {\tfrac {1}{2}}{\Biggl (}2Kx-\alpha _{k}+\displaystyle \sum \limits _{m=0,\,m\neq k}^{2K-1}x_{m}{\Biggr )}+\sum \limits _{m=-(K-1)}^{K-1}c_{k}e^{imx}}{2^{N}\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})\displaystyle \prod \limits _{m=0,\,m\neq k}^{2K-1}\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

Esto produce

${\displaystyle \alpha _{k}=\sum _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}x_{m}-2\varphi _{K}}$

y

${\displaystyle t_{k}(x)={\frac {\sin {\tfrac {1}{2}}(x-\alpha _{k})}{\sin {\tfrac {1}{2}}(x_{k}-\alpha _{k})}}\prod _{\begin{aligned}m&=0\\[-4mu]m&\neq k\end{aligned}}^{2K-1}{\frac {\sin {\tfrac {1}{2}}(x-x_{m})}{\sin {\tfrac {1}{2}}(x_{k}-x_{m})}}.}$

Téngase en cuenta que se debe tener cuidado para evitar infinitos causados ​​por ceros en los denominadores.

Nodos equidistantes

Es posible simplificar aún más el problema si los nodos son equidistantes, es decir, ${\displaystyle x_{m}}$

${\displaystyle x_{m}={\frac {2\pi m}{N}},}$

Ver Zygmund para más detalles.

Número impar de puntos

Una simplificación adicional mediante el uso de ( 4 ) sería un enfoque obvio, pero obviamente implica mucho trabajo. Un enfoque mucho más simple es considerar el núcleo de Dirichlet

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}},}$

donde es impar. Se puede ver fácilmente que es una combinación lineal de las potencias correctas de y satisface ${\displaystyle N>0}$ ${\displaystyle D(x,N)}$ ${\displaystyle e^{ix}}$

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

Dado que estas dos propiedades definen de forma única los coeficientes en ( 5 ), se deduce que ${\displaystyle t_{k}(x)}$

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\sin {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\sin {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}.\end{aligned}}}$

Aquí, la función sinc evita cualquier singularidad y está definida por

${\displaystyle \mathrm {sinc} \,x={\frac {\sin x}{x}}.}$

Número par de puntos

Para pares, definimos el núcleo de Dirichlet como ${\displaystyle N}$

${\displaystyle D(x,N)={\frac {1}{N}}+{\frac {1}{N}}\cos {\tfrac {1}{2}}Nx+{\frac {2}{N}}\sum _{k=1}^{(N-1)/2}\cos(kx)={\frac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}.}$

Nuevamente, se puede ver fácilmente que es una combinación lineal de las potencias correctas de , no contiene el término y satisface ${\displaystyle D(x,N)}$ ${\displaystyle e^{ix}}$ ${\displaystyle \sin {\tfrac {1}{2}}Nx}$

${\displaystyle D(x_{m},N)={\begin{cases}0{\text{ for }}m\neq 0\\1{\text{ for }}m=0\end{cases}}.}$

Usando estas propiedades, se deduce que los coeficientes en ( 6 ) están dados por ${\displaystyle t_{k}(x)}$

${\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\dfrac {\sin {\tfrac {1}{2}}N(x-x_{k})}{N\tan {\tfrac {1}{2}}(x-x_{k})}}{\text{ for }}x\neq x_{k}\\[10mu]\lim \limits _{x\to 0}{\dfrac {\sin {\tfrac {1}{2}}Nx}{N\tan {\tfrac {1}{2}}x}}=1{\text{ for }}x=x_{k}.\end{cases}}\\&={\frac {\mathrm {sinc} \,{\tfrac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\tfrac {1}{2}}(x-x_{k})}}\cos {\tfrac {1}{2}}(x-x_{k})\end{aligned}}}$

Nótese que no contiene el también. Finalmente, nótese que la función se anula en todos los puntos . Por lo tanto, siempre se pueden sumar múltiplos de este término, pero normalmente se omite. ${\displaystyle t_{k}(x)}$ ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ ${\displaystyle x_{m}}$

Implementación

Una implementación MATLAB de lo anterior se puede encontrar aquí y está dada por:

función  P = triginterp ( xi,x,y ) % TRIGINTERP Interpolación trigonométrica. % Entrada: % xi puntos de evaluación para el interpolante (vector) % x nodos de interpolación equiespaciados (vector, longitud N) % y valores de interpolación (vector, longitud N) % Salida: % Valores P del interpolante trigonométrico (vector) N = longitud ( x ); % Ajusta el espaciado de la variable independiente dada. h = 2 / N ; escala = ( x ( 2 ) - x ( 1 )) / h ; x = x / escala ; xi = xi / escala ; % Evaluar el interpolante. P = ceros ( tamaño ( xi )); para k = 1 : N P = P + y ( k ) * trigcardinal ( xi - x ( k ), N ); fin                         función  tau = trigcardinal ( x, N ) ws = advertencia ( 'off' , 'MATLAB:divideByZero' ); % La forma es diferente para N par e impar. si rem ( N , 2 ) == 1 % impar tau = sin ( N * pi * x / 2 ) ./ ( N * sin ( pi * x / 2 )); de lo contrario % par tau = sin ( N * pi * x / 2 ) ./ ( N * tan ( pi * x / 2 )); fin advertencia ( ws ) tau ( x == 0 ) = 1 ; % fija el valor en x=0                    

Relación con la transformada de Fourier discreta

El caso especial en el que los puntos x _n están igualmente espaciados es especialmente importante. En este caso, tenemos

${\displaystyle x_{n}=2\pi {\frac {n}{N}},\qquad 0\leq n<N.}$

La transformación que asigna los puntos de datos y _n a los coeficientes a _k , b _k se obtiene de la transformada de Fourier discreta (DFT) de orden N.

${\displaystyle Y_{k}=\sum _{n=0}^{N-1}y_{n}\ e^{-i2\pi nk/N}\,}$

${\displaystyle y_{n}=p(x_{n})={\frac {1}{N}}\sum _{k=0}^{N-1}Y_{k}\ e^{i2\pi nk/N}\,}$

(Debido a la forma en que se formuló el problema anteriormente, nos hemos restringido a números impares de puntos. Esto no es estrictamente necesario; para números pares de puntos, se incluye otro término coseno correspondiente a la frecuencia de Nyquist ).

El caso de la interpolación de solo coseno para puntos igualmente espaciados, correspondiente a una interpolación trigonométrica cuando los puntos tienen simetría par , fue tratado por Alexis Clairaut en 1754. En este caso, la solución es equivalente a una transformada de coseno discreta . La expansión de solo seno para puntos igualmente espaciados, correspondiente a simetría impar, fue resuelta por Joseph Louis Lagrange en 1762, para la cual la solución es una transformada de seno discreta . El polinomio de interpolación de coseno y seno completo, que da lugar a la DFT, fue resuelto por Carl Friedrich Gauss en un trabajo inédito alrededor de 1805, momento en el que también derivó un algoritmo de transformada rápida de Fourier para evaluarlo rápidamente. Clairaut, Lagrange y Gauss se ocuparon todos de estudiar el problema de inferir la órbita de planetas , asteroides , etc., a partir de un conjunto finito de puntos de observación; Como las órbitas son periódicas, una interpolación trigonométrica fue una elección natural. Véase también Heideman et al. (1984).

Aplicaciones en computación numérica

Chebfun , un sistema de software totalmente integrado escrito en MATLAB para realizar cálculos con funciones, utiliza interpolación trigonométrica y expansiones de Fourier para realizar cálculos con funciones periódicas. Muchos algoritmos relacionados con la interpolación trigonométrica están disponibles en Chebfun ; varios ejemplos están disponibles aquí.

Referencias

• Kendall E. Atkinson, Introducción al análisis numérico (2.ª edición), Sección 3.8. John Wiley & Sons, Nueva York, 1988. ISBN  0-471-50023-2 .
• MT Heideman, DH Johnson y CS Burrus, "Gauss y la historia de la transformada rápida de Fourier", IEEE ASSP Magazine 1 (4), 14–21 (1984).
• GB Wright, M. Javed, H. Montanelli y LN Trefethen, "Extensión de Chebfun a funciones periódicas", SIAM. J. Sci. Comput. , 37 (2015), C554-C573
• A. Zygmund, Series trigonométricas , Volumen II, Capítulo X, Cambridge University Press, 1988.

Enlaces externos

• www.chebfun.org