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

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

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

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

Тригонометрический полином степени 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 можно найти здесь :

функция   P   =   triginterp  (  xi,x,y  ) 
 % TRIGINTERP Тригонометрическая интерполяция. 
  % Входные данные: 
 % точек оценки xi для интерполянта (вектора) 
 % x равноотстоящие узлы интерполяции (вектор, длина N) 
 % значения интерполяции y (вектор, длина N) 
 % Выходные данные: 
 % P значений тригонометрического интерполянта (вектор) 
 N   =   длина  (  Икс  ); 
  % Отрегулируйте интервал данной независимой переменной. 
  ч   =   2  /  Н  ; 
  масштаб   =   (  Икс  (  2  )  -  Икс  (  1  ))   /   час  ; 
  х   =   х  /  масштаб  ;     xi   =   xi  /  масштаб  ; 
  % Оценка интерполянта. 
  P   =   нули  (  размер  (  xi  )); 
  для   k   =   1  :  N 
   P   =   P   +   y  (  k  )  *  трикардинальный  (  xi  -  x  (  k  ),  N  ); 
  конечная 

 функция   tau   =   trigcardinal  (  x,N  ) 
 ws   =   предупреждение  (  'off'  ,  'MATLAB:divideByZero'  ); 
  % Форма различна для четного и нечетного N. 
 if   rem  (  N  ,  2  )  ==  1     % нечетный 
   tau   =   sin  (  N  *  pi  *  x  /  2  )   ./   (  N  *  sin  (  pi  *  x  /  2  )); 
  else               %even 
   tau   =   sin  (  N  *  pi  *  x  /  2  )   ./   (  N  *  tan  (  pi  *  x  /  2  )); 
  конец 
 предупреждения  (  ws  ) 
 тау  (  x  ==  0  )   =   1  ;        % фиксированного значения при 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