Programování pro fyziky 2021/22 – Numerické recepty č. 2

Náš průvodní text je především v poznámkách k přednášce. Pomáhat si budeme i odkazy do knih Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR), které jsou (v Second Edition) na webu volně k dispozici. Vidět je ovšem lze jen z Adobe Readeru po instalaci pluginu FileOpen. Interní PDF prohlížeče Firefoxu a Chromu nám nepomohou.

Hádanky

  1. Je polynom 10. stupně vhodnou aproximační funkcí na ekvidistantní síti?

  2. Je polynom 10. stupně vhodnou aproximační funkcí na jakékoliv síti?

  3. Jaká je časová složitost polynomické regrese? A lineární regrese?

Polynomická aproximace

(NR kap. 3) Úlohou je nalézt aproximační funkci, která přesně nebo přibližně vystihuje tabelované hodnoty \((x_i,f_i),\ i=0,\dots,n, x_i \ne x_k\) pro \(i\ne k\), nebo která nahrazuje složitější funkci, kterou lze tabelovat. Vyčíslování aproximační funkce v intervalu \(<x_0,x_n>\) se nazývá interpolace, mimo tento interval extrapolace. První a frekventovanou volbou pro aproximační funkci jsou polynomy. Úlohu požadující přesně vystihnout tabelované hodnoty (interpolační podmínky v interpolačních uzlech) pomocí jediného polynomu řeší polynomická interpolace, jindy se aproximuje po částech polynomy nižšího stupně spojitými (až hladkými) v interpolačních uzlech, jindy se jen minimalizují odchylky aproximované a aproximační funkce a hovoří se o polynomické regresi. Tyto aproximační úlohy vedou vesměs k řešení soustav lineárních algebraických rovnic; Octave to zakrývá funkcemi polyfit a interp1.

Polynomická interpolace jedním polynomem

(NR kap. 3.1, 3.5) Aproximační funkce je \(P_n(x)\), polynom obecně \(n\)-tého stupně pro \(n+1\) interpolačních podmínek. Existuje právě jeden takový polynom, zapsat jej je možné více způsoby.

Lagrangeův interpolační polynom (NR 3.1.1) \(L_n(x)=f_0 l_n^{(0)}(x)+\dots+f_n l_n^{(n)}(x)\) je tvořen váženým součtem elementárních Lagrangeových polynomů \(l_n^{(i)}(x)\) stupně \(n\), splňujících \(l_n^{(i)}(x_i)=1\) a \(l_n^{(i)}(x_k)=0\) pro \(i\ne k\), interpolační podmínky jsou tak zjevně splněny. Polynomy \(l_n^{(i)}(x)\) lze (opět zjevně) zapsat jako normované součiny kořenových činitelů,

\[l_n^{(i)}(x)=[(x-x_0)(x-x_1)...(x-x_{i-1})(x-x_{i+1})\dots(x-x_n)] / [(x_i-x_0)(x_i-x_1)...(x_i-x_{i-1})(x_i-x_{i+1})\dots(x_i-x_n)].\]

Pro interpolační polynom ve standardním tvaru (NR 3.5.2), \(P_n(x)=c_0 + c_1x + \dots + c_nx^n\), tvoří interpolační podmínky \(P_n(x_i)=f_i\) soustavu \(n+1\) lineárních algebraických rovnic pro koeficienty \(c_i\), \(\mathbf{V}\cdot\mathbf{c}=\mathbf{f}\),

\[\begin{split}\left(\begin{array}{cccc} 1 & x_0 & \dots & x_0^n\\ 1 & x_1 & \dots & x_1^n\\ \vdots & & & \vdots\\ 1 & x_n & \dots & x_n^n \end{array}\right) \left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_n \end{array}\right)= \left(\begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_n \end{array}\right).\end{split}\]

Pro tento speciální tvar (Vandermondovy) matice soustavy existuje úsporná metoda (Octave: funkce vander) s časovou složitostí \(O(n^2)\).

n=3; x=(1:n+1)'; f=ones(n+1,1); f(2:2:end)=-1;  % to jest x=[1;2;3;4;...]; f=[1;-1;1;-1;...];
A=[ones(n+1,1),x,x.^2,x.^3]; p=flipud(A\f)      % obecný algebraický řešič
% nebo
A=vander(x); p=A\f                              % speciální řešič vander
% nebo
p=polyfit(x,f,n)                                % fitační řešič polyfit
% a kreslení...
xi=0:.1:5; plot(xi,polyval(p,xi),x,f,'o')

