Plotten der Eigenzustände

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
SIto
User
Beiträge: 8
Registriert: Donnerstag 13. April 2017, 15:36

Hallo zusammen,

ich versuche gerade einen Plot etwas in die geforderte Form zu bringen, bekomme das aber nicht wirklich hin. Erwartet wurde von mir die Eigenzustände des Henon-Heiles-Potentials zu plotten. Meine Ergebnisse sind:
Bild

Nun sollen sie aber so aussehen:
Bild

In erster Linie interessiert mich, wie ich die Farbe in den gewünschten Zustand bringen kann.
Der Code, den ich fürs plotten verwendet habe lautet:

Code: Alles auswählen

# Plot first K eigenstates
        psik = Psik[:,k].reshape((N,N))

        fig.add_subplot(2, 3, vc)
        ax = fig.gca()
        ax.set_aspect('equal')

        ax.contour(X, Y, v(X,Y), colors="gray", levels=linspace(0, 20, 15))
        ax.contourf(X, Y, abs(psik), levels=linspace(0, 0.15, 40))
Hoffe mal jemand kann mir bei dem Problem helfen, bzw. falls ihr noch mehr Infos braucht, dann könnt ihr gerne bescheid geben.
Gruss Sito
Sirius3
User
Beiträge: 17711
Registriert: Sonntag 21. Oktober 2012, 17:20

@SIto: in den unteren Bildern scheint die Skalierung eine andere zu sein. Auch die x- und y-Achse und die Ordnung scheint anders zu sein. Wenn Du für positive Zahlen eine andere Farbe, wie für negative willst, würde ich das abs weglassen.
SIto
User
Beiträge: 8
Registriert: Donnerstag 13. April 2017, 15:36

Vielen Dank für die Antwort!

Der Asisstent meinte, dass das mit der Sakllierung und Orientierung nicht so tragisch sei, bzw. das wir uns das in der Übungsstunde noch einmal zusammen genauer ansehen werden. Er wollte aber, dass ich das mit den Farben versuche zu ändern..

Ich habe einmal deinen Vorschlag versucht, Problem ist dabei, dass wohl gewisse Werte dann keinen Sinn machen, bzw. nicht richtig interpretiert werden. Der Fehler in der cmd lautet:

Code: Alles auswählen

C:\Users\-\AppData\Local\Programs\Python\Python36-32\lib\site-packages\numpy\ma\core.py:2766: ComplexWarning: Casting complex values to real discards the imaginary part order=order, subok=True, ndmin=ndmin)
Das Bild dazu sieht dann folgendermassen aus:
Bild

Ich denke ich lasse mal den Rest des Codes da, sonst wird es wohl schwer hier eine Lösung zu finden.

Code: Alles auswählen

from numpy import *
from numpy.linalg import eig, norm
from scipy.special import hermite, gamma
from matplotlib.pyplot import *

def arnoldi(A, v0, k):
    r"""Arnoldi algorithm to compute the Krylov approximation :math:`H` of a matrix :math:`A`.

    :param A: The matrix :math:`A` of shape :math:`N \times N` to approximate.
    :param v0: The initial vector :math:`v_0` of length :math:`N`.
    :param k: The number :math:`k` of Krylov steps performed.
    :return: A tuple :math:`(V, H)` where :math:`V` is the large matrix of shape
             :math:`N \times (k+1)` containing the orthogonal vectors and :math:`H` is the
             small matrix of shape :math:`k \times k` containing the Krylov approximation
             of :math:`A`.
    """
    r, c = A.shape
    V = zeros((r, k+1), dtype=complexfloating)
    H = zeros((k+1, k), dtype=complexfloating)

    V[:,0] = v0.reshape(-1) / norm(v0)

    for i in range(1, k+1):
        vi = dot(A, V[:,i-1])
        for j in range(i):
            H[j,i-1] = dot(conjugate(V[:,j]), vi)
            vi -= H[j,i-1] * V[:,j]
        H[i,i-1] = norm(vi)
        V[:,i] = vi / H[i,i-1]

    return V, H[:-1,:]

