Тригонометрическая интерполяция - 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)

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

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

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

${ displaystyle 0 leq x_ {0}$

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

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

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

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

с

${ Displaystyle д (е ^ {ix}) треугольник р (х). ,}$

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

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

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

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

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

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

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

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

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

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

Используя эти свойства, следует, что коэффициенты ${ Displaystyle т_ {к} (х)}$ в (6) даны

${ 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 можно найти здесь и определяется:

функцияп =тригинтерп(xi, x, y)% TRIGINTERP Тригонометрическая интерполяция.% Вход:% xi баллов оценки для интерполянта (вектор)% x узлов интерполяции с равным интервалом (вектор, длина N)% y значения интерполяции (вектор, длина N)% Выход:% P значений тригонометрического интерполянта (вектора)N = длина(Икс);% Отрегулируйте интервал данной независимой переменной.час = 2/N;шкала = (Икс(2)-Икс(1)) / час;Икс = Икс/шкала;  xi = xi/шкала;% Оценить интерполянт.п = нули(размер(xi));за к = 1: N  п = п + у(k)*три кардинальный(xi-Икс(k),N);конецфункция tau = trigcardinal (x, N)WS = предупреждение('выключенный','MATLAB: divByZero');% Форма различна для четного и нечетного N.если rem(N,2)==1   % странный  тау = грех(N*число Пи*Икс/2) ./ (N*грех(число Пи*Икс/2));еще % четное  тау = грех(N*число Пи*Икс/2) ./ (N*загар(число Пи*Икс/2));конецпредупреждение (ws)тау(Икс==0) = 1;     % фиксированного значения при x = 0

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

Частный случай, когда точки Икс_п равномерно распределены, что особенно важно. В этом случае мы имеем

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

Преобразование, отображающее точки данных у_п к коэффициентам а_k, б_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 и другие. (1984).

Приложения в численных вычислениях

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

Рекомендации

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

внешняя ссылка

• www.chebfun.org