第7回

ここでは，１次元波動方程式(wave equation)をコンピュータを用いて解くために差分化を行い，コンピュータを用いて波動方程式を解くアルゴリズムを示す．前回の拡散方程式に比べると、時間微分が2階になっているため、多少ややこしいが考え方は同じである。特に、初期条件の与え方がわかりにくいかも知れない。

さて，天下り的だが1次元波動方程式は弦の振動を（おおよそ）あらわす．次の動画は、1次元波動方程式をここで紹介する数値解法で解いた結果をアニメーションとしたものだ。 （こちらは、C言語＋あるグラフィックライブラリで作成したアニメーション）

下記方程式 (W) における u(x,t) は弦の（静止位置からの）変位をあらわす．ここでは1次元波動方程式を数値計算を用いて解き，その解 u(x,t) をアニメーションとして表示する事を目標とする．つまり，コンピュータ内部に弦の振動を再現し，それをアニメーションとして表示するわけだ．

さて，次の1次元波動方程式の初期値境界値問題、

但し、c >0, f(0)=0, f(L)=0, g(0)=0, g(L)=0．

を考える。境界条件，初期条件が適切に与えられれば (W) の解が弦の振動を表現していると言えるだろう．さて，数値計算する為にまず(W)を離散化する必要がある。これはこれまでの常微分方程式の場合、また拡散方程式と共通で，要は微分を差分で置き換える必要がある．分割方法は以下の２種類が考えられるが，

ここでは時間方向の離散化は (1) を，空間方向の離散化は (2) を採用します．（前回の拡散方程式と同じ）

[0,l] 区間を N 等分するとし，上記離散化方法に従い，h = l/N, τ>0 とし，x_j = (j - 1/2)h, t_k = kτ として，

と書くことにすると、(W)の第一式は次のように離散化される。

λ=cτ/h として、まとめると

という式が得られる。2階微分の差分化については，このページ最後のスライドを参照．この式は、時刻 t_k-1, t_kでの U から次の時刻 t_k+1での U が計算できるという形をしている。したがって、初期条件として、 U_j⁰,U_j¹ が必要となる。まず U_j⁰は(W)の第三式より、

となる。 U_j¹ を決める方法はいくつか考えられるが、例えば次のようにする。まず、テイラーの公式を用いると、

となる。(W)より

であり，さらに t=0 でも(W)の第一式が成り立つとすると、

となる。この式の右辺は２階差分商を用いて

と近似できる。つまり，

と求まる．

境界条件は（W）の第二式だが，まずは境界条件について触れる必要がある．ここでは，境界の処理の為に，U₀^k,U_N+1^k という架空のセルを用意することにする．

様々な境界条件があるが，今回の（W) における境界条件は Dirichlet 条件（この場合各時間 0 という値に固定されている．一般の場合，例えば a という値に固定されている場合には，次の式の右辺 0 が a となる．）と呼ばれるもの．空間離散化の方法として（１）を採用した事を思い出すと，

となればよいことが分かる．すなわち，

となる。（この辺り、拡散方程式と同じ）

まとめると、以下のようになる．

但し、 λ≦1 が安定性のための必要十分条件（安定性条件とは，数値計算がうまく行く為の条件）.

めでたく差分方程式が得られた。(2)をコンピューターを用いて解くのがここでの目標である。状況を、次の絵を使って説明する。

縦軸方向を時間の進行方向とし，横軸方向を空間方向とする．赤色の円盤で示した位置の値 U^k_j が，他の4点，U^k-1_j ,U^k-1_j-1 ,U^k-1_j+1 ,U^k-2_j ,から決定される様子．領域の端以外では上記のイメージで良いが，端においては境界条件の考慮が必要となる（前述）．

1次元波動方程式を解く為のアルゴリズムは次のようになる。なお、安定性条件を満たすため、λをパラメータとして先に与えて、τは与えられたλから求めるようにする。

アルゴリズム中の u_j, old_u_j, new_u_j はそれぞれ実数型の配列として定義すると良い．境界条件処理用に配列要素はN個よりも余分に２つ必要となり，計N+2個必要となる事に注意（new_u_j については，実際には境界処理用の要素は用いない）．例えば u_j について，Python風に配列を定義すると，

u = [0.0]*(N+2)

となる．このようにすれば，u[0], u[1], ... , u[N+1] が用意される．実際そのように用意すれば良い．

以下のアルゴリズムのイメージ．new_u という配列要素を u および old_u という配列要素から求める．

初期条件における、new_u, u, old_u の関係は、old_u が時刻0の弦の形状、すなわち f(x)であり与えられているもの、u が時刻 τ における弦の形状であり、f(x)およびg(x)から決定されるもの、new_u はそれらから求まる時刻 2τ における弦の形状であり未知であるものとなっている。

アルゴリズム

（１）初期パラメータ設定

N, T, λ, c, l を設定する (λ ≦ 1)

h = l/N, τ=(λh)/c, M = T/τ

x, new_u, u, old_u を配列（リスト）として用意

（２） 初期値設定

j := 1,....,N の順に

x[j] = (j - 0.5)h

old_u[j] = f(x[j])

を繰り返す

old_u[0]=-old_u[1], old_u[N+1]=-old_u[N]（境界条件：初期値決定のため）

j := 1,....,N の順に

u[j]=old_u[j]+τg(x[j])+(λ*λ)/2*(old_u[j-1]-2old_u[j]+old_u[j+1])

を繰り返す

（３）

k := 0,1,.....,M-1 の順に

(old_u[1] ～ old_u[N] を画面に表示)

u[0]=-u[1], u[N+1]=-u[N]（境界条件）

j := 1,....,N の順に

new_u[j]=2u[j]-old_u[j]+(λ*λ)*(u[j-1]-2u[j]+u[j+1])

を繰り返す

j := 1,....,N の順に

old_u[j] = u[j]

を繰り返す

j := 1,....,N の順に

u[j] = new_u[j]

を繰り返す

を繰り返す

課題1（演習）

l=1, c=1, f(x)=2*x*(1-x), g(x)=0, u(0,t)=u(1,t)=0 として (2) を解くプログラムを作成せよ。(Nは100程度)

課題2

l=1, c=1, u(0,t)=u(1,t)=0 はそのままで、f(x) および g(x) を他のもにかえてみよ。（もっと局在化した初期値からはじめると、波の伝播が観察される。）

f(x)=exp(-100*(x-0.4)*(x-0.4))

とし、異なる g(x) を初期条件としたものを、

(a)waveeq2

(b)waveeq3

この時の g(x) は(a),(b)それぞれどのようなものか考えよ。

実際の弦の振動現象と照らし合わせ，f(x) および g(x) の意味を再考せよ．

課題3

弦の振動ではなく，膜の振動を考えたい．すなわち，１次元波動方程式ではなく２次元波動方程式を考える．どのような方程式になるであろうか？

2次元波動方程式を数値計算で解き、その結果を OpenGL という可視化ライブラリを用いて可視化したサンプルは以下の様。

これに近いものを、おまけ回に載せたので参考にして欲しい。

差分近似について．