(32) 膜振動 - 読み書きプログラミング

円形の薄膜の振動を考えます。
極座標を使うと、運動方程式は以下になります。

$\frac{1}{v^2}\frac{\partial^2z}{\partial t^2}-\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial}{\partial r})z-\frac{1}{r^2}\frac{\partial^2z}{\partial \theta^2}=0$

ここでvは音速です。

この方程式は、時間と偏角に関して線形なので、変数分離して周波数領域で解くことができます。すなわち、

$z(r,\theta,t)=J_n(\frac{\omega r}{v})e^{in\theta}$

として一般性を失いません。
運動方程式の解となる $J_n(x)$ はBessel函数として知られています。

ドラムのように、半径aで振幅が0となる境界条件を考えます。例えば、 $J_2(x)$ の零点を探すのに、グラフを書いてみましょう。

/* 6.1.16m */
plot2d(bessel_j(2, x), [x, 0, 30]);

原点の次に5前後に最初の零点があることがわかりましたので、その値を数値的に計算しましょう。

/* 6.1.17m */
find_root(bessel_j(2, x), x, 5, 7);

$5.135622301840683$

使用した機能

find_root

求めた値を使って、ドラムのn=2の振動をアニメーション化してみましょう。

/* 6.1.18m */
load(draw);
z : 5.135622301840683^2;
for t : 0 thru 6 do
draw3d(enhanced3d=true, surface_hide=true,
  explicit(-bessel_j(2, sqrt(x^2+y^2))*x*y/(x^2+y^2) * cos(t) * (1+signum(z-x^2-y^2))-0.1*(1-signum(z-x^2-y^2)), x, -6.5, 6.5, y, -6.5, 6.5), zrange=[-0.6, 0.6]);

使用した機能

参考文献

Crandall, Mathematica―理工系ツールとしての (アジソンウェスレイ・トッパン情報科学シリーズ) p.145-p.149