Trigonometrická interpolace

Trigonometrická interpolace je v matematice interpolační metoda zvlášť vhodná pro interpolaci periodických funkcí. Cílem interpolace je nalezení funkce, která prochází zadanými datovými body. U trigonometrické interpolace se funkce interpoluje trigonometrickými polynomy, což je suma sinů a kosinů daných period.

Důležitým speciálním případem je, když dané datové body jsou stejně vzdálené; pak je řešení dáno diskrétní Fourierovou transformací.

Formulace problému interpolace

Trigonometrický polynom stupně K má tvar Šablona:Vzorec Tento výraz obsahuje 2K + 1 koeficientů, a₀, a₁, … a_K, b₁, …, b_K, které je třeba vypočítat tak, aby funkce procházela N body:

$p(x_n)=y_n,n=0,\dots ,N-1.$

Protože trigonometrický polynom je periodický s periodou 2π, lze těchto N bodů rozdělit a uspořádat do jedné periody:

$0\le x_0<x_1<x_2<\dots <x_{N-1}<2\pi .$

Přičemž obecně není požadováno, aby tyto body byly stejně vzdálené. Úkolem interpolace je pak nalézt takové koeficienty, aby trigonometrický polynom p vyhovoval interpolačním podmínkám.

Formulace v komplexní rovině

Nejpřirozenější je formulace trigonometrické interpolace v komplexní rovině. Vzorec pro trigonometrický polynom můžeme přepsat na tvar $p(x)=\sum _{k=-K}^K{c}_k{e}^{ikx},$ kde i je Imaginární jednotka. Pokud položíme z = e^ix, dostaneme tvar

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

kde

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

Tím je problém trigonometrické interpolace převeden na polynomiální interpolaci na jednotkové kružnici. Existence a jednoznačnost trigonometrické interpolace nyní okamžitě vyplývá z odpovídajících výsledků pro polynomiální interpolaci.

Více informací o formulaci trigonometrických interpolačních polynomů v komplexní rovině obsahuje text „Interpolation using Fourier Polynomials“ od strany 156.Šablona:Sfn

Řešení problému

Za výše uvedených podmínek existuje řešení problému pro jakoukoli danou množinu datových bodů {x_k, y_k}, pokud počet datových bodů N není větší než počet koeficientů v polynomu, tj. N ≤ 2K+1 (i pro N>2K+1 může existovat řešení, podle toho, jaká je konkrétní množina datových bodů). Interpolační polynom je jednoznačný právě tehdy, když počet nastavitelných koeficientů je roven počtu datových bodů, tj. N = 2K + 1. Ve zbytku tohoto článku budeme předpokládat, že je tato podmínka splněna.

Lichý počet bodů

Pokud je počet bodů N lichý (N=2K+1), použití Lagrangeva vzorce pro polynomiální interpolaci na polynomiální formulaci v komplexní rovině dostaneme řešení ve tvaru Šablona:Vzorec kde

$t_k(x)=e^{-iKx+iK{x}_k}\prod _{m=0,m\ne k}^{2K}\frac{e^{ix}-e^{i{x}_m}}{e^{i{x}_k}-e^{i{x}_m}}.$

Faktor $e^{-iKx+iK{x}_k}$ je v tomto vzorci kvůli tomu, že formulace pro komplexní rovinu obsahuje i záporné mocniny $e^{ix}$ a proto není polynomiálním výrazem pro $e^{ix}$. Korektnost tohoto výrazu lze snadno ověřit díky pozorování, že $t_k(x_k)=1$ a že $t_k(x)$ je lineární kombinací vhodných mocnin $e^{ix}$. Využitím identity Šablona:Vzorec lze koeficient $t_k(x)$ zapsat ve tvaru Šablona:Vzorec

Sudý počet bodů

Pokud počet bodů N je sudý (N=2K), použití Lagrangeova vzorce pro polynomiální interpolaci na polynomiální formulaci v komplexní rovině dostaneme řešení ve tvaru Šablona:Vzorec kde Šablona:Vzorec Přičemž konstanty ${\alpha }_k$ lze libovolně zvolit. To je způsobené faktem, že interpolační funkce (Šablona:Odkaz na vzorec) obsahuje lichý počet neznámých konstant. Běžnou volbou je vyžadovat, že nejvyšší frekvence měla tvar konstanta krát $\cos (Kx)$, tj. aby člen $\sin (Kx)$ měl nulový koeficient, ale fázi nejvyšší frekvence lze obecně zvolit tak, aby byla ${\phi }_K$. Pro výpočet ${\alpha }_k$ podle (Šablona:Odkaz na vzorec) je třeba (Šablona:Odkaz na vzorec) zapsat ve tvaru

