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

p ( x ) = a 0 + ∑ k = 1 K a k cos ⁡ ( k x ) + ∑ k = 1 K b k sin ⁡ ( k x ) . {\displaystyle 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 , n = 0 , … , N − 1. {\displaystyle p(x_{n})=y_{n},\quad n=0,\ldots ,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 ≤ x 0 < x 1 < x 2 < … < x N − 1 < 2 π . {\displaystyle 0\leq x_{0}<x_{1}<x_{2}<\ldots <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 ) = ∑ k = − K K c k e i k x , {\displaystyle p(x)=\sum _{k=-K}^{K}c_{k}e^{ikx},\,} kde i je Imaginární jednotka. Pokud položíme z = eix, dostaneme tvar

q ( z ) = ∑ k = − K K c k z k , {\displaystyle q(z)=\sum _{k=-K}^{K}c_{k}z^{k},\,}

kde

q ( e i x ) ≜ p ( x ) . {\displaystyle 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.[1]

Ř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ů). 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

p ( x ) = ∑ k = 0 2 K y k t k ( x ) , {\displaystyle p(x)=\sum _{k=0}^{2K}y_{k}\,t_{k}(x),}

 

 

 

 

(5)

kde

t k ( x ) = e − i K x + i K x k ∏ m = 0 , m ≠ k 2 K e i x − e i x m e i x k − e i x m . {\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}}}}.}

Faktor e − i K x + i K x k {\displaystyle e^{-iKx+iKx_{k}}} je v tomto vzorci kvůli tomu, že formulace pro komplexní rovinu obsahuje i záporné mocniny e i x {\displaystyle e^{ix}} a proto není polynomiálním výrazem pro e i x {\displaystyle e^{ix}} . Korektnost tohoto výrazu lze snadno ověřit díky pozorování, že t k ( x k ) = 1 {\displaystyle t_{k}(x_{k})=1} a že t k ( x ) {\displaystyle t_{k}(x)} je lineární kombinací vhodných mocnin e i x {\displaystyle e^{ix}} . Využitím identity

e i z 1 − e i z 2 = 2 i sin ⁡ ( z 1 − z 2 2 ) e i 1 2 z 1 + i 1 2 z 2 , {\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)

lze koeficient t k ( x ) {\displaystyle t_{k}(x)} zapsat ve tvaru

t k ( x ) = ∏ m = 0 , m ≠ k 2 K sin ⁡ 1 2 ( x − x m ) sin ⁡ 1 2 ( x k − x m ) . {\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)

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

p ( x ) = ∑ k = 0 2 K − 1 y k t k ( x ) , {\displaystyle p(x)=\sum _{k=0}^{2K-1}y_{k}\,t_{k}(x),}

 

 

 

 

(6)

kde

t k ( x ) = e − i K x + i K x k e i x − e i α k e i x k − e i α k ∏ m = 0 , m ≠ k 2 K − 1 e i x − e i x m e i x k − e i x m . {\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)

Přičemž konstanty α k {\displaystyle \alpha _{k}} lze libovolně zvolit. To je způsobené faktem, že interpolační funkce (1) obsahuje lichý počet neznámých konstant. Běžnou volbou je vyžadovat, že nejvyšší frekvence měla tvar konstanta krát cos ⁡ ( K x ) {\displaystyle \cos(Kx)} , tj. aby člen sin ⁡ ( K x ) {\displaystyle \sin(Kx)} měl nulový koeficient, ale fázi nejvyšší frekvence lze obecně zvolit tak, aby byla φ K {\displaystyle \varphi _{K}} . Pro výpočet α k {\displaystyle \alpha _{k}} podle (2) je třeba (3) zapsat ve tvaru

t k ( x ) = cos ⁡ ( K x − 1 2 α k + 1 2 ∑ m = 0 , m ≠ k 2 K − 1 x m ) + ∑ m = − ( K − 1 ) K − 1 c k e i m x 2 N sin ⁡ ( x k − α k 2 ) ∏ m = 0 , m ≠ k 2 K − 1 sin ⁡ ( x k − x m 2 ) . {\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}})}}.}

Což dává

α k = ∑ m = 0 , m ≠ k 2 K − 1 x m − 2 φ K {\displaystyle \alpha _{k}=\sum _{m=0,m\neq k}^{2K-1}x_{m}-2\varphi _{K}}

a

t k ( x ) = sin ⁡ 1 2 ( x − α k ) sin ⁡ 1 2 ( x k − α k ) ∏ m = 0 , m ≠ k 2 K − 1 sin ⁡ 1 2 ( x − x m ) sin ⁡ 1 2 ( x k − x m ) . {\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})}}.}

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é,[2] pokud uzly x m {\displaystyle x_{m}} jsou ekvidistantní, tj.

x m = 2 π m N . {\displaystyle x_{m}={\frac {2\pi m}{N}}.}

Lichý počet bodů

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

D ( x , N ) = 1 N + 2 N ∑ k = 1 1 2 ( N − 1 ) cos ⁡ ( k x ) = sin ⁡ 1 2 N x N sin ⁡ 1 2 x , {\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}},}

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

