基于路径积分-蒙特卡洛方法的简单一维场

发表于

相关理论

一维场即一维空间加上一维时间

基于拉氏量的路径积分

一维标量场给等效成空间内分立格点的（在空间维与时间维具有双重周期条件的）类似与固体物理的布洛赫定理的场。

核心方法：量子力学之中的变分法。稳定的经典场系统的作用量总是取最小值，量子力学之中不确定原理导致非最小作用量的状态仍然有机会出现。量子系统的具体状态由概率确定，由作用量小的状态向作用量大的状态变化（非时间演化）的概率由作用量的差值的负指数确定。

格点方法下的场的拉格朗日量。具体方法见ppt截图。

tu1

tu2

编程与计算

python与C++混合编程。C++编译计算模块加速运算性能，利用python的众多的数据处理与分析工具进行可视化输出与处理。

相关的C++程序：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <random>
#include <vector>
#include <algorithm>

using namespace std;

extern "C"
{
    int cal(double *dog, int Ndex, int Tdex, int Xdex, int Nmax, int Tmax, int Xmax, double msq, double lam, double dx);
}

double ps(double *dog, int X, int T, int N, int Xmax, int Tmax, double lam, double msq, int Nindex, int Xindex, int Tindex)
{
    int Xb = X - 1;
    int Xf = X + 1;
    int Tb = T - 1;
    int Tf = T + 1;
    if (Xb < 0)
    {
        Xb = Xmax - 1;
    }
    if (Tf >= Tmax)
    {
        Tf = 0;
    }
    if (Tb < 0)
    {
        Tb = Tmax - 1;
    }
    if (Tf >= Tmax)
    {
        Xf = 0;
    }

    return pow(dog[N * Nindex + T * Tindex + X * Xindex] - dog[N * Nindex + Tf * Tindex + X * Xindex], 2) + pow(dog[N * Nindex + T * Tindex + X * Xindex] - dog[N * Nindex + Tb * Tindex + X * Xindex], 2) + pow(dog[N * Nindex + T * Tindex + X * Xindex] - dog[N * Nindex + T * Tindex + Xf * Xindex], 2) + pow(dog[N * Nindex + T * Tindex + X * Xindex] - dog[N * Nindex + T * Tindex + Xb * Xindex], 2) + msq * pow(dog[N * Nindex + T * Tindex + X * Xindex], 2) + lam * msq * pow(dog[N * Nindex + T * Tindex + X * Xindex], 2);
}

int cal(double *dog, int Ndex, int Tdex, int Xdex, int Nmax, int Tmax, int Xmax, double msq, double lam, double dx)
{
    int acp = 0;
    int xx;
    int tt;
    std::default_random_engine seed;
    seed.seed(time(NULL));
    // seed.seed(1234);
    // std::uniform_real_distribution<double> ran1(-0.1,0.1);
    std::uniform_real_distribution<double> ran2(-1.0, 1.0);
    std::uniform_real_distribution<double> ran3(0, 1.0);

    std::vector<int> T_iter;
    std::vector<int> X_iter;
    for (int i = 0; i < Tmax; i++)
    {
        T_iter.push_back(i);
    }
    for (int i = 0; i < Xmax; i++)
    {
        X_iter.push_back(i);
    }

    for (int in = 1; in < (Nmax - 1); in++)
    {
        for (int ix = 0; ix < Xmax; ix++)
        {
            for (int it = 0; it < Tmax; it++)
            {
                dog[(in + 1) * Ndex + ix * Xdex + it * Tdex] = dog[(in)*Ndex + ix * Xdex + it * Tdex];
            }
        }
        std::random_shuffle(X_iter.begin(), X_iter.end());
        std::random_shuffle(T_iter.begin(), T_iter.end());

        for (int ix = 0; ix < Xmax; ix++)
        {
            for (int it = 0; it < Tmax; it++)
            {
                xx = X_iter[ix];
                tt = T_iter[it];

                dog[(in + 1) * Ndex + xx * Xdex + tt * Tdex] = dog[(in + 1) * Ndex + xx * Xdex + tt * Tdex] + dx * ran2(seed);
                double ds = ps(dog, xx, tt, (in + 1), Xmax, Tmax, lam, msq, Ndex, Xdex, Tdex) - ps(dog, xx, tt, (in), Xmax, Tmax, lam, msq, Ndex, Xdex, Tdex);

                double pan = ran3(seed);

                if (pan < exp(-ds))
                {
                    /* code */
                }
                else
                {

                    dog[(in + 1) * Ndex + xx * Xdex + tt * Tdex] = dog[(in)*Ndex + xx * Xdex + tt * Tdex];
                }
            }
        }
    }
    return 0;
}

编译指令：

1
2

g++ -shared -fPIC d.cpp -o d.so

python部分:

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import matplotlib.pyplot as plt
import time
#plt.rcParams['font.sans-serif']=['FangSong_GB2312'] #用来正常显示中文标签
from matplotlib.font_manager import _rebuild
_rebuild()#重新创建字体索引列表
import matplotlib
matplotlib.rcParams['font.family']=['Microsoft YaHei']
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
import math
from ctypes import *
import numpy.ctypeslib as npct

