計算物理 WS 太鼓の膜の振動

2 次元波動方程式を数値的に解く

1

今回は太鼓の膜の振動の問題（波動方程式）を数値的に C 言語を使って計算する。

1.1

波動方程式の差分化

以下の 2 次元波動方程式を与えられた境界条件で解くことを考えよう。
2
∂ 2u
∂ 2u
2 ∂ u
=c+
∂t2
∂x2 ∂y 2!

(1)

ここで u は太鼓の膜の変移であり、c は膜を振動が伝わる速さになり、膜の張力によって決まる量である。ラプラス方程式と同様に差分化を行う。前回のラプラス方程式のときと同様に、2 回微分を差分で評価すると、∂ 2u 2∆x + · · ·∂x2∆x の 4 次以上を無視すると、2 階微分が格子点の u の値から評価でき、
u(xi − ∆x, yj ) + u(xi + ∆x, yj ) = 2u(xi , yj ) +(2)∂ 2u1={u(xi − ∆x, yj ) + u(xi + ∆x, yj ) − 2u(xi , yj )}2∂x∆x2(3)1∂ 2u={u(xi , yj − ∆y) + u(xi , yj + ∆y) − 2u(xi , yj )}2∂y∆y 2(4)y 方向では、したがって右辺は ∆x = ∆y = ∆ とおけば、R.H.S. =c2(ui+1 j + ui−1 j + ui j+1 + ui j−1 − 4ui j )∆2(5)時間に関するグリッドの添え字を右肩に括弧つきで書くことにすると、左辺はL.H.S. =1  (k+1)
(k−1)
(k)
u+u−2uijijij∆t2(6)右辺の u は k で評価することにすると、

1  (k+1)
c2  (k)
(k−1)
(k)
(k)
(k)
(k)
(k)
u+
u−
2u
=u+u+u+u−
4ui ji ji−1 j
i j+1i j−1i j∆t2 i j
∆2 i+1 j(k+1)

これを ui jについて解くと、

(k+1)

ui j

(7)

(k)

(k−1)

= 2ui j − ui j+
c2 ∆t2  (k)
(k)
(k)
(k)
(k)
u+u+u+u−4ui+1
ji−1ji
j+1i
j−1ij
∆2(8)

この差分式は右辺に k + 1 の時間ステップの値を含んでいないので、Laplace 方程式のときのように反復する必要は無い。図は差分の様子を表しているが、y 方向は紙面に垂直方向なので割愛してある。ここで u0 は初期の膜の変形を表す。

66t(k+1)

ui j 66I(k)ui−1 ju(t, 0, y) ≡ 0(k)(k)ui jui+1 ju(t, L, y) ≡ 0(k−1)

ui ju(0, x, y) ≡ u0(x, y)0x -図 1: 格子点1.2

境界条件・初期条件

この問題の境界条件は明らかに膜の端の変移 u がゼロになることである。また初期条件
はたとえば初期の膜の変形を u(0, x, y) = u0(x, y) と与えればよいが、それだけでは初期条件は決まらない。図の差分の方法を見てもわかるように、2 ステップ前の時間での変移の値が必要である。
これは物理的には初期の膜の速度を与えることによって決めることができる。今たとえば簡単なケースとして膜の速度がゼロ、すなわち∂u=0∂t t=0と与え、t = 0(すなわち k = 0) を中心とする差分を作ると、
(1)

(−1)

ui j − ui j∆t(1)=0(−1)となる。したがって ui j = ui j であるので、式??に代入すると最初のステップは
c2 ∆t2  (0)
(1)
(0)
(0)
(0)
(0)
(0)
ui j = u i j +u+u+u+u−4ui+1
ji−1jij+1ij−1ij22∆とかける。第 2 ステップ以降はそのまま式??を用いて積分していけばよい。

(9)それでは具体的にコードを見ていこう。
コードを理解して動かせたら• 違うモードの初期条件で振動の様子を見てみよう。
• 単一のモードであらわせない振動（境界条件を満たしつつ sin でかけない初期条件）の振動の様子を見てみよう。
• 今回のような四角い枠の境界条件では、Laplace 方程式の時と同様、変数分離することにより、解析的に解を求めることができる。これを求め、数値解と比較せよ。

1.3

提出課題

1. 正方形の四角太鼓の問題で、対角線（一本だけ）に骨を貼り付けたときの振動の様子を適当な初期条件から調べよ。
2. 丸い枠の太鼓の振動を調べよ。

1.4

参考：丸い枠の境界条件での解析解

円形の膜の振動を取り扱うには、式 (??) のような xy のカーテシアン座標はあまり適当でない。円形の膜の中心を原点とする極座標 ρ − ϕ に座標を変換すると、右辺のラプラシアンの形が変わって、!

2
∂ 2u
1 ∂u
1 ∂ 2u2 ∂ u=c++(= c2 ∇2 u)∂t2∂ρ2 ρ ∂ρ ρ2 ∂ϕ2(10)

図 2 のような丸い境界条件の問題を考える。物理的には半径 ρ = ρ0 で膜を枠に貼り付け、動かなくした太鼓の問題である。

y

変数分離から解を求めることを考える。すなわち、

6ρu(t, ρ, ϕ) = T (t)R(ρ)Φ(ϕ)とする。これを元の波動方程式 (??) に代入して整理すると、

1 ∂ 2T∇2 (RΦ)=c2 T ∂t2RΦu(t, ρ0 , ϕ) ≡ 0

