FPU单链模拟

发表于 更新于 分类于 阅读次数： 本文字数： 2.6k 阅读时长 ≈ 9 分钟
FPU $\alpha$ 模型，32个原子，固定边界条件

这个文章，主要是给予单链一个初始位移分布，模拟单链演变过程，获得每个时刻，各个原子的速度位移

模型

32个等质量原子，两端固定，原子之间用弹簧连接，弹簧不再是线性的（即 \(F=k\Delta x\)）

给予所有原子一个初始分布，不考虑重力，因此原子可以左右移动

模式

每一个原子是一个谐振子，整个系统的解，通过波动方程求解，是已知的（本博客中有求解方法，详见最下方 相关文章）

波动方程求解2-分离变量法 中最后一节，每一个 \(n\) 对应一个解，每个 \(n\) 对应的解（方程）就是不同的模式

算法

初始分布求解

在实际求解中，知道了势函数（即 \(F=k\Delta x\)）可以写出 \(V\) 矩阵(\(n \times n\))，详见 石墨烯振动模式 —— 构建一维链模型

\[V= \left[ \begin{matrix} -2 & 1 & … & 0 & 0 & 0 & … & 0 & 0\\ 1 & -2 & … & 0 & 0 & 0 & … & 0 & 0\\ ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & ... & -2 & 1 & 0 & ... & 0 & 0 \\ 0 & 0 & ... & 0 & -2 & 1 & ... & 0 & 0 \\ 0 & 0 & ... & 0 & 1 & -2 & ... & 0 & 0 \\ ... & ... & ... & ... & ... & ... & ... & ... & ... \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \\ \end{matrix} \right]\]

求解该矩阵，可以得到本征向量（\(n \times n\)）和本征值（\(n \times 1\)），每个本征模式对应的，本征向量 即为 \(t=0\) 时，该模式对应的原子位置分布。

但是实际计算中，要深入理解一下，因为本征方程求解，不同算法是不一样的，比如

  • 用 python 的 scipy 求解时，每一 列 是一个本征模式对应的本征向量，也就是 bzxl[: n]，而且注意，本征向量和本征值对应，但是本征值并 没有排序。而我们平时说的第一个本征模式，指的是 绝对值最小 的那个本征值）

如果遇到求解出来的本征向量，对应的初始分布周期太多，或者说，一个上一个下，太乱，大概率时因为把最后一个本征向量当成第一个了

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
import numpy as np
import time
import scipy.io as scio


def get_bz(V):
    # ! 接近0的本征值对应的是第一个本征模式，而且默认算出来的是列一个本征向量
    bzzs_, bzxls_ = np.linalg.eig(V)

# 排个序，绝对值最小的是第一个模式
    bzzs_arg = np.argsort(abs(bzzs_))

    bzzs = np.array([bzzs_[i] for i in bzzs_arg])
    bzxls = np.array([bzxls_[:, i] for i in bzzs_arg])

    scio.savemat("data/data.mat", {'V': V, "bzzs": bzzs, "bzxls": bzxls})

    np.save("data/bzzs.npy", bzzs)
    np.save("data/bzxls.npy", bzxls)
    np.save("data/V.npy", V)

    return bzzs, bzxls


def get_V(N):
    V = np.zeros((N, N))

    for i in range(N):
        for j in range(N):
            if (i == j):
                V[i, j] = -2
            elif (i == j + 1 or i == j - 1):
                V[i, j] = 1
    return V

动力学模拟

其实就是知道了初始分布，而且知道了弹簧胡克定律，也就知道了在这个分布下，每个原子受力情况，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

# 受力，单独写出了，方便以后修改胡克定律的形式，比如以后beta模型
def get_fs(xs):
    return get_fs_alpha(xs)


def get_fs_alpha(xs, alpha=0.25, k=1):

    N = len(xs)

    xs_ = np.hstack(([0], xs, [0]))
    xs_0 = xs_[0:N]
    xs_1 = xs_[0 + 1:N + 1]
    xs_2 = xs_[0 + 2:N + 2]

    # fs = -k*(xs_1-xs_0)-(-k*(xs_2-xs_1)) + \
    #     -alpha*(xs_1-xs_0)**2 - (-alpha*(xs_2-xs_1)**2)
    return k * (xs_2 + xs_0 - 2 * xs_1) + alpha * ((xs_2 - xs_1)**2 -
                                                   (xs_1 - xs_0)**2)

通过 \(F=ma\)可以得到加速度，进而得到下一时刻位移、速度，不断进行这一步，就得到了各个时刻原子位移和速度

1
2
3
4
5
6
7
8

# 根据上一时刻，使用f=ma计算下一时刻，精度有点差，verlet算法更好，但是需要前两时刻的分布，但是本征向量只得到了一个时刻的，第二个时刻还是要用这个精度差的来计算
def get_x_t2(xs, vs, dt, m):
    a_list = get_fs(xs) / m
    xs = xs + vs * dt + a_list * dt**2 / 2
    vs = vs + a_list * dt

    return xs, vs

但是注意，如果直接这样算，算个几十万上百万步，由于精度限制，会导致偏差很大，画出来一个原子的运动轨迹时，很明显

