波動方程式
波動方程式は、なんらかの物理的なパターン（波動）がどのように空間を伝わっていくのかを記述した式です。３次元の場合、以下のように書かれます。
空気中を伝わる音、水面の波、弦の振動など、一見性質の違う別な物理に見えるようなものでも、要素を抽出していくと、根っこの部分でこの式を見いだすことができます。
\frac { 1 }{ { s }^{ 2 } } \frac { { \partial }^{ 2 }u }{ { \partial }^{ 2 }t } =\frac { { \partial }^{ 2 }u }{ { \partial }^{ 2 }x } +\frac { { \partial }^{ 2 }u }{ { \partial }^{ 2 }y } +\frac { { \partial }^{ 2 }u }{ { \partial }^{ 2 }z }
今回は身近な自然現象として水面の波を取り上げ、波の現象を数式で記述し、それを計算機の力を少し借りてシミュレートしてみるところまでを目標とします。
これらはすでに優秀な先人がやっていますが、それを追体験してみるのは大変面白いことだと思います。
モデル化
まずは座標系を考えます。いきなり３次元を考えると大変なので、１次元で考えることにします。
とりあえず直交座標系で、地面に対して水平な軸をx軸とします。
次に波が存在しないときの高さを0として、高さをu軸とします（なのでマイナスもありえます）。
水面には水深がありますので、適当に定数hとします。
これを図にすると以下のようになります。
ここでは波は関数u(x, t)によって記述されます。
ある時間 t のときの、x地点における波の高さ = u(x, t)
ついでに、図には書きませんでしたが、奥行きは定数 w とします。
微小空間
さて、水の運動を観察するには、大きな空間を小さな空間に細切れにするのが定石です。
なのでx方向に切り刻みましょう。x方向の微小距離を⊿xとします。
とりあえずこの微小空間の水の物理量を考えてみます。
収支式
物理量を考えるとき、基本となる考えは、収支式を考えることです。
入 - 溜 = 出
文章的に書くと、入ってきた分から、溜まる分を差し引いたものが出ていく分である、といえますね。
また、ちょっと変形して、
入 - 出 = 溜
これだと、入っていく分から、出て行く分を差し引いた分が溜まっていく。
とも解釈できます。微小空間ではつねにこのバランスを保とうとするわけです。これをもとにどんなつりあいが発生するのかを考えます。
溜まる
溜まるというのは、微小空間の物理量が増えることです。
なので微小空間の質量をまず考えます。
幅が⊿x 高さがh+u(x, t) 奥行きがwです。そしてなおかつ密度を適当にρ（ロー）として、
(h + u\left( x,t \right)) \triangle xw\rho = u\left( x,t \right)\triangle xw\rho + h\triangle xw\rho
となります。
これの増加を考える、つまり時間で微分です。微分なので、hを含む項は消滅します。
\frac { \partial u\left( x,t \right) \triangle xw\rho }{ \partial t }
これが溜です。
※いきなり偏微分を出してしまいましたが、これは時間だけをちょっとだけ進めたときの傾きを表します。
なのでほとんど通常の微分と同じです。
入、出
境界面を通して水は入ったり出たりします。つまり、水の速さに面積をかければ入ったり出たりする量が記述できそうです。
ここで、質量流速を考えます。
質量流速は、単位面積の領域を、単位時間あたりに通過する水の物理量と考えます。
数式として表現するために、質量流速もまた関数と考えることにします。
ある時間 t のときの、x地点における質量流速 = V(x, t)
ついでに単位は
kg\frac { 1 }{ { m }^{ 2 } } \frac { 1 }{ s }
となります。うっかり右辺と左辺の単位がずれないように注意しましょう。
※よく次元の記述でL^(-2)などと表記されますが個人的には割り算の書き方のほうが直観的で好きです。
ここからは便宜上、左から入ってきて、右から出ていくと考えます。
※なんでこれでいいかというと、例えば左からマイナスの物理量が入ってきて、右からマイナスの物理量が出て行くこともできるからです。これは結局のところ見る立場の解釈の問題です。
さて、高さは、思い切ってhと考えてしまいます。何事も思い切りが大事です。
すると境界面の面積はh*wとなります。
これをふまえてまず入ってくる分です。左側ですので、
\rho V\left( x,t \right) hw
次にでていく分です。右側ですので、
\rho V\left( x+\triangle x,t \right) hw
となります。
収支式Final
収支式は以下のようになるとはじめのほうで述べました。
入 - 出 = 溜
この構造に先程の式をあてはめます。
\rho V\left( x,t \right) hw - \rho V\left( x+\triangle x,t \right) hw = \frac { \partial u\left( x,t \right) \triangle xw\rho }{ \partial t }
式を変形します。
\frac { \partial u\left( x,t \right) \triangle xw\rho }{ \partial t } = \rho V\left( x,t \right) hw - \rho V\left( x+\triangle x,t \right) hw \\
\frac { \partial u\left( x,t \right) \triangle xw }{ \partial t } = V\left( x,t \right) hw - V\left( x+\triangle x,t \right) hw \\
\frac { \partial u\left( x,t \right) \triangle x }{ \partial t } = V\left( x,t \right) h - V\left( x+\triangle x,t \right) h \\
\frac { \partial u\left( x,t \right) }{ \partial t } = \frac { V\left( x,t \right) h - V\left( x+\triangle x,t \right) h} { \triangle x } \\
\frac { \partial u\left( x,t \right) }{ \partial t } = \frac { V\left( x,t \right) h - V\left( x+\triangle x,t \right) h} { \triangle x } \\
\frac { \partial u\left( x,t \right) }{ \partial t } = -\frac { V\left( x+\triangle x,t \right) h - V\left( x,t \right) h} { \triangle x } \\
\frac { \partial u\left( x,t \right) }{ \partial t } =-h\frac { \partial V\left( x,t \right) }{ \partial x }
だいぶスッキリしてきました。
ちょっと天下り的ですが、ここでtでさらに微分します。
\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } =\frac { \partial (-h\frac { \partial V\left( x,t \right) }{ \partial x } ) }{ \partial t } \\
\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } = -h\frac { \partial (\frac { \partial V\left( x,t \right) }{ \partial x } ) }{ \partial t }
右辺の∂xと∂tを交換します。（ある意味交換したいがための、tでの微分だったといえます）
\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } =-h\frac { \partial (\frac { \partial V\left( x,t \right) }{ \partial t } ) }{ \partial x }
これができる理由を考えると
\frac { \partial (\frac { \partial V\left( x,t \right) }{ \partial t } ) }{ \partial x } = \\
\frac { \frac { V\left( x+\triangle x,t+\triangle t \right) -V\left( x,t+\triangle t \right) }{ \triangle x } -\frac { V\left( x+\triangle x,t \right) -V\left( x,t \right) }{ \triangle x } }{ \triangle t }
よくよく見ると、ただの割り算なので、
この同じ色の部分を交換できます。
\frac { \frac { V\left( x+\triangle x,t+\triangle t \right) -V\left( x+\triangle x,t \right) }{ \triangle t } -\frac { V\left( x,t+\triangle t \right) -V\left( x,t \right) }{ \triangle t } }{ \triangle x } = \\
\frac { \partial (\frac { \partial V\left( x,t \right) }{ \partial t } ) }{ \partial x }
かくして、
\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } =-h\frac { \partial (\frac { \partial V\left( x,t \right) }{ \partial t } ) }{ \partial x }
が手に入りました。
右辺のV(x, t)の時間微分は何か？
これは結局のところ、質量流速の変化率なのですから、加速度と解釈できます。
なので運動方程式から考えることで糸口がつかめます。
F = ma
まずは動かそうとする対象は、今まで考えていた微小空間です。
したがって質量は、
m =\rho w(h + u(x,t))\triangle x
加速度は、
a = \frac { \partial V\left( x,t \right) }{ \partial t }
問題は外力です。
これはちょっと図をズームしてみましょう。波の上のほうです。
なにやら右肩に波の大きさの差分だけブロックが積み上がっていますね。
両肩が同じなら力は釣り合った状態であると考えると、余分な青いブロックの分だけ力が加わっていると考えることができそうです。
青ブロックの質量を考えますと、⊿yというのは、
\triangle x \frac { \partial u\left( x,t \right) }{ \partial x }
ですので、
青ブロックの質量 = \triangle x \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } w\rho
これが生み出す力は、F=maなので、
青ブロックの力 = -g\triangle x \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } w\rho
ここで青ブロックの圧力を考えます。
青ブロック圧力は下方向にかかっていますので、これを底面積で割ればいいですね
青ブロック圧力 = \frac { -g\triangle x \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } w\rho }{ w\triangle x } \\
青ブロック圧力 = -g \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } \rho
これが横方向の面積に対して力が均一にかかっていますので、
力としては、
F=-w(h + u(x,t))g \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } \rho
となります。
これにより、運動方程式が完成します。
\rho w(h + u(x,t))\triangle x \frac { \partial V\left( x,t \right) }{ \partial t } = -w(h + u(x,t))g \triangle x\frac { \partial u\left( x,t \right) }{ \partial x } \rho
打ち消し合う項を削除すると、
\frac { \partial V\left( x,t \right) }{ \partial t } =-g\frac { \partial u\left( x,t \right) }{ \partial x }
が得られます。
波動方程式
一つ前パートにて、V(x, t)の時間微分の正体が明らかになりました。これを先程のV(x, t)を含む式に代入すると、
\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } =-h\frac { \partial (-g\frac { \partial u\left( x,t \right) }{ \partial x } ) }{ \partial x } \\
\frac { 1 }{ gh } \frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }t } =\frac { { \partial }^{ 2 }u\left( x,t \right) }{ { \partial }^{ 2 }x }
と波動方程式のよく見る形式にまで持ってくることができました。
ところで、このgh、波動方程式の位相速度と対応していますので、
gh=c^{ 2 } \\ c=\sqrt { gh }
ようは水深hが深いほど、波の進行速度が早いということが言えます。
だんだん数式が重くなってきましたので、数値シミュレーションは別記事にしたいと思います。
参考資料
この記事を書くにあたって、以下の資料を大変参考に（ パクり ）させていただきました。
先人の知恵に感謝申し上げます。