Issue
I tried to solve this problem by changing a lot of variables but nothing it's working. Here is my code:
import numpy as np
import matplotlib.pyplot as plt
import math as mt
from numba import njit
N=1000
J=1
h=0.5
plt.rcParams['figure.dpi']=100
plt.xlabel('T', fontsize=14)
ns=100000000
s=np.empty([N,0])
for i in range(N):
s[i]=np.random.choice([-1,1])
E=0
M=0
for i in range(-1,N-1):
E=E-J*(s[i]*s[i+1]+s[i]*s[i-1])
M=M+s[i]
Energy=np.empty([0,0])
C=np.empty([0,0])
Magne=np.empty([0,0])
chi=np.empty([0,0])
@njit
def average(k , E, M, ns, x, Energy, C, Magne, chi, s):
Ev=0.
Ev2=0.
Mv=0.
Mv2=0.
for z in range(1,ns+1):
i=np.random.randint(-1, high=N-1)
dE=2*J*(s[i]*s[i+1]+s[i]*s[i-1])
dM=2*h*s[i]
pace=1./(1+mt.exp(k*(dE+dM)))
#pace=min(1,mt.exp(-k*dE))
if np.random.random()<pace:
s[i]=-s[i]
E=E+dE
M=M-dM/h
if z>x:
Ev=Ev+(E-h*M)
Ev2=Ev2+(E-h*M)**2
Mv=Mv +M
Mv2=Mv2 +M**2
Ev=Ev/(ns-x)
Ev2=Ev2/(ns-x)
varE=(Ev2-Ev**2)
Energy=np.append(Energy,Ev)
C=np.append(C,varE*(k**2))
Mv=Mv/(ns-x)
Mv2=Mv2/(ns-x)
varM=(Mv2-Mv**2)
Magne=np.append(Magne,Mv)
chi=np.append(chi,varM*k)
return E, M, Energy, C, s, Magne, chi
b=np.arange(0.3,1.,0.01)
x=0.75*ns
for k in b:
E, M, Energy, C, s, Magne, chi =average(k,E,M,ns, x,Energy,C, Magne, chi, s)
plt.plot(b,Energy, '.', color='r',linestyle='--')
plt.plot(b,C, '.', color='b',linestyle='--')
plt.title('Energy and heat capacity')
plt.legend(['E','C'])
plt.show()
plt.rcParams['figure.dpi']=100
plt.xlabel(r'$\beta$', fontsize=14)
plt.plot(b,Magne, '.',linestyle='--', color='r')
plt.plot(b,chi, '.', color='b',linestyle='--')
plt.title('Magnetization and susceptibility')
plt.legend(['M',r'$\chi$'])
plt.show()
I keep getting the error:
Cannot unify float64 and array(float64, 1d, C) for 'Mv2.3', defined at c:\users\usuario\documents\python scripts\ex13hw7.py (65)
File "ex13hw7.py", line 65: def average(k , E, M, ns, x, Energy, C, Magne, chi, s): <source elided> if z>x: Ev=Ev+(E-h*M) ^
During: typing of assignment at c:\users\usuario\documents\python scripts\ex13hw7.py (65)
Does anyone know what is the issue?
I tried changing the variable's name and including/deleting parameters from the function.
Solution
Here is fixed version of the code that compiles with numba:
import math as mt
import matplotlib.pyplot as plt
import numpy as np
from numba import njit
N = 1000
J = 1
h = 0.5
plt.rcParams["figure.dpi"] = 100
plt.xlabel("T", fontsize=14)
ns = 100000000
s = np.empty(N) # <-- don't use np.empty([N, 0])
for i in range(N):
s[i] = np.random.choice([-1, 1])
E = 0
M = 0
for i in range(-1, N - 1):
E = E - J * (s[i] * s[i + 1] + s[i] * s[i - 1])
M = M + s[i]
Energy = np.array([]) # <-- don't use np.empty([0, 0])
C = np.array([])
Magne = np.array([])
chi = np.array([])
@njit
def average(k, E, M, ns, x, Energy, C, Magne, chi, s):
Ev = 0.0
Ev2 = 0.0
Mv = 0.0
Mv2 = 0.0
for z in range(1, ns + 1):
i = np.random.randint(-1, high=N - 1)
dE = 2 * J * (s[i] * s[i + 1] + s[i] * s[i - 1])
dM = 2 * h * s[i]
pace = 1.0 / (1 + mt.e ** (k * (dE + dM))) # <-- use math.e ** x
if np.random.random() < pace:
s[i] = -s[i]
E = E + dE
M = M - dM / h
if z > x:
Ev = Ev + (E - h * M)
Ev2 = Ev2 + (E - h * M) ** 2
Mv = Mv + M
Mv2 = Mv2 + M**2
Ev = Ev / (ns - x)
Ev2 = Ev2 / (ns - x)
varE = Ev2 - Ev**2
Energy = np.append(Energy, Ev)
C = np.append(C, varE * (k**2))
Mv = Mv / (ns - x)
Mv2 = Mv2 / (ns - x)
varM = Mv2 - Mv**2
Magne = np.append(Magne, Mv)
chi = np.append(chi, varM * k)
return E, M, Energy, C, s, Magne, chi
b = np.arange(0.3, 1.0, 0.1)
x = 0.75 * ns
for k in b:
E, M, Energy, C, s, Magne, chi = average(k, E, M, ns, x, Energy, C, Magne, chi, s)
print("Finished!")
plt.plot(b, Energy, ".", color="r", linestyle="--")
plt.plot(b, C, ".", color="b", linestyle="--")
plt.title("Energy and heat capacity")
plt.legend(["E", "C"])
plt.show()
plt.rcParams["figure.dpi"] = 100
plt.xlabel(r"$\beta$", fontsize=14)
plt.plot(b, Magne, ".", linestyle="--", color="r")
plt.plot(b, chi, ".", color="b", linestyle="--")
plt.title("Magnetization and susceptibility")
plt.legend(["M", r"$\chi$"])
plt.show()
Shows these two graphs:
Answered By - Andrej Kesely
0 comments:
Post a Comment
Note: Only a member of this blog may post a comment.