Numerische Lösung der Schrödinger-Gleichung - Eigenwerte

Das ist meine erste Frage hier. Ich versuche, die Schrödinger-Gleichung für das Woods-Saxon-Potential numerisch zu lösen und die Energieeigenwerte und Eigenfunktionen zu finden, aber ich bin verwirrt darüber, wie genau das gemacht werden soll.

Ich habe in der Vergangenheit einige Anfangswertprobleme mit iterativen Methoden wie Runge–Kutta gelöst. Ich habe gelesen, dass Numerovs Methode der Weg ist, um Schrödingers Gleichung zu lösen, aber Wikipedia beschreibt sie auch als iterative Methode für Anfangswertprobleme.

Wie verwende ich es, um ein Eigenwertproblem zu lösen?

Das verwirrt mich aus folgenden Gründen:

  • Würde das iterative Lösen des DE nicht die Kenntnis der Energieeigenwerte erfordern, die als Eingabe für die Berechnung verwendet werden sollen? Ich kenne die Eigenwerte noch nicht; sie sind genau das, was ich zu berechnen versuche.
  • Wenn ich das täte, würde ich dann nicht einfach eine eindeutige Lösung anstelle einer Familie von Eigenfunktionen und Eigenwerten erhalten?

Ich habe einige Erwähnungen von "tridiagonalen Matrizen" gesehen, die irgendwie generiert wurden, bin mir aber nicht sicher, was die Elemente dieser Matrix sein würden oder wie das auf das Problem zutrifft. Leandro M. erwähnte, dass "die Diskretisierung ein endlichdimensionales (Matrix-) Eigenwertproblem definiert". Dies scheint der richtige Weg zu sein, den ich gehen sollte, aber ich konnte nichts finden, was diesen Prozess oder den Aufbau der Matrix explizit erklärt. Wenn dies das richtige Verfahren ist, wie wird eine solche Matrix aufgebaut?

Wäre Computational Science ein besseres Zuhause für diese Frage?
@DavidZ Ist die Formulierung der Frage jetzt zufriedenstellend?
Viel besser, ja. (und mir ist aufgefallen, dass die Frage bereits einige weitere Wiedereröffnungsstimmen gesammelt hat, seit Sie sie bearbeitet haben, also würde sie wahrscheinlich auch ohne mich wiedereröffnet werden, so wie das System funktionieren soll!) Übrigens, tut mir leid, ich habe es nicht bekommen früher anschauen; Letzte Nacht war ich viel beschäftigter, als ich dachte.
Wenn Sie Anfangswertprobleme lösen können, können Sie dann das Randwertproblem mit der Schießmethode lösen .
Ich konnte nichts finden, was diesen Prozess oder den Aufbau der Matrix explizit erklärt. Dann haben Sie sich nicht die Mühe gemacht, ein Buch über numerische PDEs aufzuschlagen.
Warum nicht stattdessen die Variationsmethode verwenden?

Antworten (2)

Ich bin froh, dass ich das endlich beantworten kann.

Numerovs Methode, wie sie auf Wikipedia beschrieben ist, ist nicht so, wie Sie hier vorgehen möchten. Um Ihnen eine Vorstellung davon zu geben, wie Sie vorgehen müssen, beginnen wir mit einer vereinfachten Version der Methode. Was ich tun werde, ist, den Differentialoperator einfach naiv zu diskretisieren, wie folgt:

d 2 d x 2 ψ ≈ ψ ( x − d ) − 2 ψ ( x ) + ψ ( x + d ) d 2 ≡ ψ n − 1 − 2 ψ n + ψ n + 1 d 2

Die letzte Gleichung ist nur eine Definition – ich behandle den Raum, als wäre er ein Gitter mit Punkten d auseinander und ich rufe ψ n der Wert der Wellenfunktion auf der n Punkt.. Nun lautet die Schrödinger-Gleichung in dieser Schreibweise:

− ℏ 2 2 m ψ n − 1 − 2 ψ n + ψ n + 1 d 2 + v n ψ n = E ψ n

Aber das ist eine Matrixgleichung ! Lassen Sie mich deutlicher werden:

