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
7. hét csütörtöki előadás
Speciális függvények értékeinek kiszámítása
Differenciálegyenletek, hatványsorok, aszimptotikus sorok
Rekurziós formulák
Pl.:
Legendre-polinomok:
$$(n+1)P_{n+1}(x)=(2n+1)xP_n(x)-nP_{n-1}(x)$$
Bessel-függvények:
$$J_{n+1}(x)=\frac{2n}{x}J_n(x)-J_{n-1}(x)$$
Exponenciális integrálok:
$$nE_{n+1}(x)={\rm e}^{-x}-xE_n(x)$$
Koszinusz- és szinuszfüggvények:
$$\cos nx=2\cos x \cos (n-1)x-\cos(n-2)x$$
$$\sin nx=2\cos x \sin (n-1)x-\sin(n-2)x$$
Általában egy
$$y_{n+1}+a_ny_{n}+b_ny_{n-1}=0$$
típusú rekurziós összefüggésnek két lineárisan független megoldása van,
jelöljük ezeket $f_n$-nel
és $g_n$-nel. Ha
$$\lim_{n\rightarrow \infty}\frac{f_n}{g_n}=0\;,$$
akkor $f_n$ a minimális, $g_n$ pedig a domináns megoldás. Növekvő $n$-ek
irányában alkalmazva a rekurziót, a domináns megoldást kapjuk, csökkenő $n$-ek
esetén pedig a minimális megoldást.
Lánctörtek
$$b_0+\frac{a_1}{b_1+\frac{a_2}{b_2+\frac{a_3}{b_3+\dots}}}$$
Kényelmesebb jelölés:
$$b_0+\frac{a_1}{b_1+}\;\frac{a_2}{b_2+}\;\frac{a_3}{b_3+}\dots$$
Pl.:
$$\frac{\sqrt{5}-1}{2}=1+\frac{1}{1+}\;\frac{1}{1+}\;\frac{1}{1+}\dots$$
$$\tan x =\frac{x}{1-}\;\frac{x^2}{3-}\;\frac{x^2}{5-}\;\frac{x^2}{7-}\dots$$
Lánctörtek kiszámítása:
$$t_n\equiv b_0+\frac{a_1}{b_1+}\;\frac{a_2}{b_2+}\dots\frac{a_n}{b_n}
=\frac{A_n}{B_n}$$
Ekkor az $A_n$, $B_n$ mennyiségek a következő rekurzióval kaphatók meg
(J.Wallis, 1655):
$$A_{-1}=1\;,\quad B_{-1}=0 \;,\quad A_{0}=b_0\;,\quad B_{0}=1$$
$$A_j=b_jA_{j-1}+a_jA_{j-2}\;,\quad B_j=b_jB_{j-1}+a_jB_{j-2}\quad (j\ge 1)$$
Ui.
$$\left(\begin{array}{c}A_{j-2}\\B_{j-2}\end{array}\right)=\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right)\left(\begin{array}{c}a_{j-2}\\b_{j-2}\end{array}\right)$$
$$\left(\begin{array}{c}A_{j-1}\\B_{j-1}\end{array}\right)=\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right)\left(\begin{array}{c}a_{j-2}b_{j-1}\\b_{j-2}b_{j-1}+a_{j-1}\end{array}\right)=\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right)\left(\begin{array}{cc}0&
a_{j-2}\\1&
b_{j-2}\end{array}\right)\left(\begin{array}{c}a_{j-1}\\b_{j-1}\end{array}\right)$$
$$\left(\begin{array}{c}A_{j}\\B_{j}\end{array}\right)=\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right)\left(\begin{array}{cc}0&
a_{j-2}\\1& b_{j-2}\end{array}\right)\left(\begin{array}{cc}0&
a_{j-1}\\1& b_{j-1}\end{array}\right)\left(\begin{array}{c}a_{j}\\b_{j}\end{array}\right)=\left(\begin{array}{cc}A_{j-2}&
A_{j-1}\\ B_{j-2}& B_{j-1}\end{array}\right)\left(\begin{array}{c}a_{j}\\b_{j}\end{array}\right)$$
Lánctörtek numerikus meghatározása (Lentz módszere):
Legyen
$$C_j=\frac{A_j}{A_{j-1}}\;,\quad D_j=\frac{B_{j-1}}{B_j}\;.$$
Ezzel
$$t_j=C_jD_jt_{j-1}\;.$$
A $C_j$, $D_j$ mennyiségekre
$$C_j=\frac{b_jA_{j-1}+a_jA_{j-2}}{A_{j-1}}=b_j+\frac{a_j}{C_{j-1}}$$
és
$$D_j=\frac{B_{j-1}}{b_jB_{j-1}+a_jB_{j-2}}=\frac{1}{b_j+a_jD_{j-1}}$$
teljesül
$$C_0=b_0\;,\quad D_0=0\;,\quad t_0=b_0$$
kezdeti feltételekkel. Ha valamelyik nevező nullává válna, akkor egy kis
számmal (pl. $10^{-30}$-cal) helyettesítjük, a beavatkozás a nullával osztást
kivédi, egy újabb ciklus múlva pedig a hatása elenyészik.
Lánctört kifejtés a rekurziós formula alapján:
$$y_{n+1}+a_ny_{n}+b_ny_{n-1}=0$$
$$\frac{y_{n+1}}{y_n}+a_n+b_n\frac{y_{n-1}}{y_n}=0$$
$$\frac{y_n}{y_{n-1}}=-\frac{b_n}{a_n+\frac{y_{n+1}}{y_n}}=-\frac{b_n}{a_n-}\;\frac{b_{n+1}}{a_{n+1}-}\dots$$
Pincherle tétele: a fenti lánctört akkor és csak akkor konvergál, ha létezik a
rekurziónak az $f_n$ minimális megoldása. Ebben az esetben a lánctört értéke
$f_n/f_{n-1}$.
Gamma-függvény, béta-függvény, faktoriális, binomiális
együtthatók
$$\Gamma(z)=\int_0^\infty {\rm e}^{-t}\;t^{z-1}dt$$
Néhány fontos összefüggés:
$$\Gamma(n+1)=n!$$
$$\Gamma(z+1)=z\Gamma(z)$$
Bizonyítás: parciális integrálással.
Reflexiós összefüggés:
$$\Gamma(1-z)=\frac{\pi}{\Gamma(z)\sin(\pi z)}$$
Közelítő formula a numerikus kiértékeléshez (az $Re(z)>0$ komplex félsíkon):
$\begin{eqnarray}\Gamma(z)=\sqrt{2\pi}\left(z+\gamma+\frac{1}{2}\right)^{z+\frac{1}{2}}{\rm
e}^{-\left(z+\gamma+\frac{1}{2}\right)}\left[c_0+\frac{c_1}{z+1}+\frac{c_2}{z+2}+\cdots +\frac{c_N}{z+N}+\epsilon\right]
\end{eqnarray}$
Itt $\gamma$ egész szám. $6$ tag már elegendő $|\epsilon|<2\times 10^{-10}$
pontossághoz
az $Re(z)>0$
komplex
félsíkon
mindenütt. Gyakorlatban
a függvény
gyors
növekedése
miatt
célszerűbb a
logaritmusát
kiszámítani.
Euler-féle béta-függvény:
$$B(z,w)=B(w,z)=\int_0^1t^{z-1}(1-t)^{w-1}dt$$
$$B(z,w)=\frac{\Gamma(z)\Gamma(w)}{\Gamma(z+w)}$$
Bizonyítás:
$$\Gamma(z)\Gamma(w)=\int_0^\infty\int_0^\infty {\rm
e}^{-(t+u)}\;t^{z-1}u^{w-1}dtdu$$
Változócsere:
$$t=xy\;,\quad u=x(1-y)$$
Belátható, hogy $x\in (0,\infty)$ és $y\in (0,1)$, és
$$dt du = x dx dy$$
Ezzel
$$\Gamma(z)\Gamma(w)=\int_0^\infty dx\int_0^1 dy\; x\; {\rm
e}^{-x} x^{z-1}y^{z-1}x^{w-1}(1-y)^{w-1}=\underbrace{\int_0^\infty dx\;{\rm
e}^{-x} \;x^{(z+w)-1}}_{\Gamma(z+w)}\underbrace{\int_0^1 dy\;y^{z-1}(1-y)^{w-1}}_{B(z,w)}$$QED.
Nem teljes gamma-függvény és alkalmazásai: hibafüggvény,
$\chi^2$-eloszlás, Poisson-eloszlás, exponenciális
integrálok
$$\gamma(a,x)=\int_0^x {\rm e}^{-t}\;t^{a-1}dt$$
$$\Gamma(a,x)=\int_x^\infty {\rm e}^{-t}\;t^{a-1}dt$$
$$P(a,x)=\frac{\gamma(x,a)}{\Gamma(a)}$$
$$Q(a,x)=\frac{\Gamma(x,a)}{\Gamma(a)}=1-P(a,x)$$
Sorfejtés ($x< a+1$ esetén gyorsan konvergál):
$$\gamma(a,x)={\rm e}^{-x}x^a\sum_{n_0}^\infty
\frac{\Gamma(a)}{\Gamma(a+1+n)}x^n$$
Lánctört-kifejtés ($x> a+1$ esetén gyorsan konvergál):
$$\Gamma(a,x)={\rm
e}^{-x}x^a\left(\frac{1}{x+}\;\frac{1-a}{1+}\;\frac{1}{x+}\;\frac{2-a}{1+}\;\cdots\right)$$
A gyorsabb konvergencia érdekében célszerű a lánctört páros részét (melynek
közelítő törtjei $t_{2n}$-ek, ha az eredetié $t_n$-ek) venni.
Hibafüggvény:
$${\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x{\rm e}^{-t^2}dt$$
$${\rm erfc}(x)=1-{\rm erf}(x)=\frac{2}{\sqrt{\pi}}\int_x^\infty{\rm
e}^{-t^2}dt$$
Kapcsolat a nem teljes gamma-függvénnyel:
$${\rm erf}(x)=P\left(\frac{1}{2},x^2\right)$$
$${\rm erfc}(x)=Q\left(\frac{1}{2},x^2\right)$$
Poisson-eloszlás és a nem teljes gamma-függvény:
$$\sum_{n=0}^{k-1}{\rm e}^{-\lambda}\; \frac{\lambda^n}{n!}=Q(k,\lambda)$$
$\chi^2$-eloszlás és a nem teljes gamma-függvény:
$$P(\chi^2|\nu)=P\left(\frac{\nu}{2},\frac{\chi^2}{2}\right)$$
Exponenciális integrálok:
$$E_n(x)=\int_1^\infty \frac{{\rm e}^{-xt}}{t^n}dt=x^{n-1}\Gamma(1-n,x)$$
Nem teljes béta-függvény és alkalmazásai: Student-eloszlás,
F-eloszlás, binomiális eloszlás
$$I_x(a,b)=\frac{B_x(a,b)}{B(a,b)}=\frac{\int_0^xt^{a-1}(1-t)^{b-1}dt}{\int_0^1t^{a-1}(1-t)^{b-1}dt}$$
Nyilván
$$I_x(a,b)=1-I_{1-x}(b,a)\;.$$
Lánctört-kifejtés:
$\begin{eqnarray}I_x(a,b)&=&\frac{x^a(1-x)^b}{aB(a,b)}\left[\frac{1}{1+}\;\frac{d_1}{1+}\;\frac{d_2}{1+}\;\cdots\right]\\
d_{2m+1}&=&-\frac{(a+m)(a+b+m)x}{(a+2m)(a+2m+1)}\\
d_{2m}&=&\frac{m(b-m)x}{(a+2m)(a+2m-1)}\end{eqnarray}$
Gyorsan konvergál, ha $x<(a+1)/(a+b+2)$. Ha ez nem teljesül, előbb a fenti reflexiós
összefüggést alkalmazzuk.
Student-eloszlás:
$$A(t|\nu)=\frac{1}{\nu^{1/2}B\left(\frac{1}{2},\frac{\nu}{2}\right)}\int_{-t}^t\left(1+\frac{x^2}{\nu}\right)^{-\frac{\nu+1}{2}}dx=1-I_{\frac{\nu}{\nu+t^2}}\left(\frac{\nu}{2},\frac{1}{2}\right)$$
F-eloszlás:
$$Q(F|\nu_1,\nu_2)=I_{\frac{\nu_2}{\nu_2+\nu_1
F}}\left(\frac{\nu_2}{2},\frac{\nu_1}{2}\right)$$
Binomiális eloszlás:
$$\sum_{j=k}^n\left(\begin{array}{c}n\\j\end{array}\right)p^j(1-p)^{n-j}=I_p(k,n-k+1)$$
Bessel-függvények
$$z^2\frac{d^2w}{dz^2}+z\frac{dw}{dz}+(z^2-\nu^2)w=0$$
$$J_\nu(z)=\left(\frac{1}{2}z\right)^\nu\sum_{k=0}^\infty\frac{\left(-\frac{1}{4}z^2\right)^k}{k!\Gamma(\nu+k+1)}$$
$$Y_\nu(z)=\frac{J_\nu(z)\cos(\nu \pi)-J_{-\nu}(z)}{\sin(\nu \pi)}$$
$$H_\nu^{(1)}(z)=J_\nu(z)+iY_\nu(z)$$
$$H_\nu^{(2)}(z)=J_\nu(z)-iY_\nu(z)$$
$$J_\nu(z)=\frac{\left(\frac{1}{2}z\right)^\nu}{\pi^{\frac{1}{2}}\Gamma\left(\nu+\frac{1}{2}\right)}\int_0^\pi\cos\left(z\cos\theta\right)\sin^{2\nu}\theta
d\theta=\frac{2\left(\frac{1}{2}z\right)^\nu}{\pi^{\frac{1}{2}}\Gamma\left(\nu+\frac{1}{2}\right)}\int_0^1(1-t^2)^{\nu-\frac{1}{2}}\cos(zt)dt$$
$$J_n(z)=\frac{1}{\pi}\int_0^\pi
\cos\left(z\sin\theta-n\theta\right)d\theta=\frac{i^{-n}}{\pi}\int_0^\pi {\rm
e}^{iz\cos\theta}\cos\left(n\theta\right)d\theta$$
$$w_{\nu-1}(z)+w_{\nu+1}(z)=\frac{2\nu}{z}w_\nu(z)$$
$$w_{\nu-1}(z)-w_{\nu+1}(z)=2w'_\nu(z)$$
Itt $w_\nu(z)$ lehet $J_\nu(z)$, $Y_\nu(z)$, $H_\nu^{(1,2)}(z)$.
$$\frac{J_\nu(z)}{J_{\nu-1}(z)}=\frac{1}{2\nu z^{-1}-}\frac{1}{2(\nu+1) z^{-1}-}\frac{1}{2(\nu+2) z^{-1}-}\cdots$$
Aszimptotikus sor $|z|\rightarrow \infty$ esetén:
$$J_\nu(z)=\sqrt{\frac{2}{\pi
z}}\left[P(\nu,z)\cos\left(z-\frac{1}{2}\nu\pi-\frac{1}{4}\pi\right)-Q(\nu,z)\sin\left(z-\frac{1}{2}\nu\pi-\frac{1}{4}\pi\right)\right]$$
Itt
$$P(\nu,z)=1-\frac{(4\nu^2-1)(4\nu^2-9)}{2!(8z)^2}+\frac{(4\nu^2-1)(4\nu^2-9)(4\nu^2-25)(4\nu^2-49)}{4!(8z)^4}-\cdots$$
$$Q(\nu,z)=\frac{4\nu^2-1}{8z}-\frac{(4\nu^2-1)(4\nu^2-9)(4\nu^2-25)}{3!(8z)^3}+\cdots$$
Gömbharmonikusok
$$Y_l^m(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta)e^{im\phi}$$
Itt $P_l^m(x)$-ek az asszociált Legendre-polinomok, a $P_l(x)$ közönséges
Legendre-polinomok deriváltjai:
$$P_l^m(x)=(-1)^m(1-x^2)^{m/2}\frac{d^m}{dx^m}P_l(x)$$
Stabil rekurzió $l$-re:
$$(l-m)P_l^m(x)=x(2l-1)P_{l-1}^m(x)-(l+m-1)P_{l-2}^m(x)$$
A rekurzió kezdeti feltétele:
$$P_m^m(x)=(-1)^m(2m-1)!!(1-x^2)^{m/2}\;,\quad P_{m-1}^m(x)=0\;.$$
Elliptikus integrálok és elliptikus
függvények
Elsőfajú elliptikus integrál:
$$F(\phi,k)=\int_0^\phi
\frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}=\int_0^{\sin\phi}\frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}}$$
vagy más definíció szerint
$$F(x,k)= \int_0^{x}\frac{dt}{\sqrt{(1-t^2)(1-k^2t^2)}}$$
Elsőfajú teljes elliptikus integrál:
$$K(k)=F(\pi/2,k)=\int_0^{\pi/2} \frac{d\theta}{\sqrt{1-k^2\sin^2\theta}}$$
Másodfajú elliptikus integrál:
$$E(x,k)=\int_0^{x}dt\sqrt{\frac{1-k^2t^2}{1-t^2}}$$
Harmadfajú elliptikus integrál:
$$\Pi(x,k,\nu)=\int_0^{x}\frac{dt}{(1-\nu t^2)\sqrt{(1-t^2)(1-k^2t^2)}}$$
Numerikus számításhoz használható alternatíva: Carlson-függvények:
$$R_F(x,y,z)= \frac{1}{2}\int_0^{\infty}\frac{dt}{\sqrt{(t+x)(t+y)(t+z)}}$$
Ezzel az elsőfajú elliptikus integrál kifejeheztő.
$$R_J(x,y,z,p)=
\frac{3}{2}\int_0^{\infty}\frac{dt}{(t+p)\sqrt{(t+x)(t+y)(t+z)}}$$
Ezzel a másod- és harmadfajú elliptikus integrálok kifejezhetők.
Degenerált esetek:
$$R_C(x,y)=R_F(x,y,y)\;,\quad R_D(x,y,z)=R_J(x,y,z,z)$$
Az elliptikus integrálok kiszámítása az ún. duplikációs tételen alapul:
$$R_F(x,y,z)=2R_F(x+\lambda,y+\lambda,z+\lambda)=R_F\left(\frac{x+\lambda}{4},
\frac{y+\lambda}{4},\frac{z+\lambda}{4} \right)$$
ahol
$$\lambda=\sqrt{yz}+\sqrt{zx}+\sqrt{xy}\;.$$
A tételt sokszor egymás után alkalmazva az argumentumok közel egyenlőkké
válnak.
Végül pedig
$$R_F(x,x,x)=x^{-1/2}\;.$$
Alkalmazások:
- Inga lengésideje:
$$T=2\sqrt{\frac{l}{g}}K\left(\sqrt{\frac{E+mgl}{2mgl}}\right)$$
- Ellipszis kerülete:
$$2\pi aE(\epsilon)=2\pi aE(1,\epsilon)$$
Itt $a$ az ellipszis fél nagytengelye, $\epsilon$ az excentricitása,
$$E(k)=E(1,k)$$
pedig a teljes másodfajú elliptikus integrál.
- Ellipszoid felszíne:
$$2\pi c^2+\frac{2\pi b}{a^2-c^2}\left(c^2E(x,k)+(a^2-c^2)F(x,k)\right)\;,$$
ahol
$$x=\frac{\sqrt{a^2-c^2}}{a}\;,\quad
k=\frac{a}{b}\frac{\sqrt{b^2-c^2}}{\sqrt{a^2-c^2}}\;,$$
$a>b>c$ pedig az ellipszoid tengelyei.
Hipergeometrikus függvény
$$\phantom{a}_{2}F_1(a,b,c,z)=1+\frac{ab}{c}z+\frac{a(a+1)b(b+1)}{c(c+1)}\frac{z^2}{2!}+\cdots
+\frac{a(a+1)\cdots (a+j-1)b(b+1)\cdots (b+j-1)}{c(c+1)\cdots
(c+j-1)}\frac{z^j}{j!}+\cdots$$
Differenciálegyenlete:
$$z(z-1)F''=abF-\left[c-(a+b+1)z\right]F'$$
Kiszámítás: a differenciálegyenlet közvetlen megoldásával.
bene@arpad.elte.hu