Programování pro fyziky 2023/24 – Numerické recepty č. 2

Při výkladu si pomáháme odkazy do knih Numerical Recipes in C++, Numerical Recipes in C a Numerical Recipes in Fortran (krátce: NR) a také Numerical Computing with MATLAB, které jsou na webu volně k dispozici. Výklad provázejí ukázky skriptů pro Python s balíčky NumPy a SciPy a také pro Matlab a GNU Octave. V těchto prostředích je k demonstraci numerických metod krátká cesta.

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; Matlab a Octave, resp. NumPy a SciPy to obalují funkcemi polyfit, resp. numpy.polyfit a interp1, resp. scipy.interpolate.interp1d.

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 (funkce vander) s časovou složitostí \(O(n^2)\).

# Python: polynomická interpolace (řešení algebraické soustavy pro standardní polynom, také solver numpy.polyfit)
import numpy as np
import matplotlib.pyplot as plt
n=3                           # stupeň interpolačního polynomu
x=np.linspace(0,n,n+1)        # vstupní data: x=[0;1;2;3;...]
f=np.ones(n+1); f[1::2]=-1    #               f=[1;-1;1;-1;...]

# obecný algebraický solver
A=np.ones([n+1,n+1])
for nn in range(1,n+1): A[:,nn]=x**nn
p=np.flipud(np.linalg.solve(A,f)); print(p)

# speciální solver pro vander
A=np.vander(x)
p=np.linalg.solve(A,f); print(p)

# fitační solver polyfit
p=np.polyfit(x,f,n); print(p)

# a kreslení...
xi=np.linspace(-1,n+1,(n+2)*100+1)
plt.plot(xi,np.polyval(p,xi),x,f,'o')
plt.show()
% Octave: polynomická interpolace
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ý solver
% nebo
A=vander(x); p=A\f                              % speciální solver pro vander
% nebo
p=polyfit(x,f,n)                                % fitační solver 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).

% Octave: polynomická interpolace (řešení algebraické soustavy pro Newtonův polynom)
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 polynomů (Čebyševových, Lagrangeových aj.). 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=0,\dots,n-1\).

# Python: polynomická aproximace vyššího stupně
import numpy as np
import matplotlib.pyplot as plt
hat=lambda x:1-np.sign(x)*x    # interpolovaná funkce
n=10                           # stupeň interpolačního polynomu
xi=np.linspace(-1,1,201)       # síť pro vykreslení

# ... na ekvidistantní síti
x=np.linspace(-1,1,n+1)
f=hat(x)
peq=np.polyfit(x,f,n)
plt.plot(xi,np.polyval(peq,xi),'-r',x,f,'dr')

# ... na síti tvořené kořeny OG Čebyševových polynomů
x=np.cos(np.pi*(np.arange(n+1)+0.5)/(n+1))
f=hat(x)
pog=np.polyfit(x,f,n)
plt.plot(xi,np.polyval(pog,xi),'--b',x,f,'ob')
plt.show()
% Octave: polynomická aproximace vyššího stupně
hat=@(x) 1-sign(x).*x;
n=10;
xi=-1:.01:1;

% ... na ekvidistantní síti
x=linspace(-1,1,n+1); f=hat(x);
peq=polyfit(x,f,n);
hold off; plot(xi,polyval(peq,xi),'-r',x,f,'dr')

% ... na síti tvořené kořeny OG Čebyševových polynomů
x=cos(pi.^((0:n)+0.5)/(n+1)); 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)\).

# Python: po částech lineární interpolace (solver scipy.interpolate.interp1d 'linear')
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

n=20; x=np.linspace(0,1,n); f=np.random.rand(n)
ni=1000
fi=sp.interpolate.interp1d(x,f,'linear')
xi=np.linspace(0,1,ni)
plt.plot(xi,fi(xi),x,f,'o')
plt.show()
% Octave: po částech lineární interpolace (solver interp1 'linear')
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 (with lines, linespoints)
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)\).

# Python: po částech kubická interpolace (solver scipy.interpolate.interp1d 'cubic')
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

n=20; x=np.linspace(0,1,n); f=np.random.rand(n)
ni=1000
fi=sp.interpolate.interp1d(x,f,'cubic')
xi=np.linspace(0,1,ni)
plt.plot(xi,fi(xi),x,f,'o')
plt.show()
% Octave: po částech kubická interpolace (solver interp1 'spline')
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 (smooth csplines)
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).

# Python: polynomická regrese (solver numpy.polyfit)
import numpy as np
import matplotlib.pyplot as plt

n=3; x=[1,2,3,4]; f=[0,3,4,2]  # vstupní data
xi=np.linspace(0,5,51)         # síť pro vykreslení

m=1                            # lineární regrese
p=np.polyfit(x,f,m); print(p)
plt.plot(xi,np.polyval(p,xi),x,f,'o')

m=2                            # polynomická regrese kvadratickým polynomem
p=np.polyfit(x,f,m); print(p)
plt.plot(xi,np.polyval(p,xi),x,f,'o')
plt.show()
% Octave: polynomická regrese (solver polyfit)
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 (fit)
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\).

Na cvičení jsme kromě pythonského polyfitu zpracovali vlastnoruční interpolační a regresní polyfit v Pythonu a Fortranu.