Mittelwert ermitteln

[Micha_uni]

Hallo liebe Freunde!
Weiß jemand vielleicht, wie man < M**4 > / < M**2 > **2 in Python konstruieren kann?
M steht für die Magnetisierung. Sie wird in einer Schleife wird durch 1000 Durchläufe ermittelt. Nun muss ich den Mittelwert bilden, weiß aber nicht, wie es geht.

Code: Alles auswählen

def initialzustand(N):
    zustand = 2*np.random.randint(2, size=(N,N))-1
    return zustand

def mcmove(konfig, beta):
    for i in range(N):
        for j in range(N):
            a = np.random.randint(0,N)
            b = np.random.randint(0,N)
            s = konfig[a,b]
            nb = konfig[(a+1)%N,b]+konfig[a,(b+1)%N]+konfig[(a-1)%N,b]+konfig[a,(b-1)%N]
            cost = 2*s*nb
            if cost < 0 :
                s*= -1
            elif rand() < np.exp(-cost*beta):
                s*= -1
            konfig[a,b]=s
    return konfig

def calcEnergie(konfig):
    Energie = 0
    for i in range(len(konfig)):
        for j in range(len(konfig)):
            S = konfig[i,j]
            nb = konfig[(i+1)%N,j]+konfig[i,(j+1)%N]+konfig[(i-1)%N,j]+konfig[i,(j-1)%N]
            Energie += -nb*S
    return Energie/4.

def calcMag(konfig):
    return np.sum(konfig)

nt = 500
N = 10
eqSteps = 1000
mcSteps = 1000

T          = np.linspace(1,4,nt)
Kumulante  = np.zeros(nt)

for m in range(len(T)):
    M1 = M2 = 0
    konfig = initialzustand (N)

    for i in range(eqSteps):
        mcmove(konfig,1.0/T[m])

    for i in range(mcSteps):
        mcmove(konfig, 1.0/T[m])
        Mag = calcMag(konfig)
        M1 = M1 + Mag
        M2 = M2 + Mag*Mag;
       [u] Kumulante[m]   =( Hier muss ich[b] < M**4 > / < M**2 > **2[/b] eingeben und anschließend plotten)[/u] 
       Der erste Term soll irgedwie         [u]M2**4/mcSteps[/u]  aussehen. oder zumindest habe ich was ähnliches in 
       einem Buch gesehen. 

plt.plot(T, Kumulante, 'o',color='red', markersize=3, label='Specific Heat');
plt.xlabel('Temperature (T)', fontsize=20);
plt.ylabel('Kumulante',fontsize=20);
plt.show()

Für eine Hilfe wäre ich dankbar
Liebe Grüße
Micha

[Sirius3]

In mcmove wird i und j gar nicht benutzt. Das Aufaddieren der nächsten Nachbarn kommt zweimal im Code vor, wäre also ein Kandidat für eine eigene Funktion. Obwohl man calcEnergie auch als

Code: Alles auswählen

def calcEnergie(konfig):
    nb = (np.roll(konfig, 1, 0) + np.roll(konfig, -1, 0) + np.roll(konfig, 1, 1) + np.roll(konfig, -1, 1))
    return (konfig * nb).sum() / 4

schreiben könnte.
M ist nicht definiert, über was soll also der Mittelwert gebildet werden? Und wie man einen Mittelwert ausrechnet, sollte Dir doch hoffentlich geläufig sein, notfalls gibt es auch np.mean.

[Micha_uni]

Hallo und danke für die schnelle Antwort.
hier ist noch mal der Code

Code: Alles auswählen

def initialstate(N):
    state = 2*np.random.randint(2, size=(N,N))-1
    return state

def mcmove(config, beta):
    for i in range(N):
        for j in range(N):
            a = np.random.randint(0,N)
            b = np.random.randint(0,N)
            s = config[a,b]
            nb = config[(a+1)%N,b]+config[a,(b+1)%N]+config[(a-1)%N,b]+config[a,(b-1)%N]
            cost = 2*s*nb
            if cost < 0 :
                s*= -1
            elif rand() < np.exp(-cost*beta):
                s*= -1
            config[a,b]=s
    return config

def calcEnergy(config):
    energy = 0
    for i in range(len(config)):
        for j in range(len(config)):
            S = config[i,j]
            nb = config[(i+1)%N,j]+config[i,(j+1)%N]+config[(i-1)%N,j]+config[i,(j-1)%N]
            energy += -nb*S
    return energy/4.

def calcMag(config):
    mag = np.sum(config)
    return mag




nt = 512
N = 16
eqSteps = 2000
mcSteps = 2000

T               = np.linspace(1,4,nt)
Energy          = np.zeros(nt)
Magnetization   = np.zeros(nt)
SpecificHeat    = np.zeros(nt)
Susceptibility  = np.zeros(nt)