ρ0ϕ

ρ0(11)

7/ρ0-xρ0

図 2: 境界条件
するとこの左辺は t のみの関数、右辺は ρ, ϕのみの関数となるので、結局ある定数 α に等しくならざるを得ない。よって1 ∂ 2T∇2 (RΦ)=α=c2 T ∂t2RΦ

(12)

さらに RΦ の部分を分解する。上の式の中辺=右辺を変形すると、

ρ2 ∂ 2 R
ρ ∂R
1 ∂ 2Φ
2+−αρ=−(= β)R ∂ρ2R ∂ρΦ ∂ϕ2(13)再び左辺は ρ のみの式、右辺は ϕ のみの関数なので、両辺ともにある定数 β に等しい。
この β が今 β < −a2 (a > 0) とする (要するに負の数) と、Φ の解は

Φ = Aeaϕ + Be−aϕ(14)であるが、この解は ϕ = 0 と ϕ = 2π で連続かつ滑らかにつなぐことができない。よって不適当な解である（連続性を仮定すると微分が連続となる条件を加える自由度がない）。したがって β は正またはゼロであるとしてよい。そこで以下 β = a2 (a ≥ 0) とおくと、

Φ = A sin aϕ + B cos aϕ(15)Φ = A cos a(ϕ + ϕ0 )(16)であるが、三角関数を合成してとしてもよい。ただし A は新たに定義しなおしてある。またここで ϕ0 は座標をどこから取るかを決めているだけなので、物理的には興味のない値である。よって ϕ0 = 0 として話を進める。さらにこの関数が ϕ = 0 と ϕ = 2π で滑らかにつながるには、周期 2π を持てばよいので、a は整数でなければならない。よって a を新たに n と書き直してやると、Φ はとかける。
すると R の方の式は、

Φ = A cos nϕ(17)∂ 2 R 1 ∂Rn2+−(α+)R = 0∂ρ2ρ ∂ρρ2(18)となる。
α = −k −2 (< 0) の時には、kρ = x と書き直せばd2 R 1 dRn2+(1−+)R = 0dx2x dxx2(19)と変形でき、α = −k −2 (> 0) の時には、やはり kρ = x と書き直せばd2 R 1 dRn2−(1++)R = 0dx2x dxx2(20)と書くことができる。(??) を Bessel の微分方程式とよぶ。また (??) を変形 Bessel の微分方程式とよぶ。前者は Bessel 関数 Jn と Neumann 関数 Nn として知られる 2 つの独立な解があることが知られている。これらがどのような関数であるかは Mathematica を用いてプロットさせてみるよい。いずれにせよ Bessel の微分方程式の一般解は Jn と Nn の線形結合で表されることを覚えておこう。また変形 Bessel の微分方程式も 2 つの変形ベッセル関数

In と Kn を独立な解として持つ。 → Mathematica で Bessel 関数をプロット
一方、α = 0 の時にはd2 R 1 dR n2+− 2R = 0dρ2ρ dρρ(21)となり、独立な 2 つの解 R ∝ ρn , ρ−n を持つ。
したがって R は α の符号に応じて、

R = AJn (kρ) + BNn (kρ)(α < 0)(22)

R = AIn (kρ) + BKn (kρ)(α > 0)(23)

R = Aρn + Bρ−n(α = 0)(24)

これらが境界条件 u(t, ρ0 , ϕ) ≡ 0 に適合するかどうか調べてみよう。まず簡単な場合で、
α = 0 の場合、原点で振幅が発散しないためには B = 0 でなければならない。また、有限の半径 ρ0 で R = 0 にするためには A = 0 でもなければならず、結局 α = 0 のケースは意味のない解となる。同様に α > 0 の場合にも、Kn が原点付近で発散するので B = 0 となり、In が単調増加の関数であるのでやはり ρ0 で R = 0 にするためには A = 0 でもなければならない。よって α > 0 も除外してよい。α < 0 の場合にはまずやはり原点で発散しないという条件から B = 0 がいえるが、Jn は有限の ρ に対してゼロ点を持つ（Mathematicaのプロット参照）ので、境界条件を満たすことができる。したがって

R = AJn (kρ)(25)であるが、境界条件を満たすために、k は、kmn ρ0 = ωmn(26)である必要がある。ただしここで ωnm は、Jn (x) の m 番目のゼロ点を与える x の値と定義する (Jn のプロットと x 軸の交点を与える x 座標)。
よって、R は!
ρR = AJn ωmn(27)ρ0の形となる。
最後に T (t) は、

d2 Tc ωmn=−2dt
ρ0!2T(28)なので、解は単振動と同じで、
!

c ωmn
c ωmn
T = Amn sint + Bmn costρ0ρ0!

(29)となる。これらの解 T, R, Φ を掛け合わせた u = T RΦ をすべての mn に関して足し合わせたものが一般的な解の表式となる。したがって与えられた境界条件を満たす一般解は

u(t, ρ, ϕ) =∞
Xm,n!

c ωmn
c ωmn
Amn sint + Bmn cost
ρ0
ρ0!!

!

Jnρωmncos nϕρ0(30)

未定の係数 Amn , Bmn は本来振動の初期条件（ある時刻で膜がどのように変形されるか、どのような速度で運動しているか）によって決まるべき量である。たとえば太鼓の場合には太鼓のどこをどのくらいの強さでたたくか、ということが初期条件を与える。