Az 1-dimenziós Schrödinger-egyenlet numerikus megoldása

Bartha Ferenc

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

Szegedi Tudományegyetem, Elméleti Fizikai Tanszék (April 16, 2002)

I. FELADAT:

Az időfüggő Schrödinger-egyenlet alkalmas (kényelmes) egységekben felírva

i\frac{∂ψ(x, t)}{∂t} = Hψ(x, t) , H = −\frac{∂^{2}}{∂x^{2}} + V (x)
kezdőfeltétel : ψ(x, t = t_{0}) = ψ_{0}(x)
peremfeltétel : ψ(±1, t) = 0

A numerikus tárgyaláshoz szükséges, hogy a peremfeltétel a végesben legyen megadva. A kitűzött feladatban (az általánosság megsértése nélkül) a koordináta-skála választásával rögzítettük, hogy a perem x = −1 és x = 1 legyen.

• Keressük meg (például a Numerov-módszer segítségével) az első néhány kötött állapotot. Néhány tetszés szerinti potenciált kiválasztva oldjuk meg a kétpontos peremértek-problémát

Hφ_{μ}(x) = E_{μ}φ_{μ}(x) , μ = 0, 1, 2, ..., NodesMax
peremfeltétel : φ_{μ}(x = ±1) = 0

Kezdjük a ”legegyszerűbb” esettel, a V (x) ≡ 0 potenciállal, hiszen erre az esetre az analitikus megoldást is ismerjük (∼részecske dobozban).

• Az így meghatározott stacionárius állapotokból építsük fel a

Ψ(x, t) = ∑_{n=0}^{NodesMax} φ_{μ}(x)e^{−iE_{μ}t}

hullámfüggvényt. Ez megoldása az időfüggő Schrödinger-egyenletnek, a numerikus módszerünket ezzel fogjuk tesztelni.

• A Crank-Nicholson sémát használva írjunk programot az időfüggő Schrödinger-egyenlet, mint kezdetiértékprobléma megoldására. Próbáljuk ki a

ψ(x, t = 0) = Ψ(x, 0) = ∑_{n=0}^{NodesMax} φ_{μ}(x)

kezdőfeltétellel a programunkat és hasonlítsuk össze a kapott numerikus megoldást az ”egzakt” Ψ(x, t)-vel.

II. MEGJEGYZÉSEK A MELLÉKELT MEGOLDÁSHOZ

A. A Numerov rekurzió a sajátérték problémához

A megoldandó differenciálegyenlet

φ′′(x) = f(x)φ(x) , f(x) = V (x) − E

Taylor-sorból látható, hogy

φ_{n+1} − 2φ_{n} + φ_{n−1} = [φ′′]_{n} h^{2} + \frac{h^{4}}{12} [φ′′′′]_{n} + O (h^{6})
= [fφ]_{n} h^{2} + \frac{h^{4}}{12} [\frac{(fφ)_{n+1} − 2 (fφ)_{n} + (fφ)_{n−1}}{h^{2}} ] + O (h^{6})

vagyis

[1 − \frac{h^{2}}{12} f_{n+1}] φ_{n+1} − [2 + \frac{10h^{2}}{12} f_{n}] φ_{n} + [1 − \frac{h^{2}}{12} f_{n−1}] φ_{n−1} ≈ 0

Tömöreb írással ez

[1 − T_{n+1}] φ_{n+1} − [2 + 10T_{n}] φ_{n} + [1 − T_{n−1}] φ_{n−1} ≈ 0 ahol T (x) = \frac{h^{2}}{12} f(x)

majd kissé átrendezve

\frac{[1 − T_{n+1}] φ_{n+1}}{[1 − T_{n}] φ_{n}} + \frac{[1 − T_{n−1}] φ_{n−1}}{[1 − T_{n}] φ_{n}} = \frac{[2 + 10T_{n}]}{[1 − T_{n}]}

Az egyenlet ebben az alakjában igen alkalmas a probléma tárgyalására, hiszen a következő rekurziókra vezet

előre: R_{n+1}^{+} = U_{n} − 1/R_{n}^{+} vagy hátra: R_{n−1}^{−} = U_{n} − 1/R_{n}^{−}

ahol