Hledání koeficientů Newtonova interpolačního polynomu \(N_n(x) = a_0 + a_1 (x-x_0) + a_2 (x-x_0)(x-x_1) + ... + a_n (x-x_0)...(x-x_{n-1})\) vede k soustavě \(n+1\) lineárních algebraických rovnic s dolní trojúhelníkovou maticí, \(\mathbf{L}\cdot\mathbf{a}=\mathbf{f}\),

\[\begin{split}\left(\begin{array}{cccc} 1 & 0 & \dots & 0\\ 1 & x_1-x_0 & \dots & 0\\ \vdots & & & \vdots\\ 1 & x_n-x_0 & \dots & (x_n-x_0)(x_n-x_1)\dots(x_n-x_{n-1}) \end{array}\right) \left(\begin{array}{c} a_0 \\ a_1 \\ \vdots \\ a_n \end{array}\right)= \left(\begin{array}{c} f_0 \\ f_1 \\ \vdots \\ f_n \end{array}\right).\end{split}\]

Soustava je snadno řešitelná dopřednou substitucí s časovou složitostí \(O(n^2)\). Přidávání interpolačních podmínek nevede ke změně dosud získaných koeficientů \(a_i\) (vlastnost permanence).

n=3; x=(1:n+1)'; f=ones(n+1,1); f(2:2:end)=-1;  % to jest x=[1;2;3;4;...]; f=[1;-1;1;-1;...];
A=zeros(n+1);
for i=0:n, for j=0:i, A(i+1,j+1)=prod(x(i+1)-x(1:j)); end, end;
newton=A\f

Interpolace polynomem vyššího stupně je na ekvidistantních sítích nevhodná pro přílišné oscilace, zvláště u krajů (Rungeův fenomenon). Aproximuje se proto buď po částech polynomy nízkého stupně nebo na speciálních nerovnoměrných sítích tvořených kořeny ortogonálních (Čebyševových, Lagrangeových aj.) polynomů. Např. pro ortogonální Čebyševovy polynomy \(T_n(x)=\cos(n\arccos x)\) plyne z uvedeného vyjádření pomocí goniometrických funkcí poloha jejich kořenů ve tvaru \(x_k=\cos[\pi(k-\frac12)/n],\ k=1,\dots,n\).

hat=@(x) 1-sign(x).*x;

% polynomická aproximace na ekvidistantní síti...
n=10; x=-1:2/n:1; f=hat(x);
peq=polyfit(x,f,n);
hold off; xi=-1:.01:1; plot(xi,polyval(peq,xi),'-r',x,f,'dr')

% polynomická aproximace na síti tvořené kořeny OG Čebyševových polynomů...
n=11; x=cos(pi*((1:n)-0.5)/n); f=hat(x);
pog=polyfit(x,f,n);
hold on; plot(xi,polyval(pog,xi),'--b',x,f,'ob')

Polynomická interpolace po částech

Po částech lineární interpolace: Aproximační funkcí jsou po částech lineární polynomy spojité v interpolačních uzlech, jako při spojování podle pravítka. Časová složitost nalezení aproximované hodnoty odpovídá časové složitosti nalezení příslušného podintervalu, ideálně \(O(1)\)\(O(\log n)\).

% Octave...
n=20; x=linspace(0,1,n); f=rand(1,n);
dx=.001; xi=0:dx:1; fi=interp1(x,f,xi,'linear');
plot(xi,fi,x,f,'o')

# Gnuplot...
file='data'
plot file with linespoints

Interpolace kubickými spliny (NR kap. 3.3): Aproximační funkcí jsou po částech kubické polynomy spojité v interpolačních uzlech až do druhých derivací, jako při spojování podle křivítka. Funkce je na \(n\) podintervalech definována pomocí \(4n\) neznámých koeficientů; požadavek splnění interpolačních podmínek, spojitosti nulté, první a druhé derivace a 2 dodatečných podmínek tvoří \(4n\) rovnic, redukovatelných na soustavu \(n\) lineárních rovnic s třídiagonální maticí. Časová složitost sestrojení kubických splinů je tak lineární, \(O(n)\).

% Octave...
n=20; x=linspace(0,1,n); f=rand(1,n);
dx=.001; xi=0:dx:1; fi=interp1(x,f,xi,'spline');
plot(xi,fi,x,f,'o')

# Gnuplot...
file='data'
plot file smooth csplines, file with points

Polynomická regrese metodou nejmenších čtverců

