Java Appletsによる振動と波動の可視化

全文

(1)Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 1. はじめに. Java Applets による振動と波動の可視化. 理工系学部の物理教育において，講義や教科書での教材をより効果的なものにする為に，双方向性を持たせることや，JavaGraphics1)–4) を取り入れて内容を可視化5)–10) すること. 菊. 池光子†1. 泉本利章†1. などが行われている．ここで注目される Java 言語の運用上の特長は，プログラム実行環境に依存しない，つまり環境ごとにプログラムを再コンパイルする必要が無いということである．ユーザは PC の OS が異なっても，標準設定のブラウザーを使ってインターネットの. 理系基礎としての古典力学の中から連続体の振動と波動をとりあげ，基準振動，固有振動モードとともに差分による数値解法にも焦点をあてた教材として Java Applets を開発した．具体的に，3 質点とばねからなる連成振動の数値解法やノーマルモードとともに 1 次元の弦の振動・波動の Fourier 級数展開法や数値計算法を取り上げた. そこで，Fourier 展開の有限項近似による誤差や数値積分の不安定条件なども取りあげた. また, 波動方程式のダランベール解も取り上げるとともに，両端が固定された弦の振動も方向の異なる２つの波動の重ね合わせで表現されることを示した。2 次元膜の振動についても，同様な Fourier 級数展開や数値積分による解法を議論した. これらの Java Applets を使った可視化により，連続媒体の振動と波動の現象を統一的に波動方程式を通して理解することができるようになると期待される．. Web ページ上の Applets を呼び出し，いつでも自由に教材として利用できるのである．また，物理を学ぶ過程で基本の方程式に支配される系を扱うが，それらを組み込んだ Java. Applets を教材としてその理解に役立てることができる．さらに偏微分方程式の初期値や境界値問題についても，解析解と数値解を Applets で確認することができるという利点がある. 私達は，この報告で古典力学の振動と波動11)–19) を取り上げる．これらは量子力学の学習準備としても一般に大切であると考えられている．第２章では，Java 言語の特長をいかして，１次元の連成振動としてばねと３つの重りが交互に連なった系を扱い，その運動方程式に基づいたシミュレーションを可能にする Java Applet を制作した．その際，基準振動を求めるとともに，Euler 法や Runge-Kutta 法による差分. Simulations of Classical Vibrations and Waves in terms of the Java Applets. 方程式の数値解法をとりあげ，その時間刻みの粗さによって計算精度が変化する様子を可視化した．. Mitsuko Kikuchi. †1. and Toshiaki Izumoto. †1. 次に，第３章と第４章では，多体連成系の質点の数を無限と考えた連続体でもある弦と膜の振動を取り扱った．第３章では，弦の振動を記述する 1 次元波動方程式の解析解・数値解を求め，Java Applets を制作し，弦の垂直変位を時間と位置の関数として可視化した．初期. We have developed some Java applets to simulate the motion of three coupled oscillators, and the vibrational displacements for a string and a membrane: The motion of three coupled-oscillators is briefly studied by employing numerical methods of the Euler and the Runge-Kutta methods, which are compared with the normal mode methods. The vibration of a string is also studied by employing two methods, i.e. the Fourier expansion and the numerical integration in solving the wave equation. The convergence of the Fourier series and the numerical instability condition are studied. As demonstrated by the d’Alembert solution, the vibration of a string bounded on both ends are well described as the superposition of the waves propagating to the opposite ways. The vibrational displacement for 2-dimensional case of a square membrane is studied in the same way as in the 1-dimensional case. These applets may facilitate students to learn the basics of vibrations and waves on an equal basis of solving the wave equation.. 条件として, 弦を平衡状態から摘み静かに離す場合の他に，一直線に張った弦にゼロでない初速度を与える場合も取りあげた．Fourier 級数による解法では，基準振動の有限項数の和として解が記述される．また差分による偏微分方程式の数値解の不安定性にも言及した．波動の現象を理解するため，無限区間に対する D’Alembert 解をとりあげ有限区間における反射を含む解析解・数値解との比較も行った．第４章では，2 次元膜の振動の時間変化の記述についてとりあげ，１次元の場合と同じく，波動方程式の Fourier 展開を用いた解法とともに差分法による数値解法を採用して，シミュレーションを可能にする Java Applet を制. †1 立教大学理学部物理学科 Department of Physics, Rikkyo University. 1. c 2011 Information Processing Society of Japan °.

