FPU单链模拟

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

模型

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]，而且注意，本征向量和本征值对应，但是本征值并 没有排序。而我们平时说的第一个本征模式，指的是 绝对值最小 的那个本征值）

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

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

动力学模拟

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

# 受力，单独写出了，方便以后修改胡克定律的形式，比如以后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\)可以得到加速度，进而得到下一时刻位移、速度，不断进行这一步，就得到了各个时刻原子位移和速度

# 根据上一时刻，使用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) \]

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

# 好处是，三阶导数那里约掉了，比直接用牛顿方程，精度高
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 都在内存里）

# 使用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脚本方便

mkdir cache;
mkdir cache/data;
mkdir cache/imgs;
mkdir cache/imgs/bz;
mkdir test
mkdir data;
mkdir imgs;

第三方库

pip3 install numpy scipy

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

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
版权声明：转载请注明出处！