Számítógépes Fizika (Computational Physics) :Numerikus módszerek összefoglaló.

Bartha Ferenc

e-mail: barthaf@physx.u-szeged.hu, http : www.jate.u-szeged.hu/˜ barthaf/

Szegedi Tudományegyetem, Elméleti Fizikai Tanszék (készültségi fok: April 24, 2002)

Folyamatosan bővítve, javítva, ahogy az anyagban haladunk. Kérem hogy az elgépeléseket, helyesírási hibákat, hibás és téves képleteket, kijelentéseket keressétek és segítsetek kijavítani azokat!

Contents

I  Közönséges differenciálegyenletek (ODE: Ordinary Differential Equations)  2

A Numerikus differenciálás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 B Numerikus integrálás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 C ODE redukálása egyenletrendszerre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 D Euler módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 E Pontoság . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 F Stabilitás . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 G Többlépéses módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 H Implicit módszerek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 I Runge-Kutta módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 J Kiegészítések . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

1 Gauss-kvadratúra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Adaptív lépések a differenciálegyenletek integrálásakor . . . . . . . . . . . . . . . . . . . . . . . 9

II  Parciális differenciálegyenletek  11

A A lineáris egyenletek osztályozása . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 B Hiperbolikus egyenletek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 C Parabolikus egyenletek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 D Kiegészítés: Tridiagonális egyenletrendszer megoldása rekurzióval . . . . . . . . . . . . . . . . . . . 17 E Kiegészítés: a diffúziós probléma több dimenzióban . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 F Kiegészítés: Av=b lineáris egyenletrendszer megoldása . . . . . . . . . . . . . . . . . . . . . . . . . 19 G Elliptikus egyenletek: peremérték probléma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

1 Relaxáció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2 Váltakozó irányú implicit módszer: ADI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Ciklikus redukció . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Fourier módszer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