$t_k(x)=\frac{\cos (Kx-\frac{1}{2}{\alpha }_k+\frac{1}{2}\sum _{m=0,m\ne k}^{2K-1}{x}_m)+\sum _{m=-(K-1)}^{K-1}{c}_k{e}^{imx}}{2^N\sin (\frac{x_k-{\alpha }_k}{2})\prod _{m=0,m\ne k}^{2K-1}\sin (\frac{x_k-x_m}{2})}.$

Což dává

${\alpha }_k=\sum _{m=0,m\ne k}^{2K-1}{x}_m-2{\phi }_K$

a

$t_k(x)=\frac{\sin \frac{1}{2}(x-{\alpha }_k)}{\sin \frac{1}{2}(x_k-{\alpha }_k)}\prod _{m=0,m\ne k}^{2K-1}\frac{\sin \frac{1}{2}(x-x_m)}{\sin \frac{1}{2}(x_k-x_m)}.$

Přitom je třeba se vyhnout nekonečnům způsobeným nulovými hodnotami ve jmenovateli.

Ekvidistantní uzly

Další zjednodušení problému je možné,Šablona:Sfn pokud uzly $x_m$ jsou ekvidistantní, tj.

$x_m=\frac{2\pi m}{N}.$

Lichý počet bodů

Přestože je možné další zjednodušení pomocí (Šablona:Odkaz na vzorec), mnohem jednodušší přístup je uvažovat Dirichletovo jádro

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

kde $N>0$ je liché. Je zřejmé, že $D(x,N)$ je lineární kombinací vhodných mocnin $e^{ix}$ a platí

$D(x_m,N)=\{\begin{array}{c}0\text{ pro }m\ne 0\\ 1\text{ pro }m=0\end{array}.$

Protože tyto dvě vlastnosti definují koeficienty $t_k(x)$ v (Šablona:Odkaz na vzorec) jednoznačně, z toho plyne, že

$\begin{array}{cc}{t}_k(x)& =D(x-x_k,N)=\{\begin{array}{c}\frac{\sin \frac{1}{2}N(x-x_k)}{N\sin \frac{1}{2}(x-x_k)}\text{ pro }x\ne x_k\\ \lim {\mathrm{\backslash limits}}_{x\to 0}\frac{\sin \frac{1}{2}Nx}{N\sin \frac{1}{2}x}=1\text{ pro }x=x_k\end{array}\\ & =\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\frac{1}{2}N(x-x_k)}{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\frac{1}{2}(x-x_k)}.\end{array}$

Funkce Sinc definovaná vztahem

$\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}x=\frac{\sin x}{x}.$

zabrání singularitě.

Sudý počet bodů

Pro $N$ sudé definujeme Dirichletovo jádro jako

$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\mathrm{tg}\frac{1}{2}x}.$

Opět je zřejmé, že $D(x,N)$ je lineární kombinací vhodných mocnin $e^{ix}$, neobsahuje člen $\sin \frac{1}{2}Nx$ a vyhovuje

