FPU单链选择定则

发表于 更新于 分类于 阅读次数： 本文字数： 1.6k 阅读时长 ≈ 6 分钟
FPU $\beta$ 模型，32个原子，固定边界条件，各个模式能量被激活规律

前两篇，已经模拟计算了单链 αβ 两种模型，并且复现了 α模型的FPU经典图，接下来，在 β 模型下推导选择定则，并且做相关模拟验证。

推导

N + 1 个原子，两端原子固定，中间 N 个非线性弹簧连接，

x0(t) = ẋ0(t) = xN(t) = ẋN(t) = 0

第j个原子受力（j = 1, 2, ..., N − 1）

$$ \begin{aligned} F_j =& [(x_{j+1}-x_{j})-(x_j-x_{j-1})] \\ &+ \alpha [(x_{j+1}-x_{j})^2-(x_j-x_{j-1})^2] \\ &+ \beta[(x_{j+1}-x_{j})^3-(x_j-x_{j-1})^3] \end{aligned} $$

哈密顿量为：

$$ H=\sum_{j=1}^{N}\left[\dot{x}_{j}^{2} / 2+\left(x_{j}-x_{j-1}\right)^{2} / 2+\alpha\left(x_{j}-x_{j-1}\right)^{3} / 3+\beta\left(x_{j}-x_{j-1}\right)^{4} / 4\right] $$

上面的推导可参见 相关文章

β模型下，3次项可出现正负，令α = 0

对于固定边界条件，谐振子来说解为

$$x_j = i\sum_{k=1}^{N}\frac{a_k}{\omega_k} e^{\frac{-i\pi j}{N}}$$ 其中

$$\omega = 2 sin(\frac{\pi k}{2N})$$

这里，因为是固定边界条件，只看实数部分，即为初始状态振幅，$x_j=\cos{\frac{\pi}{N}j}$，J从0取到N，正好初始状态时，“拱形”固定边界条件

$$ x_{j}-x_{j-1}=i\sum_{k=1}^{N}\frac{a_k}{\omega_k} e^{\frac{-i\pi j}{N}} - i\sum_{k=1}^{N}\frac{a_k}{\omega_k} e^{\frac{-i\pi (j-1)}{N}}= i\sum_{k=1}^{N}\frac{a_k}{\omega_k} e^{\frac{-i\pi j}{N}}(1-e^{\frac{-i\pi}{N}}) $$

$$ \begin{aligned} \sum_{j=1}^{N}\left(x_{j}-x_{j-1}\right)^{4}=& \sum \frac{a_{k} a_{k'} a_{k''} a_{k'''}}{\omega_{k} \omega_{k'} \omega_{k} \omega_{k'''}} \exp \left[-i \pi j\left(k+k'+k''+k'''\right) / N\right] \\ & \times[1-\exp (i \pi k / N)]\left[1-\exp \left(i \pi k' / N\right)\right] \\ & \times\left[1-\exp \left(i \pi k'' / N\right)\right]\left[1-\exp \left(i \pi k''' / N\right)\right] \\ =& \sum a_{k} a_{k'} a_{k''} a_{k''} \exp \left[-i \pi\left(j-\frac{1}{2}\right)\left(k+k'+k''+k'''\right) / N\right] \end{aligned} $$

上面方程最后一步化简：

已知$\omega_k = 2sin(\frac{\pi k}{2N})$，令 $\frac{\pi k}{2N}=y$

可得 $$ \begin{aligned} &1-\exp (i \pi k / N)\\ &=1-[cos(2y)+isin(2y)]\\ &=1-[1-2sin(y)^2+i2sin(y)cos(y)]\\ &=2sin(y)[sin(y)-icos(y)]\\ &=\omega_k (-i)[cos(y)+isin(y)]\\ &= -i \omega_k e^{iy}\\ &=-i \omega_k e^{i\frac{\pi k}{2N}} \end{aligned} $$ 则

$$ \begin{aligned} &\frac{[1-\exp (i \pi k / N)]\left[1-\exp \left(i \pi k' / N\right)\right] \times\left[1-\exp \left(i \pi k'' / N\right)\right]\left[1-\exp \left(i \pi k''' / N\right)\right] }{\omega_k \omega_k' \omega_k'' \omega_k'''}\\ &=(-i)^4 e^{i\pi \frac{1}{2N}(k+k' + k'' + k''')} \end{aligned} $$

实数部分只有余弦项

$$ \sum_{j=1}^{N} \cos \left[\frac{(j-\frac{1}{2})(k+k'+k''+k''') \pi}{N}\right]=N \cdot D\left(k+k'+k''+k'''\right) $$

因此 k为整数时，有

D(0) = 1  D(±2N) = −1  D(*) = 0

$$ \sum_{j=1}^{N}\left(x_{j}-x_{j-1}\right)^{4}=N \sum a_{k} a_{k'} a_{k''} a_{k'''}D\left(k+k'+k''+k'''\right) $$