Aproximační funkcí nechť je lineární kombinace \(m+1\) bázových funkcí, např. polynomů \(x^0,\dots,x^m\). Úmyslem není splnit \(n+1\) (obvykle \(n \gg m\)) interpolačních podmínek přesně (třeba pro jejich nepřesnost), ale minimalizovat odchylku hodnot aproximované a aproximační funkce. Odchylka se ve smyslu metody nejmenších čtverců definuje jako (vážený) součet druhých mocnin („čtverců”) rozdílů hodnot aproximované a aproximační funkce (\(L_2\) aproximace). Minimalizace odchylky vzhledem k \(m+1\) koeficientům aproximační funkce vyžaduje podle nich odchylku derivovat, což vede k soustavě \(m+1\) lineárních (tzv. normálních) rovnic pro koeficienty, \(\mathbf{G}\cdot\mathbf{c}=\mathbf{b}\). Při volbě lineárně nezávislých bázových funkcí je (Gramova) matice soustavy symetrická a pozitivně definitní, pro ortogonální bázové funkce je dokonce diagonální. Případu s polynomickými bázovými funkcemi se říká polynomická regrese, při \(m=1\) lineární regrese (proložení přímky).

% Octave...
m=1; n=3; x=[1;2;3;4]; f=[0;3;4;2];
p=polyfit(x,f,m)
xi=0:.1:5; plot(xi,polyval(p,xi),x,f,'o')

m=2;
p=polyfit(x,f,m)
plot(xi,polyval(p,xi),x,f,'o')

# Gnuplot...
f(x)=a*x+b
fit f(x) 'xf.dat' via a,b
plot 'xf.dat',f(x)

Odvození pro aproximační funkci \(\varphi(x)=\sum_{j=0}^m c_j\varphi_j(x)\). Minimalizuje se výraz \(\sum_{k=0}^n\omega_k\left[\varphi(x_k)-f_k\right]^2\) vzhledem k \(c_i\):

\[\begin{split}\begin{array}{l} \displaystyle\frac\partial{\partial c_i}\sum_{k=0}^n\omega_k\left[\sum_{j=0}^mc_j\varphi_j(x_k)-f_k\right]^2 = \displaystyle\sum_{k=0}^n2\omega_k\varphi_i(x_k)\left[\sum_{j=0}^mc_j\varphi_j(x_k)-f_k\right] = 0,\quad i=0,\dots,m,\\ \displaystyle\sum_{j=0}^m\left[\sum_{k=0}^n\omega_k\varphi_i(x_k)\varphi_j(x_k)\right]c_j = \displaystyle\sum_{k=0}^n\omega_kf_k\varphi_i(x_k),\quad i=0,\dots,m. \end{array}\end{split}\]

Pro polynomy \(\varphi_j(x)=x^j\) a \(\omega_k=1\):

\[\begin{split}\left(\begin{array}{llll} \sum_k 1 & \sum_k x_k & \dots & \sum_k x_k^m\\ \sum_k x_k & \sum_k x_k^2 & \dots & \sum_k x_k^{m+1}\\ \vdots & & & \vdots\\ \sum_k x_k^m & \sum_k x_k^{m+1} & \dots & \sum_k x_k^{2m} \end{array}\right) \left(\begin{array}{l} c_0 \\ c_1 \\ \vdots \\ c_m \end{array}\right)= \left(\begin{array}{l} \sum_k f_k \\ \sum_k f_kx_k \\ \vdots \\ \sum_k f_kx_k^m \end{array}\right) .\end{split}\]

Časová složitost řešení této soustavy s plnou maticí je kubická vzhledem k počtu hledaných koeficientů \(m\), který však má být mnohem menší než počet dat \(n\). Časová složitost sumací v prvcích matice i pravé strany je přitom lineární vzhledem k \(n\), a při pevném \(m\) tak lze ocenit časovou složitost polynomické regrese lineárním odhadem \(O(n)\).

Pro ortogonální Čebyševovy polynomy \(\varphi_j(x)=T_j(x)=\cos(j\arccos x)\), \(\omega_k=1\) a uzly \(x_k\) v kořenech \(\vphantom{\varphi_j(x)=T_j}T_{n+1}(x)\):

\[\begin{split}\left(\begin{array}{cccc} m+1 & 0 & \dots & 0\\ 0 & \frac{m+1}2 & \dots & 0\\ \vdots & & & \vdots\\ 0 & 0 & \dots & \frac{m+1}2 \end{array}\right) \left(\begin{array}{l} c_0 \\ c_1 \\ \vdots \\ c_m \end{array}\right)= \left(\begin{array}{l} \sum_k f_kT_0(x_k) \\ \sum_k f_kT_1(x_k) \\ \vdots \\ \sum_k f_kT_m(x_k) \end{array}\right),\quad \textstyle x_k=\cos[\pi(k+\frac12)/(n+1)],\ k=0,\dots,n,\end{split}\]

neboť ortogonalita zde znamená \(\sum_k T_i(x_k)T_j(x_k)=0\) pro \(i\ne j\).