(2) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 2.1 連成振動の解析ここでは，図１の上段に示されているような連成振動システムを取り上げる．つまり，3 個の質点と重さが無視できる 4 個のばねが連なり，両端を壁に固定してある系を仮定する．これらの３個の質点の質量をそれぞれ m1 , m2 , m3 とし，それらの平衡の位置からのずれを. (x1 , x2 , x3 ) とする．４個のばね定数は全て等しく k であるとすると，各質点の運動を記述する運動方程式は，フックの法則より，. md2 x1 /dt2 = −kx1 + k(x2 − x1 ). = −2kx1 + kx2. md x2 /dt = −k(x2 − x1 ) + k(x3 − x2 ) = −2kx2 + kx3 2. 2. md2 x3 /dt2 = −k(x3 − x2 ) − kx3. (1). = kx2 − 2kx3. と書き下すことができる．さらに，3 個の振動子からなる連成振動の場合，解は基準振動数. ω1 , ω2 , ω3 の 3 個の単振動の合成として表わされることがわかる．基準振動の角振動数は， ω0 =. √. k/m とすると，次のように表される：. √ ω1 =. 2−. √. 2 ω0 , ω 2 =. √. √ 2 ω0 , ω 3 =. 2+. √. 2 ω0. (2). 各質点 (x1 , x2 , x3 ) の運動は，次の 3 つの基準振動の振幅を A1 , A2 , A3 とするとき，. (x1 (t), x2 (t), x3 (t)) = ( 図 1 3 個のばね振動子の連成振動をシミュレーションする Applet の一場面．上段には，３質点の運動を，中段にはそこから選択した１質点の振動を，下段には３つの基準振動を時間 t(単位は 1/ω0 ) の関数として表示している． Fig. 1 The horizontal displacements of the three coupled spring oscillators are simulated in the upper part. The displacement of a oscillator chosen among them can be displayed as the vertical displacement in the middle, which is decomposed into three normal modes in the bottom.. √. 1/2,. + ( 2/2, +(. √ √. 2/2,. 1/2) A1 cos(ω1 t + α1 ) √ 0, − 2/2) A2 cos(ω2 t + α2 ). 1/2, − 2/2,. (3). 1/2) A3 cos(ω3 t + α3 ). と表される．初期条件として，3 質点の変位を (x1 (0), x2 (0), x3 (0)) とし，その速度をゼロ（静止）と設定すると，右辺の固有振動の振幅 A1 , A2 , A3 が求められ，またその位相は,α1 = α2 = α3 = 0 となる．ここで，逆に連立微分方程式 (1) の解 (3) を用いて，固有振. 作した．Fourier 級数による基準振動による記述や差分化による数値解の不安定性もシミュ. 動を次のように表すこともできのは面白いことである:. レーションで確認できるようにした．. A1 cos(ω1 t + α1 ) = (. 2. 振動子の連成振動. √. x1 (t) +. A2 cos(ω2 t + α2 ) = ( 2x1 (t) A3 cos(ω3 t + α3 ) = (. 2 重振り子の連成振動20) もしばしば取り上げられているが，ここでは，3 個のばね振動子の. x1 (t) −. √ √. 2x2 (t) + − 2x2 (t) +. √. x3 (t))/2, 2x3 (t))/2,. (4). x3 (t))/2. 連成振動12) を取り上げる．その運動方程式を，Euler 法16),18),19) と Runge-Kutta 法16),18). また，Euler 法と Runge-Kutta 法を採用して, 微分方程式を差分化し，それを数値積分す. を用いて，数値的に解くプログラムを作成する．また，重りの運動と同時に基準モードを見. ることでも解を求めた．その解が, 正しければ, 上の（4）式の右辺に代入することによって. ることができるようにする．. 固有振動が求められることになる. これらの処方箋は, 微分方程式の教科書などにも詳しい．以下では，数値計算における時間 t を計る単位として，1/ω0 =. 2. √. m/k を採用する．. c 2011 Information Processing Society of Japan °.