H = (N/2)∑[(ȧk/ωk)2 + ak2] + (μN/4)∑akak′ak″ak‴(k + k′ + k″ + k‴)

求导得

$$\ddot{a}_{k}=-\omega_{k}^{2} a_{k}-\mu \omega_{k}^{2} \sum_{k' k'' k'''} a_{k'} a_{k''} a_{k'''} D\left(k+k'+k''+k'''\right) $$

欧拉公式可推出 $\cos x = \frac{e^{ix}+e^{-ix}}{2}$ $$ \begin{aligned} &=\sum_{j=1}^{N} \cos \left[\left(j-\frac{1}{2}\right)\left(k+k'+k''+k'''\right) \pi / N\right] \\ &=\sum_{j=1}^{N} \cos[\frac{(2j-1)k \pi}{2N}] = \sum_{j=1}^{N} \frac{e^{i\frac{(2j-1)k \pi}{2N}} + e^{-i\frac{(2j-1)k \pi}{2N}}}{2}\\ &=a_1 \frac{(1-q^n)}{1-q} + a_1'\frac{(1-q'^n)}{1-q'} \end{aligned} $$

其中，$q=e^{\pm \frac{ik\pi}{N}}$， q = 1的时候即为2N  0

等比数列以后为（忽略分母2）： $$ \begin{aligned} &=e^{i\frac{k \pi}{2N}}\frac{(1- e^{ik\pi})}{1-e^\frac{ik\pi}{N}}+e^{-i\frac{k \pi}{2N}}\frac{(1- e^{-ik\pi})}{1-e^\frac{-ik\pi}{N}} =e^{i\frac{k \pi}{2N}}(\frac{(1- e^{ik\pi})}{1-e^\frac{ik\pi}{N}})+e^{i\frac{k \pi}{2N}}\frac{(1- e^{-ik\pi})}{e^\frac{ik\pi}{N}-1}\\ &=e^{i\frac{k \pi}{2N}}(\frac{(1- e^{ik\pi})-(1- e^{-ik\pi})}{1-e^\frac{ik\pi}{N}})=e^{i\frac{k \pi}{2N}}(\frac{e^{-ik\pi}- e^{ik\pi}}{1-e^\frac{ik\pi}{N}})\\ &=e^{i\frac{k \pi}{2N}}(\frac{\cos{(-k\pi)}+i\sin(-k\pi)- (\cos{(k\pi)}+i\sin(k\pi))}{1-e^\frac{ik\pi}{N}})\\ &=e^{i\frac{k \pi}{2N}}\frac{-2i\sin{k\pi}}{1-e^\frac{ik\pi}{N}} \end{aligned} $$

∴ 只要 k 为整数，求和即为0

这个方程的含义是，整个势能项，是由多个模式（kk′k″k‴）组合而成，其中（k′k″k‴）是原有模式，k是被原有模式激发起来的新模式，模式的取值为 1、2、……、N

$$ \sum_{j=1}^{N}\left(x_{j}-x_{j-1}\right)^{4}=N \sum a_{k} a_{k'} a_{k''} a_{k'''}D(k+k'+k''+k''') $$

假设激发了 k0 = 1 模式，也就是初始状态是按照第一个模式分布的，因此初始各个模式的振幅只有第一个不为0： ak0 ≠ 0

如果上面这个方程不为0，由于 akak′ak″ak‴ 的存在，其中每一项都不能是0，因此k′ = k″ = k‴ = k0，k被原有模式激发，也不为0 D(0) = 1  D(±2N) = −1  D(*) = 0 由于上面这个式子，D(k + k′ + k″ + k‴) 这个函数不为 0，只能是 k + k′ + k″ + k‴ 为 2N 的倍数，注意 k = −k

如果，初始激发模式为 k0 = 3 = k′ = k″ = k‴，那么，下一个被激发的模式只能是 9 − 3 − 3 − 3，即 k = 9，再接下来 15 − 3 − 3 − 9，即 k = 15，接下来：21  27

模拟

β 模型，其实前两篇blog已经做了模拟，改一下参数就行，现在验证 β 模型下，激发第3个模式，看看各个模式能量分布。

代码又重新改了一下，整体上和前两个bolg的差不多，一共三个文件，直接点击打开可能乱码，其实打开以后右键保存，是没有乱码的，编码用的 utf-8，在 ubuntu2004、python3.8 下测试

\beta模型，奇数模式能量

可以看到，与推测的一致：初始激发3模式，然后9、15、21、27被激发

参考文献

  • Nonlinear coupled oscillators: Modal equation approach

    10.1016/0021-9991(73)90169-1

  • q -Breathers and the Fermi-Pasta-Ulam Problem

    10.1103/PhysRevLett.95.064102

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