Cesko.wiki Trigonometrická interpolace
Author
Albert FloresTrigonometrická 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 {{Vzorec| p(x) = a_0 + \sum_{k=1}^K a_k \cos(kx) + \sum_{k=1}^K b_k \sin(kx). \, |1}} Tento výraz obsahuje 2K + 1 koeficientů, a0, a1, … aK, b1, …, bK, které je třeba vypočítat tak, aby funkce procházela N body: : p(x_n) = y_n, \quad n=0, \ldots, N-1. +more \, Protože trigonometrický polynom je periodický s periodou 2π, lze těchto N bodů rozdělit a uspořádat do jedné periody: : 0 \leq x_0 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. +more Pokud položíme z = eix, 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.
Řešení problému
Za výše uvedených podmínek existuje řešení problému pro jakoukoli danou množinu datových bodů {xk, yk}, 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ů). +more 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 {{Vzorec| p(x) = \sum_{k=0}^{2K} y_k\,t_k(x),|5}} kde : t_k(x) = e^{-iKx+iKx_k} \prod_{m=0,m\ne k}^{2K} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}. Faktor e^{-iKx+iKx_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}. +more 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 {{Vzorec|e^{iz_1}-e^{iz_2}=2i\sin\left(\frac{z_1-z_2}{2}\right)e^{i\frac12 z_1+i\frac12 z_2},|2}} lze koeficient t_k(x) zapsat ve tvaru {{Vzorec| t_k(x) = \prod_{m=0,m\ne k}^{2K} \frac{\sin\frac12(x-x_m)}{\sin\frac12(x_k-x_m)}. |4}}.
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 {{Vzorec| p(x) = \sum_{k=0}^{2K-1} y_k\,t_k(x),|6}} kde {{Vzorec| 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\ne k}^{2K-1} \frac{e^{ix}-e^{ix_m}}{e^{ix_k}-e^{ix_m}}. |3}} Přičemž konstanty \alpha_k lze libovolně zvolit. +more To je způsobené faktem, že interpolační funkce 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 \varphi_K. Pro výpočet \alpha_k podle je třeba zapsat ve tvaru : t_k(x) = \frac{\cos\left(Kx-\frac12\alpha_k+\frac12\sum\limits_{m=0,m\ne 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\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 \varphi_K a : t_k(x) = \frac{\sin\frac12(x-\alpha_k)}{\sin\frac12(x_k-\alpha_k)}\prod_{m=0,m\ne k}^{2K-1} \frac{\sin\frac12(x-x_m)}{\sin\frac12(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é, 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í , mnohem jednodušší přístup je uvažovat Dirichletovo jádro :D(x,N)=\frac{1}{N} +\frac{2}{N} \sum_{k=1}^{\frac12 (N-1)}\cos(kx) = \frac{\sin\frac12 Nx}{N\sin\frac12 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{cases}0\text{ pro } m\neq0 \\1\text{ pro } m=0\end{cases}. +more Protože tyto dvě vlastnosti definují koeficienty t_k(x) v jednoznačně, z toho plyne, že :\begin{align} t_k(x) &= D(x-x_k,N)=\begin{cases} \frac{\sin\frac12 N(x-x_k)}{N\sin\frac12 (x-x_k)} \text{ pro } x\neq x_k\\ \lim\limits_{x\to 0} \frac{\sin\frac12 Nx}{N\sin\frac12 x}=1 \text{ pro } x= x_k \end{cases}\\&= \frac{\mathrm{sinc}\,\frac12 N(x-x_k)}{\mathrm{sinc}\,\frac12 (x-x_k)}. \end{align} Funkce Sinc definovaná vztahem : \mathrm{sinc}\,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 \frac12 Nx + \frac{2}{N} \sum_{k=1}^{\frac12 N-1}\cos(kx) = \frac{\sin\frac12 Nx}{N\operatorname{tg}\frac12 x}. Opět je zřejmé, že D(x,N) je lineární kombinací vhodných mocnin e^{ix}, neobsahuje člen \sin \frac12 Nx a vyhovuje :D(x_m,N)=\begin{cases}0\text{ pro } m\neq0 \\1\text{ pro } m=0\end{cases}. +more Z těchto vlastností vyplývá, že koeficienty t_k(x) v jsou dány vzorcem :\begin{align} t_k(x) &= D(x-x_k,N)=\begin{cases} \frac{\sin\frac12 N(x-x_k)}{N\operatorname{tg}\frac12 (x-x_k)}\text{ pro } x\neq x_k\\ \lim\limits_{x\to 0} \frac{\sin\frac12 Nx}{N\operatorname{tg}\frac12 x}=1 \text{ pro } x= x_k. \end{cases}\\&= \frac{\mathrm{sinc}\,\frac12 N(x-x_k)}{ \mathrm{sinc}\,\frac12 (x-x_k)}\cos\frac12 (x-x_k) \end{align} Přitom t_k(x) neobsahuje \sin \frac12 Nx . Funkce \sin \frac12 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:
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é. +more 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) . +more/ (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 xn jsou stejně vzdálené. V tomto případě platí : x_n = 2 \pi \frac{n}{N}, \qquad 0 \leq n
Transformace převádějící datové body yn na koeficienty ak, bk je získaná z diskrétní Fourierovy transformace (DFT) řádu N.
: Y_k = \sum_{n=0}^{N-1} y_n \ e^{-i 2 \pi \frac{nk}{N}} \,
: y_n = p(x_n) = \frac{1}{N} \sum_{k=0}^{N-1} Y_k \ e^{i 2 \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. +more).
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í. +more 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.
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.
Odkazy
Reference
Literatura
Externí odkazy
[url=http://www.chebfun.org]www.chebfun.org[/url]