Stabilitás

Az elõzõ fejezetben felírt (6) diszkretizált egyenlet nem mûködik megfelelõen mint az a 2. feladat megoldásából kiderült. Ennek az az oka, hogy a kontinuitási egyenletre az FTCS módszer numerikusan instabil. Ezt könnyen láthatjuk, ha elvégezzük a Neumann féle stabilitási vizsgálatot. Bontsuk fel a kezdeti feltételt Fourier komponensekre, és vizsgáljuk meg, hogy az egyes komponensek amplitúdója hogyan változik az FTCS módszer alkalmazása során. Mivel a (6) egyenlet lineáris $\rho$-ban, valamint a határfeltételeket periodikusnak választottuk, a stabilitás analízis egzakt eredményt fog adni. A $k$ hullámszámú Fourier komponens a

$$\rho^n_j = G_n \exp(i k x_j)~~~~~~~~~~j=1,\ldots,N$$ (8)

alakban írható fel, ahol $G_n$ az amplitúdó, és $i=\sqrt{-1}$ az imaginárius egység. Ha behelyettesítjük a (8) formulát az (6) egyenletbe, akkor a $G=G_{n+1}/G_n$ amplitúdó növekedési faktorra a következõ megoldást kapjuk:

$$\vert G\vert=\left\vert 1-v\frac{\Delta t}{2\Delta x} \lef... ...vert 1-v\frac{\Delta t}{\Delta x}i\sin(k\Delta x)\right\vert>1$$ (9)

Mint látható tetszõlegesen kis $\Delta t$ idõlépés mellett is $\vert G\vert>1$, azaz a Fourier komponens amplitúdója minden lépésben nõ. Tehát az FTCS módszer a kontinuitási egyenletre nézve instabil.

Vizsgáljuk meg, hogyan javítható a diszkretizált egyenlet viselkedése. Az egyik megoldás az, ha a térbeli deriváltat féloldalasra változtatjuk:

$$\rho_j^{n+1}=\rho_j^n-v \Delta t \frac{\rho^n_{j}-\rho^n_{j-1}} {\Delta x}$$ (10)

Természetesen nem mindegy, hogy a bal vagy a jobb oldali differenciát használjuk. Mivel itt $v>0$, azaz az áramlás balról jobbra halad, logikus, hogy a bal oldali $\rho^n_{j}-\rho^n_{j-1}$ differenciát használjuk. Ezt áramlásfelöli (angolul ,,upwind'') deriváltnak nevezzük.

Végezzük el a stabilitás vizsgálatot a (10) módszerre vonatkozóan. A (8) Fourier komponsnek az upwind módszer (10) képletébe való helyettesítésével az erõsítési faktorra a

$$G=1-v\frac{\Delta t}{\Delta x} \left(1-e^{-ik\Delta x}\right) = 1-C + C e^{-ik\Delta x}$$ (11)

megoldást kapjuk, ahol

$$C=v\frac{\Delta t}{\Delta x}$$ (12)

a Courant szám, ami azt adja meg, hogy $\Delta t$ idõ alatt az információ hány rácstávolságnyira terjed. A (11) egyenlet szerint az erõsítési faktor a komplex számsíkon egy $C$ sugarú körön helyezkedik el, melynek középpontja a valós tengelyen $1-C$ -nél van amint ez a 2. ábrán látható.

Ábra 2.: A G komplex erõsítési faktor az upwind módszerre. G a folytonos vonallal rajzolt körön helyezkedik el. A stabil tartományt a szaggatott kör jelöli.

A stabilitás feltétele

$$0 \le C \le 1$$ (13)

azaz a Courant szám nem lehet 1-nél nagyobb. A $C=1$, azaz $\Delta x=v\Delta t$ választás rendkívül érdekes eredményt ad, ugyanis ekkor $G=\exp(ikx)$, azaz a Fourier komponens értéke

$$\rho^n_j = (G)^n \rho^0_j = e^{i k x_j -i k n \Delta x} = e^{i k (x_j - v t_n)}$$ (14)

ami éppen az egzakt megoldása a kontinuitási egyenletnek. Mivel minden Fourier komponensre egzakt a diszkrét megoldás, általában is igaz, hogy az upwind módszer a konstans sebességû kontinuitási egyenletre $C=1$ Courant szám mellett egzakt megoldást ad. Ez másképpen is belátható: $C=1$ esetén az upwind módszer minden idõlépésben éppen egy rácstávolságnyival viszi arrébb a megoldást. Valójában persze ez nem túl hasznos, hiszen a konstans sebességû kontinuitási egyenlet megoldásához nincs szükség számítógépes modellre. Reális, általánosabb egyenletre is jellemzõ eredményt a $C<1$ Courant számmal végzett szimulációk mutatnak.

Az is látható, hogy ha $v$ elõjele negatív lenne, akkor $\Delta t$ -tõl függetlenül instabillá válik a módszer. Éppen ezért az általános upwind módszer $v$ elõjelének függvényében választja a bal vagy jobboldali differencia formulát:

$$\rho_j^{n+1}=\left\{ \begin{array}{ll} \rho_j^n-C (\rho^n_... ...ho^n_{j+1}-\rho^n_{j}) & \mbox{ha $v\le0$} \end{array}\right.$$ (15)

Az így definiált upwind módszer feltételesen stabil $\vert C\vert<1$-re, most már $v$ elõjelétõl függetlenül.