III Fourier transzformáció 24 A A Fourier transzformációról általában . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 B Mintavételezés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 C Diszkrét Fourier transzformáció (DFT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 D Gyors Fourier Transzformáció (FFT) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 E FFT valós függvényekre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 F sinFFT szinusz és cosFFT koszinusz Fourier transzformáltak . . . . . . . . . . . . . . . . . . . . . . 28 G Többdimenziós FFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 H Konvolúciós egyenletek . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

I. KÖZÖNSÉGES DIFFERENCIÁLEGYENLETEK (ODE: ORDINARY DIFFERENTIAL EQUATIONS)

A. Numerikus differenciálás

Keressük f(x) deriváltját az x = 0 pontban, f′(0)-t miközben ismerjük f(x) értékét az x = 0 pont körüli ekvidisztáns beosztáson

f_{k} = f(x_{k}) ; x_{k} = kh ; k = 0, ±1, ±2, ...  (1.1)

A függvényt x = 0 körül Taylor sorral közelítjük:

f(x) = f(0) + \frac{x}{1!}f′(0) + \frac{x^{2}}{2!} f′′(0) + \frac{x^{3}}{3!} f′′′(0) + .... ⟹ f_{k} = ∑_{i=0}^{n} \frac{(kh)^{i}}{i!} f^{(i)}(0) + O (h^{n+1})  (1.2)

Ebből különböző kifejezéseket származtathatunk f′(0)-ra. Így kétpontos formula

f_{0} = f(0) ; f_{±1} = f(0) ± hf′(0) + O (h^{2}) ⟹ f′(0) = ±\frac{f_{±1} − f_{0}}{h} + O (h)  (1.3)

illetve hárompontos

f_{±1} = f(0) ± hf′(0) + \frac{h^{2}}{2} f′′(0) + O (h^{3}) ⟹ f′(0) = \frac{f_{+1} − f_{−1}}{2h} + O (h^{2})  (1.4)

vagy 5-pontos

f′(0) = \frac{1}{12h} [f_{−2} − 8f_{−1} + 8f_{+1} + f_{+2}] + O (h^{4})  (1.5)

Magasabb rendű deriváltakhoz jutunk például

f_{±1} = f(0) ± hf′(0) + \frac{h^{2}}{2} f′′(0) ± \frac{h^{3}}{6} f′′′(0) + O (h^{4})  (1.6)
f_{+1} + f_{−1} = 2f(0) + \frac{h^{2}}{2} f′′(0) + O (h^{4}) ⟹ f′′(0) = \frac{f_{+1} − 2f_{0} + f_{−1}}{h^{2}} + O (h^{2})  (1.7)

három ponttal. További hasznos formulák

4-pont  5-pont
hf^{(1)} ±\frac{1}{6} [−2f_{∓1} − 3f_{0} + 6f_{±1} − f_{±}]  \frac{1}{12} [f_{−2} − 8f_{−1} + 8f_{+1} + f_{+2}]
h^{2}f^{(2)}  [4f_{−1} − 2f_{0} + f_{+1}] \frac{1}{12} [−f_{−2} + 16f_{−1} − 30f_{0} + 16f_{+1} − f_{+2}]
h^{3}f^{(3)}  ± [−f_{∓1} + 3f_{0} − 3f_{±1} + f_{±}]  \frac{1}{2} [−f_{−2} + 2f_{−1} − 2f_{+1} + f_{+2}]
h^{4}f^{(4)}  ...........  [f_{−2} − 4f_{−1} + 6f_{0} − 4f_{+1} + f_{+2}]

(1.8)

B. Numerikus integrálás

Szükség lehet arra, hogy közelítőleg meghatározzuk az

∫_{a}^{b} f(t)dt , a > b  (1.9)

határozott integrált, amikor f(t) primitív függvénye nem ismert. (Esetleg nem is létezik , mint pl. f(t) = \frac{sin(x)}{x}, vagy f(t) = \sqrt{1 + x^{2}} esetben. Vagy az integrandus nem is adott analitikus formában, csak pl. tabulálva, stb.)

Felosztjuk az integrációs tartományt részintervallumokra, az egyes részintervallumokon valamilyen egyszerű függvénnyel (pl. polinommal) illesztjük (helyettesítjük) a függvényt és integráljuk az illesztő függvényt. Minél kisebb részintervallumokon illesztünk és integrálunk, annál pontosabb értéket várunk. Trapéz szabály: Az [a, b] intervallumot n egyenlő részre vágjuk. Ezek hossza h = \frac{b−a}{n}. A függvényt az egyes részintervallumkonon első fokú polinommal közelítjük

f(t) ≈ f_{0} + t · f′(0) = f_{0} + t · \frac{f_{1} − f_{0}}{h}  (1.10)
∫_{0}^{h} f(t)dt ≈ h · f_{0} + \frac{h^{2}}{2} · \frac{f_{1} − f_{0}}{h} = \frac{h}{2} [f_{0} + f_{1}] + O(h^{3})  (1.11)

illetve

∫_{a}^{b} f(t)dt ≈ \frac{h}{2} [f_{0} + 2f_{1} + ... + 2f_{n−1} + f_{n}]  (1.12)

A közelítés hibája: Ha f^{(2)} (második derivált) folytonos az [a, b]-n és létezik olyan M_{2} felső korlát, hogy |f^{(2)}| ≤ M_{2} az egész intervallumon, akkor

E_{T} ≤ \frac{b − a}{12}h^{2}M_{2}

Azt látjuk, hogy ez esetben valóban nagyon kis h-t választva nagy pontosságot érhetnénk el. A gyakorlatban azonban nem lehet tetszőlegesen kis h-t használni a numerikus bizonytalanságok miatt. Ha a trapéz szabály ”értelmes” h-val nem ad kielégítő pontosságot, akkor más módszert kell használni. Simpson szabály: Az [a, b] intervallumot páros számú, egyenlő hosszúságú részintervallumra osztjuk. Minden részintervallum páron másodfokú polinommal közelítve a függvényt

f(t) ≈ f_{0} + t · f′(0) + \frac{t^{2}}{2} · f′′(0) = f_{0} + t · \frac{f_{1} − f_{−1}}{2h} + \frac{t^{2}}{2} · \frac{f_{+1} − 2f_{0} + f_{−1}}{h^{2}} + O (h^{3})  (1.13)
∫_{−h}^{h} f(t)dt ≈ 2h · f_{0} + \frac{h^{3}}{6} · \frac{f_{+1} − 2f_{0} + f_{−1}}{h^{2}} = \frac{h}{3} [f_{−1} + 4f_{0} + f_{1}] + O(h^{5})  (1.14)

Vegyük észre, hogy a t^{3} tag integrálja nulla lenne, így a vártnál jobb közelítést kaptunk, ezért az O(h^{5}) hiba! A teljes integrál tehát

∫_{a}^{b} f(t)dt ≈ \frac{h}{3} [f_{0} + 4f_{1} + 2f_{2} + ... + 2f_{n−2} + 4f_{n−1} + f_{n}]  (1.15)

A közelítés hibája: Ha f^{(4)} (negyedik derivált) folytonos az [a, b]-n és létezik olyan M_{4} felső korlát, hogy |f^{(4)}| ≤ M_{4} az egész intervallumon, akkor

E_{S} ≤ \frac{b − a}{180}h^{4}M_{4}

Simpson 3/8 szabály: harmadfokú Taylor polinommal és megfelelő numerikus deriváltakkal

∫_{x_{0}}^{x_{3}} f(t)dt ≈ \frac{3h}{8} [f_{0} + 3f_{1} + 3f_{2} + f_{3}] + O(h^{5})  (1.16)

Bode szabály: negyedfokú polinommal

∫_{x_{0}}^{x_{4}} f(t)dt ≈ \frac{2h}{45} [7f_{0} + 32f_{1} + 12f_{2} + 32f_{3} + 7f_{4}] + O(h^{7})  (1.17)

C. ODE redukálása egyenletrendszerre

A közönséges N-edrendű differenciálegyenletek elsőrendű differenciélegyenlet rendszerré alakíthatók.

Példa: A másodrendű

\frac{d^{2}y}{dx^{2}} + q(x)\frac{dy}{dx} = r(x)  (1.18)

egyenlet elsőrendű egyenletek rendszerévé átírható, mint

\frac{dy}{dx} = z(x) és \frac{dz}{dx} = r(x) − q(x)z(x)  (1.19)

ahol z(x) egy új változó.

Szokásos és nyilvánvaló választás olyan új változókat bevezetni, melyek egymás deriváljai, ezzel minden ODE átalakítható a kívánt alakra. Néha hasznos, hogy egyéb faktorokat beépítsünk az új változókba, ezzel nemcsak egyszerűbb egyenletekre juthatunk, de a megoldások várható szinguláris viselkedéséből származó numerikus problémákat is megelőzhetjük. Ha azt látjuk, hogy az eredeti egyenlet megoldása sima, de az önkényesen bevezett segédfüggvények őrült dolgokat csinálnak, akkor célszerű a dolognak utánajárni és más változókat választani.

Az általános feladat tehát N csatolt elsőrendű differenciálegyenlet,

\frac{dy}{dx} = f(x, y)  (1.20)

megoldása. Itt a ’vektor’ jelöléseben y(x) = {y_{1}(x), ..., y_{N}(x)} az ismeretlenek, míg f = {f_{1}, ..., f_{N}} adott függvények f_{i} = f_{i}(x, y_{1}(x), ..., y_{N}(x)).

Egy jól határozott probléma nemcsak a differenciálegyenlet megadását jelenti, hanem ú.n. kezdeti feltételek (vagy határfeltételek) kijelölését is. A határfeltételek algebrai feltételek lehetnek y_{i}(x_{a}) értékeire bizonyos x_{a} pontokban. A legegyszerűbb esetben megkövetelhetjük, hogy egy változó adott pontban adott értékrt vegyen fel, y_{i}(x_{a}) = c_{i,a} , de igen bonyolultak is lehetnek, mint például nemlineáris egyenlet rendszer a változók valamely pontban felvett értéke között.

A továbbiakban egyetlen (skalár)

\frac{dy}{dx} = f(x, y)  (1.21)

differenciálegyenlet megoldásával foglalkozunk. A csatolt egyenletrendszer hasonlón oldható meg természetes mátrixalgebra segítségével.

Az y(x = x_{0}) = y_{0} típusú kezdőfeltételt fogjuk használni. Az egyszerűség kedvéért vizsgáljuk a feladatot az x ∈ [0, 1] intervallumon, a kezdőfeltétel legyen y(0) = y_{0}. Osszuk be az intervallumot N egyenlő részere, legyen h = 1/N , x_{n} = nh és y_{n} = y(x_{n}); n = 0, 1, ..., N.

D. Euler módszer

Az n-edik osztópontban a deriváltat a kétpontos előre formulával helyettesítve

\frac{dy}{dx}|_{x=x_{n}} ≈ \frac{y_{n+1} − y_{n}}{h} = f(x_{n}, y_{n})  (1.22)

kapjuk a rekurzióra alkalmas

y_{n+1} = y_{n} + hf(x_{n}, y_{n}) + O(h^{2})  (1.23)

egyenletet. Az aszimmetrikus formula h intervallumot léptet előre x_{n}-ről x_{n+1} = x_{n} + h-ra. , miközben csak az intervallum kezdetén levő derivált információt használja fel. A lokális hiba csak egy nagyságrenddel kisebb a dy = y_{n+1} − y_{n} ∼ O(h) korrekciónál.

Javíthat a helyzeten, ha a Taylor sorból több tagot veszünk figyelembe, pl.

y_{n+1} = y_{n} + hy′_{n} + \frac{h^{2}}{2} y′′_{n} + O(h^{3})  (1.24)

Behelyettesítve a

y′ = \frac{dy}{dx} = f és y′′ = \frac{df}{dx} = \frac{∂f}{∂x} + \frac{∂f}{∂y} \frac{dy}{dx} = \frac{∂f}{∂x} + f \frac{∂f}{∂y}  (1.25)

deriváltakat

y_{n+1} = y_{n} + hf_{n} + \frac{h^{2}}{2} [\frac{∂f}{∂x} + f \frac{∂f}{∂y} ]_{n} + O(h^{3})  (1.26)

jobb rekurziós egyenletet kaptunk. Persze ez csak akkor hasznos, ha f(x, y) parciális deriváltjai könnyen (analitikusan) számolhatók.

E. Pontoság

A közelítő eljárásaink csonkolási lépéshibái a különböző sorfejtésekből tipikusan O(h^{m+1}) vislkedésűek. A globális hiba, miközben N ∼ O(1/h) lépést tettünk NO(h^{m+1}) ∼ O(h^{m}) lesz. A gyakorlati használat során valójában a pontatlanságnak csak egyik oka származik a csonkolásból. A numerikus hibák másik forrása a számítógépek pontatlan aritmetikája. A véges számábrázolásból eredően az elemi aritmetikai műveletek eredménye egy - a géptől és a programnyelvtől függő - tipikus hibával rendelkezik. Ha egy művelet elvégzésekor az aritmetikai hiba O(ε), akkor az O(1/h) lépés során felhalmozódó összes aritmetikai hiba O(ε/h). A teljes eljárás hibája tehát

ϵ = O(ε/h) + O(h^{m})  (1.27)

ami akkor a legkisebb, ha

h ∼ ε^{\frac{1}{m+1}} ϵ ∼ ε^{\frac{m}{m+1}}  (1.28)

Ennél nagyobb lépésköznél a sorfejtés hibája, kisebbnél a számítógép hibája dominál. Tipikus értékek: 8 byte-on ábrázolt lebegőpontos számoknál ε ∼ 10^{−16} míg 4 byte egyszeres pontosságnál ez csak ε ∼ 10^{−7} Tájékozódásul egy táblázat ε ∼ 10^{−16} gépi adattal

\begin{matrix}m & h & ϵ \\ 1 & 1.0·10^{−8} & 1.0·10^{−8} \\ 2 & 4.6·10^{−6} & 2.2·10^{−11} \\ 3 & 1.0·10^{−4} & 1.0·10^{−12} \\ 4 & 6.3·10^{−4} & 1.6·10^{−13} \\ 5 & 2.2·10^{−3} & 4.6·10^{−14} \end{matrix}  (1.29)

F. Stabilitás

Az Euler módszeren bemutatva.....

A numerikus eljárás során y_{n} értékét több okból is hibával kapjuk meg. Fontos kérdés, hogy valamely korábbi pontban begyüjtött hiba hogyan ’fejlődik’ a további lépések során. Tegyük fel, hogy y_{n} helyett y_{n} + δy_{n} hibával terhelt értékről kívánunk továbblépni az

y_{n+1} + δy_{n+1} = y_{n} + δy_{n} + h [f(x_{n}, y_{n}) + \frac{∂f}{∂y} |_{n} δy_{n}]  (1.30)

Euler módszerrel. A hiba ebből

δy_{n+1} = [1 + h \frac{∂f}{∂y} |_{n}] δy_{n}  (1.31)

ami csökkenő görgetett hibát jelent, ha

|1 + h\frac{∂f}{∂y} |^{2} < 1  (1.32)

egyébként a hiba nő. Az f(x, y) parciális deriváltja (általában komplex) ismeretében eldönthetjük, hogy lehet-e stabil az eljárás, illetve, hogy mekkora lépésközzel dolgozzunk ennek elérésére.

Példa: α > 0

\frac{dy}{dx} = −αy \frac{∂f}{∂y} = −α stabil h < \frac{2}{α}  (1.33)
\frac{dy}{dx} = +αy ; \frac{dy}{dx} = iαy instabil  (1.34)

G. Többlépéses módszerek

Az Euler módszerrel az n-edik pontbeli értékekből léptünk az x_{n+1}-be. Több korábbi pont figyelembel vételével további rekurziós formulát készíthetünk. Integráljuk formálisan a differenciálegyenletet

y_{n+1} = y_{n} + ∫_{x_{n}}^{x_{n+1}} f(x, y)dx  (1.35)

A probléma, hogy nem ismerjük y(x)-et az integrációs tartományon, így f(x, y)−t sem. Ismerjük viszont a korábbi lépéseinkből az y_{n−1} és y_{n} értékét. Ezen két érték segítségével.használjunk lineáris extrapolációt az integrandusra:

f ≈ \frac{x − x_{n−1}}{h} f_{n} − \frac{x − x_{n}}{h} f_{n−1} + O(h^{2})  (1.36)

Beírva és integrálva kapjuk az Adam-Bashforth két lépéses (prediktor) módszert

y_{n+1} = y_{n} + h [\frac{3}{2}f_{n} − \frac{1}{2}f_{n−1}] + O(h^{3})  (1.37)

Több lépéses formulákat kaphatunk, ha magasabbrendű extrapolációt használunk. Ilyen pl. az Adam-Bashforth négy lépéses (prediktor) módszer, amikor harmadfokú polinómmal illesztjük a függvény f_{n}, f_{n−1}, f_{n−2} és f_{n−3} értékét. Az eredmény

y_{n+1} = y_{n} + \frac{h}{24} [55f_{n} − 59f_{n−1} − 37f_{n−2} − 9f_{n−3}] + O(h^{4})  (1.38)

A több lépéses módszerek indításához nem elég egyetlen pont (kezdőfeltétel). Egyéb módszerrel, pl. az Euler módszerrel előbb le kell gyártani a szükséges számú induló pontot.

H. Implicit módszerek

Eleddig y_{n+1} értékét explicite kifejeztük. Az implicit módszerekben y_{n+1} értékére egy egyenletet állítunk fel, ezt az egyenletet még külön meg kell oldani y_{n+1} explicit előállításához.

Tekintsük egy x_{n+ \frac{1}{2}} = (n + \frac{1}{2}) h pontot az osztópontok között félúton, ahol is

\frac{dy}{dx}|_{x_{n+ \frac{1}{2}}} = f(x_{n+ \frac{1}{2}} , y_{n+ \frac{1}{2}} )  (1.39)

Használjunk szimmetrikus formulát az x_{n+ \frac{1}{2}} pontban a deriváltra a jobboldalon, a baloldalon pedig helyettesítsük f_{n+ \frac{1}{2}} értékét az osztópontbeli értékek számtani közepével:

\frac{dy}{dx}|_{x_{n+ \frac{1}{2}}} ≈ \frac{y_{n+1} − y_{n}}{h} + O(h^{2}) = \frac{f_{n} + f_{n+1}}{2} + O(h^{2})  (1.40)

Mindkét oldalon közelítés hibája O(h^{2}). A Heun rekurziós összefüggés tehát

y_{n+1} = y_{n} + \frac{h}{2} [f(x_{n} , y_{n}) + f(x_{n+1} , y_{n+1})] + O(h^{3})  (1.41)

ami az ismeretlen y_{n+1} értéket mindkét oldalon tartalmazza. Konkrét differenciálegyenletekre (adott f mellett) az egyenlet analitikusan, vagy numerikusan (gyökkereső algoritmusokkal) megoldható.

Az Adams-Bashforth-Moulton módszer egyszerre többlépéses és implicit. Másodfokú polinomot illesztünk a függvény f_{n−1}, f_{n} és f_{n+1} értékeire majd kapjuk az

y_{n+1} = y_{n} + \frac{h}{12} [5f_{n+1} + 8f_{n} − f_{n−1}] + O(h^{4})  (1.42)

kétlépéses formulát. Harmadfokú illesztéssel a három lépéses formula

y_{n+1} = y_{n} + \frac{h}{24} [9f_{n+1} + 19f_{n} − 5f_{n−1} + f_{n−2}] + O(h^{5})  (1.43)

Az implicit egyenletek megoldásának egyik módszere a prediktor-korrektor eljárás. Előbb valamely ’prediktor’ eljárással megbecsüljük y_{n+1} értékét, majd ezt az értéket helyettesítjük az egyenlet jobb oldalán f_{n+1} = f(x_{n+1} , y_{n+1})-be, hogy ezzel a bal oldalon ’korrigált’ y_{n+1} -et állítsunk elő.

I. Runge-Kutta módszer

Az implicit módszer mintájára írjuk fel, hogy

\frac{dy}{dx}|_{x_{n+ \frac{1}{2}}} ≈ \frac{y_{n+1} − y_{n}}{h} + O(h^{2}) = f(x_{n+ \frac{1}{2}} , y_{n+ \frac{1}{2}} )  (1.44)

majd a baloldalon helyettesítsük f_{n+ \frac{1}{2}} értékét

y_{n+ \frac{1}{2}} ≈ y_{n} + \frac{h}{2}f(x_{n}, y_{n}) + O(h^{2}) ⟹ f(x_{n+ \frac{1}{2}} , y_{n+ \frac{1}{2}} ) ≈ f(x_{n+ \frac{1}{2}} , y_{n} + \frac{h}{2}f(x_{n}, y_{n})) + O(h^{2})  (1.45)

és így jutottunk el a másodrendű RK eljáráshoz:

k_{1} = hf(x_{n}, y_{n}) , k_{2} = hf(x_{n} + \frac{h}{2} , y_{n} + \frac{k_{1}}{2} )  (1.46)
y_{n+1} = y_{n} + k_{2} + O(h^{3})  (1.47)

Hasonlóan magasabb rendű formulák készíthetők. Példaként a negyedrendű RK megoldás formulái:

k_{1} = hf(x_{n}, y_{n}) , k_{2} = hf(x_{n} + \frac{h}{2} , y_{n} + \frac{k_{1}}{2} )  (1.48)
k_{3} = hf(x_{n} + \frac{h}{2} , y_{n} + \frac{k_{2}}{2} ) , k_{4} = hf(x_{n} + h , y_{n} + k_{3})  (1.49)
y_{n+1} = y_{n} + \frac{k_{1}}{6} + \frac{k_{2}}{3} + \frac{k_{3}}{3} + \frac{k_{4}}{6} + O(h^{5})  (1.50)

J. Kiegészítések

1. Gauss-kvadratúra

A korábbiakban megismert eljárások során a határozott integrálokat úgy számoltuk ki, hogy az integrálási tartományt egyenlő hosszúságú részintervallumokra bontottuk, majd azokon a függvényt valamilyen (alacsony) fokszámú polinommal helyettesítettük. A trapéz szabálytól (elsőfokú) a Bode (negyedfokú) integrálásig megismert formulák magasabb rendű polinomokra elvileg könnyen felírhatók, mégis ritkán használják azokat. Két oka is van a mellőzésnek:

• A magasabbrendű képletekben az együtthatók váltakozó előjelűek, a numerikus hibák nagyon felnőhetnek a sok taggal végzett műveletekben.

• A Gauss-Legendre integrálás a polinommal jól közelíthető függvényekre viszonylag kevés pontban elvégzett függvényszámolással is nagy pontosságú eredményt ad, ami különösen fontos több dimenziós integrálok számolásánál.

A módszer bemutatásához tekintsük az

I = ∫_{−1}^{1} f(x)dx  (1.51)

integrált. Alkalmas koordináta-transzformációval minden integrálunk ilyen alakra hozható. Az elemi kvadratúra formulák szerint az integrál

I ≈ ∑_{i=1}^{N} w_{i}f(x_{i}) , x_{i} = −1 + 2 · \frac{i − 1}{N − 1}  (1.52)

közelítő összeggel számolható, ahol csak a w_{i} súlyok változnak a különböző eljárásokban. Nyilvánvaló, hogy ha az f(x) függvény N − 1-edfokú polinóm, akkor alkalmas súlyokkal egzaktul integrálhatjuk ezzel a formulával, hiszen az

∫_{−1}^{1} x^{p}dx = ∑_{i=1}^{N} w_{i}x_{i}^{p} , p = 0, 1, ..., N − 1  (1.53)

egyenletrendszer w_{i}-re megoldható. Ha N páratlan, akkor ez igaz N-edfokú polinomra is.

Ha elengedjük a megszorítást, hogy x_{i}-k egy ekvidisztáns beosztás pontjai, akkor az N darab súly mellett az N darab x_{i} osztópont is szabadon variálható paraméter lesz. A súlyok és az osztópontok optimális megválasztásával elérhetjük, hogy 2N − 1-nél nem magasabb fokszámú polinomokra

I = ∑_{i=1}^{N} w_{i}f(x_{i})  (1.54)

egzakt legyen. A feladat megoldásához a P_{i}(x) Legendre polinomokat használjuk. P_{i}(x) egy i-edfokú polinom, i darab zéróhelye van a [−1, 1] intervallumon. A Legendre polinomok teljes otrogonális rendszert alkotnak a [−1, 1] intervallumon, azaz

∫_{−1}^{1} P_{i}(x)P_{j}(x)dx = \frac{2}{2i + 1}δ_{ij}  (1.55)

Ha f(x) = p_{2N−1}(x) egy 2N − 1-nél nem magasabb fokszámú polinom (illetve ezzel jól közelíthető), akkor polinomosztással felírhatjuk, hogy

p_{2N−1}(x) = p_{N−1}(x) · P_{N}(x) + q_{N−1}(x)  (1.56)

ahol p_{N−1}(x) és q_{N−1}(x) polinomok fokszáma maximum N − 1. Az integrál ekkor

∫_{−1}^{1} f(x)dx = ∫_{−1}^{1} {p_{N−1}(x) · P_{N}(x) + q_{N−1}(x)} dx = ∫_{−1}^{1} q_{N−1}(x)dx  (1.57)

mert P_{N}(x) ortogonális minden N-nél kisebb fokszámú polinomra, így p_{N−1}(x)-re is. Ezzel a feladatot egy durván feleakkora fokszámú polinom integrálására redukáltuk, ami egy rögzített beosztás mellett a súlyok alkalmas megválasztásával egzaktul elvégezhető. A probléma csak az, hogy el kellene elvileg végeznünk a polinomosztást, hogy q_{N−1}(x) előálljon. Ha azonban előírjuk, hogy az x_{i}-k P_{N}(x) zéróhelyeire essenek, akkor

∫_{−1}^{1} f(x)dx = ∑_{i=1}^{N} w_{i}q_{N−1}(x_{i}) = ∑_{i=1}^{N} w_{i}f(x_{i})  (1.58)

tehát ekkor nem kell ismernünk q_{N−1}(x)-et.

A Legendre polinomok zéróhelyei megadhatók, az ehhez tartozó sulyok pedig megmutathatóan

w_{i} = \frac{2}{(1 − x_{i}^{2}) [P′_{N}(x)]^{2}}  (1.59)

módon választandók. Sima függvények, melyek véges intervallumon hatványsorba fejthetők jól közelíthetők polinomokkal, így azokra ez az integrálási séma nagy pontosságot ad.

Hasonló kvadratúra formulák vezethetők le más típusú integrálokra egyéb ortogonális függvényekkel is, ilyenek például

∫_{0}^{∞} e^{−x}f(x)dx ≈ ∑_{i=1}^{N} w_{i}f(x_{i}) , L_{N}(x_{i}) = 0 ahol L_{N}(x) Laguerre polinom  (1.60)

∫_{−∞}^{∞} e^{−x^{2}}f(x)dx ≈ ∑_{i=1}^{N} w_{i}f(x_{i}) , H_{N}(x_{i}) = 0 ahol H_{N}(x) Hermite polinom  (1.61)

2. Adaptív lépések a differenciálegyenletek integrálásakor

Eddigi módszereinkben a megoldást a teljes megcélzott intervallumon egyetlen előre kiválasztott lépésközzel állítottuk elő. A meghatározandó függvény viselkedése viszont nagyon különbözhet az egyes résztartományokon. Olyan tartományban, ahol a függvény sima a kis lépésköz feleslegesen sok számolási munkát, indokolatlanul halmozódó kerekítési hibát okoz. Ahol a függvény gyorsan változik, finomabb beosztást kellene használni.

Kézenfekvő megoldás minden lépésnél megvizsgálni az f(x_{n}, y_{n}) elsőrendű (és esetleg magasabb) deriváltjainak a nagyságát és ennek megfelelően menet közben változtatni lépésközt. Kvalitatíve jó lehet ez az eljárás, de nehéz általános kvantitatív szabályt a derivált alapján felállítani, arról nem is szólva, hogy ezen deriváltaknak az előállítása (ha egyáltalán lehetséges) nagyon munkaigényes feladat.

Univerzálisabb megoldás, ha a lépésközt valamely tervezett hibához igazítjuk. Tegyük fel, hogy a lépésenkénti hibát szeretnénk valamilyen optimális ϵ_{opt} értéken tartani. Ha sikerülne az egyes lépések során kis többletmunka árán becslést kapni arra, hogy mekkora a lépés ϵ_{h} hibája, akkor h-t menet közben úgy növelhetnénk/csökkenthetnénk, hogy ϵ_{h} ≈ ϵ_{opt} legyen. Példa: Lépésduplázós negyedrendű Runge-Kutta eljárás

Az aktuális RK lépésben x_{n}-ről indulva kétféleképpen is ellépünk x_{n} + 2h-ra:

• a) az aktuális h-val kétszer lépünk