− ℏ 2 2 m d 2 ( − 2 1 1 − 2 1 ⋱ ⋱ ⋱ 1 − 2 1 1 − 2 ) ( ψ 1 ψ 2 ⋮ ψ N − 1 ψ N ) + ( v 1 v 2 ⋱ v N − 1 v N ) ( ψ 1 ψ 2 ⋮ ψ N − 1 ψ N ) = E ( ψ 1 ψ 2 ⋮ ψ N − 1 ψ N )

Natürlich musste ich eine Ganzzahl auswählen N das entspricht einfach dem Platzieren des Systems in einer Box, die "groß genug" ist, um ein endliches System zu haben.

Es ist jetzt klar, dass Sie ein Matrix-Eigenwertproblem der Form haben

EIN ψ = E ψ

und Sie können fortfahren, es auf beliebige Weise zu diagonalisieren. Beachten Sie, dass wir anrufen EIN aus offensichtlichen Gründen eine tridiagonale Matrix. Vielleicht möchten Sie sich zuerst um die Randbedingungen kümmern – das tun Sie, indem Sie festlegen ψ 1 = ψ n = 0 bevor Sie die Matrix konstruieren, was dem Setzen der ersten und letzten Spalte auf Null entspricht. Dadurch erhalten Sie einige falsche Eigenfunktionen mit einem Eigenwert von Null, die Sie einfach verwerfen können. Wenn Sie unterschiedliche Randbedingungen haben, haben Sie Pech - ich kenne keinen Weg, der es allgemein zum Laufen bringt.

Jetzt müssen Sie das Obige nur noch mit der vollständigen Numerov-Methode wiederholen, was etwas komplizierter sein wird, und schon sind Sie fertig.

Dies scheint der Weg zu sein, nach dem Sie suchen, aber natürlich ist dies nicht der einzige Weg, dies zu tun. Griffiths beschreibt einen namens „Wag den Hund“, der aus Folgendem besteht: Sie schätzen einen Anfangswert für eine Energie und berechnen die Wellenfunktion, wie Sie wollen (zum Beispiel RK). Es besteht die Möglichkeit, dass es kein Eigenwert ist, sodass die Wellenfunktion im Unendlichen explodiert. Es könnte gehen + ∞ für groß x oder es könnte gehen − ∞ . Sie variieren nun langsam die Energie, bis das Verhalten im Unendlichen „kippt“, also bis der Schwanz „wedelt“. Auf diese Weise können Sie den Wert eines einzelnen Energieeigenwerts einschränken und die Form der Wellenfunktion auf eine maximale Größe ausgeben, die mit der gewählten Größe zunimmt E nähert sich dem exakten Eigenwert.

Du solltest besser einstellen ψ 0 = ψ N + 1 = 0 anstatt ψ 1 = ψ N = 0 (dh Punkt setzen mit r = 0 bei n = 0 ). Auf diese Weise schließen Sie gezielt die Punkte aus, die eigentlich keinen Nutzen haben, und tatsächlich repräsentiert Ihre Matrix bereits eine solche Auswahl. Dann gibt es keine falschen Lösungen, und als Bonus könnten Sie auch Coulomb-Potenzialprobleme lösen (weil die Singularität von den Punkten ausgeschlossen ist, nach denen Sie lösen). Beachten Sie auch, dass, da das Problem ein 3D-Problem ist und Sie nach dem radialen Teil der Wellenfunktion lösen, die Lösung, die Sie finden, dividiert werden sollte r .

Sie können mit der Bisektionsmethode nach Eigenwerten suchen.

Priliminaries: Um die Eigenwerte aus der Numerov-Methode zu erhalten, müssen Sie die Wellenfunktion an den Grenzen kennen. Im Allgemeinen würde dies bedeuten, dass Sie das Potenzial an den Grenzen auf unendlich setzen müssen, wodurch die Wellenfunktion an diesen Punkten auf Null gesetzt wird. Ändern Sie es für Ihr Potenzial wie folgt:

V = unendlich, wenn -50fm>r>50fm V = 0, wenn |r|>8,5fm V = Woods-Potential sonst

Der eigentliche Deal: Schreiben Sie jetzt ein Programm, um die Wellenfunktion im obigen Potential unter Verwendung einer beliebigen Energie im Bereich -50 fm zu berechnen