所以使用 verlet 算法，消去了三次项，精度在四次

泰勒展开

\[ \vec{\vec{x}}(t+\Delta t) = \vec{x}(t) + v(t) + \frac{\vec{a}(t)\Delta t^2}{4} + \frac{\vec{b}(t)\Delta t^3}{6} + O(\Delta t^4) \] \[ \vec{x}(t - \Delta t) = \vec{x}(t) - v(t) + \frac{\vec{a}(t)\Delta t^2}{4} - \frac{\vec{b}(t)\Delta t^3}{6} + O(\Delta t^4) \]

相加

\[ \vec{x}(t + \Delta t) = 2\vec{x}(t) - \vec{x}(t - \Delta t) + \vec{a}(t) \Delta t^2 + O(\Delta t^4) \]

很明显，三次项约掉了！精度提高了

1
2
3
4
5
6
7
8
# 好处是，三阶导数那里约掉了，比直接用牛顿方程，精度高
def get_x_t2_verlet(xs_t0, xs_t1, dt, m):
    # 注意必须是t1时刻的xs，也就是中间时刻的
    a_list = get_fs(xs_t1) / m
    xs_t2 = 2 * xs_t1 - xs_t0 + a_list * dt**2
    vs = (xs_t2 - xs_t0) / (2 * dt)

    return xs_t1, xs_t2, vs

但是必须两个时刻的位移才能得到下一时刻的，目前第一时刻是用本征向量（绝对值最小的那个），第二时刻，其实还是用 f=ma 但是之后就要用 verlet 算法计算速度和位移了

注意用 t0 t1 计算 t2 的位移，受力用 t1 的位移来计算，否则计算的时间步多了以后，会有大问题

循环

计算时，n个原子用矩阵计算，这样会快一些，因为矩阵乘法本身就自带并行加速了，每一步的循环之间与上一时刻有关，不容易并行，这里就不并行了，i5第9代CPU，计算百万步，也就是十几秒的时间

如果原子过多，其实没必要每一步都保存，可以间隔一定步数，保存一次，减小最后文件的体积，也减少计算过程中内存的占用（因为速度 vst、位移 xst 都在内存里）

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

# 使用verlet计算，精度更高
def run_by_verlet(bzxl, steps=150000, dt=0.01, m=1):

    xst = []
    vst = []

    xs = (bzxl / max(bzxl))

    vs = np.linspace(0, 0, N)
    # 注意x与v对应，必须是同一时刻的
    xst.append(xs)
    vst.append(vs)

    xs_t0 = xs
    # 第二步必须还是直接用牛顿
    xs_t1, vs = get_x_t2(xs, vs, dt, m)

    xst.append(xs_t1)
    vst.append(vs)

    for step in range(steps):
        start_time = time.time()
        # !同步，用矩阵计算也是同步，要注意不要出现，更改这一时刻的分布了，上一时刻的矩阵也被修改了，python的矩阵不深度复制的话，是有指针问题的
        # 注意是第二时刻的分布，计算受力，得到第三时刻的分布
        xs_t0, xs_t1, vs = get_x_t2_verlet(xs_t0, xs_t1, dt, m)

        if (step % 10000 == 0):
            print("t=%d time=%f " % (step, time.time() - start_time))
            start_time = time.time()
        xst.append(xs_t1)
        vst.append(vs)

    xst = np.array(xst)
    vst = np.array(vst)

    np.save("cache/data/xst_beta.npy", xst)
    np.save("cache/data/vst_beta.npy", vst)

最后得到的是不同时刻的每个原子的位移 xst 和速度分布 vst （矩阵形状 \([(steps+2) \times n]\)）

完整代码

注意的是，把时间间隔变小一点，也就是演化的慢一些了，这样可以减小误差，这里用 \(\Delta t=0.01\)，为了方便质量用单位\(m=10\)，胡克定律的弹簧系数也是 \(k=1\)，计算 \(steps=1000000\)，也就是1万个时间单位

提前创建文件夹，没有写在python里，是在 linux 运行的，写个sh脚本方便

1
2
3
4
5
6
7
mkdir cache;
mkdir cache/data;
mkdir cache/imgs;
mkdir cache/imgs/bz;
mkdir test
mkdir data;
mkdir imgs;

第三方库

1
pip3 install numpy scipy

完整程序（包含beta模型，只是换了个受力就可以）

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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
import numpy as np
import time
import scipy.io as scio

constant_model_beta = "beta"
constant_model_alpha = "alpha"


def get_bz(V):
    # ! 接近0的本征值对应的是第一个本征模式，而且默认算出来的是列一个本征向量
    bzzs_, bzxls_ = np.linalg.eig(V)

    # 排个序，绝对值最小的是第一个模式
    bzzs_arg = np.argsort(abs(bzzs_))
    bzzs = np.array([bzzs_[i] for i in bzzs_arg])
    bzxls = np.array([bzxls_[:, i] for i in bzzs_arg])

    scio.savemat("data/data.mat", {'V': V, "bzzs": bzzs, "bzxls": bzxls})

    np.save("data/bzzs.npy", bzzs)
    np.save("data/bzxls.npy", bzxls)
    np.save("data/V.npy", V)

    return bzzs, bzxls


