Рунге-Кутта для связанного векторного поля в Python

Question

Я хочу решить нестационарное уравнение движения Гейзенберга для спин-решеточной системы, используя Рунге-Кутту в Python. Но я открыт для любого предложения, как решить уравнение.

Начиная с системы из N * N классических векторов на решетке tri angular в 120-градусной структуре (каждый вектор повернут на 120 ° к своим соседям) с взаимодействиями только ближайших соседей, определяя матрицу с (N , N, 3) записей.

Обычно можно было бы ожидать, что при некотором возбуждении на спинах (в виде зависящего от времени переворота спина) появилось бы некоторое распространение по решетке, но я вижу, что пока нет движения векторов ( или просто немного) Уравнение выглядит как , но в моем случае немного проще (источник: https://arxiv.org/pdf/1711.03778.pdf). Я использую только один J, и у меня нет поля magneti c B. Местное поле - это первое слагаемое в этом уравнении:

def localField(J,S,i,j):
        n = len(S)
        h = []
        hx = - J*(S[(i - 1) % n,j][0]+S[(i + 1) % n,j][0] +S[i,(j - 1) % n][0] +S[(i - 1) % n,(j + 1) % n][0]+S[(i + 1) % n,(j - 1) % n][0] +S[i,(j + 1) %n][0]) 
        hy = - J*(S[(i - 1) % n,j][1]+S[(i + 1) % n,j][1] +S[i,(j - 1) % n][1] +S[(i - 1) % n,(j + 1) % n][1]+S[(i + 1) % n,(j - 1) % n][1] +S[i,(j + 1) %n][1])
        hz =- J*(S[(i - 1) % n,j][2]+S[(i + 1) % n,j][2] +S[i,(j - 1) % n][2] +S[(i - 1) % n,(j + 1) % n][2]+S[(i + 1) % n,(j - 1) % n][2] +S[i,(j + 1) %n][2]) 

        h.append(hx)
        h.append(hy)
        h.append(hz)
        return(h)
    def HeisenEqM2(conf,J1):
        S=np.copy(conf)
        n=len(S)
        conf_sum = np.zeros(3)
        Snew = np.zeros((n,n,3))
        locF = np.zeros(3)
        for i in range(n):
            for j in range(n):

                localFie = localField(J1,S,i,j)



                Snew[i,j] += np.cross(S[i,j],localFie)

        return(Snew)

Затем я начинаю вращать один спин в середине решетки и хочу Посмотрим, распространяется ли волна через решетку (согласно теории спиновых волн).

def circ_ex_2(omega,t,config):
    n = len(config)
    S = np.copy(config)

    i = int(n/2)
    j = int(n/2)

    S[i,j][0] += np.cos(omega*t)
    S[i,j][1] += np.sin(omega*t)
    S[i,j][2] += 0
    S = norm_rand(S)

    return(S)

Моя идея: Это мое определение метода Рунге-Кутты. В конце я нормализую свою новую конфигурацию вращения, чтобы гарантировать, что размер векторов остается постоянным, и вычисляю энергию, чтобы увидеть, что происходит:

def runge_kutta_4th(S,omega,J1,Tstart,Tmax,tsteps):
    n=len(S)

    S2 = np.copy(S)
    tt = np.linspace(Tstart,Tmax,tsteps)
    dt = (Tmax-Tstart)/tsteps

    en = []
    for i in range(tsteps):

        S1 = np.copy(circ_ex(omega,tt[i],S2))
        k1 = HeisenEqM2(S1,J1)
        k2 = HeisenEqM2(circ_ex(omega,tt[i]+dt/2,S1+ dt/2*k1) ,J1) 
        k3 = HeisenEqM2(circ_ex(omega,tt[i]+dt/2,S1+ dt/2*k2)  ,J1) 
        k4 = HeisenEqM2(circ_ex(omega,tt[i]+dt,S1+ dt*k3)  ,J1) 
        S2 =S1 +  dt*(1/6*k1 + 1/3*k2 + 1/3*k3 + 1/6*k4)

        S2 = norm_rand(S2)

        en.append(JEnergy(J1,S2)/n**2)
    return(S2,en)

Для воспроизведения кода: Если Вы хотите воспроизвести код, который вы можете начать со следующей системы спин-решеток:

#defining a neel lattice
def neelM(nsize):
    X=np.ones((nsize,nsize))*np.pi/2
    x= np.pi/2
    for i in range(nsize):
        for j in range(nsize):
            X[i,j] += 2*np.pi/3*(i+2*j)
    P = polar(X)
    return(P)

#returning neel lattice as XY components 
def neelMXY(M):
    n = len(M[0])
    v = np.zeros((int(n**2),3))
    X = M[0].flatten()
    Y = M[1].flatten()
    for i in range(n**2):
        v[i,0] += X[i]
        v[i,1] += Y[i]
        v[i,2] = 0
    return(v)

Она называется Neel-State и является основным состоянием три angular Heisenberg Antiferroma gnet. Тогда вам также потребуется нормализация спиновых векторов, определенных очень просто:

def norm_rand(M):
    n = len(M)
    Mnew = np.zeros(np.shape(M))
    for i in range(n):
        for j in range(n):
            norm = np.linalg.norm(M[i,j])
            Mnew[i,j] += M[i,j]/norm
    return(Mnew)

Чтобы увидеть, где находится энергия вашей системы, я определил энергию Гейзенберга, которая является просто скалярным произведением одного вектора со всеми его соседи.

def JEnergy(J,S):
    n=np.shape(S[1])[0]
    H=0
    counter = 0
    for i in range(n):
        for j in range(n):

            H += 1/2*J*(np.dot(S[i,j],S[(i - 1) % n,j])+np.dot(S[i,j],S[(i + 1) % n,j])+np.dot(S[i,j],S[i,(j - 1) % n])+np.dot(S[i,j],S[(i - 1) % n,(j + 1) % n])+np.dot(S[i,j],S[(i + 1) % n,(j - 1) % n])+np.dot(S[i,j],S[i,(j + 1) %n]))
            counter+=1
    return(H) 

В конце концов вы можете выполнить следующее:

neel24 = np.reshape(neelMXY(neelM(24)),(24,24,3))

rk24 = runge_kutta_4th(neel24,5*np.pi,1,0,10,10000)

Tt = np.linspace(0,10,10000)
list0 = []
for i in range(10000):
    list0.append(JEnergy(1,circ_ex(5*np.pi,Tt[i],iter24)))

plt.plot(np.linspace(0,10,10000),rk24[1])
plt.plot(np.linspace(0,10,10000),np.dot(list0,1/24**2))

Это показывает, что простое переключение одного вращения во времени имеет больший эффект, чем метод Рунге-Кутта. Но это не может иметь место, так как в случае перекрестного произведения в уравнении движения изменение одного спина должно влиять и на спины и, следовательно, приводит к более высокому изменению энергии. Я также нанес на карту векторов, используя quiver, и это показывает, что векторы не сильно меняются со временем. Этот график показывает энергию:

Было бы здорово, если бы кто-то мог мне помочь. Эту проблему следует решить, поскольку в приведенном выше документе проделана аналогичная работа, но с более сложной системой.