$D(x_m,N)=\{\begin{array}{c}0\text{ pro }m\ne 0\\ 1\text{ pro }m=0\end{array}.$

Z těchto vlastností vyplývá, že koeficienty $t_k(x)$ v (Šablona:Odkaz na vzorec) jsou dány vzorcem

$\begin{array}{cc}{t}_k(x)& =D(x-x_k,N)=\{\begin{array}{c}\frac{\sin \frac{1}{2}N(x-x_k)}{N\mathrm{tg}\frac{1}{2}(x-x_k)}\text{ pro }x\ne x_k\\ \lim {\mathrm{\backslash limits}}_{x\to 0}\frac{\sin \frac{1}{2}Nx}{N\mathrm{tg}\frac{1}{2}x}=1\text{ pro }x=x_k.\end{array}\\ & =\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\frac{1}{2}N(x-x_k)}{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\frac{1}{2}(x-x_k)}\cos \frac{1}{2}(x-x_k)\end{array}$

Přitom $t_k(x)$ neobsahuje $\sin \frac{1}{2}Nx$. Funkce $\sin \frac{1}{2}Nx$ bude mít nulovou hodnotu také pro všechny body $x_m$, proto se násobky tohoto členu obvykle vynechávají.

Implementace

Na webu University of Delaware byla zveřejněna tato implementace trigonometrické interpolace pro MATLAB:Šablona:Sfn

function P = triginterp(xi,x,y)
% TRIGINTERP Trigonometrická interpolace.
% Vstup:
%   xi  vektor bodů, v nichž se vyhodnocuje funkce
%   x   ekvidistantní interpolační body (vektor délky N)
%   y   interpolační hodnoty (vektor délky N)
% Výstup:
%   P   vektor hodnot trigonometrického interpolantu
N = length(x);
% Nastavení odstup nezávislé proměnné.
h = 2/N;
scale = (x(2)-x(1)) / h;
x = x/scale;  xi = xi/scale;
% Vyhodnocení interpolantu.
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');
% Výpočet je různý pro sudá a lichá N.
if rem(N,2)==1   % liché
  tau = sin(N*pi*x/2) ./ (N*sin(pi*x/2));
else             % sudé
  tau = sin(N*pi*x/2) ./ (N*tan(pi*x/2));
end
warning(ws)
tau(x==0) = 1;     % oprava hodnotu pro x=0

Spojitost s diskrétní Fourierovou transformací

Zvláště důležitý je případ, kdy body x_n jsou stejně vzdálené. V tomto případě platí

$x_n=2\pi \frac{n}{N},0\le n<N.$

Transformace převádějící datové body y_n na koeficienty a_k, b_k je získaná z diskrétní Fourierovy transformace (DFT) řádu N.

$Y_k=\sum _{n=0}^{N-1}{y}_n{e}^{-i2\pi \frac{nk}{N}}$

$y_n=p(x_n)=\frac{1}{N}\sum _{k=0}^{N-1}{Y}_k{e}^{i2\pi \frac{nk}{N}}$

(Kvůli tomu, jak byl výše problém formulován, je třeba se omezit na lichý počet bodů. To není nezbytně nutné, ale pro sudý počet bodů je třeba zahrnout další kosinový člen odpovídající Nyquistově frekvenci.)

Případ interpolace používající pouze kosiny pro stejně vzdálené body odpovídající trigonometrické interpolaci se sudou symetrie bodů zkoumal Alexis Clairaut v roce 1754. Řešení je v tomto případě ekvivalentní s diskrétní kosinovou transformací. Rozvoj pro stejně vzdálené body obsahující pouze siny odpovídající liché symetrii vyřešil v roce 1762 Joseph-Louis Lagrange; řešením je v tomto případě diskrétní sinová transformace. Úplný interpolační polynom s kosiny i siny, který odpovídá DFT, vyřešil okolo roku 1805 Carl Friedrich Gauss v nepublikované práci, v níž také odvodil algoritmus rychlé Fourierovy transformace. Clairaut, Lagrange a Gauss pracovali na výpočtech parametrů oběžných drah planet, asteroidů, atd. z konečného počtu pozorování; protože oběžné dráhy jsou periodické, byla trigonometrická interpolace přirozenou volbou.Šablona:Sfn

Použití v numerických výpočtech

Chebfun je plně integrovaný softwarový systém pro výpočty s funkcemi napsaný v MATLABu používající trigonometrickou interpolaci a Fourierovy rozvoje pro výpočty s periodickými funkcemi. V Chebfun je okamžitě dostupné množství algoritmů týkajících se trigonometrické interpolace; několik příkladů je dostupných v popisu výpočtů s periodickými funkcemi.Šablona:Sfn

Odkazy

Literatura

• Šablona:Citace monografie
• Šablona:Citace periodika
• Šablona:Citace periodika Šablona:Wayback
• Šablona:Citace monografie
• Šablona:Citace elektronické monografie
• Šablona:Citace elektronické monografie
• Šablona:Citace elektronické monografie

Externí odkazy

• www.chebfun.org