• b) majd 2h-val egyetlen lépéssel.

Az egzakt y(x + 2h) értékről és az iménti lépésekben kapott két RK O(h^{5}) rendben pontos eredményről feltehetjük, hogy

y(x + 2h) = y_{a} + 2 (h)^{5} · Φ + O(h^{6})  (1.62)
y(x + 2h) = y_{b} + (2h)^{5} · Φ + O(h^{6})  (1.63)

ahol Φ ≈ y^{(5)}(x)/5! nagyjából ugyannak vehető a két egyenletben. A két hiba vezető rendben tehát

ϵ_{a} ≈ 2h^{5} · Φ ≈ 2ϵ_{h} és ϵ_{b} ≈ 32h^{5} · Φ ≈ 32ϵ_{h}  (1.64)

Ezekből az aktuális lépésköz melletti tipikus hiba hatodrendben előáll:

ϵ_{h} ≈ \frac{1}{30} |ϵ_{b} − ϵ_{a}| ≈ \frac{1}{30} |y_{a} − y_{b}| + O(h^{6})  (1.65)

Ugyan a hatodrendű tag nem feltétlenül jelent elhanyagolható mennyiséget, mégis általában elfogadhatjuk, hogy az O(h^{6}) korrekciótól eltekintsünk. Mostmár megvizsgálhatjuk, hogy mekkora h_{opt} lépésközre érdemes áttérni, hiszen ötödrendű hibatagról lévén szó:

(\frac{h_{opt}}{h} )^{5} = \frac{ϵ_{opt}}{ϵ_{h}} h_{opt} = h (\frac{ϵ_{opt}}{ϵ_{h}} )^{0.2}  (1.66)

Ha h_{opt} nem tér el lényegesen h-tól, akkor megtartjuk a kiszámolt y_{n+2} = y_{a} számot, ha lényegesen kisebb, akkor újra számolunk kisebb lépésközzel, ha lényegesen nagyobb, akkor a következő lépésekben növelt lépéshosszal próbálkozunk. Természetesen vigyázzunk arra, hogy h_{opt} ne legyen sem túl nagy, sem túl kicsiny.

A vázolt eljárás túlmunkával jár, hiszen egy RK lépés során 4 alkalommal kell függvényértéket számolni. A két h nagyságú lépéshez szükséges 8 függvényhívás helyett most a 2h hosszú extra lépés miatt további három közbülső pontban, azaz összesen 11 alkalommal számoltuk ki f(x, y)-t. Ez mintegy 1.3-szoros számolási túlmunkát jelent. A ’rázós’ intervallumokon azonban nem nő meg ellenőrizetlenül a hiba, míg a ’sima’ szakaszokon olyan nagy lépésekkel haladhatunk, hogy összességében nagyságrendekkel csökkenhet a számolási munka. Példa: Beágyazott Runge-Kutta

Egy M-edrendű Runge Kutta eljárásban a számolási séma

k_{1} = hf(x_{n}, y_{n})  (1.67)
k_{2} = hf(x_{n} + a_{2}h , y_{n} + b_{21}k_{1})  (1.68)
⋮  (1.69)
k_{M} = hf(x_{n} + a_{M}h , y_{n} + b_{M1}k_{1} + ... + b_{M,M−1}k_{M−1})  (1.70)
y_{n+1} = y_{n} + ∑_{i=1}^{M} c_{i}k_{i} + O(h^{M+1})  (1.71)

Ha M > 4, akkor lehetőség van arra, hogy másképp is elvégezzük a lépést ugyanazokkal a k_{i} számokkal de más együtthatókkal úgy, hogy

ỹ_{n+1} = y_{n} + ∑_{i=1}^{M} c̃_{i}k_{i} + O(h^{M})  (1.72)

legyen. Megfelelő együtthatókat (táblázatosan adott) választva a lépés tipikus hibája megbecsülhető:

ϵ_{h} = |y_{n+1} − ỹ_{n+1}|  (1.73)

Az egyik leggyakoribb ilyen beágyazott RK eljárás a Fehlberg 4-5 rendű Runge-Kutta módszer, ahol az ötödrendű RK mellett ugyanazon függvényértékekkel egy negyedredű lépést is teszünk és ebből becsüljük a hibát. Ilyen a MAPLE rkf45 rutinja.

II. PARCIÁLIS DIFFERENCIÁLEGYENLETEK

A. A lineáris egyenletek osztályozása

A fizikában leggyakrabban másodrendű parciális differenciálegyenletek megoldása válik szükségessé. A lineáris 2 dimenziós egyenletek általános alakja

[a\frac{∂^{2}}{∂x^{2}} + 2b\frac{∂^{2}}{∂x∂y} + c\frac{∂^{2}}{∂y^{2}} + d\frac{∂}{∂x} + e\frac{∂}{∂y} + f] V (x, y) = g  (2.1)

Az a, b, c, d, e, f, g együtthatók maguk is (x, y) adott függvényei lehetnek. Ha az együtthatók függenek az ismeretlen V (x, y)−tól is, akkor kvázilineáris egyenletről beszélünk. Konstans együtthatók esetén a három másodrendű parciális derivált együtthatói alapján a következő táblázat szerint csoportosítjuk az egyenleteket:

feltétel  egyenlet típusa  jellegzetes fizikai problémakör  tipikus egyenlet tipikus határfeltétel
b^{2} < ac  elliptikus  Laplace/(P oisson)-egyenlet  \frac{∂^{2}V}{∂x^{2}} + \frac{∂^{2}V}{∂y^{2}} = 0  peremérték
b^{2} > ac  hiperbolikus  Hullámegyenlet  \frac{∂^{2}V}{∂x^{2}} − \frac{1}{v^{2}} \frac{∂^{2}V}{∂t^{2}} = 0 kezdetiérték
b^{2} = ac  parabolikus  Diffúzió/(Schrödinger)-egyenlet λ\frac{∂^{2}V}{∂x^{2}} − \frac{∂V}{∂t} = 0 kezdetiérték

(2.2)

A csoportosítás nem konstans együtthatók mellett is tartható, de akkor csak lokálisan értendő. A fizikában gyakori egyenletcsoportok jól osztályozhatók a táblázat alapján, ezek típusfeladatok a megoldási módszer szempontjából is nagyon elkülönülnek.

• A hiperbolikus és a parabolikus egyenleteink egyik változója általában az idő. Többnyire kezdetiértékproblémával állunk szemben: adott valamely t = t_{0}-kor a V (x, t_{0}) kezdőmegoldás és ki kell számolnunk, hogy egy későbbi időpontban mi lesz.V (x, t). Lényegében ez nem különbözik a már megismert integrálási eljárásoktól, az időbeli lépéseket diszkretizálva integrálhatunk előre.

• Az elliptikus egyenletek többnyire csak térszerű változókat tartalmaznak. Valamely zárt határfelületen adott peremfeltételek mellett kell megoldanunk az egyenletet az ismeretlen V (x, y) előállításása végett. A peremértékfeladatoknál nem lehet egyszerűen ’beintegrálni’ a peremről, a megoldást a tér minden pontját egyformán figyelve globális módszerrel kell megközelíteni.

B. Hiperbolikus egyenletek

A (homogén) hullámegyenlet (konstans sebesség mellet) átírható

[\frac{∂^{2}}{∂t^{2}} − c^{2} \frac{∂^{2}}{∂x^{2}} ] u(x, t) = 0 ∼ [\frac{∂}{∂t} + c\frac{∂}{∂x}] [\frac{∂}{∂t} − c\frac{∂}{∂x}] u(x, t) = 0  (2.3)

alakra, ahonnan látható, hogy a feladat az

[\frac{∂}{∂t} + c\frac{∂}{∂x}] z(x, t) = 0 , [\frac{∂}{∂t} − c\frac{∂}{∂x}] u(x, t) = z(x, t)  (2.4)

elsőrendű egyenletek rendszerére redukálható. Másik lehetséges felbontás

s(x, t) = c\frac{∂u}{∂x} , r(x, t) = \frac{∂u}{∂t} , \frac{∂r}{∂t} = c\frac{∂s}{∂x} , \frac{∂s}{∂t} = c\frac{∂r}{∂x}  (2.5)

Az általános hullámegyenlet (inhomogén, nem-konstans együtthatók) tárgyalása helyett a numerikus megoldási módszerek alapötletét és a felmerülő problémákat egy egyszerűbb modellen vizsgáljuk meg. A következőkben a

[\frac{∂}{∂t} + v \frac{∂}{∂x}] z(x, t) = 0 , z(x, t_{0}) = z_{0}(x)  (2.6)

kezdetiérték problémát fogjuk tanulmányozni.

A véges differenciákra való áttérés végett ekvidisztáns

x_{j} = x_{0} + jδx , t_{n} = t_{0} + nδt , z_{j}^{n} ≡ z(x_{j}, t_{n})  (2.7)

beosztást veszünk mindkét (x és t) tengelyen. Miután a t_{n}-ről való időben előre történő lépésünkkör az összes z_{1}^{n}, z_{2}^{n}, ... rendelkezésre áll, az x szerinti deriváltat tetszőlegesen sok pontból felírhatjuk. Az idő szerinti deriváltra azonban csak korábbi időpontokról van információnk, így explicit rekurzióhoz kézenfekvő egyszerű előre deriváltat használni. A kétpontos időbeli első derivált mellett nincs szükség arra, hogy a hárompontos szimmetrikus formulánál pontosabbat használjunk az x változóban. Kompromisszumként legyen a disztkretizált egyenlet tehát

|[\frac{∂}{∂t} + v \frac{∂}{∂x}] z(x, t)|_{j,n} = \frac{z_{j}^{n+1} − z_{j}^{n}}{δt} + O(δt) + v_{j}^{n} · \frac{z_{j+1}^{n} − z_{j−1}^{n}}{2δx} + O(δx^{2}) = 0  (2.8)

Ebből az FTCS (Forward Time, Centered Space) rekurzió

z_{j}^{n+1} = z_{j}^{n} − \frac{1}{2}v_{j}^{n} (z_{j+1}^{n} − z_{j−1}^{n}) \frac{δt}{δx}  (2.9)

Sajnos ez az egyszerű formula nemcsak pontatlan, de használhatatlanul instabil, nagyon ’rövid’ idő alatt a megoldás ’elszáll’. A von Neumann stabilitásvizsgálthoz tekintsük a v(x, t) ≡ c speciális hullámegyenlet

z(x, t) = e^{ik(x−ct)}  (2.10)

síkhullám alap-megoldásait. Ekkor

z_{j}^{n} = w(k)^{n} · e^{ikx_{j}}  (2.11)

ahol a w(k) egységnyi modulusú komplex szám n-edik hatványa szerepel. A rekurziós egyenletbe beírva

w(k)^{n+1} = [1 − i\frac{cδt}{δx}sin(kδx)] w(k)^{n} ⟹ w(k) = [1 − i\frac{cδt}{δx}sin(kδx)] ⟹ ∥w(k)∥ > 1  (2.12)

mégis azt látjuk, hogy minden k-val vett alapmegoldásra a formulánk feltételenül instabil. Az iteráció a parciális hullámok amplitudóját különböző mértékben felnöveli az idő haladtával. Megfontolásunkat a homogén hullámegyenletek közül is csak a konstans sebességgel felírtakra végeztük el. Mégis, lassan változó v(x, t) sebesség mellett és akár inhomogén egyenletek esetén is elfogadhatjuk, hogy ez az eljárás nem stabil.

Viszonylag könnyen segíthetünk az imént feltárt instabilitáson a Lax módszer segítségével. A rekuziónkban a jobboldalon végezzük el a

z_{j}^{n} ⟶ \frac{z_{j+1}^{n} + z_{j−1}^{n}}{2}  (2.13)

helyettesítést, azaz z_{j}^{n}-t két térbeli szomszédja átlagára cseréljük le. Ezzel kapjuk a Lax formulát

z_{j}^{n+1} = \frac{z_{j+1}^{n} + z_{j−1}^{n}}{2} − \frac{v}{2} \frac{δt}{δx} (z_{j+1}^{n} − z_{j−1}^{n})  (2.14)

A stabilitás vizsgálatát ismét a von Neumann módszerrel vizsgálva most azt kapjuk, hogy

w(k) = [cos(kδx) − i\frac{cδt}{δx}sin(kδx)]  (2.15)

∥w(k)∥^{2} = cos^{2}(kδx) + (c\frac{δt}{δx})^{2} sin^{2}(kδx) = 1 − [1 − (c\frac{δt}{δx})^{2}] sin^{2}(kδx)  (2.16)

Minden k-ra stabil az eljárásunk, ha

|v| \frac{δt}{δx} ≤ 1 vagy |v| δt ≤ δx  (2.17)

teljesül. Ez a Courant-Friedrichs-Lewy stabilitási feltétel. Fizikai magyarázat: a hullámnak (x_{j±1}, t_{n})-ből oda kell érnie x_{j}, t_{n+1}-be. Matematikai magyarázat: A Lax formulát kicsit átírva

\frac{z_{j}^{n+1} − z_{j}^{n}}{δt} = \frac{(δx)^{2}}{2δt} · \frac{z_{j+1}^{n} − 2z_{j}^{n} + z_{j−1}^{n}}{(δx)^{2}} − v (\frac{z_{j+1}^{n} − z_{j−1}^{n}}{2δx} )  (2.18)

észrevehetjük, hogy ez a

[\frac{∂}{∂t} + v \frac{∂}{∂x}] z(x, t) = λ\frac{∂^{2}}{∂x^{2}} z(x, t) , λ = \frac{(δx)^{2}}{2δt}  (2.19)

differenciálegyenlet diszkretizációjából is jöhetett volna.A jobboldali diffúziós tagot jelentő második derivált például súrlódó folyadékok áramlásánal jelenik meg a fizikai egyenletekben, lényegében disszipációs tagot jelent. A stabiltás feltétele

