Modern numerikus módszerek

Bene Gyula
Eötvös Loránd Tudományegyetem, Elméleti Fizikai Tanszék
1117 Budapest, Pázmány Péter sétány 1/A
5. hét
Interpoláció egy dimenzióban: kifejtés ortogonális polinomok szerint

Ortogonális polinomok
Skalárszorzat $W(x)\ge 0$ súlyfüggvénnyel: $$< f |g>=\int_a^b W(x)f(x)g(x)dx$$ Ortogonális polinomok rekurziós relációja: $$\begin{eqnarray} p_{-1}(x)&=&0\\ p_{0}(x)&=&1\\ p_{j+1}(x)&=&(x-a_j)p_j(x)-b_jp_{j-1}(x) \end{eqnarray}$$ Itt $$\begin{eqnarray} a_j&=&\frac{< xp_j |p_j>}{< p_j |p_j>}\\ b_j&=&\frac{< p_j |p_j>}{< p_{j-1} |p_{j-1} >} \end{eqnarray}$$ Christoffel-Darboux-azonosság:
```
               A rekurziós formulából:
```
$$\begin{eqnarray} xp_{j}(x)&=& p_{j+1}(x)+a_jp_j(x)+b_jp_{j-1}(x)\\ yp_{j}(y)&=& p_{j+1}(y)+a_jp_j(y)+b_jp_{j-1}(y) \end{eqnarray}$$
```
               Ebből következik:
```
$$(x-y)p_{j}(x)p_{j}(y)=p_{j+1}(x)p_{j}(y)-p_{j}(x)p_{j+1}(y)+b_j\left(p_{j-1}(x)p_{j}(y)-p_{j}(x)p_{j-1}(y)\right)$$
```
               Másképpen:
```
$$\frac{(x-y)p_{j}(x)p_{j}(y)}{< p_j |p_j>}=\frac{p_{j+1}(x)p_{j}(y)-p_{j}(x)p_{j+1}(y)}{< p_j |p_j>}-\frac{p_{j}(x)p_{j-1}(y)-p_{j-1}(x)p_{j}(y)}{< p_{j-1} |p_{j-1} >}$$
```
               Összegzés j-re 0-tól n-ig:
```
$$\sum_{j=0}^n\frac{p_{j}(x)p_{j}(y)}{< p_j |p_j>}=\frac{p_{n+1}(x)p_{n}(y)-p_{n}(x)p_{n+1}(y)}{< p_n |p_n>(x-y)}$$ Zérushelyek:
T.f.h $p_n(x)$ az $[a,b]$ intervallumon a $k< n$ számú $x_1, x_2,... x_k$ helyeken vált előjelet. Ekkor $p_n(x)(x-x_1)(x-x_2)\cdots (x-x_k)$ az $[a,b]$ intervallumon határozott előjelű, így $$\int_a^b W(x) p_n(x)(x-x_1)(x-x_2)\cdots (x-x_k)dx \ne 0\;.$$ Másfelől a $k$-adfokú $(x-x_1)(x-x_2)\cdots (x-x_k)$ felírható a $p_0(x), p_1(x),...p_k(x)$ polinomok lineárkombinációjaként, amelyek mindegyike ortogonális $p_n(x)$-re, ami ellentmondásra vezet. Következésképpen $p_n(x)$-nek az $[a,b]$ intervallumon pontosan $n$ db különböző zérushelye van.

Oszcillációs tétel: Christoffel-Darboux-azonosság a $\lim_{y\rightarrow x}$ határesetben $$\frac{p_{n+1}'(x)p_{n}(x)-p_{n}'(x)p_{n+1}(x)}{< p_n |p_n>}=\sum_{j=0}^n\frac{p_{j}(x)^2}{< p_j |p_j>} > 0$$ Legyen $x_k$ és $x_{k+1}$ a $p_{n+1}(x)$ polinom két egymást követő valós zérushelye. Ekkor $$p_{n+1}'(x_k)p_{n}(x_k)> 0$$ és $$p_{n+1}'(x_{k+1})p_{n}(x_{k+1})>0$$ következik. Mivel $p_{n+1}'(x_k)$ és $p_{n+1}'(x_{k+1})$ ellentétes előjelűek, az $(x_k,x_{k+1})$ intervallumban $p_{n}(x)$-nek van (pontosan egy) zérushelye.

