画固定方程曲线图像csdn_用 Mathematica 作振子的色散图像

在上固体物理的时候，大家一定接触过最简单的模型：一维振子链模型。在这个简单模型中，振子被看做是在一个固定点附近做小振动，我们一般用晶体中的格点来表示这个固定点，在对势能与耦合做一定的近似之后我们会得到一个比较简单的运动方程，进而得到简单的色散关系。然而这个关系是否可以通过软件计算然后画出来呢？我在网上找了很久，没有看到有相关的文献研究这个问题，大概是因为太trivial了吧。然而作为十分无聊的大学生，我觉得我可以稍微填补以下这里的空白，给我以后画各种色散图铺一下路，当然如果你有幸看到这篇文章，又对振子色散关系有一些疑问，不妨仔细阅读一下。

一、书本上的内容

首先我们先把课本上的内容简单重复一下：选取偏移量为自变量，假设有多个振子，则坐标可以写作

，然后考虑两个振子之间的相互作用势能

于是取原点为零势能点，然后由

，可以得到最终的势能表达式

所以相邻格点作用力为

，这些都是老生常谈的显然的事实。于是我们可以得到以下振子图像：

振子图像

根据这个图像，我们可以得到振子运动方程

接下来假设振子具有平面波解

代入方程，于是可以得到色散关系

接下来我们需要通过编程计算得出这个色散关系啦，想想都有点小激动。

二、基本方程的设定

首先我们必须理解一下这个色散关系到底是什么意思。按道理来说，这样的一个振子链中传播一个波，那么这个波的频率和两振子之间的相位之间必然满足色散关系，所以我们应该先做出一个振子链，然后让一个波在里面传播。振子链比较好构造，如下

n = 120; m = 1; tm = 300; k = 0.5; [Omega] = 1;
x[0][t_] := Sin[[Omega] t];
xvar = Array[x, n];
x0 = Table[x[i][0] == 0, {i, n}];
dx0 = Table[x[i]'[0] == 0, {i, n}];
equs0 = Table[
   m x[i]''[t] == 
    k ((x[i - 1][t] - x[i][t]) + (x[i + 1][t] - x[i][t]) UnitStep[n - i - 1]), {i, n}];
equs = Join[equs0, x0, dx0];

这里取了120个振子，相关参数如代码所示，这里要注意的是为了产生振源，我令

作为驱动力，给出一个频率为

的驱动，这样就可以产生一个传播的振动，然后解这个微分方程组

s = NDSolve[equs, xvar, {t, 0, tm}, MaxSteps -> [Infinity]];
xvar = xvar /. s[[1]];

就可以解方程了，首先看一下运行结果：

第二个振子的运动方程

横坐标是t，纵坐标是振幅，首先振子经历了一个短暂的驱动过程，然后就进入了频率恒定的周期运动，我们知道这里的周期运动图像是正弦函数，在图上看确实如此，然而凑近了看会发现相邻的峰还是有些许的不同，主要原因是后面的振子还需要启动，所以需要一定的能量传输。

下面画了两个gif图，展示了不同振子的运动图像，然后会发现行波确实是存在的，而且到后面会发生反射，毕竟振子是有界的。

不同振子的运动图像

频率较低的运动图像

所以到这里为止，我们的工作还是比较顺利的。

三、两个物理量的提取

弄好了运动方程之后，我们就可以想办法弄出来相位差和频率了。首先我们看一下频率，有人可能会以为频率与驱动力频率是相同的，事实上当频率过快的时候系统运动频率不可能那么快，固体物理里面有一个截止的频率（在色散关系里面也可以看出来），而且除了单振子链之外其他很多情况都需要自己找出振动频率。

利用快速傅里叶变换，我们可以找出振子振动的频率：

sample = 2000;
signal = Table[xvar[[2]][t], {t, 40, 320, 280/sample}];
xps = Abs[Fourier[signal]]^2;
len = Length[xps];
xps = Table[{(j - 1)/280 2 [Pi], xps[[j]]}, {j, len}];
y = Interpolation[xps];

在频率比较小的情况下，它与驱动频率类似。

然后是相位差的提取，如果直接在图像上面截取点，然后用三角函数拟合，会发现各种各样的问题，然而用如下的办法，可以轻松得到相位差。

我们取相邻的两个振子，将它们的运动方程相加，于是有

这样只需要取测量和函数的振幅就可以了。

代码如下：

F[t_] := xvar[[1]][t] + xvar[[2]][t];
[Phi] = t /. NMaximize[{Abs[F[t]], 40 < t < 100}, t][[2]];
[Phi] = 2 ArcCos[Abs[F[[Phi]]]/2];
AppendTo[T, {[Phi], [Omega]}];

四、结果

接下来只需要写一个循环，不断输出即可。整体如下：

n = 150; m = 1; tm = 350; k = 0.5; T = {};

Do[
 x[0][t_] := Sin[[Omega] t];
 xvar = Array[x, n];
 
 x0 = Table[x[i][0] == 0, {i, n}];
 dx0 = Table[x[i]'[0] == 0, {i, n}];
 equs0 = Table[
   m x[i]''[t] == 
    k ((x[i - 1][t] - x[i][t]) + (x[i + 1][t] - x[i][t]) UnitStep[
         n - i - 1]), {i, n}];
 equs = Join[equs0, x0, dx0];
 
 s = NDSolve[equs, xvar, {t, 0, tm}, MaxSteps -> [Infinity]];
 xvar = xvar /. s[[1]];
 
 sample = 2000;
 signal = Table[xvar[[1]][t], {t, 40, 320, 280/sample}];
 xps = Abs[Fourier[signal]]^2;
 len = Length[xps];
 xps = Table[{(j - 1)/280 2 [Pi], xps[[j]]}, {j, len}];
 y = Interpolation[xps];
 [Omega]0 = t /. NMaximize[{Abs[y[t]], 0 < t < [Pi]}, t][[2]];
 
 F[t_] := xvar[[1]][t] + xvar[[2]][t];
 [Phi] = t /. NMaximize[{Abs[F[t]], 40 < t < 100}, t][[2]];
 [Phi] = 2 ArcCos[Abs[F[[Phi]]]/2];
 AppendTo[T, {[Phi], [Omega]}];
 
 Clear[s, xvar, i, y, T1, T2, T3, [Phi]], {[Omega], 0.1, 1.4, 0.01}]

输出结果就是我的封面了

输出色散关系

这里的线当然是用理论计算得到的了，到这里我们基本上是搞定了这个色散关系的作图，写下来的话我发现还是挺少的，不知道为什么自己搞了挺久的才搞明白这些。

直接用驱动力频率

直接用驱动力频率，其实效果差不多，到后面基本上就是驻波了。