λ = \frac{δx}{δt} \frac{δx}{2} ≥ \frac{δx}{2 |v|}  (2.20)

tehát azt jelenti, hogy elegendően erős numerikus disszipációt, numerikus viszkozitást imitáltunk. A súrlódásra szükségünk volt, hogy stabil legyen az eljárás, ugyanakkor a nagy súrlódás ’elemészti’ a megoldásainkat. Ha olyan problémákat vizsgálunk, ahol a lényegesnek vélt komponensekre teljesül, hogy kδx ≫ 1, akkor ezekre a komponensekre mindkét eljárás használható: ∥w(k)∥^{2} ≈ 1. A rövid hullámhosszú, azaz kδx ∼ 1 komponensek azonban az első eljárásban felnövekszenek és uralkodni kezdenek, míg az utóbbiban kihalnak.

Amiatt, hogy az időbeli deriváltat kétpontos formulával írtuk fel, míg a térbelit hárompontossal, a pontosság miatt a δt beosztás finomságára jobban kell ügyelnünk, mint δx-re. Ez alkalmasint azzal járhat, hogy |v| δt ≪ δx lesz, ami a Courant feltételnél indokolatlanul erősebb megszorítás. A leapfrog (bakugrásos) módszerrel segíthetünk ezen, olyan egyenletet írunk fel, hogy az mindkét változóban másodrendű hibájú legyen:

|[\frac{∂}{∂t} + v \frac{∂}{∂x}] z(x, t)|_{j,n} = \frac{z_{j}^{n+1} − z_{j}^{n−1}}{2δt} + O(δt^{2}) + v_{j}^{n} · \frac{z_{j+1}^{n} − z_{j−1}^{n}}{2δx} + O(δx^{2}) = 0  (2.21)

z_{j}^{n+1} = z_{j}^{n−1} − \frac{δt}{δx}v_{j}^{n} (z_{j+1}^{n} − z_{j−1}^{n})  (2.22)

Ebben a rekurzióban szükség van arra, hogy a korábbi időpillantra kiszámolt z^{n−1} értékeket megint elővegyük. A stabilitásvizsgálatot a korábbiakhoz hasonlóan elvégezve kapjuk, hogy

w^{2} = 1 − w2i\frac{cδt}{δx}sin(kδx)  (2.23)

Kis kézimunkával megmutathatjuk, hogy

\frac{cδt}{δx} ≤ 1 ⟹ ∥w(k)∥ = 1  (2.24)

tehát a Courant feltétel teljesülése szükséges, de ilyenkor sem erősítés, sem disszipáció nincs. Az egyenlet annál is inkább érdekes, mert az eredeti másodrendű hullámegyenletünkre visszavezetve annak közvetlen diszkretizálásával ekvivalens

[\frac{∂^{2}u}{∂t^{2}} − v^{2} \frac{∂^{2}u}{∂x^{2}}]_{j,n} = \frac{u_{j}^{n+1} − 2u_{j}^{n} + u_{j}^{n−1}}{(δt)^{2}} − v^{2} \frac{u_{j+1}^{n} − 2u_{j}^{n} + u_{j−1}^{n}}{(δx)^{2}} = 0  (2.25)

A bakugrásos módszerrel bakot lőhetünk bonyolultabb egyenletekre. A hiba egyik forrása éppen a bakugrás, hiszen a páros és a páratlan osztópontok (akár a sakktábla világos és sötét mezői) tejesen szétcsatoltak az egyenletben. A megoldás a kétlépéses Lax-Wendroff módszer Az első lépésben Lax módszerrel csak fele akkor beosztású rácson lépünk a centrális pontba

z_{j+1/2}^{n+1/2} = \frac{z_{j+1}^{n} + z_{j}^{n}}{2} − \frac{v}{2} \frac{δt}{δx} (z_{j+1}^{n} − z_{j}^{n})  (2.26)

majd a centrális értékkel bakugrás következik

z_{j}^{n+1} = z_{j}^{n} − \frac{δt}{δx}v_{j+1/2}^{n+1/2} (z_{j+1/2}^{n+1/2} − z_{j−1/2}^{n+1/2})  (2.27)

Az eljárás stabilitásának feltétele megintcsak a |v| δt ≤ δx Courant szabály. Csillapítás itt is előfordul a nagy hullámszámokra, de ez a disszipáció kisebb, mint az eredeti Lax módszernél.

C. Parabolikus egyenletek

A diffúziós egyenlet

\frac{∂V}{∂t} = λ\frac{∂^{2}V}{∂x^{2}}  (2.28)

megoldásával foglalkozunk először, miközben a λ > 0 diffúziós együtthatót állandónak vesszük. Diszkretizálva a feladatot V_{j}^{n} ≡ V (x_{j}, t_{n}) jelöléssel az FTCS egyenletünk

\frac{V_{j}^{n+1} − V_{j}^{n}}{δt} = λ\frac{V_{j+1}^{n} − 2V_{j}^{n} + V_{j−1}^{n}}{(δx)^{2}}  (2.29)
V_{j}^{n+1} = V_{j}^{n} + α [V_{j+1}^{n} − 2V_{j}^{n} + V_{j−1}^{n}]  ,  α = \frac{λδt}{(δx)^{2}}  (2.30)

A stabilitásvizsgálathoz vegyük a diffúziós egyenletnek egy a térben síkhullám alapmegoldását

V (x, t) = w(t) · e^{ikx} , w(t) = e^{−λk^{2}t}  (2.31)

és vizsgáljuk meg, hogy a

V_{j}^{n} = w^{n} · e^{ikx_{j}}  (2.32)

helyettesítéssel w(t)-re a differencia egyenlet milyen viselkedést eredményezne. Azt kapjuk, hogy

w = 1 + α [e^{ikδx} − 2 + e^{−ikδx}] = 1 + 2α [cos(kδx) − 1] = 1 − 4αsin^{2}(\frac{kδx}{2})  (2.33)

amiből

4α sin^{2}(\frac{kδx}{2}) ≤ 2  (2.34)

A stabilitás feltétele ezek szerint, hogy

α ≤ 1/2 azaz δt ≤ \frac{1}{2} \frac{(δx)^{2}}{λ} ∼ τ  (2.35)

ahol τ a δx távolságra való diffúzió karakterisztikus ideje. A feltételes stabilitás ugyan biztató,.mégis a módszer alkalmazása a gyakorlatban nem cálszerű. A reális problémákra jellemző L tipikus hosszúság mellett a diffúziós karakterisztikus idő

τ ∼ \frac{L^{2}}{λ} = \frac{L^{2}}{(δx)^{2}} \frac{(δx)^{2}}{λ}  (2.36)

Miközben értelemszerűen a pontosabb számoláshoz L/δx ≫ 1 térbeli beosztás szükséges, a τ ideig való szimulációhoz L^{2}/ (δx)^{2} nagyságrendben kell időbeli lépést tenni.

Lépéshossztól függetlenül stabil egyenletekhez jutunk, ha a bal oldalon a térváltozó szerinti második deriváltat nem az iménti módon a t_{n}-ben, hanem a t_{n+1}-ben írjuk fel:

\frac{V_{j}^{n+1} − V_{j}^{n}}{δt} = λ\frac{V_{j+1}^{n+1} − 2V_{j}^{n+1} + V_{j−1}^{n+1}}{(δx)^{2}}  (2.37)

Így egy teljesen implicit egyenletrendszert kapunk:

−αV_{j+1}^{n+1} + (1 + 2α)V_{j}^{n+1} − αV_{j−1}^{n+1} = V_{j}^{n} , j = 1, 2, ...., J − 1  (2.38)

A stabilitás vizsgálatánál kapjuk, hogy

w [1 + 2α − 2αcos(kδx)] = 1 w = \frac{1}{1 + 4αsin^{2}(\frac{kδx}{2})}  (2.39)

azaz az eljárás minden lépésközzel stabil, |w| ≤ 1. Az implicit egyeletrendszer vektor jelölésben

A · V^{n+1} = V^{n}  (2.40)

ahol A tridiagonális mátrix:

A = ( \begin{matrix}1 & 0 & · & · & · & 0 \\ −α & 1 + 2α & −α & 0 & · & 0 \\ 0 & −α & 1 + 2α & −α & 0 & · \\ · & 0 & · & · & · & · \\ · & · & · & −α & 1 + 2α & −α \\ · & · & · & 0 & 0 & 1 \end{matrix} )  (2.41)

Az ilyen egyenletet például rekurzióval minden időbeli lépésre meg tudjuk oldani.

Az eddigi két eljárás elsőrendben volt pontos. Figyelembe véve, hogy a kétpontos első derivált másodrendben pontos szimmetrikus formulaként is felfogható a köztes t_{n+1/2} pillanatra, másodrendben pontos egyenletet kapunk, ha a jobboldalt is így fogjuk fel

\frac{∂V}{∂t} |_{j,n+1/2} = λ \frac{∂^{2}V}{∂x^{2}} |_{j,n+1/2} ≈ \frac{λ}{2} [ \frac{∂^{2}V}{∂x^{2}} |_{j,n+1} + \frac{∂^{2}V}{∂x^{2}} |_{j,n}]  (2.42)

A deriváltakat diszkretizálva kapjuk

\frac{V_{j}^{n+1} − V_{j}^{n}}{δt} = \frac{λ}{2 (δx)^{2}} [(V_{j+1}^{n+1} − 2V_{j}^{n+1} + V_{j−1}^{n+1}) + (V_{j+1}^{n} − 2V_{j}^{n} + V_{j−1}^{n})]  (2.43)

ami átrendezve a Crank-Nicholson formula

V_{j}^{n+1} − \frac{α}{2} (V_{j+1}^{n+1} − 2V_{j}^{n+1} + V_{j−1}^{n+1}) = V_{j}^{n} + \frac{α}{2} (V_{j+1}^{n} − 2V_{j}^{n} + V_{j−1}^{n})  (2.44)
−αV_{j+1}^{n+1} + (2 + 2α)V_{j}^{n+1} − αV_{j−1}^{n+1} = αV_{j+1}^{n} + (2 − 2α)V_{j}^{n} + αV_{j−1}^{n}  (2.45)

vagy mátrix alakban

A · V^{n+1} = B · V^{n}  (2.46)

rősítése

w = \frac{1 − α + αcos(kδx)}{1 + α − αcos(kδx)} = \frac{1 − 2αsin^{2}(\frac{kδx}{2})}{1 + 2αsin^{2}(\frac{kδx}{2})} ≤ 1  (2.47)

tehát a módszer stabil.

Az implicit módszerekben a tridiagonális mátrixokkal való munka időigényes lehet. Kis módosítással olyan egyenletet is felírhatunk, ahol hasonló pontossággal (és stabilitással) explicit formulát kapunk. Ilyen a Dufort-Frankel módszer. Továbbra is másodrendben korrekt idő szerinti deriváltat írunk fel, a jobb oldalon a tér szerinti második deriváltban pedig a

2V_{j}^{n} ⟶ V_{j}^{n+1} + V_{j}^{n−1}  (2.48)

helyettesítéssel élünk. Ekkor a

\frac{V_{j}^{n+1} − V_{j}^{n−1}}{2δt} = λ\frac{V_{j+1}^{n} − V_{j}^{n+1} − V_{j}^{n−1} + V_{j−1}^{n}}{(δx)^{2}}  (2.49)
V_{j}^{n+1} = V_{j}^{n−1} + 2α [V_{j+1}^{n} − V_{j}^{n+1} − V_{j}^{n−1} + V_{j−1}^{n}]  (2.50)

látszólag implicit egyenletet kapjuk. Az ismeretlen V_{j}^{n+1} azonban triviális módon explicite kifejezhető:

V_{j}^{n+1} = (\frac{1 − 2α}{1 + 2α}) V_{j}^{n−1} + (\frac{2α}{1 + 2α}) [V_{j+1}^{n} + V_{j−1}^{n}]  (2.51)

A stabilitásvizsgálakor másodfokú egyenletet kapunk w-re

w^{2} = (\frac{1 − 2α}{1 + 2α}) + 2w (\frac{2α}{1 + 2α}) cos(kδx  (2.52)

amiből

w = \frac{1}{1 + 2α} (2αcos(kδx) ± \sqrt{1 − 4α^{2}sin^{2}(kδx)})  (2.53)

Ha 4α^{2}sin^{2}(kδx) ≤ 1, akkor a gyök valós, és

− (2α + 1) ≤ cos(kδx) ± \sqrt{1 − 4α^{2}sin^{2}(kδx)} ≤ (2α + 1) miatt |w| ≤ 1  (2.54)

Komplex gyökre, amikor 4α^{2}sin^{2}(kδx) > 1

∥w∥^{2} = \frac{1}{(1 + 2α)^{2}} (cos^{2}(kδx) + 4α^{2}sin^{2}(kδx) − 1) = \frac{4α^{2} − 1}{(1 + 2α)^{2}} = \frac{2α − 1}{2α + 1} miatt ∥w∥^{2} ≤ 1  (2.55)

Tehát az eljárás minden beosztás mellett stabil. A Dufort-Frankel módszernek is vannak hátrányai. Egyrészt az eljárás során két megelőző időpontban kell ismernünk a megoldást a térbeli beosztáson, ami a V_{j}^{n−1} értékek tárolását kívánja meg. Ez jelentős tárolókapcitást igényelhet. Továbbá kellemetlen, hogy kezdőfeltételként is két időpontban is meg kell adnunk a V_{j}^{n−1} értékeket, ami indokolatlan egy időben elsőredű differenciálegyenletnél. Példa: Időfüggő Scrödinger egyenlet 1 dimenzióban

Egy ψ(x, t) ’hullámcsomag’ szóródását a U(x) potenciálon a

i\frac{∂ψ}{∂t} = −\frac{∂^{2}ψ}{∂x^{2}} + U(x)ψ , ψ(x, t = t_{0}) = ψ_{0}(x) , ψ(±∞, t) = 0  (2.56)

kezdetiérték probléma írja le. Az előzőekben tárgyalt egyenletekhez képest ez az egyenlet alapvetően két dologban új:

a) a λ együttható imaginárius, b) a potenciállal kapcsolatos tag eddig nem szerepelt. Az egyenletek természetesen általánosíthatók erre az esetre is. Tekintsük a Crank-Nicholson eljárást, miszerint

\frac{ψ_{j}^{n+1} − ψ_{j}^{n}}{δt} = \frac{i}{2} [\frac{ψ_{j+1}^{n+1} − 2ψ_{j}^{n+1} + ψ_{j−1}^{n+1}}{(δx)^{2}} − U_{j}(x)ψ_{j}^{n+1} + \frac{ψ_{j+1}^{n} − 2ψ_{j}^{n} + ψ_{j−1}^{n}}{(δx)^{2}} − U_{j}(x)ψ_{j}^{n}]  (2.57)

