Тригонометрическая интерполяция

В математике тригонометрическая интерполяция — это интерполяция с помощью тригонометрических полиномов . Интерполяция — это процесс нахождения функции, которая проходит через некоторые заданные точки данных . Для тригонометрической интерполяции эта функция должна быть тригонометрическим полиномом, то есть суммой синусов и косинусов заданных периодов. Эта форма особенно подходит для интерполяции периодических функций .

Важным особым случаем является случай, когда заданные точки данных расположены на одинаковом расстоянии друг от друга; в этом случае решение дается дискретным преобразованием Фурье .

Тригонометрический полином степени К имеет вид

${\displaystyle p(x)=a_{0}+\sum _{k=1}^{K}a_{k} \cos(kx)+\sum _{k=1}^{K}b_{k}\sin(kx).\,}$

1

Это выражение содержит 2 K + 1 коэффициентов, a ₀ , a ₁ , … a _K , b ₁ , …, b _K , и мы хотим вычислить эти коэффициенты так, чтобы функция проходила через N точек:

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

Поскольку тригонометрический полином является периодическим с периодом 2π, N точек можно распределить и упорядочить в одном периоде следующим образом:

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

(Обратите внимание, что мы , как правило, не требуем, чтобы эти точки были расположены на одинаковом расстоянии друг от друга.) Задача интерполяции теперь состоит в том, чтобы найти коэффициенты, такие, чтобы тригонометрический полином p удовлетворял условиям интерполяции.

Задача становится более естественной, если мы сформулируем ее в комплексной плоскости . Мы можем переписать формулу для тригонометрического полинома как где i — мнимая единица . Если мы положим z = e ^ix , то это станет${\displaystyle p(x)=\sum _{k=-K}^{K}c_{k}e^{ikx},\,}$

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

с

${\displaystyle q(e^{ix})\треугольникq p(x).\,}$

Это сводит задачу тригонометрической интерполяции к задаче полиномиальной интерполяции на единичной окружности . Существование и единственность для тригонометрической интерполяции теперь немедленно следует из соответствующих результатов для полиномиальной интерполяции.

Более подробную информацию о формулировке тригонометрических интерполяционных полиномов в комплексной плоскости см. на стр. 156 книги «Интерполяция с использованием полиномов Фурье».

При указанных выше условиях существует решение задачи для любого заданного набора точек данных { x _k , y _k }, пока N , число точек данных, не больше числа коэффициентов в полиноме, т. е. N  ≤ 2 K +1 (решение может существовать или не существовать, если N >2 K +1, в зависимости от конкретного набора точек данных). Более того, интерполирующий полином является уникальным тогда и только тогда, когда число регулируемых коэффициентов равно числу точек данных, т. е. N  = 2 K  + 1. В оставшейся части этой статьи мы будем предполагать, что это условие выполняется.

Если число точек N нечетно, скажем, N=2K+1 , то применение формулы Лагранжа для полиномиальной интерполяции к полиномиальной формулировке в комплексной плоскости приводит к тому, что решение можно записать в виде

${\displaystyle p(x)=\sum _{k=0}^{2K}y_{k}\,t_{k}(x),}$

5

где

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

Множитель в этой формуле компенсирует тот факт, что формулировка комплексной плоскости содержит также отрицательные степени и, следовательно, не является полиномиальным выражением по . Правильность этого выражения можно легко проверить, заметив, что и что является линейной комбинацией правильных степеней . При использовании тождества${\displaystyle e^{-iKx+iKx_{k}}}$${\displaystyle e^{ix}}$${\displaystyle e^{ix}}$${\displaystyle t_{k}(x_{k})=1}$${\displaystyle t_{k}(x)}$${\displaystyle e^{ix}}$

${\displaystyle e^{iz_{1}}-e^{iz_{2}}=2i\sin {\tfrac {1}{2}}(z_{1}-z_{2})\,e^{(z_{1}+z_{2})i/2},}$

2

коэффициент можно записать в виде${\displaystyle t_{k}(x)}$

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

4

Если число точек N четное, скажем, N=2K , то применение формулы Лагранжа для полиномиальной интерполяции к полиномиальной формулировке в комплексной плоскости приводит к тому, что решение можно записать в виде

${\displaystyle p(x)=\sum _{k=0}^{2K-1}y_{k}\,t_{k}(x),}$

6

где

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

3

Здесь константы могут быть выбраны свободно. Это вызвано тем фактом, что интерполирующая функция ( 1 ) содержит нечетное число неизвестных констант. Обычный выбор заключается в требовании, чтобы наивысшая частота имела форму константы , умноженной на , т.е. член исчезает, но в общем случае фаза наивысшей частоты может быть выбрана как . Чтобы получить выражение для , мы получаем, используя ( 2 ), что ( 3 ) можно записать в виде${\displaystyle \альфа _{k}}$${\displaystyle \cos(Kx)}$${\displaystyle \sin(Kx)}$${\displaystyle \varphi _{K}}$${\displaystyle \альфа _{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})}}.}$

Это дает

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

и

${\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})}}.}$

Обратите внимание, что необходимо соблюдать осторожность, чтобы избежать бесконечностей, вызванных нулями в знаменателях.

