0%

Mathematica-绘制瑞利波

数学模型

瑞利 (Lord Rayleigh) 研究了弹性力学线性化方程的如下形式解答:

$u, \; v, \; w$ 表示地球内质点的位移
假设波产生于地球内部,地球表面自由:即作用于地表处的应力矢量为零。在考察运动方程和边界条件后,瑞利找到了常数 A, B, b 和 c,并得到以下解:

常数 $c_R$ 称为瑞利波速。

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

初步观察:$u, \; v$ 分别为 $\cos,\; \sin$ ,初步判断出质点的运动轨迹为一个椭圆:

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]