vagy rendezve

−aψ_{j+1}^{n+1} + (1 + 2a + b_{j}) ψ_{j}^{n+1} − aψ_{j−1}^{n+1} = aψ_{j+1}^{n} + (1 − 2a − b_{j}) ψ_{j}^{n} + aψ_{j−1}^{n} ; a = \frac{iδt}{2 (δx)^{2}} ; b_{j} = \frac{iδt}{2}U(x_{j}) (2.58)

A tridiagonális egyenletrendszer numerikusan hatékonyan invertálható.

Tanulságos az egyenletet egy másik származtatását is megvizsgálni. Az időfüggő Schrödinger egyenlet formális megoldásával

e^{−iHδt}ψ(x, t_{n}) = ψ(x, t_{n+1})  (2.59)

Ha ebben az egyenletben az

e^{−iHδt} ≈ 1 − iHδt ⟹ ψ_{j}^{n+1} = (1 − iHδt) ψ_{j}^{n}  (2.60)

sorfejtést helyettesítjük, majd a térváltozóban második deriváltat centrált véges differenciával írjuk fel, akkor az FTCS egyenletet kapjuk. Imaginárius λ-val az FTCS azonban instabil. A teljesen implicit egyenletet akkor kapjuk, ha az alternatív

(1 − iHδt)^{−1} ψ_{j}^{n+1} = ψ_{j}^{n}  (2.61)

felírást választjuk. Ez az eljárás stabil imaginárius λ-val is, de nem unitér (ahogy az FTCS sem). Ez azzal jár együtt, hogy az eredetileg normált hullámfüggvény normája elromlik a szimuláció során. Ha a Cayley sorfejtést alkalmazzuk:

e^{−iHδt} ≈ \frac{1 − \frac{1}{2}iHδt}{1 + \frac{1}{2}iHδt} ⟹ (1 + \frac{1}{2}iHδt) ψ_{j}^{n+1} = (1 − \frac{1}{2}iHδt) ψ_{j}^{n}  (2.62)

akkor másodrendben pontos, stabil és unitér eljárást kapunk. A véges differenciákat behelyettesítve ez az egyenlet éppen a Crank-Nicholson formulára vezet.

D. Kiegészítés: Tridiagonális egyenletrendszer megoldása rekurzióval

Tekintsük az

A · x = d  (2.63)

egyenletrendszert, ahol A egy tridiagonális mátrix, d adott vektor és x a meghatározandó ismeretlen. Komponensekben:

( \begin{matrix}β_{1} & γ_{1} & 0 & · & · & · \\ α_{2} & β_{2} & γ_{2} & 0 & · & · \\ 0 & α_{3} & β_{3} & γ_{3} & 0 & · \\ · & · & · & · & · & · \\ · & · & 0 & α_{n−1} & β_{n−1} & γ_{n−1} \\ · & · & · & 0 & α_{n} & β_{n} \end{matrix} ) ( \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \\ · \\ x_{n−1} \\ x_{n} \end{matrix} ) = ( \begin{matrix}d_{1} \\ d_{2} \\ d_{3} \\ · \\ d_{n−1} \\ d_{n} \end{matrix} )  (2.64)

vagy

α_{i}x_{i−1} + β_{i}x_{i} + γ_{i}x_{i+1} = d_{i} , i = 1, ..., n ; α_{1} = γ_{n} = 0  (2.65)

Vezessük be a g_{i} és h_{i} segédváltozókat úgy, hogy

x_{i+1} = g_{i}x_{i} + h_{i} , i = 1, ..., n − 1  (2.66)

legyen. Ezekkel

d_{i} = α_{i}x_{i−1} + β_{i}x_{i} + γ_{i}x_{i+1}  (2.67)
= α_{i}x_{i−1} + β_{i}x_{i} + γ_{i} (g_{i}x_{i} + h_{i})  (2.68)
= α_{i}x_{i−1} + (β_{i} + γ_{i}g_{i}) x_{i} + γ_{i}h_{i}  (2.69)
= α_{i}x_{i−1} + (β_{i} + γ_{i}g_{i}) (g_{i−1}x_{i−1} + h_{i−1}) + γ_{i}h_{i}  (2.70)
= (α_{i} + g_{i−1} (β_{i} + γ_{i}g_{i})) x_{i−1} + (β_{i} + γ_{i}g_{i}) h_{i−1} + γ_{i}h_{i}  (2.71)

ami kielégíthető, ha

0 = α_{i} + g_{i−1} (β_{i} + γ_{i}g_{i}) ⟹ g_{i−1} = −\frac{α_{i}}{β_{i} + g_{i}}  (2.72)
d_{i} = (β_{i} + γ_{i}g_{i}) h_{i−1} + γ_{i}h_{i} ⟹ h_{i−1} = \frac{d_{i} − γ_{i}h_{i}}{β_{i} + γ_{i}g_{i}}  (2.73)

Láthatóan g_{n−1} és h_{n−1} ismeretében ’lelfelé’ rekurzív módón előállítható az összes g és h. Ellenőrízhető közvetlen kiírással, hogy az induló értékek:

g_{n−1} = −\frac{α_{n}}{β_{n}} és h_{n−1} = \frac{d_{n}}{β_{n}}

Az eredeti egyenletrendszer első egyenletéből

β_{1}x_{1} + γ_{1}x_{2} = d_{1} ⟹ β_{1}x_{1} + γ_{1} (g_{1}x_{1} + h_{1}) = d_{1}  (2.74)

azaz g_{1} és h_{1} segitségével x_{1} kifejezhető:

x_{1} = \frac{d_{1} − γ_{1}h_{1}}{β_{1} + γ_{1}g_{1}}  (2.75)

Erről az értékről most ’felfele’ használva a

x_{i+1} = g_{i}x_{i} + h_{i}  (2.76)

rekurziót minden keresett x_{i} kiszámolható.

E. Kiegészítés: a diffúziós probléma több dimenzióban

A diffúziós egyenlet 2-dimenzióban

\frac{∂V}{∂t} = λ (\frac{∂^{2}V}{∂x^{2}} + \frac{∂^{2}V}{∂y^{2}} )  (2.77)

Crank-Nicholson szerint diszkretizálva a V_{j,l}^{n} ≡ V (x_{j}, y_{l}, t_{n}) jelöléssel

V_{j,l}^{n+1} − V_{j,l}^{n} = \frac{α}{2} (L_{x}V_{j,l}^{n+1} + L_{x}V_{j,l}^{n} + L_{y}V_{j,l}^{n+1} + L_{y}V_{j,l}^{n})

vagy

(2 − α [L_{x} + L_{y}]) V_{j,l}^{n+1} = (2 + α [L_{x} + L_{y}]) V_{j,l}^{n}

ahol

δx = δy = Δ , α = \frac{λδt}{Δ^{2}} , L_{x}V_{j,l} = [V_{j+1,l} − 2V_{j,l} + V_{j−1,l}] , L_{y}V_{j,l} = [V_{j,l+1} − 2V_{j,l} + V_{j,l−1}]  (2.78)

Az eljárás térben és időben egyaránt másodrendben pontos továbbá feltétel nélkül stabil. Rendezzük a kétindexes mennyiségeket egyetlen vektorba valamely j, l ⟶ s megfeleltetéssel. Így az ismeretlen v⃗ ∼ {V_{j,l}^{n+1}} vektorra az egyenletrendszer

A · v⃗ = b⃗  (2.79)

alakú. Sajnos, az 1-dimenziós problémával szemben a A most nem tridiagonális, szerencsére azonban ritka mátrix. Az ilyen típusú egyenletek numerikusan standard módszerekkel (Jacobi, Gauss-Seidel és SOR) megoldhatók, erről külön kiegészítésben szólunk.

Kissé átalakítva a CN egyenletet alternatív rokonszenvesebb eljáráshoz jutunk. A váltakozó irányú implicit: ADI (alternating-direction implicit) módszerben a teljes időlépést két részre bontjuk és mindkét fél lépésben más térbeli dimenziót kezelünk implicit módszerrel

V_{j,l}^{n+1/2} − V_{j,l}^{n} = \frac{α}{2} (L_{x}V_{j,l}^{n+1/2} + L_{y}V_{j,l}^{n})  (2.80)
V_{j,l}^{n+1} − V_{j,l}^{n+1/2} = \frac{α}{2} (L_{x}V_{j,l}^{n+1/2} + L_{y}V_{j,l}^{n+1})

Az eljárás pontossága és stabilitása az előző direkt formulával azonos, de az egyenletek megoldása könyebb. Látható ugyanis, hogy ’csak’ két tridiagonális egyenletrendszert kell minden lépésben megoldanunk.

Az ADI módszerrel rokon az OS ’operator splitting’. Tegyük fel, hogy a kezdetiérték probléma jobb oldalán álló differenciáloperátor valamely operátorok összegeként írható fel:

\frac{∂V}{∂t} = λDV = λ (D_{1}V + ... + D_{m}V)  (2.81)

Tegyük fel továbbá, hogy az itt szereplő mindegyik D_{i} olyan, hogy a

\frac{∂V}{∂t} = λD_{i}V  (2.82)

egyenleből viszonylag egyszerű lenne diszkretizált eljárással

v^{[n+1]} = F_{i}(v^{[n]}, δt)  (2.83)

időbeli lépéseket tenni. Ilyenkor egy teljes időbeli lépést feloszthatunk m darab kisebb lépésre és a

v^{[n+ \frac{1}{m}]} = F_{1}(v^{[n]}, \frac{δt}{m})  (2.84)
v^{[n+ \frac{2}{m}]} = F_{2}(v^{[n+ \frac{1}{m}]}, \frac{δt}{m})  (2.85)
⋮  (2.86)
v^{[n+1]} = F_{m}(v^{[n+ \frac{m−1}{m}]}, \frac{δt}{m})  (2.87)

módon állítjuk elő az új értékeket.

Az ADI és az OS rokonsága abban áll, hogy az iménti egyenletben az F_{i} léptetési utasítást másképp is felfoghatjuk. Például lehet F_{1} egy diszkretizált lépés a teljes D operátorral, de pl. úgy felírva, hogy az csak D_{1}-ben stabil, míg F_{2} az D_{2}-ben stabil lépés, és így tovább. Így az előbbi OS egyenlet a korábbi ADI sémát adja.

F. Kiegészítés: Av=b lineáris egyenletrendszer megoldása

Bontsuk fel az A mátrixot a következőképpen

A = L + D + U  (2.88)

ahol D az A diagonális része, L az alsó triangulárisa és U pedig a felső trianguláris. Jacobi iteráció: Az egyenlet átrendezve

D · v = − [L + U] ·v + B = (D − A) ·v + B ⟹ v = G · v + D^{−1}B , G_{J}= I − D^{−1}A  (2.89)

szucesszív módszerrel oldjuk meg: induljunk ki valamely V^{(0)}-ból és iteráljunk

v^{(n+1)} = G_{J}·v^{(n)}+D^{−1}B  (2.90)

Az eljárás konvergenciája a G_{J} mátrix sajátértékeitől függ. A legnagyobb sajátérték abszolut értékének egynél kisebbnek kell lennie a konvergenciához. Ha G_{J} legnagyobb sajátértéke (spektrál sugara) λ_{J}, akkor asszimptotikusan

q_{J} = \frac{|v^{(n+1)} − v^{(n)}|}{|v^{(n)} − v|} ≈ |1 − λ_{J}|  (2.91)

Az eljárás nagyon lassan kovergál, ha a sajátértékek között van egyhez közeli. Gauss-Seidel iteráció: A baloldalra átvihetjük az alsó trianguláris részt, ha az algoritmust úgy szervezzük, hogy v_{1}^{(n+1)}, v_{2}^{(n+1)}, .... sorrendben határozzuk meg a vektor elemeit. Az így kapott iterációs séma

[L + D] ·v^{(n+1)} = −U · v^{(n)}+B  (2.92)

A konvergencia most a

G_{GS}= − [L + D]^{−1} U  (2.93)

mátrix sajátértékeitől függ. Megmutatható, hogy

λ_{GS} = λ_{J}^{2}  (2.94)

így az eljárás az előbbinél gyorsabban konvergál

q_{GS} ≈ |1 − λ_{J}^{2} |  (2.95)

Túlrelaxálás (Successive Over-Relaxation, SOR): Írjuk át a GS egyenletet

[L + D] ·v^{(n+1)} = [A − U] ·v^{(n)} − [A · v^{(n)} − B]  (2.96)

majd ’rontsuk el’ a jobboldalt egy ω relaxációs paraméterrel az alábbi módon:

[A − U] ·v^{(n)} − ω [A · v^{(n)} − B]  (2.97)

Az így kapott

[L + D] ·v^{(n+1)} = − [U−(1 − ω)A] ·v^{(n)}+ωB  (2.98)

SOR formula konvergál, ha 0 ≤ ω ≤ 2. Ha ω < 1, akkor alul relaxálásról, ha ω > 1 akkor túlrelaxálásról beszélünk. A Gauss-Seidel iterációt kapjuk vissza, ha ω = 1. Az SOR formula végsősoron megfelel egy olyan GS eljárásnak, ahol az egyes lépésekben kapott vektorokat összekeverjük az előzővel

ωv_{GS}^{(n+1)} + (1 − ω)v^{(n)} ⟶ v_{SOR}^{(n+1)}  (2.99)

Az eljárás konvergenciáját most a

G_{SOR}= − [L + D]^{−1} [U − (1 − ω)A]  (2.100)

mátrix határozza meg. Az optimális ω

ω_{SOR} = \frac{2}{1 + \sqrt{1 − λ_{J}^{2}}}  (2.101)

amikor is a spektrálsugár

λ_{SOR} = (\frac{λ_{J}}{1 + \sqrt{1 − λ_{J}^{2}}} )^{2}  (2.102)

és az asszimptotikus konvergencia hányados

q_{SOR} ≈ |1 − λ_{SOR}|  (2.103)

Iteratív javítás: Bármelyik előző után jól jön.

G. Elliptikus egyenletek: peremérték probléma

A kezdetiérték feladatok tárgyalásánál a fő probléma azzal volt, hogy a kézenfekvő megoldási sémák stabilitását biztosítani kellett. A peremérték problémák megoldásánál a stabilitás viszonylag könnyen elérhető, a fő szempont most az lesz, hogy mennyire hatékony, gazdaságos algoritmust vagyunk képesek felállítani.

Alapfeladatként a kétdimenziós Laplace egyelet

\frac{∂^{2}V}{∂x^{2}} + \frac{∂^{2}V}{∂y^{2}} = ρ(x, y)  (2.104)

véges differenciákkal való megoldását vizsgáljuk az ismeretlen V (x, y) meghatározására adott ρ(x, y) mellett. Peremfeltételként előírjuk a vizsgált kétdimenziós tartományt határoló görbén az ismeretlen V (x, y) függvény vislekedését (pl. Dirichlet, vagy Neumann határfeltételeket).