def Part3():
    # Teil 3: Unteraufgaben g), h), i)

    # Harmonischer Osillator
    v = lambda x, y: 0.5*(x**2 + y**2)

    # Henon-Heiles
    a = 2.0
    b = 0.4
    v = lambda x, y: 0.5*a*(x**2 + y**2) + b*(x**2*y - y**3/3.0)

    N = 32
    x, h = linspace(-3, 3, N, retstep=True)
    X,Y = meshgrid(x,x)
    x = X.reshape(-1)
    y = Y.reshape(-1)

    # Just K for later use
    K = 6


    # Unteraufgabe g)
    
    # Laplace Part
    D = zeros((N*N, N*N), dtype=floating)
    for r in range(0, N):
        for c in range(0, N):
            # X direction stencil
            if c == 0:
                # Forward
                D[N*r+c, N*r+c] += 1.0
                D[N*r+c, N*r+(c+1)] += -2.0
                D[N*r+c, N*r+(c+2)] += 1.0
            elif c == N-1:
                # Backward
                D[N*r+c, N*r+(c-2)] += 1.0
                D[N*r+c, N*r+(c-1)] += -2.0
                D[N*r+c, N*r+c] += 1.0
            else:
                # Central
                D[N*r+c, N*r+(c-1)] += 1.0
                D[N*r+c, N*r+c] += -2.0
                D[N*r+c, N*r+(c+1)] += 1.0

            # Y direction stencil
            if r == 0:
                # Forward
                D[N*r+c, N*r+c] += 1.0
                D[N*r+c, N*(r+1)+c] += -2.0
                D[N*r+c, N*(r+2)+c] += 1.0
            elif r == N-1:
                # Backward
                D[N*r+c, N*(r-2)+c] += 1.0
                D[N*r+c, N*(r-1)+c] += -2.0
                D[N*r+c, N*r+c] += 1.0
            else:
                # Central
                D[N*r+c, N*(r-1)+c] += 1.0
                D[N*r+c, N*r+c] += -2.0
                D[N*r+c, N*(r+1)+c] += 1.0

    D /= h**2
    # Potential Part
    V = diag(v(x,y))
    # Hamiltonian
    H = -0.5*D + V

    figure()
    matshow(H[:50,:50])
    colorbar()
    savefig("discrete_hamiltonian.pdf")
    
    
    # Unteraufgabe h)

    ew, ev = eig(H)
    # Sort the eigenvalues
    I = argsort(ew)
    ew = ew[I]
    Psi = ev[:,I]


    fig = figure(figsize=(12,8))
    for k in range(K):
        # Plot first K eigenstates
        psik = Psi[:,k].reshape((N,N))

        fig.add_subplot(2, 3, k+1)
        ax = fig.gca()
        ax.set_aspect('equal')

        ax.contour(X, Y, v(X,Y), colors="gray", levels=linspace(0, 20, 15))
        ax.contourf(X, Y, abs(psik), levels=linspace(0, 0.15, 40))

        #ax.grid(True)
    fig.savefig("henon_eigenfunctions.pdf")


    # Unteraufgabe i)

    # Startvektor fuer Arnoldi
    v0 = ones(N*N)/N
    VK, HK = arnoldi(H, v0, 220)

    ewk, evk = eig(HK)
    I = argsort(ewk)
    ewk = ewk[I]
    evk = evk[:,I]

    Psik = dot(VK[:,:-1], evk)


    fig = figure(figsize=(12,8))
    vc = 0
    for k in range(0, 2*K):
        # Skip some numerical artefacts
        if abs(ewk[k]) < 1.0:
            continue
        else:
            vc += 1
            if vc > K:
                break

        # Plot first K eigenstates
        psik = Psik[:,k].reshape((N,N))

        fig.add_subplot(2, 3, vc)
        ax = fig.gca()
        ax.set_aspect('equal')

        ax.contour(X, Y, v(X,Y), colors="gray", levels=linspace(0, 20, 15))
        ax.contourf(X, Y, abs(psik), levels=linspace(0, 0.15, 40))

        #ax.grid(True)
    fig.savefig("henon_eigenfunctions_arnoldi.png")




if __name__ == "__main__":
    Part3()
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

SIto hat geschrieben: In erster Linie interessiert mich, wie ich die Farbe in den gewünschten Zustand bringen kann.
Die Farbe änderst Du mit der Colormap:
Bild
http://matplotlib.org/examples/color/co ... rence.html
SIto hat geschrieben:ComplexWarning: Casting complex values to real discards the imaginary part order=order, subok=True, ndmin=ndmin)
Du plottest eine Array mit komplexen Zahlen. Matplotlib warnt Dich, dass es nur den Realteil plottet und nicht den Imaginärteil. Der Fehler tritt aber nicht in dem Code auf, den Du gepostet hast.

Die weißen Gebiete in Deinen Plots kommen von NaN-Werten in den geplotteten Arrays. Auch dieser Effekt tritt nicht in dem geposteten Code auf.

Auf meiner Seite habe ich ein bisschen Beispielcode zum Plotten von Eigenzuständen:
Bild
Bild
http://www.magben.de/?h1=mathematik_fue ... orbital_2d

Bild
http://www.magben.de/?h1=mathematik_fue ... funktionen
a fool with a tool is still a fool, www.magben.de, YouTube
Antworten