Дальнейшее упрощение задачи возможно, если узлы равноудалены, т.е.${\displaystyle x_{m}}$

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

более подробную информацию см. в статье Зигмунда.

Дальнейшее упрощение с использованием ( 4 ) было бы очевидным подходом, но, очевидно, это запутано. Гораздо более простой подход заключается в рассмотрении ядра Дирихле

${\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}},}$

где нечетно. Легко видеть, что является линейной комбинацией правильных степеней и удовлетворяет${\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}}.}$

Поскольку эти два свойства однозначно определяют коэффициенты в ( 5 ), то отсюда следует, что${\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}}}$

Здесь sinc -функция предотвращает любые сингулярности и определяется как

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

Для четных мы определяем ядро ​​Дирихле как${\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}}.}$

Опять же, легко видеть, что является линейной комбинацией правых степеней , не содержит члена и удовлетворяет условию${\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}}.}$

Используя эти свойства, следует, что коэффициенты в ( 6 ) определяются как${\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}}}$

Обратите внимание, что не содержит также. Наконец, обратите внимание, что функция обращается в нуль во всех точках . Поэтому кратные этого члена всегда можно добавлять, но обычно его опускают.${\displaystyle t_{k}(x)}$${\displaystyle \sin {\tfrac {1}{2}}Nx}$${\displaystyle \sin {\tfrac {1}{2}}Nx}$${\displaystyle x_{m}}$

Реализацию вышеизложенного в MATLAB можно найти здесь, она выглядит следующим образом:

функция  P = triginterp ( xi,x,y ) % TRIGINTERP Тригонометрическая интерполяция. % Вход: % xi точки оценки для интерполянта (вектор) % x равноотстоящие узлы интерполяции (вектор, длина N) % y значения интерполяции (вектор, длина N) % Выход: % P значения тригонометрического интерполянта (вектор) N = длина ( x ); % Отрегулировать интервал заданной независимой переменной. h = 2 / N ; масштаб = ( x ( 2 ) - x ( 1 )) / h ; x = x / масштаб ; xi = xi / масштаб ; % Вычислить интерполянт. P = нули ( размер ( xi )); для k = 1 : N P = P + y ( k ) * trigcardinal ( xi - x ( k ), N ); конец                         function  tau = trigcardinal ( x,N ) ws = warning ( 'off' , 'MATLAB:divideByZero' ); % Форма отличается для четных и нечетных N. if rem ( N , 2 ) == 1 % odd tau = sin ( N * pi * x / 2 ) ./ ( N * sin ( pi * x / 2 )); else % even tau = sin ( N * pi * x / 2 ) ./ ( N * tan ( pi * x / 2 )); end warning ( ws ) tau ( x == 0 ) = 1 ; % fix value at x=0                    

Особый случай, когда точки x _n расположены на одинаковом расстоянии, особенно важен. В этом случае мы имеем

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

Преобразование, которое отображает точки данных y _n в коэффициенты a _k , b _{k ,} получается из дискретного преобразования Фурье (ДПФ) порядка 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}\,}$

(Из-за вышеизложенной формулировки задачи мы ограничились нечетным числом точек. Это не является строго необходимым; для четного числа точек включается еще один косинусный член, соответствующий частоте Найквиста .)

Случай интерполяции только косинусов для равномерно разнесенных точек, соответствующий тригонометрической интерполяции, когда точки имеют четную симметрию , был рассмотрен Алексисом Клеро в 1754 году. В этом случае решение эквивалентно дискретному косинусному преобразованию . Разложение только синусов для равномерно разнесенных точек, соответствующее нечетной симметрии, было решено Жозефом Луи Лагранжем в 1762 году, для которого решение представляет собой дискретное синусное преобразование . Полный косинусный и синусный интерполяционный полином, который дает начало ДПФ, был решен Карлом Фридрихом Гауссом в неопубликованной работе около 1805 года, и в этот момент он также вывел алгоритм быстрого преобразования Фурье для его быстрой оценки. Клеро, Лагранж и Гаусс были заняты изучением проблемы вывода орбиты планет , астероидов и т. д. из конечного набора точек наблюдения; поскольку орбиты являются периодическими, тригонометрическая интерполяция была естественным выбором. См. также Heideman et al. (1984).

Chebfun , полностью интегрированная программная система, написанная на MATLAB для вычислений с функциями, использует тригонометрическую интерполяцию и разложения Фурье для вычислений с периодическими функциями. Многие алгоритмы, связанные с тригонометрической интерполяцией, легко доступны в Chebfun ; несколько примеров доступны здесь.

• Кендалл Э. Аткинсон, Введение в численный анализ (2-е издание), Раздел 3.8. John Wiley & Sons, Нью-Йорк, 1988. ISBN  0-471-50023-2 .
• М. Т. Хайдеман, Д. Х. Джонсон и К. С. Беррус, «Гаусс и история быстрого преобразования Фурье», Журнал IEEE ASSP 1 (4), 14–21 (1984).
• GB Wright, M. Javed, H. Montanelli и LN Trefethen, "Расширение Chebfun до периодических функций", SIAM. J. Sci. Comput. , 37 (2015), C554-C573
• А. Зигмунд, Тригонометрические ряды , том II, глава X, Cambridge University Press, 1988.

• www.chebfun.org