A síktartományban felvett diszkrét pontokon

(x_{j} = x_{0} + jδx , j = 0, 1, ..., J) ; (y_{l} = y_{0} + lδy , l = 0, 1, ..., L)  (2.105)

V (x, y)-t a gridpontokon felvett

u_{j,l} = V (x_{j}, y_{l})  (2.106)

értékeivel adjuk meg. Az egyszerűség kedvéért most válasszuk a beosztást úgy, hogy legyen

δx = δy = Δ  (2.107)

A differencia egyenlet ekkor (a belső pontokban (j = 1, ..., J − 1) , (l = 1, ..., L − 1))

\frac{u_{j+1,l} − 2u_{j,l} + u_{j−1,l}}{Δ^{2}} + \frac{u_{j,l+1} − 2u_{j,l} + u_{j,l−1}}{Δ^{2}} = ρ_{j,l}  (2.108)

Használjunk most Dirichlet peremfeltételeket, azaz legyenek az {u_{j,0}, u_{0,l}, u_{j,L}, u_{J,l}|j = 0..J, l = 0..L} adottak.

1. Relaxáció

Az u_{j} = (u_{j,1}, u_{j,2}, ..., u_{j,L−1}) és az r_{j} = Δ^{2}(ρ_{j,1}, ρ_{j,2}, ..., ρ_{j,L−1}) − (u_{j,0}, 0, ..., 0, u_{j,L}) vektorok bevezetésével, az ismert mennyiségeket a jobboldalra rendezve az egyenletrendszer

T · u_{1} + I · u_{2} = r_{1} − u_{0}  (2.109)
I · u_{j−1} + T · u_{j} + I · u_{j+1} = r_{j} j = 2, ..., J − 2  (2.110)
I · u_{J−2} + T · u_{J−1} = r_{J−1} − u_{J}  (2.111)

alakba írható, ahol

T = ( \begin{matrix}−4 & 1 \\ 1 & −4 & 1 \\ & · & · & · \\ & & · & · & · \\ & & & 1 & −4 & 1 \\ & & & & 1 & −4 \end{matrix} ) és I = ( \begin{matrix}1 \\ & 1 \\ & & · \\ & & & · \\ & & & & 1 \\ & & & & & 1 \end{matrix} ) : (L − 1) × (L − 1)  (2.112)

A (J − 1) · (L − 1) dimenziós v⃗ = {u_{1} ⊕ u_{2}⊕, ..., ⊕u_{J−1}} és b⃗ = {(r_{1} − u_{0}) ⊕ r_{2}⊕, ..., ⊕ (r_{J−1} − u_{J})} szupervektorokkal a meghatározandó

v_{r} = u_{j,l} , r = (j − 1)L + l  (2.113)

mennyiségek az

A · v⃗ = b⃗ ahol A = ( \begin{matrix}T & I \\ I & T & I \\ & I & T & I \\ & & · & · & · \\ & & & I & T & I \\ & & & & I & T \end{matrix} )  (2.114)

lineáris egyenletrendszer megoldásával megkaphatók. Az ilyen pentadiagonális (következésképp ritka) mátrixokkal felírt lineáris egyenletek a már megismert (Jacobi, Gauss-Seidel és SQR) módszerekkel megoldhatók.

2. Váltakozó irányú implicit módszer: ADI

Az DV = ρ peremértékproblémával kapcsolatba hozható a

\frac{∂V}{∂t} = DV − ρ  (2.115)

diffúziós probléma. A diffúziós egyenlet t → ∞ megoldásai stacionáriusak, ezek éppen a peremértékfeladat keresett megoldásai. Esetünkben kézenfekvő az ADI eljárást felírni

\frac{v^{[n+1/2]} − v^{[n]}}{δt/2} = \frac{L_{x}v^{[n+1/2]} + L_{y}v^{[n]}}{Δ^{2}} − ρ⃗ és \frac{v^{[n+1]} − v^{[n+1/2]}}{δt/2} = \frac{L_{x}v^{[n+1/2]} + L_{y}v^{[n+1]}}{Δ^{2}} − ρ⃗

Átrendezve

[qI − L_{x}] v^{[n+1/2]} = [qI + L_{y}] v^{[n]} − Δ^{2}ρ⃗
[qI − L_{y}] v^{[n+1]} = [qI+L_{x}] v^{[n+1/2]} − Δ^{2}ρ⃗ , ahol q = 2\frac{Δ^{2}}{δt}

a bal oldalon álló [qI±L_{i}] mátrixok tridiagonálisak, így az egyenlet közvetlenül megoldható. Kiindulunk valamely v^{[0]}-ból, ezzel előállítjuk az első sor szerint v^{[1/2]}-et, majd a második sorral v^{[1]}-et. Ezt az első egyenletbe helyettesítjük és addig folytatjuk az eljárást, amíg |v^{[n+1]} − v^{[n]}| < ϵ nem lesz. Az eljárás konvergenciája a képzeletbeli időtengelyen felvett δt lépéshossztól és így q−tól függ. Ezen értékek dinamikus (iterációnként különböző) megválasztásával elérhető, hogy lényegesen gyorsabban jussunk célhoz, mint pl. az SOR módszerrel.

3. Ciklikus redukció

Válasszuk a beosztást úgy, hogy J = 2^{p} legyen. Tekintsünk három egymás utáni

I · u_{j−2} + T · u_{j−1} + I · u_{j} = r_{j−1}  (2.116)
I · u_{j−1} + T · u_{j} + I · u_{j+1} = r_{j}  (2.117)
I · u_{j} + T · u_{j+1} + I · u_{j+2} = r_{j+1}  (2.118)

egyenletet. A középsőt megszorozva −T-vel és összeadva a hármat

I · u_{j−2} + T^{(1)}·u_{j} + I · u_{j+2} = r_{j}^{(1)}  (2.119)

ahol

T^{(1)} = [2I − T^{2}] és r_{j}^{(1)} = r_{j−1} − T · r_{j} + r_{j+1}  (2.120)

Ezzel a centrális u_{j} meghatározható első szomszédai helyett csupán a második szomszédaiból, azaz a meghatározandó.u-k számát a felére csökkentettük (ebben a lépésben kiszórtunk minden páratlan indexű u_{j}-t). Az így kapott egyenlet szerkezetében azonos az eredetivel, így ezt újból és újból redukálhatjuk. A végén

⇓  (2.121)
u_{0} + T^{(p−1)}·u_{J/4} + u_{J/2} = r_{J/4}^{(p−1)}  (2.122)
u_{J/4} + T^{(p−1)}·u_{J/2} + u_{3J/4} = r_{J/2}^{(p−1)}  (2.123)
u_{J/2} + T^{(p−1)}·u_{3J/4} + u_{J} = r_{3J/4}^{(p−1)}  (2.124)
⇓  (2.125)
u_{0} + T^{(p)}·u_{J/2} + u_{J} = r_{J/2}^{(p)}  (2.126)

Az utolsó egyenletben T^{(p)} és r_{J/2}^{(p)} ismert (előállítható) az u_{0} és u_{J} pedig adott a peremfeltételből. Következésképpen u_{J/2} kiszámítható. Ebből az (p − 1) lépés egyenleteire visszalépve u_{J/4} és u_{3J/4} számolható, és így tovább számolva sorra előállíthatjuk a megoldást a többi pontban is.

4. Fourier módszer

A Fourier transzformáltakkal fontosságuk miatt külön fejezetben foglalkozunk részletesen. Megelőlegezve most néhány formulát a diszkrét Fourier transzformáció témaköréből, megvizsgáljuk, hogy a konstans együtthatós parciális differenciálegyenleteket hogyan oldhatjuk meg ezek segítségével. Példaként a Poisson egyenletet vesszük elő, ezt véges differenciákkal az (2.108) egyenletben

u_{j+1,l} + u_{j−1,l} + u_{j,l+1} + u_{j,l−1} − 4u_{j,l} = Δ^{2}ρ_{j,l}  (2.127)

alakra hoztuk. Az {u_{j,l}} és {ρ_{j,l}} mennyiségek kétdimenziós diszkrét Fourier transzformáltja:

ũ_{m,n} = ∑_{j=0}^{J−1} ∑_{l=0}^{L−1} u_{j,l} · e^{i2πjm/J}e^{i2πln/L} u_{j,l} = \frac{1}{JL} ∑_{m=0}^{J−1} ∑_{n=0}^{L−1} ũ_{m,n} · e^{−i2πjm/J}e^{−i2πln/L}  (2.128)

ρ̃_{m,n} = ∑_{j=0}^{J−1} ∑_{l=0}^{L−1} ρ_{j,l} · e^{i2πjm/J}e^{i2πln/L} ρ_{j,l} = \frac{1}{JL} ∑_{m=0}^{J−1} ∑_{n=0}^{L−1} ρ̃_{m,n} · e^{−i2πjm/J}e^{−i2πln/L}  (2.129)

Ezeket a differencia egyenletbe helyettesítve

ũ_{m,n} [e^{−i2πm/J} + e^{i2πm/J} + e^{−i2πn/L} + e^{i2πn/L} − 4] = Δ^{2}ρ̃_{m,n}  (2.130)

vagy

ũ_{m,n} = \frac{Δ^{2}}{2} \frac{ρ̃_{m,n}}{cos(2πm/J) + cos(2πn/L) − 2}  (2.131)

A probléma megoldása ezek után a következő lépésekben történik • meghatározzuk a ρ̃_{m,n} Fourier transzformáltakat

• az (2.131) egyenlet alapján kiszámítjuk az ũ_{m,n} mátrixot

• az inverz formulákkal u_{j,l} (2.128) alapján előáll

A most vázolt eljárás bármely periodikus határfeltétel mellett használható, azaz a megoldásunkra u_{j,l} = u_{j+J,l} = u_{j,l+L} és különösen a peremen u_{0,l} = u_{J,l} u_{j,0} = u_{j,L}. Ha Dirichlet peremfeltételünk van, azaz a megoldás a peremen eltűnik: u_{0,l} = u_{j,0} = u_{J,l} = u_{j,L} = 0, akkor a szinusz transzformáltakat célszerű felírni:

ũ_{m,n} = ∑_{j=1}^{J−1} ∑_{l=1}^{L−1} u_{j,l} · sin (\frac{πjm}{J}) sin (\frac{πln}{L}) u_{j,l} = \frac{4}{JL} ∑_{m=1}^{J−1} ∑_{n=1}^{L−1} ũ_{m,n} · sin (\frac{πjm}{J}) sin (\frac{πln}{L})  (2.132)

ρ̃_{m,n} = ∑_{j=1}^{J−1} ∑_{l=1}^{L−1} ρ_{j,l} · sin (\frac{πjm}{J}) sin (\frac{πln}{L}) ρ_{j,l} = \frac{4}{JL} ∑_{m=1}^{J−1} ∑_{n=1}^{L−1} ρ̃_{m,n} · sin (\frac{πjm}{J}) sin (\frac{πln}{L})  (2.133)

Ilyen felírással a peremfeltétel automatikusan kielégül. A transzformált egyenlet

ũ_{m,n} = \frac{Δ^{2}}{2} \frac{ρ̃_{m,n}}{cos(πm/J) + cos(πn/L) − 2}  (2.134)

aminek a megoldási folyamata az előzőével azonos. Neumann peremfeltételnél ∇u = 0 a peremen. Ez esetben a koszinusz transzformáltat használjuk

ũ_{m,n} = ∑_{j=0}^{J(!)} ∑_{l=0}^{L(!)} u_{j,l} · cos (\frac{πjm}{J}) cos (\frac{πln}{L}) u_{j,l} = \frac{4}{JL} ∑_{m=0}^{J(!)} ∑_{n=0}^{L(!)} ũ_{m,n} · cos (\frac{πjm}{J}) cos (\frac{πln}{L})  (2.135)

ρ̃_{m,n} = ∑_{j=0}^{J(!)} ∑_{l=0}^{L(!)} ρ_{j,l} · cos (\frac{πjm}{J}) cos (\frac{πln}{L}) ρ_{j,l} = \frac{4}{JL} ∑_{m=0}^{J(!)} ∑_{n=0}^{L(!)} ρ̃_{m,n} · cos (\frac{πjm}{J}) cos (\frac{πln}{L})  (2.136)

ahol a (!) jelzés arra figyelmeztet, hogy az összegzéseben, ha az index a határon van, akkor a megfelelő tagnak csak a felét kell számítani.

III. FOURIER TRANSZFORMÁCIÓ

A. A Fourier transzformációról általában

Az f(x) függvény Fourier transzformáltját a fizikában a

H(ω) = ∫_{−∞}^{∞} h(t)e^{−iωt}dt vagy H(ω) = \frac{1}{\sqrt{2π}} ∫_{−∞}^{∞} h(t)e^{−iωt}dt  (3.1)

módon szokás definiálni. Az inverz Fourier transzformáció szerint (amennyiben az alábbi integrálok léteznek) rendre

h(t) ≃ \frac{1}{2π} ∫_{−∞}^{∞} H(ω)e^{iωt}dω vagy h(t) ≃ \frac{1}{\sqrt{2π}} ∫_{−∞}^{∞} H(ω)e^{iωt}dω  (3.2)

A direkt és az inverz Fourier transzformációban az exponens előjele ellentétes, de néha szokás fordított sorrendben érteni. A most tárgyalni kívánt numerikus Fourier transzformációk összefüggésében az (1, 2π, \sqrt{2π}) faktorokban megjelenő diverzitást a

f = \frac{ω}{2π}  (3.3)

változócserével oldjuk fel, az exponens előjelét pedig úgy választjuk, hogy Fourier transzformált {h(x) ⟺ H(f)} párról beszélünk, ha

H(f) = ∫_{−∞}^{∞} h(t)e^{i2πft}dt és h(t) = ∫_{−∞}^{∞} H(f)e^{−i2πft}df  (3.4)

Emlékeztetőül néhány összefüggés a Fourier párokra

\begin{aligned}h(at) & ⟺ & \frac{1}{|a|}H(\frac{f}{a}) & \frac{1}{|a|}h(\frac{t}{a}) & ⟺ & H(af) \\ h(t − a) & ⟺ & H(f)e^{i2πfa} & h(t)e^{−i2πta} & ⟺ & H(f − a) \\ h^{(n)}(t) & ⟺ & (−i2πf)^{n} H(f) & (i2πt)^{n} h(t) & ⟺ & H^{(n)}(f) \end{aligned}  (3.5)

További hasznos fogalmak, tételek:

• Konvolúciós tétel

g ∗ h ≡ ∫_{−∞}^{∞} h(t − τ)g(τ)dτ = h ∗ g  g ∗ h ⟺ G · H  (3.6)

• Korreláció

h|g ≡ ∫_{−∞}^{∞} h(t + τ)g(τ)dτ  h|g ⟺ H · G^{∗} (^{∗} komplex konjugált)  (3.7)

• Autokorreláció