(3) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 2.2 Applet による連成振動の数値シミュレーション図１に示す Applet では，初期条件として，赤色系，黄（または緑）色系，青色系の各重りが左壁（位置 x = 0 cm）から，20 cm 間隔で並べられている．これらのうち，右端の重りだけを，右方向に 10 cm だけ引っ張った後，静かに放したときの連成振動を取り扱う．この図の最上段部では，各質点とばねの運動を時間の経過に従って表示している．中段には，3 質点の変位の時間変化 (x1 (t), x2 (t), x3 (t)) のうち, １つずつ右の制御パネルで選択表示できる．最下段には，先に示した（4）式にしたがって，数値計算の解から 3 つの基準振動 ω1 , ω2 , ω3 を持った基準振動のグラフを求めて表示している．ここに表示されている基準振動成分は, 初期条件によって定まる振幅を含んでいる．数値計算の精度が悪くなっている場合は，表示されたものは基準振動から外れたものになる可能性があることに注意する．右側の制御パネルでは，運動方程式の数値解法を Euler 法または Runge-Kutta 法の内から選択できる．両者を同時に選択表示することもできる．なお，時間ステップ ∆t の可変入力によって，運動方程式を差分法によって解いた時の計算精度を，選択した質点の振動表示か. 図 2 フーリェ展開の有限項近似において，両端を固定した弦を左から x = 75 cm の位置で摘みあげたときの垂直変位とそれを構成する主要なフーリェ成分． Fig. 2 The vertical displacement for a string bounded on both sides, which is initially plucked at a position of x = 75 cm, is shown with those for a few leading terms of its Fourier expansion.. ら確認することができる．＜表示速度＞のスクロールバーや，アニメーションの開始/再開，停止/終了，条件設定用のボタンも用意されている．. 3. 弦の一次元振動と波動 3.1 両端を固定された弦の振動. Xn (x) = An sin(nπx/`) + Bn cos(nπx/`). 3.1.1 Fourier 級数による解の表示. (9). Tn (t) = Cn cos(nπct/`) + Dn sin(nπct/`). ここでは，両端が固定された長さ ` の弦を考える．弦は一様な密度 ρ と一定の張力τを. ここで，弦の両端を固定いう境界条件から，上の展開係数 Bn の値はゼロとなる.. 持ち，摩擦も重力も受けないとする．このとき，弦の垂直方向の変位を，水平方向の位置 x. まず，上の線形波動方程式は，後でその複数の解をあわせて 1 つにまとめることもできるこ. と時刻 t の関数 y(x, t) で表わす．弦の 1 次元波動方程式は，次のように表される．. 13) とに注意する．では，具体的な初期条件としていくつかのケースを考察する：. ∂ 2 y(x, t)/∂x2 = (1/c2 )∂ 2 y(x, t)/∂t2 ここで波動の伝搬速度は c =. y(x, 0) = φ(x),. √. まず，CASE 1 として, t = 0 において x = a (但し，0 ≤ a ≤ `）の位置で弦を持ち上げ，. (5). τ /ρ である．この方程式を解くための初期条件は，. ∂y(x, t)/∂t|t=0 = ψ(x). 初速度は与えないで静かに離す，つまり次の場合を考える：. (6). y(x, 0) = φ1 (x),. と書ける．また，境界条件は，次のように表される：. y(0, t) = y(`, t) ≡ 0. (10). り，展開係数 Dn = 0 となる．故に，両端が固定され，初期に静止状態にある弦の一般解. (7). 解析解は，変数分離法11)–19) を適用して求める．. y(x, t) = X(x)T (t). ∂y(x, t)/∂t|t=0 = ψ1 (x) = 0. この弦の初速度がゼロという条件から，Tn (t) を t で微分したものは t = 0 において 0 となは，基準振動モード n の線形重ね合わせとして次のように表わされる：. (8) y(x, t) =. これを波動方程式に代入して，X(x) および T (t) に対する方程式に分離して n 番目の基準. ∞ ∑ n=1. 振動を求めると，. 3. Cn sin(. nπct 2 nπx ) cos( ), Cn = ` ` `. ∫. `. φ1 (x) sin( 0. nπx )dx `. (11). c 2011 Information Processing Society of Japan °.

(4) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 図 3 両端が固定されている弦を左から x = 75 cm の位置でつまみあげ静かに離したときの弦の垂直変位 y(x, t) を，時間経過を追って 1 周期 τ = 2`/c 分表示している． Fig. 3 The vertical displacement for a string which is initially plucked at a position of x = 75 cm, is shown for the period of vibration τ = 2`/c.. 図 4 両端が固定されている弦を静止状態で左から x = 75 cm の位置を叩いて初速度を与えたときの弦の垂直変位 y を，時間経過を追って 1 周期 τ = 2`/c 分表示している Fig. 4 The vertical displacement y for an string which is initially struck at a position of x = 75 cm, is shown for the period of vibration τ = 2`/c.. {. ここで, 基準振動 n = 1 の振動数が 1 番低く，他の基準振動ではその倍数になっているの. φ1 (x) =. で、その周期 τ = 2`/c が系の振動の周期を代表することに注意する．さて，このタイプの初期条件として，具体的に次の２つのケースを考えることにする：. CASE 1a として，t = 0 において x = a (但し，0 ≤ a ≤ `）の位置で弦を持ち上げ，静. x/a,. (0 ≤ x ≤ a). (` − x)/(` − a),. (a ≤ x ≤ `). (12). この初期条件を満たすように一般解における展開係数 Cn を求めると：. Cn = 2`2 /{n2 π 2 a(` − a)} sin(nπa/`). かに離す場合を考える，つまり，. (13) 2. となる．この係数の振る舞いは，Cn /c1 = sin(nπa/`)/{n sin(πa/`)} から分かるように，位置 a の値にもよるが n について 1/n2 の依存性をもっている．もし，a = `/2 の時は，n. 4. c 2011 Information Processing Society of Japan °.

(5) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. が偶数の基準振動の発生が抑えられることに気がつく. こうして得られた弦の垂直変位 y は，図３に，位置 x の関数として 1 周期（τ = 2`/c）分プロットされている. CASE 1b として，弦上の孤立波として三角波を考える．つまり，. {. φ1 (x) =. 1 − |x − a|/∆, 0,. |x − a| ≤ ∆. (14). otherwise. という場合を考える．先と同様に Cn を求めると. Cn = 4`/(n2 π 2 ∆) sin(nπa/`) {1 − cos(nπ∆/`)}. (15). となる. この係数は，n が大きくなってゆくと 1/n2 に比例して減少する．次に，CASE 2 として，弦の初期変位がゼロの状態で，初速度分布を与える，つまり. y(x, 0) = φ2 (x) = 0,. ∂y(x, t)/∂t|t=0 = ψ2 (x). (16). を考える. t = 0 に変位がゼロである弦の一般的な解の形は，次のように表される：. y(x, t) =. ∞ ∑. Dn sin(. n=1. nπx nπct 2 ) sin( ), Dn = ` ` mπc. ∫. `. ψ2 (x) sin( 0. nπx )dx `. (17). ここで，CASE 2a として，次のような δ 関数タイプの初速度分布を考える：. ψ2 (x) = chδ(x − a),. (0 ≤ x ≤ `). 図 5 左から 70 cm の位置に初期設定された三角形の孤立波による有限長の弦上の振動と無限長の弦上の波動を，時間の経過とともに図の上段および下段に比較表示している． Fig. 5 The vertical displacements y(x,t)’s are simultaneouly shown as a function of t for strings with finite and infinite lengths, which are both initially plucked at a position of x = 70 cm forming a triagle shape.. (18). この初期条件に対応した展開係数 Dn は次のように求められる：. Dn = (2h/nπ) sin(nπa/`). (19). この係数 Dn の振る舞いは，Dn /D1 = sin(nπa/`)/{n sin(πa/`)} から分かるように，n について 1/n の依存性があり，緩やかに減少する．つまり，Furier 級数の有限項近似をする. 3.1.2 弦の振動の数値解法（差分法）. と結果の収束性が悪いと思われる．こうして得られた弦の垂直変位 y は，図４に，位置 x. 2 次元格子点を考え，その格子点上で先の境界条件と初期条件 CASE 1 を満たす波動方程式の解を探す．式 (5) において，x, t を離散的な変数に直し xi = i∆x, tj = j∆t とし，変. の関数として 1 周期（τ = 2`/c）分プロットされている. 余談ですが，ここで取り扱った初期条件 CASE 1a は，弦を引っ張って静かに離す弦楽器. 位 y は y(i, j) で表す．. y(i, j + 1) = 2y(i, j) − y(i, j − 1). 2. ハープの場合に似通っており，振動に含まれる高調波の割合が 1/n で減小することが分か. (20). +(c2 ∆t2 /∆x2 ){y(i + 1, j) − 2y(i, j) + y(i − 1, j)}. る．また初期条件 CASE 2 は，鋼鉄線をハンマーで叩く楽器ピアノの場合に似通っている．この場合，（19 ）式から，高調波の振動（n > 1）の割合がハープの場合より多く含まれ，音. この式を用いて偏微分方程式を数値的に解くことができる．しかし式の右辺を見てみると，. 12) 色が硬く感じられる．実際の楽器では，共鳴板がついていること，また空気との摩擦によ. 時間に関して 2 ステップ前の変位 y を必要としているにもかかわらず，初期条件は 1 つの. る音の減衰などの違いがある．ここでの解析において，弦の波動方程式にそのような減衰効. 時間についてしか与えられていない．そのうちの初期変位は，. 果. 16)–19). は取り入れていない．. y(i, 0) = φ(i∆x). (21). と与えられている．次に，初速度の式 ∂y(x, t)/∂t|t=0 = ψ(x) を中心差分を用いて表すと，. 5. c 2011 Information Processing Society of Japan °.

(6) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 図 7 両端を固定した弦の振動/波動の現象を進行波成分（赤色）と退行波成分（青色）の２つに分離して表示することができる． Fig. 7 The vertical displacement for a string bounded on both sides can be described in terms of two waves in red and blue, respectively traveling to the right and to the left. 図 6 静止した弦の左から 70 cm の位置を叩いて δ 関数型の初期速度分布を与えたとき，有限長の弦上に発生する振動と無限長の弦上に発生する波動を，時間の経過とともに図の上段および下段に比較表示している． Fig. 6 The vertical displacements y(x,t)’s are simultaneouly shown as a fuction of t for strings with finite and infinite lengths, which are both initially struck at a position of x = 70 cm.. 離れて隣り合う格子点に達するまでの時間より小さく設定されなければならないことを示している．制作した Applet では ∆x, ∆t の値の変更が可能であるので，計算の破綻の様子も観察することができる．. {y(x, ∆t) − y(x, −∆t)}/(2∆t) = ψ(x) となる．これから，y(i, −1) = y(i, 1) − 2ψ(i∆x)∆t. 3.2 弦を伝わる波. として良いことがわかる．よって，. 3.2.1 無限の長さの弦の波動ここでは，無限の長さの弦を考える．一次元の波動方程式に対して，初期条件が，y(x, 0) =. y(i, 1) = y(i, 0) + ψ(i∆x)∆t +(c2 ∆t2 /2∆x2 ){y(i + 1, 0) − 2y(i, 0) + y(i − 1, 0)}. φ(x) および ∂y(x, 0)/∂t = ψ(x) と与えられた時，d’Alembert の波動公式として. (22). ∫. が得られる．初期条件を指定する関数 φ(x), ψ(x) が，δ 関数のように特異なものでなく，差. x+ct. y(x, t) = {φ(x + ct) + φ(x − ct)}/2 + (1/2c). 分化に対して緩やかに変化する場合は，これで計算を行うことができる．. ψ(s)ds. (24). x−ct. これと同じような数値解法17)–19) を利用する際の注意点は，∆x, ∆t の値の関係である．数. 12)–14) が求められている．. 値不安定性の条件を満たすと，数ステップ進んだところで計算が破綻してしまう事態を引き. 初期条件として先の CASE 1b の区間 [0, `] を無限区間 [−∞, ∞] に拡大して CASE 3 と. 17)–19). 起こす．数値安定性を保証するクーラン条件. c∆t/∆x ≤ 1. すると，弦の垂直変位は，容易に,d’Alembert 解として. （Courant-Friedrichs-Lewy 条件）:. y(x, t) = {φ(x + ct) + φ(x − ct)}/2. (23). (25). と求められる．これは，先の初期条件として設定した１つの三角波が高さ半分の 2 つの三角. を満たすように決めなければならない．この条件は，計算に用いる ∆t は，波動が ∆x だけ. 6. c 2011 Information Processing Society of Japan °.

(7) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 波に分かれて，それぞれ左右に進行していく解になっている．. 図 2 は，その時の弦の垂直変位の Fourier 展開有限項近似とその初めの数項の初期値を表示. 先の CASE 2 の区間を同様に拡大して，新しい初期条件 CASE 4 とする．φ(x) =. する Applet の 1 場面を示している．右側の操作パネルでは，チェックボックスをマウスに. 0, ψ(x) = hcδ(x − a) を (24) 式に代入して，d’Alembert 解を求めると，. {. a − ct < x < a + ct. h/2,. y(x, t) =. 0,. より選択することにより波動方程式の解法として，Fourier 級数による解と差分法による解が選択できる．（これらは同時に選択することも可能．）同様にその下のチェックボックスを. (26). otherwise. 選択することにより，基準振動 1∼6 番目の時間変動を表示することもできる．＜速度表示＞のスクロールバーではアニメーションの表示スピードを調節することができる. また，「再. となる．これは，高さ h/2 のステップ関数状の波が，左右に広がる様を表している.. 開」ボタンでアニメーションの開始/再開，「停止」ボタンで一時停止，「設定」ボタンでアニ. 3.2.2 両端固定の弦の振動. メーションは，既に指定済みの条件が設定される．．. 子供の時に糸電話で遊んだ事を思い出せば，弦上を波が伝わるのは全く不思議な現象では. 図 3 には，この弦の初期条件 CASE 1a に対する垂直変位の Fourier 級数展開による解の. ない．先の 3.1.1 章で波動方程式の Fourier 級数による振動の一般解を求めたが，よく見る. 周期振動（周期 τ = 2`/c）の様子を一覧できる. これは，図 2 の Applet のアニメーション. 14),15) と，それらも実は d’Alembert 解になっていることが分かる．このことを，3 角関数の. でも確認できるものである．図 4 には, 異なる初期条件 CASE 2a に対する弦の垂直変位. 積・和の公式を用いて，簡単化のため２つの初期条件それぞれの場合に分けて示す：. の周期振動（周期 τ = 2`/c, c は弦を伝わる波動の伝搬速度とする．）がまとめられている.. まず，初期条件 CASE 1 に対する Fourier 級数解（11）式について検討すると，. 初期条件 CASE 1a, 1b に対しては，それぞれ先の図 2 および後に出てくる図 5 の Applet. y(x, t) =. ∞ ∑ n=1. =. ∞ ∑. Cn. n=1. 1 2. で，差分法による解を表示することができる．またその差分のステップサイズ c∆t, ∆x を指. nπ c x cos nπt ` `. Cn sin. { sin. 定することにより，数値計算の不安定性（クーラン条件）を確認することもできる．Fourier. nπ(x + ct) nπ(x − ct) + sin ` `. }. 級数による解は，何番目の項までの和を表示するかを，テキストフィールドに数字 n ( ただし n ≤ 20) を入力することによって指定/変更できる．図 5，図 6 に示す Applet 表示画面は，長さ ` cm あるいは長さが無限 ∞ の一様な弦に. 1 {φ1 (x + ct) + φ1 (x − ct)} (27) 2 さらに初期条件 CASE 2 に対する Fourier 級数解（17）式について検討すると， =. y(x, t) =. ∞ ∑. Dn sin. n=1 ∞. =. ∑. (. Dn −. n=1. 1 = 2c. ∫. 放す．この場合，右側のパネルでは，チェックボックスの項目を選択することにより波動. nπ nπ x sin ct ` ` 1 2. ){. ∑. x+ct ∞. x−ct. おける波の伝搬を，それぞれ，上段および下段に表示することができる．図 5 の初期条件. CASE 1b および CASE 3 では，弦の途中の位置 (x = a cm) で弦を持ち上げ，静かに. n=1. cos(. 方程式の解法として，差分法または Fourier 級数による解を選択できる．図 6 の初期条件. nπ(x + ct) nπ(x − ct) ) − cos( ) ` `. nπc nπξ 1 sin( )dξ = Dn ` ` 2c. ∫. }. CASE 2a および CASE 4 では，弦の途中を叩いて δ 関数的な初期速度分部を設定している．この場合の解法は，CASE 2a では，Fourier 級数展開による解そして，CASE 4 では d’Alembert 解が表示される．その下のチェックボックスを選択することにより，基準. x+ct. ψ2 (ξ)dξ. 振動モード 1∼6 番目を見ることもできる．. (28). x−ct. 図 7 には，初期条件 CASE 1a および CASE 2a の場合の Furier 展開有限項近似における解を確認することができる．また，進行波/退行波の分別表示をクリックすることにより，. どちらの場合も速度 c で進行する波と退行する波の和を表していることがわかった．. 3.3 Applet による振動・波動の数値シミュレーション. 両端固定の振動解を 2 つの成分に分けて表示することもできる．また「多重表示」をチェッ. ここでは，両端が固定された長さ `(= 100) cm の一様な弦を考える．初期条件として，左. クすることによりいわゆる写真における多重露出の結果を見ることができる.. 壁（x = 0）から x(< `) cm の位置で弦を 1mm 持ち上げ，静かに放す．. 7. c 2011 Information Processing Society of Japan °.

