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 determinados. 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 discreta de Fourier .

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 queremos calcular esos coeficientes para que la función pase por N puntos:

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

Dado que 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 .\,}$

(Obsérvese que, en general, no requerimos que estos puntos estén equidistantes.) El problema de interpolación ahora consiste en 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 de 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 ahora se derivan inmediatamente de los resultados correspondientes de la interpolación polinomial.

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

Solución del problema

Bajo 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 sólo 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 la interpolación polinomial 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 la 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 opción común es exigir que la frecuencia más alta tenga la forma a veces constante , es decir, el término desaparece, pero en general se puede elegir que la fase de la frecuencia más alta sea . Para obtener una expresión para , obtenemos usando ( 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})}}.}$

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

Nodos equidistantes

Es posible una mayor simplificación del 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 mayor simplificación mediante el uso de ( 4 ) sería un enfoque obvio, pero obviamente está involucrado. 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 extraño. Se puede ver fácilmente que es una combinación lineal de los poderes correctos 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 previene cualquier singularidad y está definida por

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

Número par de puntos

Incluso , 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}}}$

Tenga en cuenta que tampoco contiene el . Finalmente, observe que la función desaparece en todos los puntos . Por lo tanto, siempre se pueden agregar múltiplos de este término, pero comúnmente se omite. ${\displaystyle t_{k}(x)}$ ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ ${\displaystyle \sin {\tfrac {1}{2}}Nx}$ ${\displaystyle x_{m}}$

Implementación

Puede encontrar una implementación de MATLAB de lo anterior aquí y viene dada por:

función  P = triginterp ( xi,x,y ) % TRIGINTERP Interpolación trigonométrica. % Entrada: % puntos de evaluación xi para el interpolador (vector) % x nodos de interpolación equiespaciados (vector, longitud N) % valores de interpolación y (vector, longitud N) % Salida: % Valores P del interpolador trigonométrico (vector) N = longitud ( X ); % Ajusta el espaciado de la variable independiente dada. h = 2 / norte ; escala = ( x ( 2 ) - x ( 1 )) / h ; x = x / escala ; xi = xi / escala ; % Evaluar 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. if rem ( N , 2 ) == 1 % impar tau = sin ( N * pi * x / 2 ) ./ ( N * sin ( pi * x / 2 )); else % par tau = sin ( N * pi * x / 2 ) ./ ( N * tan ( pi * x / 2 )); advertencia final ( ws ) tau ( x == 0 ) = 1 ; % valor fijo en x=0                    

Relación con la transformada discreta de Fourier

Es especialmente importante el caso especial en el que los puntos x _n están igualmente espaciados. 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 discreta de Fourier (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 sólo coseno para puntos equidistantes, 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 equivale a una transformada discreta de coseno . La expansión solo seno para puntos equidistantes, correspondiente a simetría impar, fue resuelta por Joseph Louis Lagrange en 1762, cuya solución es una transformada 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 estaban todos preocupados por estudiar el problema de inferir la órbita de planetas , asteroides , etc., a partir de un conjunto finito de puntos de observación; Dado que las órbitas son periódicas, una interpolación trigonométrica fue la 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 calcular con funciones, utiliza interpolación trigonométrica y expansiones de Fourier para calcular 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 (segunda 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. Ciencias. Computadora. , 37 (2015), C554-C573
• A. Zygmund, Serie trigonométrica , volumen II, capítulo X, Cambridge University Press, 1988.

enlaces externos

• www.chebfun.org