h|h ≡ ∫_{−∞}^{∞} h(t + τ)h(τ)dτ  h|h ⟺ |H|^{2}  (3.8)

A h(t) függvényhez tartozó P (f) ≡ |H(f)|^{2} más néven a spektrális teljesítménysűrűség.

B. Mintavételezés

Vizsgáljuk azt az esetet, amikor a transzformálandó függvényünk ekvidisztáns beosztáson az

h_{n} = h(nΔ) n = ..., −2, −1, 0, 1, 2, ....  (3.9)

értékekkel adott. Ehhez a mintához tartozik egy úgynevezett ’kritikus frekvencia’, a Nyquist frekvencia

f_{c} = \frac{1}{2Δ}  (3.10)

A kritikus frekvencia jelentősége megmutatkozik az alábbi tételben: Mintavételezési tétel: Ha a h(t) folytonos függvény sávhatárolt, azaz létezik olyan f_{M}, hogy

H(f) = 0 ha |f| > f_{M}  (3.11)

akkor h(t) teljesen megadható a

h_{n} = h(\frac{n}{2f_{M}})  (3.12)

mintavételezéssel és ekkor

h(t) = Δ ∑_{n=−∞}^{∞} h_{n} \frac{sin [2πf_{M} (nΔ − t)]}{π (nΔ − t)}  (3.13)

A tétel nagy jelentőségű sok fizikai/műszaki problémánál, amikor biztosak lehetünk abban, hogy a jel sávhatárolt. Ugyanakkor felhívja a figyelmet arra is, hogy a jelek nem megfelelő mintavételezése (f_{c} < f_{M}) alapján kiszámolt Fourier transzformáltak hamisak lesznek. Különösen biztosak lehetünk ebben, ha a függvény nem sávhatárolt, azaz f_{M} ∼ ∞. Sajnos ez is tipikus eset a gyakorlati problémákban. Ilyenkor a P (f) tejesítménysűrűség azon része, amire |f| > f_{M} hozzáadódik az |f| < f_{M} tartományon felvett értékekhez, a mintából készített Fourier transzformált hamis lesz.

Nem mindig tudjuk előre, hogy a függvényünk sávhatárolt-e, vagy mennyire jó közelítéssel az. Egy mintavételezési gyakoriságról akkor mondhatjuk el, hogy praktikusan kielégítő, ha az eredményül kapott Fourier transzformáltra teljesül, hogy az kellő mértékben eltűnik miközben a frekvencia felülről közelít valamely −f_{c} vagy alulról egy f_{c} értéket.

C. Diszkrét Fourier transzformáció (DFT)

Tegyük most fel, hogy h(t) tartója (vagy legalábbis az a tartomány, ahol a függvény viselkedése érdekes) valamely véges intervallumba esik. A vizsgált tartomány a függvény alkalmas eltolásával választható úgy, hogy a [0, T] intervallum legyen. Valamely

h_{k} = h(kΔ) , k = 0, 1, ..., N − 1  (3.14)

beosztást véve a Fourier transzformáltat az integrál közelítésével kiszámolhatjuk

H(f) = ∫_{0}^{T} h(t)e^{i2πtf}dt ≈ Δ ∑_{k=0}^{N−1} h_{k}e^{i2πfkΔ}  (3.15)

Valójában nincs értelme tetszőleges f frekvenciára ezt kiszámolni, hiszen az N darab h_{k} számból csak N darab független transzformált szám állítható elő. Legyen N páros szám és keressük csak a

H(f_{n}) = H(\frac{n}{NΔ}) ≈ Δ ∑_{k=0}^{N−1} h_{k}e^{i2πkn/N} , f_{n} = \frac{n}{NΔ} = \frac{n}{T} , n = −\frac{N}{2} , ...., +\frac{N}{2}  (3.16)

függvényértékeket. Az itt szereplő

H_{n} ≡ ∑_{k=0}^{N−1} h_{k}e^{i2πkn/N}  (3.17)

számokat az N darab h_{k} pont diszkrét Fourier transzformáltjának nevezzük. Vegyük észre, hogy H_{−N/2} = H_{N/2}, így pontosan N független értéket kapunk. Ezekkel a függvény Fourier transzformáltja

H(f_{n}) ≈ ΔH_{n}  (3.18)

A diszkrét Fourier transzformáltban a negatív indexek kényelmetlenné válhatnak. Látható a definícióból, hogy az indexhatárokat tetszőlegesre kiterjesztve a kifejezés periódikus az indexben, azaz H_{−n} = H_{N−n}, így kézenfekvő, hogy a negatív indexekkel jellemzett tartományt áthelyezzük egy index-periódussal feljebb. Válasszuk tehát az indexeléshez az n = 0, 1, ..., N − 1 számokat. Így n = 0 a nulla frekvenciához tartozó szám, a pozitív frekvenciás 0 < f < f_{c} transzformáltak 1 ≤ n ≤ N/2 − 1 indexűek, míg a negatív frekvenciákhoz tartozó értékek az N/2 + 1 ≤ n ≤ N − 1 helyen vannak. Az n = N/2 indexű transzformált f_{c} és −f_{c}-hez egyaránt hozzátartozik.

Az inverz diszkrét Fourier transzformáció

h_{k} ≡ \frac{1}{N} ∑_{n=0}^{N−1} H_{n}e^{−i2πkn/N}  (3.19)

hasonló módon számolható, mint a direkt DFT, a rutinban az exponens előjelét kell csak változtatni, a normálás pedig a rutinon kívül elvégezhető. A formula igazolására helyettesítsünk be

h_{k} ≡ \frac{1}{N} ∑_{n=0}^{N−1} ∑_{l=0}^{N−1} h_{l}e^{i2πln/N}e^{−i2πkn/N} = ∑_{l=0}^{N−1} h_{l} ∑_{n=0}^{N−1} (e^{i2π(l−k)/N})^{n}  (3.20)

és vegyük észre, hogy

∑_{n=0}^{N−1} (e^{i2π(l−k)/N})^{n} = N , ha l = k  (3.21)

egyébként pedig

∑_{n=0}^{N−1} (e^{i2π(l−k)/N})^{n} = \frac{1 − (e^{i2π(l−k)/N})^{N}}{1 − e^{i2π(l−k)/N}} = \frac{1 − 1}{1 − e^{i2π(l−k)/N}} = 0  (3.22)

D. Gyors Fourier Transzformáció (FFT)

A diszkrét Fourier transzformáltak

H_{n} = ∑_{k=0}^{N−1} h_{k}w_{n,k} ahol w_{n,k} = e^{i2πkn/N}  (3.23)

kiszámítása egy komplex mátrixszorzás munkigényének felel meg. NxN-es mátrixokra ez O(N^{2}) művelet. A gyors Fourier transzformáció (FFT: Fast Fourier Transform) speciális estekben ennél lényegesen kevesebb számolási munkával teszi lehetővé a transzformáció elvégzését. Válasszuk ugyanis a mintavételezésben használt pontok számát úgy, hogy az kettőnek valamely egész hatványa legyen

N = 2^{p}  (3.24)

Ha valamilyen okból nem áll módunkban így megválasztani a minták számát, akkor egészítsük ki a mintát fiktív nulla értékekkel amíg kettő valamely hatványa nem lesz a kibővített minta. Ugyanis olyan eljárást fogunk adni, ami ilyen esetben csak O(pN) számítási műveletet igényel. Az FFT algoritmus lényege, hogy az N pont diszkrét Fourier transzformációját visszavezetjük két fele akkora méretű, azaz N/2 méretű pontsereg transzformációjára. A dolog menetét a Danielson és Lánczos nevével jegyzett alábbi lemma mutatja.

∑_{k=0}^{N−1} h_{k}e^{i2πkn/N} = ∑_{k=0}^{N/2−1} h_{2k}e^{i2π(2k)n/N} + ∑_{k=0}^{N/2−1} h_{2k+1}e^{i2π(2k+1)n/N}  (3.25)

= ∑_{k=0}^{N/2−1} h_{2k}e^{i2πkn/(N/2)} + W^{n} ∑_{k=0}^{N/2−1} h_{2k+1}e^{i2πn/(N/2)}  (3.26)
H_{n} = F_{n}^{e} + W^{n}F_{n}^{o} , W = e^{i2π/N}  (3.27)

ahol az utóbbi egyenletben F_{n}^{e} a minta páros sorszámú, F_{n}^{o} pedig csak a páratlan sorszámú pontjaiból készült DFT. Figyeljünk arra, hogy az egyenlet jobb oldalán n = 0, 1, ..., N, míg a fele akkora méretű ponthalmazból készült F_{n}^{e} és F_{n}^{o} indexei n = 0, 1, ..., N/2 Az indexek szerinti periodikusság F_{n} = F_{n±N/2} miatt azonban egyszerűen kétszer egymás után kell leírnunk az F_{n} számokat.

Ha az eredeti pontjaink száma N = 2^{p}, akkor ezt a ’felezéses’ ejárást p-szer rekurzív módon használhatjuk, míg végül csak egyetlen pontból álló halmaz triviális diszkrét Fourier transzformációját kell végrehajtanunk.

Az FFT rutinok segítségével elvégezhetjük a diszkrét Fourier transzformációt akkor is, ha a pontjaink száma nem kettő hatványa. Ekkor azonban az előbbi redukció nem, vagy csak részben végezhető el és nem biztos, hogy az ”FFT” gyorsabb lesz, mint a közönséges direkt számolás. Néhány szerencsés eset:

• A Danielson-Lánczos teljes faktorizációhoz hasonlóan részleges faktorizációt alkalmazhatunk, ha a pontjaink számának van egész osztója. Jó esetben N néhány kis prímszám szorzataként előáll. Az iménti módon akkora méretű tömbök transzformációjára vezethetjük vissza a feladatot, mint N legnagyobb prím osztója.

• Néhány speciális tömbméret mellett lehetőség van arra, hogy nagyon effektív számítógépes rutint írjunk. Így például N = 2, 4, 8 esetén a transzformációban szereplő trigonometrikus függvények értéke speciálisan egyszerű és hatékonyan elvégezhető a transzformáció. Az ekkor a méretű tömbökre jól megírt résztranszformációkat használva nem kell a Lánczos redukciót végigvinni, ami 20-30% nyereséget hozhat a számításban.

• Hatékony transzformációs eljárások léteznek nemcsak az előbbi tömbméretekre, de pl. N = 3, 5, 7, 11, 16.... számú adathalmazra is. Ilyen például a Winograd algoritmus. Előfordulhat, hogy ezzel az eljárással kétszer gyorsabban tudjuk transzformálni az adatainkat, mint a legközelebb eső 2^{p}-re kiegészített adatsoron FFT-vel. Nagyméretű adatsorokra az algoritmus nagy hátránya, hogy extra ideiglenes tárolási igénye van, szemben a ”közönséges” FFT-vel, ahol a számolás nem igényel munkatömböket.

Részletesebben az FFT programozásával nem foglalkozunk, a megfelelő procedurák, szubrutinok általában készen rendelkezésre állnak a különböző programnyelvekhez. Az előzőekből láthatóan indokolt a tanács: Győződjünk meg arról, hogy milyen rutinok állnak rendelkezésre és készítőik melyik problémához melyik rutint ajánlják.

E. FFT valós függvényekre

Igen gyakori eset, hogy valós függvényt kell Fourier transzformálni. Nyilván pocséklás lenne, ha a DFT-t úgy végeznénk el, hogy a komplex bemenő függvényre megírt rutinnak átadnánk a valós részben a tényleges mintákat, míg a minta képzetes részét kitöltenénk nullákkal.

• Ha a probléma olyan, hogy több azonos méretű tisztán valós (vagy tisztán képzetes) adatsort kell transzformálni, akkor két mintát összefoglalva egyetlen komplex mintába egy lépésben két transzformációt is elvégezhetünk amivel teljes mértékben kihasználjuk az általános FFT rutinunkat. Legyenek ugyanis f(t) és g(t) valós függvények. A h(t) = f(t) + i · g(t) komplex függvény Fourier transzformáltjára a linearitás miatt igaz, hogy

H(f) = ∫_{−∞}^{∞} h(t)e^{i2πft}dt = ∫_{−∞}^{∞} f(t)e^{i2πft}dt + i ∫_{−∞}^{∞} g(t)e^{i2πft}dt = F (f) + iG(f)  (3.28)

Vegyük figyelembe, hogy valós függvényekre

(F (f))^{∗} = (∫_{−∞}^{∞} f(t)e^{i2πft}dt)^{∗} = ∫_{−∞}^{∞} f(t)e^{−i2πft}dt = F (−f)  (3.29)

így

(H (f))^{∗} = F (−f) − iG(−f) illetve (H (−f))^{∗} = F (f) − iG(f)  (3.30)

A H komplex Fourier transzformáltból tehát a keresett két transzformált egyszerűen visszafejthető:

F (f) = \frac{(H (−f))^{∗} + H(f)}{2} és G (f) = i\frac{(H (−f))^{∗} − H(f)}{2}  (3.31)

Ne felejtsük, hogy a DFT-ben szokásos módon a negatív frekvenciás együtthatókat a tömbben átmásoltuk és a pozitív frekvenciásak után tároljuk, azaz

H (−f_{n}) = H_{N−n} , f_{n} = \frac{1}{NΔ} \{0, 1, 2, ..., N/2\}  (3.32)

• Ha csak egyetlen adatsorunkat kell transzformálni, akkor az N elemű valós mintát alakítsuk át N/2 elemű komplex mintává a redundancia elkerülésére. Válasszuk két részre a mintát, célszerűen a páros és a páratlan sorszámú pontokra és legyen

h_{k} = f_{2k} + i · f_{2k+1} , k = 0, 1, ..., N/2 − 1  (3.33)

Az DFT-t elvégzzük erre az N/2 méretű tömbre. Eredeményül olyan koefficienseket kapunk amiben a páros és a páratlan pontok külön-külön vett transzformáltjai az alábbi módon vannak kombinálva:

H_{n} = F_{n}^{e} + iF_{n}^{o}  (3.34)

A Lánczos lemma kapcsán leírtak szerint

F_{n} = F_{n}^{e} + e^{i2πn/N}F_{n}^{o}  (3.35)

lenne az eredetileg keresett transzformált. Ezt az iménti H_{n}-ekből kifejezhetjük

F_{n} = \frac{1}{2} [H_{n} + H_{−n}^{∗}] − \frac{i}{2} [H_{n} − H_{−n}^{∗}] e^{i2πn/N} , n = 0, 1, ..., N − 1  (3.36)

ahol felhasznátuk, hogy a valós bemenő adatok miatt (F_{−n}^{x})^{∗} = F_{n}^{x} x = e/o. A formulában a negatív indexekre továbbra is érvényes a (3.32) tárolási konvenció.

F. sinFFT szinusz és cosFFT koszinusz Fourier transzformáltak