D ( x m , N ) = { 0  pro  m ≠ 0 1  pro  m = 0 . {\displaystyle D(x_{m},N)={\begin{cases}0{\text{ pro }}m\neq 0\\1{\text{ pro }}m=0\end{cases}}.}

Protože tyto dvě vlastnosti definují koeficienty t k ( x ) {\displaystyle t_{k}(x)} v (5) jednoznačně, z toho plyne, že

t k ( x ) = D ( x − x k , N ) = { sin ⁡ 1 2 N ( x − x k ) N sin ⁡ 1 2 ( x − x k )  pro  x ≠ x k lim x → 0 sin ⁡ 1 2 N x N sin ⁡ 1 2 x = 1  pro  x = x k = s i n c 1 2 N ( x − x k ) s i n c 1 2 ( x − x k ) . {\displaystyle {\begin{aligned}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{ pro }}x\neq x_{k}\\\lim \limits _{x\to 0}{\frac {\sin {\frac {1}{2}}Nx}{N\sin {\frac {1}{2}}x}}=1{\text{ pro }}x=x_{k}\end{cases}}\\&={\frac {\mathrm {sinc} \,{\frac {1}{2}}N(x-x_{k})}{\mathrm {sinc} \,{\frac {1}{2}}(x-x_{k})}}.\end{aligned}}}

Funkce Sinc definovaná vztahem

s i n c x = sin ⁡ x x . {\displaystyle \mathrm {sinc} \,x={\frac {\sin x}{x}}.}

zabrání singularitě.

Sudý počet bodů

Pro N {\displaystyle N} sudé definujeme Dirichletovo jádro jako

D ( x , N ) = 1 N + 1 N cos ⁡ 1 2 N x + 2 N ∑ k = 1 1 2 N − 1 cos ⁡ ( k x ) = sin ⁡ 1 2 N x N tg ⁡ 1 2 x . {\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\operatorname {tg} {\frac {1}{2}}x}}.}

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

D ( x m , N ) = { 0  pro  m ≠ 0 1  pro  m = 0 . {\displaystyle D(x_{m},N)={\begin{cases}0{\text{ pro }}m\neq 0\\1{\text{ pro }}m=0\end{cases}}.}

Z těchto vlastností vyplývá, že koeficienty t k ( x ) {\displaystyle t_{k}(x)} v (6) jsou dány vzorcem

t k ( x ) = D ( x − x k , N ) = { sin ⁡ 1 2 N ( x − x k ) N tg ⁡ 1 2 ( x − x k )  pro  x ≠ x k lim x → 0 sin ⁡ 1 2 N x N tg ⁡ 1 2 x = 1  pro  x = x k . = s i n c 1 2 N ( x − x k ) s i n c 1 2 ( x − x k ) cos ⁡ 1 2 ( x − x k ) {\displaystyle {\begin{aligned}t_{k}(x)&=D(x-x_{k},N)={\begin{cases}{\frac {\sin {\frac {1}{2}}N(x-x_{k})}{N\operatorname {tg} {\frac {1}{2}}(x-x_{k})}}{\text{ pro }}x\neq x_{k}\\\lim \limits _{x\to 0}{\frac {\sin {\frac {1}{2}}Nx}{N\operatorname {tg} {\frac {1}{2}}x}}=1{\text{ pro }}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{aligned}}}

Přitom t k ( x ) {\displaystyle t_{k}(x)} neobsahuje sin ⁡ 1 2 N x {\displaystyle \sin {\frac {1}{2}}Nx} . Funkce sin ⁡ 1 2 N x {\displaystyle \sin {\frac {1}{2}}Nx} bude mít nulovou hodnotu také pro všechny body x m {\displaystyle 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:[3] 

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 xn jsou stejně vzdálené. V tomto případě platí

x n = 2 π n N , 0 ≤ n < N . {\displaystyle x_{n}=2\pi {\frac {n}{N}},\qquad 0\leq n<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 = ∑ n = 0 N − 1 y n   e − i 2 π n k N {\displaystyle Y_{k}=\sum _{n=0}^{N-1}y_{n}\ e^{-i2\pi {\frac {nk}{N}}}\,}
y n = p ( x n ) = 1 N ∑ k = 0 N − 1 Y k   e i 2 π n k N {\displaystyle 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.[4]

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.[5]

Odkazy

Reference

V tomto článku byl použit překlad textu z článku Trigonometric interpolation na anglické Wikipedii.

Literatura

  • ATKINSON, Kendall E., 1988. An Introduction to Numerical Analysis. 2. vyd. New York: John Wiley & Sons. Dostupné online. ISBN 0-471-50023-2. 
  • HEIDEMAN, M. T.; JOHNSON, D. H.; BURRUS, C. S., 1984. Gauss and the history of the fast Fourier transform. IEEE ASSP Magazine. Roč. 1, čís. 4, s. 14–21. Dostupné online. 
  • WRIGHT, G.B.; JAVED, M.; MONTANELLI, H.; TREFETHEN, L.N., 2015. Extension of Chebfun to periodic functions. SIAM. J. Sci. Comput.. S. C554-C573. Dostupné online.  Archivováno 15. 9. 2015 na Wayback Machine.
  • ZYGMUND, A., 1988. Trigonometric Series. Svazek II. [s.l.]: Cambridge University Press. Kapitola X. 
  • Interpolation using Fourier Polynomials [online]. University of Arizona. Dostupné online. 
  • WRIGHT, Grady B.; TREFETHEN, Lloyd N., 2014. Periodic Chebfuns [online]. Prosinec 2014, červenec 2019. Dostupné online. 
  • Dostupné online. 

Externí odkazy

  • www.chebfun.org