(8) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 4. 膜の 2 次元振動と波動弦と同様，正方形の膜を平衡位置から摘み，静かに放すときの波動のパターンを予測する. Applet を制作する．膜の振動の 2 次元の波動方程式に基づいて，2 重 Fourier 級数で表される解を解析的に求め，基準モードを取り出して表示できるようにする．また，3.1.2 で述べた差分法を 2 次元に拡張することにより，運動方程式を，初期条件を加えて数値積分する定式化を示す．. 4.1 Fourier 展開による膜の振動の解析ここでは，4 辺が固定された長さ ` の正方形膜を考える．膜は一定の面密度 σ と張力. τ を持ち，摩擦も重力も受けないとする．膜の安定位置からの垂直変位を，3 変数の関数 u(x, y, t) で表わすと，x, y はそれぞれ静止した膜の辺に沿って測った水平位置，t は時間である．. 2 次元の膜の振動/波動を記述する波動方程式は, 基本的に 1 次元の場合と同様である：. (. ). ∂ 2 u/∂x2 + ∂ 2 u/∂y 2 = (1/c2 )∂ 2 u(x, y, t)/∂t2. 波の伝搬速度は，c =. (29). √. τ /σ である．4 辺固定の膜の境界条件は次のように表される：. u(0, y, t) = u(x, 0, t) = u(`, y, t) = u(x, `, t) ≡ 0. (30). 図 8 １辺の長さ ` = 60 cm の正方形の膜を (x, y) = (30, 30) の位置で Gauss 型 (分散 d2 = 9 と設定) に摘み上げた場合の差分法による振動シミュレーション. この数値計算における摘み上げの位置や x, y 方向の差分ステップの大きさは可変入力である． Fig. 8 The simulation of the vertical displacement of a membrane of ` cm ×` cm, initially plucked at a position (a, b) forming a Gaussian shape with the variance d2 = 9. Here the length ` = 30 cm and the position (a, b) = (`/2, `/2) are assumed.. 初期条件として，膜の垂直初期変位と, その初速度がある：. u(x, y, 0) = φ(x, y),. ∂u(x, y, t)/∂t|t=0 = ψ(x, y). (31). 膜を摘み上げ静かに離す初期条件を，CASE 5 として考える：. [. ]. φ(x, y) = exp −{(x − a)2 + (y − b)2 }/(2d2 ) ,. ψ(x, y) = 0. (32). この場合，波動方程式の解は 2 重 Fourier 級数によって次のように表される：. ∑∑ ∞. u(x, y, t) =. ∞. Cm,n sin(mπx/`) sin(nπy/`) cos. m=1 n=1 2 2. [. ]. [√. いないことにも気がつく．. ]. m2 + n2 πct/` ,. Cm,n = (8πd /` ) exp −π 2 d2 (m2 + n2 )/(2`2 ) sin(mπa/`) sin(nπb/`). また，参考までにここでは数値解を与えないが，膜を叩いて初速度を与える初期条件を. (33). CASE 6 として考えると： φ(x, y) = 0,. (34). ψ(x, y) = hc δ(x − a)δ(y − b). (35). この場合，波動方程式の解は 2 重 Fourier 級数によって次のように表される：. この解は，重み関数 Cm,n で，基準振動モード (m, n) を線形に重ね合わせたものになって. u(x, y, t) =. いる．初期条件において，膜の摘み上げ幅 d が狭い場合は，より大きな m > 1 または n > 1. ∞ ∞ ∑ ∑. √. Dm,n sin(mπx/`) sin(nπy/`) sin. [√. ] m2 + n2 πct/` ,. (36). m=1 n=1. をもった高調波の基準振動モード (m, n) の振幅も少なからず寄与することが分かる. また，. Dm,n = (2h/. 弦の場合のように，振動が m = n = 1 の時の最低振動数の整数倍の組み合わせにはなって. 8. m2 + n2 π`) sin(mπa/`) sin(nπb/`). (37). c 2011 Information Processing Society of Japan °.