for m in range(len(T)):
    E1 = M1 = E2 = M2 = 0
    config = initialstate (N)

    for i in range(eqSteps):
        mcmove(config,1.0/T[m])

    for i in range(mcSteps):
        mcmove(config, 1.0/T[m])
        Ene = calcEnergy(config)
        Mag = calcMag(config)

        E1 = E1 + Ene
        M1 = M1 + Mag
        M2 = M2 + Mag*Mag;
        E2 = E2 + Ene*Ene;

        Energy[m]           = E1/(mcSteps*N*N)
        Magnetization[m]    = M1/(mcSteps*N*N)
        SpecificHeat[m]     =(E2/mcSteps - E1*E1/(mcSteps*mcSteps))/(N*T[m]*T[m]);
        Susceptibility[m]   =(M2/mcSteps - M1*M1/(mcSteps*mcSteps))/(N*T[m]);


f = plt.figure(figsize=(18,10), dpi=80, facecolor='w', edgecolor='k');

sp = f.add_subplot(2,2,1);
plt.plot(T, Energy, 'o', color='#A60628',label='Energy');
plt.xlabel('Temperature (T)', fontsize=20);
plt.ylabel('Energy', fontsize=20);


sp = f.add_subplot(2,2,2);
plt.plot(T, abs(Magnetization),'*', color='#348ABD', label='Magnetization');
plt.xlabel('Temperature (T)', fontsize=20);
plt.ylabel('Magnetization', fontsize=20);


sp = f.add_subplot(2,2,2);
plt.plot(T, abs(Magnetization),'*', color='#348ABD', label='Magnetization');
plt.xlabel('Temperature (T)', fontsize=20);
plt.ylabel('Magnetization', fontsize=20);


sp = f.add_subplot(2,2,3);
plt.plot(T, SpecificHeat, 'd', color='black', label='Specific Heat');
plt.xlabel('Temperature (T)', fontsize=20)
plt.ylabel('Susceptibility', fontsize=20);

sp = f.add_subplot(2,2,4);
plt.plot(T, Susceptibility, '+',color='green', label='Specific Heat');
plt.xlabel('Temperature (T)', fontsize=20);
plt.ylabel('Susceptibility',fontsize=20);
plt.show()

Dass der Code einige Ungeschicktheiten hat, weiß ich. Die Fehler sind aber für die Berechnung nicht relevant, weil der Code funktioniert.
Hier ist Mag die Magnetisierung. M1=M1+Mag wird dann 1000 mal Wiederholt. (mcSteps steht für die Monte Carlo Schritte)
Für die Suszebtibilität haben Wir folgende Formel Suscebtitility = (< M^2 > - < M >^2)/T
Das ist hier in der Zeile
Susceptibility[m] =(M2/mcSteps - M1*M1/(mcSteps*mcSteps))/(N*T[m])
als Code geschrieben worden, was auch funktioniert.
Der Unterschied zwischen M1 und M2 liegt daran, dass M1 Vorzeichen behaftet ist, das Vorzeichen von M2 wird direkt durch das Quadrat abgefangen.
die Mag-Werte sind für die 1000 Durchläufe unterschiedlich. und nun muss darüber gemittelt werden.

Ich möchte einfach Kumulante = < M^4 > / < M^2 >^2 eingeben. Es klappt aber nicht. Ich erhalte immer die Konstante Null.
Die Funktion aber soll eine monoton steigende Kurve geben, von der Form her ähnlich wie die e-Funktion.
Der erster Term in Suscebtibility erntspricht, meines Wissens, dem Nenner in Kumulante. Wobei es muss noch quadriert werden.
Ich mache irgendwo einen Denkfehler, weiß aber nicht wo.

Liebe Grüße
Micha

[Sirius3]

@Micha_uni: neben der Tatsache, dass die Einrückung stellenweise keinen Sinn macht, wo ist nun das konkrete Problem, höhere Momente auszurechnen, außer vielleicht dass Du immer mehr durchnummerierte Variablen bekommst, die Du günstigerweise in ein Array wandelst.

[Micha_uni]

Mein Kumpel und ich arbeiten an unseren Abschlussarbeiten.
Was Physik und Mathe. angeht, kennen wir uns mit dem Thema aus.
Was Python angeht, sind wir blutige Anfänger.
Diesen Code haben wir gerade mal mit ach und krach gebastelt.

Für jeden Tipp für die Verbesserung sind wir euch dankbar.
Wir können <M^4 > / < M^2>^2 nicht bilden.

Über array und Liste gibt es einiges formales im Internet, wir haben aber keine Ahnung, wie genau man diese Werkzeuge für unser Anliegen benutzen kann.

Liebe Grüße
Micha

[Sirius3]

@Micha_uni: dass Du <M^2> bilden kannst, hast Du schon gezeigt. Also wo ist das Problem bei <M^4>?

[Siegfried]

@Micha_uni: Könnte es an der Division liegen? Bei Python 2.x ergibt 5 / 4 solange 1, bis man

[codebox=pys60 file=Unbenannt.txt]
from __future__ import division
[/code]
eingibt, danach ergibt 5 / 4 wie erwartet 1.25