Mathematica-绘制瑞利波

发表于 更新于

数学模型

瑞利 (Lord Rayleigh) 研究了弹性力学线性化方程的如下形式解答： \[ \begin{equation} \begin{split} u &= A e^{-by} exp[ik(x-ct)] \\ v &= B e^{-by} exp[ik(x-ct)] \\ w &= 0 \end{split} \end{equation} \]

\(u, \; v, \; w\) 表示地球内质点的位移 假设波产生于地球内部，地球表面自由：即作用于地表处的应力矢量为零。在考察运动方程和边界条件后，瑞利找到了常数 A, B, b 和 c，并得到以下解： \[ \begin{split} u &= A (e^{-0.8475ky}-0.5773e^{-3933ky}) \cos k(x-c_Rt) \\ v &= B (-0.8475e^{-0.8475ky}+1.46793e^{-3933ky}) \sin k(x-c_Rt) \\ w &= 0 \end{split} \] 常数 \(c_R\) 称为瑞利波速。

对应于数学模型的波形长什么样子？

初步观察：\(u, \; v\) 分别为 \(\cos,\; \sin\) ，初步判断出质点的运动轨迹为一个椭圆： \[ \frac{u^2}{a^2} + \frac{v^2}{b^2} = 1 \]

Mathematica go！

输入数学模型

1
2
3
4
5
6
c1 = -0.8475;
c2 = -0.5773;
c3 = -0.3933;
c4 = 1.4679;
u[x_, y_, t_] = ( Exp[c1 y] + c2 Exp[c3 y]) Cos[x - t];
v[x_, y_, t_] = ( c1 Exp[c1 y] + c4 Exp[c3 y]) Sin[x - t];

为便于做定性演示，\(k,\; c_R\) 的值取为1。

坐标原点的运动

循序渐进，先来看看一个点的运动（当然预期的是一个椭圆）

简单地画个图

1
ParametricPlot[{u[x, y, t], v[x, y, t]} /. {x -> 0, y -> 0}, {t, 0, 2 Pi}, PlotRange -> {{-1, 1}, {-1, 1}}]

这是初始时刻位于原点的质点的运动轨迹；图中 y 轴与题中 y 轴方向相反，最后再进行处理

加个函数动起来

1
Manipulate[ParametricPlot[{u[x, y, t], v[x, y, t]} /. {x -> 0, y -> 0}, {t, 0, a}, PlotRange -> {{-1, 1}, {-1, 1}}], {a, 0.001, 2 Pi}]

所有质点的运动轨迹

使用 Table 函数列出要观察的质点，(t = 0 时刻)

1
ListPlot@Table[{x + u[x, y, 0], y + v[x, y, 0]}, {x, -4, 4, 0.5}, {y, 0, 5, 0.5}]

观察所有质点随时间的运动，把 u[x, y, 0] 和 v[x, y, 0] 中的参数改为 [x, y, t]，动起来

1
2
3
4
Manipulate[
 ListPlot[Table[{x + u[x, y, t], y + v[x, y, t]}, {x, -4, 4, 0.5}, {y,
      0, 5, 0.5}] /. {t -> a}, PlotRange -> {{-5, 5}, {-1, 6}}], {a, 
  Pi/20, 2 Pi}]

在图中，接近地表的质点顺时针运动，远离地表的质点逆时针运动

更进一步地，在动画中显示完整的轨迹 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
fig1 = Table[
   ParametricPlot[{x + u[x, y, t], y + v[x, y, t]}, {t, 0, 2 Pi}, 
    PlotRange -> {{-5, 5}, {-1, 6}}, Axes -> False, 
    PlotStyle -> {Black, Thin}], {x, -4, 4, 1}, {y, 0, 5, 1}] // Show;
Manipulate[
 fig2 = ListPlot[
   Table[{x + u[x, y, t], y + v[x, y, t]}, {y, 0, 5, 1}, {x, -4, 4, 
      0.2}] /. {t -> a}, PlotRange -> {{-5, 5}, {-1, 6}}, 
   Axes -> False, PlotStyle -> {{Red, Thin, Dashed}}, Joined -> True];
 fig3 = ListPlot[
   Table[{x + u[x, y, t], y + v[x, y, t]}, {y, 0, 5, 1}, {x, -4, 4, 
      1}] /. {t -> a}, PlotRange -> {{-5, 5}, {-1, 6}}, Axes -> False,
    PlotStyle -> {{Red, Thin}}]; 
 Show[fig1, fig2, fig3] // Rotate[#, 180 Degree] &, {a, Pi/10., 2 Pi, 
  Pi/20}]

运动轨迹分析

从哪一个深度开始，质点椭圆运动方向开始发生变化呢？观察模型，可以发现，在初始时刻，\(\cos, \sin\)的系数同号；当两个系数异号，质点椭圆运动方向相反。 绘图估计一下临界点在哪个位置，进行求解

1
2
3
4
Plot[( Exp[c1 y] + c2 Exp[c3 y]) ( c1 Exp[c1 y] + c4 Exp[c3 y]), {y, 
  0, 2}]
Solve[{( Exp[c1 y] + c2 Exp[c3 y]) ( c1 Exp[c1 y] + c4 Exp[c3 y]) == 
   0, 0 < y < 2}, y, Reals]