(9) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 図 9 １辺の長さ ` = 30cm の正方形の膜を (x, y) = (15, 15) の位置で Gauss 型 (分散 d2 = 9 と設定) に摘み上げた場合の 2 重フーリェ級数展開法による振動シミュレーション. この数値計算における摘み上げの位置は可変入力である. Fig. 9 The simulation of the vertical displacement for a membrane, initially plucked at a position (x, y) = (15, 15) forming a Gaussian shape with the variance d2 = 9.. 図 10 正方形膜の振動計算におけるフーリェ展開成分の１つ，(m, n) = (3, 3) を表示 Fig. 10 A vibrational mode (m, n) = (3, 3) of the square membrane.. u(i, j, 0) = φ(i∆x, j∆y). (39). と与えられている．また，膜の垂直変位の初速度に着目して，それを中心差分を用いて表すと，{u(x, y, ∆t) − u(x, y, −∆t)}/(2∆t) = ψ(x, y) となる．つまり，u(i, j, −1) =. u(i, j, 1) − 2ψ(i∆x, j∆y)∆t として良いことがわかる．これから，. 4.2 波動方程式の数値解法先の差分による数値積分法を 2 次元に拡張する. 格子点の座標 x = i∆x, y = j∆y, t =. u(i, j, 1) = u(i, j, 0) + ψ(i∆x, j∆y). k∆t を格子点番号 i, j, k を用いて表現すると，変位 u(i, j, k) に関する差分方程式は次のよ. + (c2 ∆t2 /2∆x2 ) {u(i + 1, j, 0) − 2u(i, j, 00) + u(i − 1, j, 0)}. うに表される：. + (c2 ∆t2 /2∆y 2 ) {u(i, j + 1, 0) − 2u(i, j, 0) + u(i, j − 1, 0)}. u(i, j, k + 1) = 2u(i, j, k) − u(i, j, k − 1). が得られる．初期条件を指定する関数 φ(x, y), ψ(x, y) が，δ 関数のように特異なものでな. + (c2 ∆t2 /∆x2 ) {u(i + 1, j, k) − 2u(i, j, k) + u(i − 1, j, k)} + (c ∆t /∆y ) {u(i, j + 1, k) − 2u(i, j, k) + u(i, j − 1, k)} 2. 2. 2. (40). く，差分化に対して緩やかに変化する場合は，これで計算を行うことができる．. (38). しかし式の右辺では時間に関して 2 ステップ前の変位 y を必要としているが，初期条件は. 4.3 Applet による膜の振動シミュレーション. 1 つの時間についてしか与えられていない．初期条件（31）式の内，膜の初期変位は，. ここでは，4 辺が固定された正方形膜を考える．膜は一定の面密度σと張力τを持ち，. 9. c 2011 Information Processing Society of Japan °.

(10) Vol.2011-CE-109 No.12 2011/3/19. 情報処理学会研究報告 IPSJ SIG Technical Report. 摩擦も重力も受けないとする．初期条件として，中央で膜を 1mm 持ち上げ，静かに放す．. 2) C.S. Lindsey, J.S. Tolliver, and T. Lindblad, JavaTech: An Introduction to Scientific and Technical Computing with Java, Cambridge University Press, 2005. 3) 峯村吉泰:Java によるコンピュータグラフィックス（2003, 森北出版株式会社）. 4) 中山茂，『Java2 グラフィックスプログラミング入門』，技報堂出版，1999 年. 5) 安藤高志, 高田壮起, 武田裕太, 津金澤雅人, 畑中涼, 泉本利章：JAVA を用いた物理現象の数値シミュレーションと可視化，日本物理学会講演概要集（2009）28pYE-3. 6) 泉本利章，畑中涼, 武田裕太, 柳町朋樹：Java Graphics を用いた光波の干渉実験のシミュレーション，情報処理学会研究報告 Vol.2009-CE-99, No.2,1-8, 2009. 7) 泉本利章，色川由季：Java Applets による古典力学のシミュレーションとその可視化，情報処理学会研究報告 Vol.2011-CE-108, No.20,1-9, 2011. 8) W. Fendt, http://www.walter-fendt.de/ph14e/． 9) PhET, Interactive Simulations, univ. of Colorado at Boulder, http://phet.colorado.edu/． 10) 覧具博義, 能動的学習を支援する PheET シミュレーション教材, 大学の物理教育 16 （2010）34-37． 11) 大槻義彦，『セメスター物理波動』，学術図書出版社，1998 年. 12) 吉岡大二郎，『振動と波動』，東京大学出版会，2005 年. 13) 竹之内脩，『微分方程式とその応用』，サイエンス社，2004 年. 14) 須藤彰三，『波動方程式の解き方』，共立出版，1994 年. 15) 登坂宣好・大西和榮，『偏微分方程式の数値シミュレーション』，東京大学出版会，2003 年. 16) 小柳義夫，『計算物理学基礎編』，朝倉書店，2001 年. 17) 小柳義夫，『計算物理学応用編』，朝倉書店，2001 年. 18) R.H. Landau and M.J.P. Mejia,”Computational Physics - Problem Solving with Computers -” (John Wiley & Sons, Inc, 1997). 19) R.H. Landau, M.J. Paez and C.C. Bordeianu,”Computational Physics - Problem Solving with Computers -” (Wiley-Vch, 2007). 20) 土井正男，滝本淳一編，「物理仮想実験室」（名古屋大学出版会，2004）. 図 8 に示す Applet の中央画面では，` = 60cm，∆x = 2cm の場合の膜の振動を差分法で求めて表示している．右側のコントロールパネルでは，チェックボックスの項目選択により波動方程式の解法が選択できる．（差分法または 2 重 Fourier 級数による解．同時に表示することも可能．）差分法での解を表示する際，クーラン条件を確認しながら，数値不安定性に陥らないように c∆t, ∆x の値を選択すればよい．図 9 では，` = 30cm の場合の 2 重. Fourier 級数による解を m ≤ 7, n ≤ 7 の条件を満たす基準振動 (m, n) の和として表示するが，テキストフィールドに別の数字を入力することによって，m と n の値を変更することができる．（ただし m, n ≤ 30．）図 10 では，その下のチェックボックスから特定の基準モード (m, n) = (3, 3) を選択表示している．アニメーションの「開始」，「停止」，「設定」ボタンや表示速度調節も，前と同様である．. 5. まとめこの報告では，３質点の連成振動の運動，および連続体としての弦の振動と２次元膜の振動を，それらを支配する運動方程式や２階偏微分波動方程式に基づいて，その解析解・数値解を求め，その振る舞いをシミュレーションで可視化する作業を行った．まとめとして，ここで制作された Applets 教材の有用性について述べてみる：まず，Java 言語を用いて Applets を制作する際に，多くの学習者に馴染みのあるマイクロソフトのウィンドウ部品（AWT）を利用しているので，完成した作品は簡単操作で動かす事ができる．また，これらの Applets をこの報告と合わせて利用することにより，連成振動とともに，１次元波動方程式の Fourier 級数による解や数値積分による解，そして d’Alembert 解の意味するところを物理的に理解するだけでなくシミュレーションにより視覚的にとらえることも可能である．２次元膜の振動も同様に，Fourier 展開法による基準振動（固有モード）の重ね合わせや差分による数値計算法によりシミュレーションが可能である．今回扱った運動方程式や波動方程式以外の偏微分方程式の初期値や境界値問題においても，こうして経験した振動・波動のプログラムの一部分を変更するだけで，解析解と数値解の比較をし，視覚的に確認することができるようになることが期待される．謝辞貴重な議論をして頂いた物理教室の皆様に，謹んで感謝の意を表します．. 参. 考. 文. 献. 1) Oracle Corporation, Java Platform, Standard Edition 6 JDK.. 10. c 2011 Information Processing Society of Japan °.

(11)