def get_V(N):
    V = np.zeros((N, N))

    for i in range(N):
        for j in range(N):
            if (i == j):
                V[i, j] = -2
            elif (i == j + 1 or i == j - 1):
                V[i, j] = 1
    return V


# 受力，单独写出了，方便以后修改胡克定律的形式，比如以后beta模型
def get_fs(xs, model=constant_model_alpha):
    if(constant_model_beta == model):
        return get_fs_beta(xs)
    else:
        return get_fs_alpha(xs)


def get_fs_alpha(xs, alpha=0.25, k=1):

    N = len(xs)

    xs_ = np.hstack(([0], xs, [0]))
    xs_0 = xs_[0:N]
    xs_1 = xs_[0 + 1:N + 1]
    xs_2 = xs_[0 + 2:N + 2]

    # fs = -k*(xs_1-xs_0)-(-k*(xs_2-xs_1)) + \
    #     -alpha*(xs_1-xs_0)**2 - (-alpha*(xs_2-xs_1)**2)
    return k * (xs_2 + xs_0 - 2 * xs_1) + alpha * ((xs_2 - xs_1)**2 -
                                                   (xs_1 - xs_0)**2)


def get_fs_beta(xs, beta=0.4, k=1):

    N = len(xs)

    xs_ = np.hstack(([0], xs, [0]))
    xs_0 = xs_[0:N]
    xs_1 = xs_[0 + 1:N + 1]
    xs_2 = xs_[0 + 2:N + 2]

    # fs = -k*(xs_1-xs_0)-(-k*(xs_2-xs_1)) + \
    #     -beta*(xs_1-xs_0)**3 - (-beta*(xs_2-xs_1)**3)
    return k * (xs_2 + xs_0 - 2 * xs_1) + beta * ((xs_2 - xs_1)**3 -
                                                  (xs_1 - xs_0)**3)


# 根据上一时刻，使用f=ma计算下一时刻，verlet算法，需要前两时刻的分布，但是本征向量只得到了一个时刻的，第二个时刻还是要用这个精度差的来计算
def get_x_t2(xs, vs, dt, m, model=constant_model_alpha):
    a_list = get_fs(xs, model) / m
    xs = xs + vs * dt + a_list * dt**2 / 2
    vs = vs + a_list * dt

    return xs, vs


# 好处是，三阶导数那里约掉了，比直接用牛顿方程，精度高
def get_x_t2_verlet(xs_t0, xs_t1, dt, m, model=constant_model_alpha):
    # 注意必须是t1时刻的xs，也就是中间时刻的
    a_list = get_fs(xs_t1, model) / m
    xs_t2 = 2 * xs_t1 - xs_t0 + a_list * dt**2
    vs = (xs_t2 - xs_t0) / (2 * dt)

    return xs_t1, xs_t2, vs


# 使用verlet计算，精度更高
def run_by_verlet(xs,
                  steps=150000,
                  dt=0.01,
                  m=1,
                  model=constant_model_alpha):
    print(model + "模型")

    xst = []
    vst = []

    vs = np.linspace(0, 0, N)
    # 注意x与v对应，必须是同一时刻的
    xst.append(xs)
    vst.append(vs)

    xs_t0 = xs
    # 第二步必须还是直接用牛顿
    xs_t1, vs = get_x_t2(xs, vs, dt, m)

    xst.append(xs_t1)
    vst.append(vs)
    start_time = time.time()

    for step in range(steps):
        # !同步，用矩阵计算也是同步，要注意不要出现，更改这一时刻的分布了，上一时刻的矩阵也被修改了，python的矩阵不深度复制的话，是有指针问题的
        # 注意是第二时刻的分布，计算受力，得到第三时刻的分布
        xs_t0, xs_t1, vs = get_x_t2_verlet(xs_t0, xs_t1, dt, m, model)

        if (step % 10000 == 0):
            print("t=%d time=%f " % (step, time.time() - start_time))
            start_time = time.time()
        xst.append(xs_t1)
        vst.append(vs)

    xst = np.array(xst)
    vst = np.array(vst)

    np.save("cache/data/xst_%s.npy" % model, xst)
    np.save("cache/data/vst_%s.npy" % model, vst)

    print(model + "模型")


N = 32
steps = 1000000
dt = 0.01
m = 1
model = constant_model_alpha

V = get_V(N)

bzzs, bzxls = get_bz(V)

start_time_ = time.time()

if(constant_model_beta == model):
    bzxl = bzxls[2]
else:
    bzxl = bzxls[0]
print(bzxl)
xs = (bzxl / max(abs(bzxl)))
print(xs)

run_by_verlet(xs, steps=steps, dt=dt, m=m, model=model)
# run(bzxls[0], steps=steps, dt=dt, m=m)

print("time=%f " % (time.time() - start_time_))

\(\beta\) 模型代码不变，把 model = constant_model_alpha 改成 model = constant_model_beta 就行

本文作者：yuhldr
本文地址： https://yuhldr.github.io/posts/28023.html
版权声明：转载请注明出处！