Hier ist ein Python-Code, der einen Eigenwert im angegebenen Intervall (E1, E2) findet, falls vorhanden.

import numpy as np

rmin = 0 rmax = 0

def vradial(r): global rmin global rmax if r < rmin or r > rmax: #This is important. >= and <= won't work #as they interfere with the bc on psi return np.inf elif r > rmin and r < rmin + 2: return (r - rmin - 2)**2 elif r > rmax - 2 and r < rmax: return (r - rmax + 2)**2 else: return 0

def f(l,E,r): """Calculate the f(r) in the Schrodinger equation of the form D2(Psi(r)) = f(r)Psi(r)""" if r == 0: return np.inf return l*(l+1)/(r**2) + vradial(r) - E

def psitophi(psi, l, E, r, delta): if psi == 0: return 0 #To handle the 0*infinity case of boundary return psi*(1-(f(l, E, r)*delta**2)/12)

def phitopsi(phi, l, E, r, delta): #if phi == 0: return 0 return phi/(1-(f(l, E, r)*delta**2)/12)

def calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, E):

psiarray = []
r = np.linspace(rmin, rmax, N)
delta = r[1]-r[0]

#Assume psi(rmin) = psibcmin and psi(rmin+delta) = 1 and then
#calculate phi0 and phi1
psiarray.append(psibcmin)
psiarray.append(1)
phi0 = psitophi(psiarray[0], l, E, r[0], delta) #r[0]=rmin
phi1 = psitophi(psiarray[1], l, E, r[1], delta)

#Now populate the psiarray for each value of r
for i in range(2, len(r)):
    phi2 = 2*phi1 - phi0 + (delta**2)*f(l, E, r[i-1])*psiarray[i-1]
    psi = phitopsi(phi2, l, E, r[i], delta)
    psiarray.append(psi)
    phi0 = phi1
    phi1 = phi2

return r, psiarray

def normalize(delta, psi): area = 0 for i in range(1,len(psi)-1): area = area + abs(psi[i])*delta

for i in range(len(psi)):
    psi[i] = psi[i]/area

return psi

def locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax, l, e1, e2,\ tol): #Any value of Psi smaller than psi is while abs(e2-e1) > tol:
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1) psi1 = psi.pop() r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e2) psi2 = psi.pop()

    if psi1*psi2 < 0:
        emid = e1 + (e2 - e1) * 0.5
        r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, emid)
        psimid = psi.pop()

        if psimid*psi1 < 0:
            e2 = emid

        elif psimid*psi2 < 0:
            e1 = emid
    elif psi1*psi2 > 0:
        print "There are either no eigenvalues or too many of them in"+\
        " the given energy interval. For e1={} psi={} and e2={} psi=".format(e1, psi1, e2)+\
        "{}".format(psi2)
        return e1,e2
return e1,e2

def main(): import sys import matplotlib.pyplot as plt global rmin global rmax

e1 = float(sys.argv[1])
e2 = float(sys.argv[2])
l = int(sys.argv[3])
if l == 0:
    rmin = -1
else:
    rmin = 0
tol =1e-5
rmax = 1
psibcmin = 0
psibcmax = 0
N = 100

e1, e2 = locateEvalueInBracket(rmin, rmax, N, psibcmin, psibcmax,\
        l, e1, e2, tol)

fig = plt.figure()
ax = fig.add_subplot(111)
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l, e1)
ax.plot(r,psi, 'o')
r, psi = calcpsiradial(rmin, rmax, N, psibcmin, psibcmax, l ,e2)
ax.plot(r, psi, 'g')
ax.set_title("l = {} and E_blue = {}, E_green={}".format(l,e1,e2))
plt.show()

if name=="main": main()

Bitte verwenden Sie beim Kopieren und Einfügen des obigen Codes sinnvolle Einzüge. Ich kann nicht herausfinden, wie ich es hier schön rendern kann. Ich sende Ihnen eine Kopie, wenn Sie mir eine E-Mail senden.
Numerische Lösung der Schrödinger-Gleichung - Eigenwerte
AntwortenHierDatenschutz-Bestimmungen1y, 5v