Differenciálegyenletekkel kapcsolatos problémákban gyakran teljesül, hogy a peremfeltételek szerint a megoldásnak el kell tűnnie a vizsgált tartomány határain (a peremen). A [0, T] intervallum határán eltűnő függvények Fourierszinusz sorba fejthetők

h(x) = ∑_{n=1}^{∞} b_{n} · sin (x\frac{nπ}{T} ) b_{n} = \frac{2}{T} ∫_{0}^{T} h(x) sin (x\frac{nπ}{T} ) dt  (3.37)

A függvényekről a mintavételezéssel kapott {h_{k} | k = 0, 1, ..., N − 1} értékekre áttérve a b_{n} Fourier együtthatók helyett használjuk a diszkrét szinusz Fourier transzformáltakat, definíció szerint a

B_{n} = ∑_{k=1}^{N−1} h_{k} sin (\frac{πnk}{N} ) ≈ \frac{T}{2Δ}b_{n} = \frac{N}{2} b_{n} n = (0), 1, 2, ..., N − 1, (N)  (3.38)

számokat . (Feltevésünk szerint h_{0}, B_{0}, h_{N}, B_{N} = 0 .) A diszkrét Fourier szinusz transzformáció faktortól eltekintve ”önmaga” inverze, pontosabban

h_{k} = \frac{2}{N} ∑_{n=1}^{N−1} B_{n} sin (\frac{πnk}{N} )  (3.39)

Ha nem áll rendelkezésre speciálisan az erre a célra használható gyors szinusz Fourier transzformáció (sinFFT) rutin, akkor az áltlános FFT eljárást használhatjuk a következő megfontolások alapján. Vezessük be a

h̃(x) = \{ \begin{matrix}h(x) & , x > 0 \\ −h(x) & , x < 0 \end{matrix} \}  (3.40)

páratlan segédfüggvényt. A [0, T] intervallum helyett a [−T, T] duplázott intervallumon

H̃(f) = ∫_{−T}^{T} h̃(t)e^{i2πtf}dt = 2i ∫_{0}^{T} h(t) sin(2πtf)dt ≈ Δ2i ∑_{k=1}^{N−1} h_{k} sin(2πkΔf)  (3.41)

Ha ezek után a 2T intervallum 2N darab {h̃_{k} | k = 0, 1, ..., 2N − 1} mintájához a szokásos

f_{n} = \frac{n}{2ΔN} = \frac{n}{2T}  (3.42)

frekvenciákat választjuk, akkor H̃(f_{n}) = ΔH̃_{n}, ahol {H̃_{n} | n = 0, 1, .., 2N − 1} a h̃ diszkrét Fourier transzformáltja. A megoldás tehát az, hogy az N számmal adott {h_{k} | k = 0, 1, ..., N − 1} mintát inverzióval megkétszerezzük úgy, hogy h̃_{k} = h_{k} ha k < N és h̃_{2N−k} = −h_{k}, ha k > N és ezen végrehajtunk egy diszkrét Fourier transzformációt. Utóbbi eredményéből pedig

H̃_{n} = 2iB_{n}  (3.43)

alapján a szinusz transzformáltak az FFT együtthatókból számothatók.

• Azt vártuk volna, hogy az N pont szinusz transzformációja fele akkora munkával végezhető el, mint ugyanannyi pont teljes komplex diszkrét transzformációja, ehelyett olyan eljárást adtunk, hogy az egy kétszer akkora méretű komplex transzformációval lett egyenértékű.

• A pazarlás mindaddig többnyire elfogadható, amíg egyszer, vagy csak néhányszor kell elkövetni.

• Súlyos problémává dagadthat a többletszámolás, ha sokszori transzformációban vagyunk kénytelenek lenyelni, mint például egy többdimenziós Fourier transzformáció számolásakor. Ebben az esetben - és, ha nem áll rendelkezésre megfelelő könyvtári rutin -, szükséges és lehetséges egyéb trükkökkel minimális szinten tartani a számolási munkát. A részletekkel kapcsolatban utalunk a szakirodalomra.

A koszinus transzformáltak szintén gyakran szükségesek. Tipikus eset, hogy a peremértékproblémában Neumann határfeltét adunk, azaz a megoldásról előírjuk, hogy annak deriváltja legyen nulla a határon. A peremfeltételeket kielégítő megoldások ilyenkor koszinusz sorba fejthetők:

h(x) = \frac{a_{0}}{2} + ∑_{n=1}^{∞} a_{n} · cos (x\frac{nπ}{T} ) a_{n} = \frac{2}{T} ∫_{0}^{T} h(x) cos (x\frac{nπ}{T} ) dt  (3.44)

A megfelelő diszkrét koszinusz transzformáltak

A_{n} = \frac{1}{2} \{h_{0} + (−1)^{n} h_{N}\} + ∑_{k=1}^{N−1} h_{k} cos (\frac{πnk}{N} )  (3.45)

Éppúgy, mint a szinusz transzformáltaknál megmutathajuk, hogy a párosan kiterjesztett

h̃(x) = \{ \begin{matrix}h(x) & , x > 0 \\ h(x) & , x < 0 \end{matrix} \}  (3.46)

függvényből vett {h̃_{k} | k = 0, 1, ..., 2N − 1} minta diszkrét {H̃_{n} | n = 0, 1, .., 2N − 1} Fourier transzformáltjából

H̃_{n} = 2A_{n}  (3.47)

A koszinusz transzformáció megint (faktortól eltekintve) önmaga inverze

h_{k} = \frac{2}{N} [\frac{1}{2} \{A_{0} + (−1)^{n} A_{N}\} + ∑_{n=1}^{N−1} A_{n} cos (\frac{πnk}{N} )]  (3.48)

G. Többdimenziós FFT

Elsőként a fizikában igen gyakori 3-dimenziós Fourier transzformálttal foglalkozunk. Egy f(r): R^{3} → C függvény Fourier transzformáltja a konvencionális definíciókkal

F (q) = ∫ f(r)e^{−iqr}d^{3}r ⟺ f(r) = \frac{1}{(2π)^{3}} ∫ F (q)e^{iqr}d^{3}q  (3.49)

Ha a transzformálandó függvénynek valamilyen szimmetriája van, azt a transzformáció során célszerű kihasználni. A legfontosabb eset, ha a függvény gömbszimmetrikus, azaz f(r) ≡ f(r), r = |r|. Ekkor ugyanis

∫ f(r)e^{−iqr}d^{3}r =2π ∫_{0}^{∞} r^{2}dr ∫_{−1}^{1} du f(r)e^{−iqru} = 2π ∫_{0}^{∞} r^{2}dr f(r)\frac{2 sin(qr)}{qr} = \frac{4π}{q} ∫_{0}^{∞} r · f(r) · sin(qr)dr (3.50)

és így a 3-dimenziós Fourier transzformált helyett a g(r) = r · f(r) függvény 1-dimenziós szinusz Fourier transzformáltját elég kiszámolni. Ha hasonló elvi megfontolásokkal nem sikerül a feladatot alacsonyabb dimenzióra redukálni, akkor kénytelenek leszünk ’gyalog’ megoldással élni. Tekintsük például a kétdimenziós

H_{n,m} = ∑_{k=0}^{Nk−1} ∑_{l=0}^{Nl−1} h_{k,l} e^{i2πkn}e^{i2πlm}  (3.51)

diszkrét transzformáltat, amit a következő módon tudunk kiszámolni. Előbb az egyik dimenzióban, esetünkben minden rögzített k index mellett Nk-alkalommal elvégzett DFT-vel kiszámoljuk az összes ideiglenes

(k)H_{m} = ∑_{l=0}^{Nl−1} h_{k,l} e^{i2πlm} k = 0, 1, .., Nk − 1  (3.52)

diszkrét transzformáltat. Az így kapott tömböket átrendezzük

(k)H_{m} →^{(m)} H_{k}  (3.53)

alakra, majd Nl alkalommal használjuk ezekre újból az egydimenziós

H_{n,m} = ∑_{k=0}^{Nk−1} ^{(m)}H_{k} e^{i2πkn} l = 0, 1, .., Nl − 1  (3.54)

transzformációt. Összesen tehát Nk × Nl alkalommal kellett az FFT rutint meghívni, jó lesz ügyelni arra, hogy az FFT optimálisan legyen kihasználva. Némi programozási gyakorlattal azt is láthatjuk, hogy a két egydimenziós transzformáció közötti átrendezés (oszlop-sor csere) nagy mennyiségű komplex számra igen forrásigényes lehet.

Célszerű tehát meggyőződni, hogy nincs-e a használt programnyelvhez olyan kész rutin, amelyik a többdimenziós transzformáltakat hatékonyabban számolja.

Magasabb dimenziós Fourier transzformáltak számításakor különösen fontos, hogy minden speciális körülményt figyelembe vegyünk és ne számoljunk redundáns dolgokat. Ilyen gyakori eset az, amikor valós függvényeket kell transzformálni. A fizikában a Poisson egyenletre vezető problémák többsége ilyen valós függvényekre felírt differenciálegyenlet. Az igencsak nagyon fontos 2 vagy 3 dimenziós képfeldolgozási problémák is valós számokkal adott információk Fourier analízisét igénylik. Nagyon nagy számú transzformálandó valós adat esetén ne használjuk az előbbi eljárást, hanem keressünk erre a célra írt (remélhetőleg) effektív rutint. Szükség esetén írjunk magunk, segítséget ehhez a szakirodalomban találunk.

H. Konvolúciós egyenletek

A g(t) és h(t) függvények konvolúcióját a

g ∗ h ≡ ∫_{−∞}^{∞} h(t − τ)g(τ)dτ = h ∗ g  (3.55)

formulával definiáltuk. A konvolúció kommutatív, a két függvény matematikailag azonos szerepet játszik. A gyakorlati problémákban azonban eredetét és természetét illetően lényeges különbség szokott lenni a konvolúciós partnerek között. Míg az egyik függvény valamilyen jel folyam, többnyire végtelen adathalmaz, addig a másik egy többnyire véges tartójú válasz függvény. Ebben az összefüggésben a válaszfüggvénynek szemléletes jelentést adhatunk. Ha a bemenő jel egyetlen tű-impulzus, amit egy Dirac-delta általánosított függvénnyel írhatunk le, akkor a konvolúció eredményeként

g ∗ δ = g  (3.56)

alapján magát a válaszfüggvényt kapjuk. Másképp mondva: a válaszfüggvény az a mért hatás, amit a készülékbe bemenő rövid impulzus generál.

A két különböző jelentésű függvényt célszerű jelölésben is megkülönböztetni, így a továbbiakban s(t) jelről (signal) és r(t) válaszfüggvényről (response) értekezünk. Tegyük fel, hogy a válaszfüggvény véges intervallumon különbözik nullától, azaz valamely {r_{k} | − M/2 < k ≤ M/2} mintával megadható. A szokásos elnevezés szerint

(s ∗ r)_{j} = ∑_{k=−M/2+1}^{M/2} s_{j−k} · r_{k}  (3.57)

számsor s és az r diszkrét konvolúciója. Az integrállal adott definícióból láthatjuk, hogy a numerikus diszkrét konvolúcióra

(s ∗ r) (t_{j}) ≈ Δ · (s ∗ r)_{j}  (3.58)

A konvolúciós tétel s ∗ r ⟺ S · R diszkrét változata a következő: Ha a jel N-periódikus akkor

(s ∗ r)_{j} = ∑_{k=−N/2+1}^{N/2} s_{j−k} · r_{k} ⟺ S_{n} · R_{n}  (3.59)

ahol {s_{k} | 0 ≤ k ≤ N − 1} ⟺ {S_{k} | 0 ≤ k ≤ N − 1} és {r_{k} | 0 ≤ k ≤ N − 1} ⟺ {R_{k} | 0 ≤ k ≤ N − 1}. Értelemszerűen az összegzésben előforduló negatív indexek kezeléséhez használjuk ki, hogy mindegyik számhalmaz indexében periódikus, azaz x_{−n} = x_{N−n} , x = s, r, S, R.

A tétel kellemetlennek tűnő megszorítása az, hogy mindkét mintáról feltételezi, hogy azok periódikusak, ráadásul ugyanazzal a periódussal.

• Ha a vizsgált jel történetesen periódikus és csak az a probléma, hogy a válaszfüggvény tartója rövidebb, mint a jel periódusa, akkor a válasz függvényből vett mintát nullákkal kiegészítjük, azaz

r_{k} = 0 ha M/2 ≤ k ≤ N/2 vagy − N/2 + 1 ≤ k ≤ −M/2 + 1  (3.60)

• Óvatosan kell eljárnunk, ha a jel a valóságban nem periódikus. A konvolúciós eljárás a vizsgált jel-intervallumba ’behozza’ az intervallumot megelőző és az azt követő tartományból a jel ottani feltételezett értékét. A diszkrét konvolúcióban implicite a jel periódikusan megismételt, így a jelet nullákkal ki kell egészítenünk annyira, hogy a kibővített és periódikus mintában a válaszfüggvény ’hatósugaránál’ jobban szeparáltak legyenek a tényleges jel adatai. Ha a véges tartójú válasz olyan, hogy r_{k} = 0, ha |k| > K, akkor a jel mintájának elejére vagy a végére K darab nullát be kell írnunk.

Két függvény konvolúcióját az előbbiek alapján viszonylag veszélytelenül kiszámolhatjuk. A fordított feladat, a dekonvolúció. Ilyenkor a konvolúciót ismerjük, és szeretnénk az egyik partner ismeretében a másikat kiszámolni, például f ∗ g = h konvolúciós egyenletben h és g ismeretében meghatározni, hogy mi lehetett az f. A megoldás formálisan egyszerűnek látszik, hiszen

F[h] = F[f] · F[g] ⟹ F[f] = F[h]/F[g] ⟹ \{ \begin{matrix}f = F^{−1} [F[h]/F[g]] \\ f = h ∗ F^{−1} [1/F[g]] \end{matrix}  (3.61)

A probléma ezekkel azonban az, hogy a jelzett műveletek általában nem végezhetők el, vagy megbízhatatlan eredményt adnak. Gondoljunk arra, hogy az F[g]-vel való osztás csaknem biztosan problémához vezet amikor a Fourier transzformált kicsi, vagy egyenesen nulla. De, még ha ettől el is tekintünk, a további lépésekben jelzett inverz Fourier transzformáltak létezése, numerikus értékük értelmezése további megfontolásokat igényel.

Két függvény korrelációját korábban a

h|g ≡ ∫_{−∞}^{∞} h(t + τ)g(τ)dτ  (3.62)

módon definiáltuk Az N-periódikus jelekre a diszkrét korreláció és a diszkrét Fourier transzformáció képletei:

(h|g)_{j} = ∑_{k=0}^{N−1} h_{j+k} · g_{k} ⟺ H_{n} · G_{n}^{∗}  (3.63)

A diszkrét korreláció numerikus számítására a konvolúció kapcsán elmondottak egyszerűen átvihetők.