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

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

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

Постановка задачи интерполяции.

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

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

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

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

Проблема становится более естественной, если мы сформулируем ее в комплексной плоскости . Мы можем переписать формулу тригонометрического полинома в виде где 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 д (е ^ {ix}) \ треугольник р (х). \,}$

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

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

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

При вышеуказанных условиях существует решение проблемы для любого заданного набора точек данных { 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 _ {m = 0, m \ neq k} ^ {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 т_ {к} (х_ {к}) = 1}$${\ Displaystyle т_ {к} (х)}$${\ displaystyle e ^ {ix}}$

${\ displaystyle e ^ {iz_ {1}} - e ^ {iz_ {2}} = 2i \ sin \ left ({\ frac {z_ {1} -z_ {2}} {2}} \ right) e ^ {i {\ frac {1} {2}} z_ {1} + i {\ frac {1} {2}} z_ {2}},}$

( 2 )

коэффициент можно записать в виде ${\ Displaystyle т_ {к} (х)}$

${\ displaystyle t_ {k} (x) = \ prod _ {m = 0, m \ neq k} ^ {2K} {\ frac {\ sin {\ frac {1} {2}} (x-x_ {m })} {\ sin {\ frac {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 _ {m = 0, m \ neq k} ^ {2K-1} {\ frac {e ^ {ix} -e ^ {ix_ {m}}} {e ^ {ix_ {k}} - e ^ {ix_ {m}}}}.}$

( 3 )

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

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

Это дает

${\ displaystyle \ alpha _ {k} = \ sum _ {m = 0, m \ neq k} ^ {2K-1} x_ {m} -2 \ varphi _ {K}}$

а также

${\ displaystyle t_ {k} (x) = {\ frac {\ sin {\ frac {1} {2}} (x- \ alpha _ {k})} {\ sin {\ frac {1} {2} } (x_ {k} - \ alpha _ {k})}} \ prod _ {m = 0, m \ neq k} ^ {2K-1} {\ frac {\ sin {\ frac {1} {2} } (x-x_ {m})} {\ sin {\ frac {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} ^ {{\ frac {1} {2} } (N-1)} \ cos (kx) = {\ frac {\ sin {\ frac {1} {2}} Nx} {N \ sin {\ frac {1} {2}} x}},}$

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

${\ displaystyle {\ begin {align} t_ {k} (x) & = D (x-x_ {k}, N) = {\ begin {cases} {\ frac {\ sin {\ frac {1} {2} }} N (x-x_ {k})} {N \ sin {\ frac {1} {2}} (x-x_ {k})}} {\ text {for}} x \ neq x_ {k} \\\ lim \ limits _ {x \ to 0} {\ frac {\ sin {\ frac {1} {2}} Nx} {N \ sin {\ frac {1} {2}} x}} = 1 {\ text {for}} x = x_ {k} \ end {cases}} \\ & = {\ frac {\ mathrm {sinc} \, {\ frac {1} {2}} N (x-x_ { k})} {\ mathrm {sinc} \, {\ frac {1} {2}} (x-x_ {k})}}. \ end {выравнивается}}}$

Здесь sinc -функция предотвращает возникновение сингулярностей и определяется формулой

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

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

Для четных определим ядро ​​Дирихле как ${\ displaystyle N}$

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

Опять же, легко увидеть, что это линейная комбинация правильных степеней , не содержит члена и удовлетворяет ${\ Displaystyle D (х, N)}$${\ displaystyle e ^ {ix}}$${\ displaystyle \ sin {\ frac {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 т_ {к} (х)}$

${\ displaystyle {\ begin {align} t_ {k} (x) & = D (x-x_ {k}, N) = {\ begin {cases} {\ frac {\ sin {\ frac {1} {2} }} N (x-x_ {k})} {N \ tan {\ frac {1} {2}} (x-x_ {k})}} {\ text {for}} x \ neq x_ {k} \\\ lim \ limits _ {x \ to 0} {\ frac {\ sin {\ frac {1} {2}} Nx} {N \ tan {\ frac {1} {2}} x}} = 1 {\ text {for}} x = x_ {k}. \ end {cases}} \\ & = {\ frac {\ mathrm {sinc} \, {\ frac {1} {2}} N (x-x_ {k})} {\ mathrm {sinc} \, {\ frac {1} {2}} (x-x_ {k})}} \ cos {\ frac {1} {2}} (x-x_ { k}) \ end {выровненный}}}$

Обратите внимание, что он также не содержит . Наконец, обратите внимание, что функция обращается в нуль во всех точках . Следовательно, всегда можно добавить кратные этого члена, но обычно его опускают. ${\ Displaystyle т_ {к} (х)}$${\ displaystyle \ sin {\ frac {1} {2}} Nx}$${\ displaystyle \ sin {\ frac {1} {2}} Nx}$${\ displaystyle 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

Связь с дискретным преобразованием Фурье

Особый случай, когда точки 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 {\ frac {nk} {N}}} \,}$

${\ displaystyle y_ {n} = p (x_ {n}) = {\ frac {1} {N}} \ sum _ {k = 0} ^ {N-1} Y_ {k} \ e ^ {i2 \ пи {\ frac {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 Magazine 1 (4), 14–21 (1984).
• Г.Б. Райт, М. Джавед, Х. Монтанелли и Л.Н. Трефетен, " Расширение Chebfun на периодические функции ", SIAM. J. Sci. Comput. , 37 (2015), C554-C573
• А. Зигмунд, Тригонометрические ряды , Том II, Глава X, Cambridge University Press, 1988.

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

• www.chebfun.org