面向对象的编程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
class zutai(object):
    def __init__(self,a=0.1,N=6000,la=0,m=1,ls=1):
        self.mass=m
        self.a=a
        self.N=N
        self.nx=int(5/a)
        self.nt=int(10/a)
        self.msq=math.pow(m*a,2)
        self.lam=la*a*a
        self.ls=ls
        self.dog=np.zeros([self.N,self.nt,self.nx],dtype='float64')
        self.Ndex=self.dog.shape[1]*self.dog.shape[2]
        self.Tdex=self.dog.shape[2]
        self.Xdex=1
        self.Nmax=self.dog.shape[0]
        self.Tmax=self.dog.shape[1]
        self.Xmax=self.dog.shape[2]
        self.ldoglib=npct.load_library('d','.')
        array_1d_double = npct.ndpointer(dtype=np.float64, ndim=3, flags='CONTIGUOUS')
        self.ldoglib.cal.argtypes=[array_1d_double,c_int,c_int,c_int,c_int,c_int,c_int,c_double,c_double,c_double]
        self.ldoglib.cal.restype=c_int
        pass
    def cal(self):
        self. ldoglib.cal(self.dog,self.Ndex,self.Tdex,self.Xdex,self.Nmax,self.Tmax,self.Xmax,self.msq,self.lam,self.ls)
        pass
    def tu(self):
        av=np.zeros(self.dog.shape[0])
        for lu in range(0,self.dog.shape[0]):
            av[lu]=np.mean(self.dog[lu,:,:]**2)
            pass
        plt.figure()
        plt.plot(av)
        tit='$\phi^{2}$平均值变化图'+' a='+str(self.a)+'$\lambda = $ '+str(self.lam)+' m='+str(self.mass)
        plt.title(tit)
        plt.xlabel('组态数')
        plt.savefig('phi'+str((time.time()))+'.png')
        plt.show()
#     def jt(self,NUM=70,Lang=25,St=1000,jian=50,zhi=21):
#         fu=np.zeros([NUM,Lang])
#         llk=np.zeros(NUM)
#         for kl in range(0,NUM):
#             for ti in range(0,Lang):
#                 fu[kl,ti]=np.mean((self.dog[St+jian*kl,0,:])*(self.dog[St+jian*kl,ti,:]))
#                 pass
#             pass
       
        
#         llu=fu.mean(axis=0)
#         llue=fu.var(axis=0)
# #         print('误差：'+str(llue))
#         lluu=(np.log(np.abs(llu)))[:zhi]
#         x=np.linspace(0,(zhi-1)/10,zhi)
#         ooo=np.polyfit(x,lluu,1)
#         yyy=ooo[0]*x+ooo[1]
#         print(ooo[0])
#         self.eng=ooo[0]
#         plt.figure()
#         plt.plot(x,lluu,'rx')
#         plt.plot(x,yyy,'-')
# #         for nnn in range(0,NUM):
# #             ddog=(np.log(np.abs(llu)))[nnn,:zhi]
# #             llk[nnn]=np.polyfit(x,ddog,1)
# #             pass
# #         self.res=llk
        pass

    def ek(self,St=3000,jiange=98,dianshu=30,shijian=10,qv=6):
        lp=np.zeros([dianshu,self.nx])
#         tfrt=np.zeros([50,10,50])
        for up in range(dianshu):

            hu=self.dog[St+jiange*up,:,:]

            ss=np.zeros_like(hu,dtype='complex')
            for ii in range(self.nt):
                ss[ii,:]=np.fft.fft(hu[ii,:])
            frt=np.zeros([shijian,self.nx],dtype='complex')
            for aa in range(shijian):
                er=np.zeros(self.nx)
                for bb in range(self.nt):
                    er=er+(ss[bb,:].conjugate()*ss[(bb+aa)%self.nt,:])
        #             if aa==9 and up==1:
        #                 print((bb+aa)%100)
        #             pass
                frt[aa,:]=er
            for down in range(self.nx):
                sls=frt[:,down]
                wew=np.log(sls)
        #         wew.shape
                xcx=np.linspace(0,1-self.a,shijian)
                xcx.shape
                llll=np.polyfit(xcx[:qv],wew[:qv],1)[0]
                #print()
                lp[up,down]=(float(llll))
        mme=np.zeros(self.nx)
        meer=np.zeros(self.nx)
        xday=np.zeros(self.nx)
        for kku in range(self.nx):
            mme[kku]=-lp[:,kku].mean()
            meer[kku]=lp[:,kku].var()
            xday[kku]=kku
        plt.figure()
        plt.errorbar(xday[:int(self.nx/2)],mme[:int(self.nx/2)],yerr=meer[:int(self.nx/2)],fmt='rx:',ecolor='b',capthick=2,capsize=4)
        tit='色散关系图 a='+str(self.a)+'$\lambda = $ '+str(self.lam)+' m='+str(self.mass)
        plt.title(tit)
        plt.xlabel("k")
        plt.ylabel('m')
#         plt.savefig()
        plt.savefig('san'+str((time.time()))+'.png')
        plt.show()
        return [mme,meer]
1
2
3
4
dd=zutai()#初始化
dd.cal()#进行计算
dd.tu()#画出平均数据
[a,b]=dd.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(m=2)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
[a,b]=d2.ek(jiange=74,St=3040)#数据不够理想，可以更改采样点
#算了不调了，就这样了
/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(m=3)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
[a,b]=d2.ek(jiange=27)#生成色散关系
/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

m越大越离谱，减小m试试看

1
2
3
4
d2=zutai(m=0.5)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(m=0.25)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=0.5)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=10)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=0.5)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=1,a=0.1)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=1,a=0.05)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=1,a=0.2)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

1
2
3
4
d2=zutai(la=1,a=0.125)#初始化
d2.cal()#进行计算
d2.tu()#画出平均数据
[a,b]=d2.ek()#生成色散关系

png

/home/yu/.local/lib/python3.6/site-packages/ipykernel_launcher.py:94: ComplexWarning: Casting complex values to real discards the imaginary part

png

数据点过少，结果看看就行，不要当真。