Diszkrét ortogonalitási reláció: Legyenek $x_1,x_2,...x_{n+1}$ a $p_{n+1}(x)$ polinom zérushelyei. Ekkor $$\sum_{j=0}^n\frac{p_{j}(x_k)p_{j}(x_m)}{< p_j |p_j>}=0\;,$$ ha $k\ne m$, és $$\sum_{j=0}^n\frac{p_{j}^2(x_k)}{< p_j |p_j>}=\frac{p_{n+1}'(x_k)p_{n}(x_k)}{< p_n |p_n>}$$ Ez azt jelenti, hogy a $$p_{j}(x_k)\sqrt{\frac{< p_n |p_n>}{< p_j |p_j>p_{n+1}'(x_k)p_{n}(x_k)}}$$ mátrix ortogonális. Emiatt $$\sum_{k=1}^{n+1}\underbrace{\frac{< p_n|p_n>}{p_{n+1}'(x_k)p_{n}(x_k)}}_{w_k}p_{j}(x_k)p_{l}(x_k)=< p_j |p_j>\delta_{jl}$$ Függvényinterpoláció ortogonális polinomokkal: $$f(x)\approx \sum_{j=0}^n c_j p_j(x)$$ A diszkrét ortogonalitási relációt felhasználva $$c_j=\frac{1}{< p_j|p_j>}\sum_{k=1}^{n+1} w_k p_j(x_k) f(x_k)$$ írható. A Christoffel-Darboux-formula alapján könnyű belátni, hogy ekkor $$f(x_k)=\sum_{j=0}^n c_j p_j(x_k)$$ teljesül, tehát valóban interpolációról van szó.

Gauss-integrálás:
Kiszámítandó $$\int_a^b W(x)f(x)dx$$ Ebből a célból $f(x)$-et a $W(x)$ súlyfüggvényhez és az $(a, b)$ intervallumhoz tartozó $p_j(x)$ ortogonális polinomokkal közelítjük: $$\int_a^b W(x)f(x)dx=\sum_{j=0}^nc_j\int_a^b W(x)p_j(x)p_0(x)dx=c_0< p_0|p_0 >=\sum_{k=1}^{n+1} w_k f(x_k)$$ (Vigyázat! $w_k=\frac{< p_n|p_n>}{p_{n+1}'(x_k)p_{n}(x_k)}\ne W(x_k)\;$)

Összeg kiértékelése: Clenshaw-formula:
Keressük az $$s_N=\sum_{j=0}^{N-1}c_jf_j(x)$$ összeget (adott $x$ pontban), ahol a $c_j$ együtthatók adottak, az $f_j(x)$ függvények pedig az $$f_{j+1}(x)=a_j(x)f_j(x)+b_j(x)f_{j-1}(x)$$ rekurziós összefüggésnek tesznek eleget.
Bevezetjük az $y_k$ mennyiségeket ($x$ függvényei) úgy, hogy $$y_{N+1}=y_N=0$$ és $$y_k=c_k+a_k(x) y_{k+1}+b_{k+1}(x)y_{k+2}\;.$$ Ekkor $$ \begin{eqnarray} s_N & = & [y_{N-1}\underbrace{- a_{N-1}(x) y_{N}-b_{N}(x)y_{N+1}}_{=0}]f_{N-1}(x)\\ &+& \dots \\ &+& \left[{\color{red} {y_{j+1}}}- \right.\left. a_{j+1}(x) y_{j+2}-b_{j+2}(x)y_{j+3}\right]{\color{red} {f_{j+1}(x)}}\\ &+&\left[y_{j}{\color{red} {- a_{j}(x) y_{j+1}}}-b_{j+1}(x)y_{j+2}\right]{\color{red} {f_{j}(x)}}\\ &+&\left[y_{j-1}- a_{j-1}(x) y_{j}{\color{red} {-b_{j}(x)y_{j+1}}}\right]{\color{red} {f_{j-1}(x)}}\\ &+& \dots \\ &+& \left[{\color{red} {y_{2}}}- \right.\left. a_{2}(x) y_{3}-b_{3}(x)y_{4}\right]{\color{red} {f_{2}(x)}}\\ &+&\left[{\color{blue} {y_{1}}}{\color{red} {- a_{1}(x) y_{2}}}-b_{2}(x)y_{3}\right]{\color{red} {f_{1}(x)}}\\ &+&\left[{\color{blue} {y_{0}- a_{0}(x) y_{1}}} {\color{red} { -b_{1}(x)y_{2} }}\right]{\color{red} {f_{0}(x)}}\\ &=& \left[y_{0}- a_{0}(x) y_{1}\right]f_{0}(x)+y_1f_1(x)=\left[c_{0}+ b_{1}(x) y_{2}\right]f_{0}(x)+y_1f_1(x) \end{eqnarray} $$
Közelítés Csebisev-polinomokkal
Minimax-polinom: adott intervallumon adott függvényt minimális maximális hibával közelítő polinom. A Csebisev-polinomok a $[-1,1]$ intervallumon "majdnem" minimax polinomok, ugyanis ezen az intervallumon $-1$ és $1$ között oszcillálnak, a közelítés hibája így közel egyenletes.
Csebisev-polinomok: $$T_n(x)=\cos\left(n\arccos(x)\right)$$ Az első néhány Csebisev-polinom: $$T_0(x)=1$$ $$T_1(x)=x$$ $$T_2(x)=2x^2-1$$ $$T_3(x)=4x^3-3x$$ Általában teljesül a $$T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)$$ rekurziós összefüggés, ugyanis $y=\arccos(x)$ jelöléssel $$\cos((n+1)y)=\cos(y)\cos(ny)-\sin(y)\sin(ny)$$ $$\cos((n-1)y)=\cos(y)\cos(ny)+\sin(y)\sin(ny)$$ A két egyenlet összegéből pedig $$\cos((n+1)y)=2\cos(y)\cos(ny)-\cos((n-1)y)$$ Ortogonalitási reláció: $$\int_{-1}^1\frac{T_j(x)T_k(x)}{\sqrt{1-x^2}}dx=\left\{\begin{array}{lcl} 0\;,& {\rm ha} & j\ne k\\ \pi/2\;, & {\rm ha} & j=k\ne 0 \\ \pi\;, & {\rm ha} & j=k= 0 \\ \end{array} \right.$$ (Az $y=\arccos(x)$ változóra áttérve közvetlen számolással bizonyítható.)
Az $N$-edik Csebisev-polinom gyökei: $$x_k=\cos\left(\pi\frac{k-\frac{1}{2}}{N}\right)\;,\quad k=1,2,\dots , N$$
Diszkrét ortogonalitási reláció: $$\sum_{k=1}^N T_m(x_k)T_j(x_k)=\left\{\begin{array}{lcl} 0\;,& {\rm ha} & j\ne m\\ N/2\;, & {\rm ha} & j=m\ne 0 \\ N\;, & {\rm ha} & j=m= 0 \\ \end{array} \right.$$ ugyanis $$\begin{eqnarray} &&\sum_{k=1}^N T_m(x_k)T_j(x_k)=\sum_{k=1}^N \cos\left(\pi \frac{m\left(k-\frac{1}{2}\right)}{N}\right)\cos\left(\pi \frac{j\left(k-\frac{1}{2}\right)}{N}\right)\\ &&=\frac{1}{2}\left[\sum_{k=1}^N \cos\left(\pi \frac{(m+j)\left(k-\frac{1}{2}\right)}{N}\right)+\cos\left(\pi \frac{(m-j)\left(k-\frac{1}{2}\right)}{N}\right)\right]\\ &&=\frac{1}{2}\Re\left[\sum_{k=1}^N \exp\left(i\pi \frac{(m+j)\left(k-\frac{1}{2}\right)}{N}\right)+\exp\left(i\pi \frac{(m-j)\left(k-\frac{1}{2}\right)}{N}\right)\right] \end{eqnarray} $$ Legyen $$\alpha=\pi \frac{m\pm j}{N}\;,$$ ezzel, ha $\alpha\ne 0$, $$\begin{eqnarray} &&\sum_{k=1}^N \exp\left(i\alpha\left(k-\frac{1}{2}\right)\right)=\exp\left(i\frac{\alpha}{2}\right)\sum_{p=0}^{N-1} \exp\left(i\alpha p\right)\\ &&=\exp\left(i\frac{\alpha}{2}\right)\frac{1-\exp\left(i\alpha N\right)}{1-\exp\left(i\alpha \right)}=\frac{1-\left(-1\right)^{j\pm m}}{\exp\left(-i\frac{\alpha}{2}\right)-\exp\left(i\frac{\alpha}{2}\right)}\\ &&=i\frac{1-\left(-1\right)^{j\pm m}}{2\sin\left(\frac{\alpha}{2}\right)} \end{eqnarray} $$ Ez tisztán képzetes, tehát a valós rész képzésekor $0$-t kapunk. Ha $\alpha=0$, akkor az összeg $N$.
Csebisev-approximáció: $$f(x)\approx \sum_{j=0}^{N-1}c_jT_j(x)\;,$$ ahol $$c_j=\left\{\begin{array}{lcl} \frac{2}{N}\sum_{k=1}^N f(x_k)T_j(x_k)\;,& {\rm ha} & j\ne 0\\ \frac{1}{N}\sum_{k=1}^N f(x_k)\;, & {\rm ha} & j= 0 \\ \end{array} \right.$$ Deriválás és integrálás: $$c'_{j-1}=c'_{j+1}+2jc_j$$ $$C_j=\frac{c_{j-1}-c_{j+1}}{2j}$$

bene@arpad.elte.hu