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

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

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

Формулировка задачи интерполяции

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

$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 точек:

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

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

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

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

Формулировка в комплексной плоскости

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

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

с

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

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

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

Решение проблемы

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

Нечетное количество очков

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

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

где

$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}}}}.$

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

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

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

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

Четное количество очков

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

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

где

$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 )

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

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

Это дает

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

и

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

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

Эквидистантные узлы

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

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

см. Зигмунд для более подробной информации.

Нечетное количество очков

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

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

где $N>0$ странно. Легко можно увидеть, что $D(x,N)$ представляет собой линейную комбинацию правых степеней $e^{ix}$ и удовлетворяет

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

Поскольку эти два свойства однозначно определяют коэффициенты $t_{k}(x)$ в ( 5 ) следует, что

${\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 предотвращает любые сингулярности и определяется формулой

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

Четное количество очков

Для $N$ даже мы определяем ядро ​​Дирихле как

$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}}.$

Опять же, легко увидеть, что $D(x,N)$ представляет собой линейную комбинацию правых степеней $e^{ix}$, не содержит термина $\sin {\tfrac {1}{2}}Nx$ и удовлетворяет

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

Используя эти свойства, следует, что коэффициенты $t_{k}(x)$ в ( 6 ) имеют вид

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

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

Реализация

Реализацию вышеизложенного в MATLAB можно найти здесь :

function P = triginterp(xi,x,y)
% TRIGINTERP Trigonometric interpolation.
% Input:
%   xi  evaluation points for the interpolant (vector)
%   x   equispaced interpolation nodes (vector, length N)
%   y   interpolation values (vector, length N)
% Output:
%   P   values of the trigonometric interpolant (vector)
N = length(x);
% Adjust the spacing of the given independent variable.
h = 2/N;
scale = (x(2)-x(1)) / h;
x = x/scale;  xi = xi/scale;
% Evaluate interpolant.
P = zeros(size(xi));
for k = 1:N
  P = P + y(k)*trigcardinal(xi-x(k),N);
end

function tau = trigcardinal(x,N)
ws = warning('off','MATLAB:divideByZero');
% Form is different for even and odd 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

с дискретным

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

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

Преобразование, которое отображает точки данных y _n в коэффициенты a _k , b _k, получается из дискретного преобразования Фурье (ДПФ) порядка N.

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

$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. Джон Уайли и сыновья, Нью-Йорк, 1988 год. ISBN   0-471-50023-2 .
• М.Т. Хайдеман, Д.Х. Джонсон и К.С. Буррус, « Гаусс и история быстрого преобразования Фурье », журнал IEEE ASSP Magazine 1 (4), 14–21 (1984).
• ГБ Райт, М. Джавед, Х. Монтанелли и Л. Н. Трефетен, « Расширение Чебфуна на периодические функции », SIAM. Дж. Наук. Вычислить. , 37 (2015), С554-С573
• А. Зигмунд, Тригонометрический ряд , Том II, Глава X, Cambridge University Press, 1988.

Внешние ссылки

• www.chebfun.org