R_{n+1}^{+} = \frac{[1 − T_{n+1}] φ_{n+1}}{[1 − T_{n}] φ_{n}} , R_{n−1}^{−} = \frac{1}{R_{n}^{+}} = \frac{[1 − T_{n−1}] φ_{n−1}}{[1 − T_{n}] φ_{n}} , U = \frac{[2 + 10T]}{[1 − T]}

Ha a beosztást úgy választottuk, hogy x_{1} = −1 és x_{N} = 1, akkor

φ_{1} = 0 → R_{2}^{+} = ∞ és φ_{N} = 0 → R_{N−1}^{−} = ∞

Balról indulva az előre rekurziót tudjuk hajtani ([R_{2}^{+}] → R_{3}^{+} → R_{4}^{+} → ...) és a jobb oldali peremtől pedig a hátra rekurziót ([R_{N−1}^{+}] → R_{N−2}^{+} → R_{N−3}^{+} → ...) használhatjuk. A programozott eljárásban balról haladunk mindaddig, míg R_{n}^{+} ≤ 1. Az R_{M}^{+} ≈ 1 osztópontnál megállunk az előre rekurzióval és a másik szélről belépegetünk ehhez a találkozási ponthoz. A találkozásnál kiszámolt

|R_{M⇀}^{+} − 1/R_{M↼}^{−}| = Δ (E)

mennyiség alkalmas annak az eldöntésére, hogy az aktuális energia sajátenergia-e, hiszen sajátállapotban

Δ (E_{μ}) ≈ 0

• Az alapállapot (és így minden sajátállapot) energiájára alsó korlátot adhatunk, hiszen (miért?)

E_{0} ≥ V_{min} \frac{π^{2}}{4} ahol V_{min} = min \{V (x) | x ∈ [−1, 1]\}

• Az 1-dimenziós Schrödinger-egyenlet kötött állapotairól tudjuk, hogy az energiájuk szerint növekvő sorba rendezett {φ_{μ} (x) , μ = 0, 1, 2, ...} függvényekre igaz, hogy φ_{μ} (x)-nek pontosan μ darab belső zéróhelye van. A Numerov integrálás során megszámoljuk, hogy az aktuális E érték mellett R hány osztópontnál negatív, ez adja a zéróhelyek számát. Tekintettel arra, hogy ”normális” körülmények között

T (x) = \frac{h^{2}}{12} [V (x) − E] < 1

R_{n−1}^{−} negatív értéke valóban azzal kapcsolatos, hogy φ_{n−1} és φ_{n} ellenkező előjelű.

• Az m = 0, 1, 2, ..., NodesMax adott számú zéróhellyel rendelkező megoldás keresésekor az energiát növeljük, ha a próba során kevesebb zéróhelyet és csökkentjük, ha többet találtunk. Ha a zéróhelyek száma megfelelő (és bízunk abban, hogy közel járunk az igazi sajátértékhez) lineáris extrapolációval keressük a Δ (E) = 0 egyenlet gyökét. Ez az eljárás alaposan eltévedhet, de nem volt most cél, hogy univerzális (idiótabiztos) eljárást adjunk.

B. Crank-Nicholson eljárás az időfejlődéshez

A beprogramozott eljárás a CN formulát használja:

−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}V (x_{j})

vagy mátrix alakban

A · ψ^{n+1} = d ahol d_{j} = aψ_{j+1}^{n} + (1 − 2a − b_{j}) ψ_{j}^{n} + aψ_{j−1}^{n}

és A tridiagonális mátrix:

A = ( \begin{matrix}1 & 0 & · & · & · & · & · \\ AL_{2} & AD_{2} & AU_{2} & · & · & · & · \\ · & AL_{3} & AD_{3} & AU_{3} & · & · & · \\ · & · & ... & ... & ... & · & · \\ · & · & · & ... & ... & ... & · \\ · & · & · & · & AL_{N−1} & AD_{N−1} & AU_{N−1} \\ · & · & · & · & · & 0 & 1 \end{matrix} ) és \{ \begin{matrix}AL_{j} = AU_{j} = −a \\ AD_{j} = 1 + 2a + b_{j} \end{matrix} j = 2, ..., N − 1\}