Part.17

１次元波動方程式の数値解法

ここでは，１次元波動方程式(wave equation)をコンピュータを用いて解くために差分化を行い，コンピュータを用いて波動方程式を解くアルゴリズムを示します．

さて，天下り的ですが１次元波動方程式は弦の振動をあらわします．次の動画は、1次元波動方程式をここで紹介する数値解法で解いた結果をアニメーションとしたものです。 （なお、GLSCで画面を静止画に保存して、アニメーションを作成する方法は「macOSでシミュレーション（C言語+GLSC） 」にある。）

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

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

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

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

ここでは時間方向の離散化は(1)を，空間方向の離散化は(2)を採用します．（これまでの演習で常微分方程式の数値解法を扱ったが，そこでは時間方向の離散化として(1)を採用してきた）

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

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

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

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

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

となります。(W)より

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

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

と近似できます。つまり，

と求まります．

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

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

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

となります．

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

但し、 λ≦1 が安定性のための必要十分条件（安定性条件とは，数値計算がうまく行く為の条件です．なぜにこういった条件が必要であるかは，後に説明します）.

めでたく差分方程式が得られました。(2)をコンピューターを用いて解くのがここでの目標です。

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

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

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

double u[N+2];

となる．このようにすれば，u[0], u[1], ... , u[N+1] が用意される．実際そのように用意すれば良い．または、配列のところで説明した動的確保(malloc を用いる)でもよい。

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

アルゴリズム

（１）初期パラメータ設定
 N, T, λ, c, l を設定する (λ≦ 1)

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

（２） 初期値設定
 j := 1,....,N の順に
  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]) 
 を繰り返す
 u[0]=-u[1], u[N+1]=-u[N]（境界条件）
（３）
 k := 0,1,.....,M-1 の順に
  u[1] ～ u[N] を画面に表示 (g_move + g_plot, u[0], u[N+1] 
  は境界処理用なので表示しない)
 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]
 を繰り返す
 old_u[0]=-old_u[1],old_u[N+1]=-old_u[N]（境界条件）
 j := 1,....,N の順に
  u[j] = new_u[j]
 を繰り返す
 u[0]=-u[1], u[N+1]=-u[N]（境界条件）
を繰り返す

上記アルゴリズム中のグラフの画面表示は，全ての k について行う必要は無い．例えば，計算100回に一度グラフを表示するには，次のようにする．(100回に一度表示するのが適当かどうかは試してみないとわからない．)

if(k%100 == 0)
{
  u[1] ～ u[N] を画面に表示 (g_move + g_plot, u[0], u[N+1] 
  は境界処理用なので表示しない)
}

また，アニメーションを作成するには，画面の消去が必要であるが，画面全体を g_cls() で消去するのでは無く，書き換えが必要な部分を必要最小限に消去（背景色と同じ色で上書きする等）すればなめらかなアニメーションとなる．

課題１

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

上記パラメータでのサンプルは以下。

課題２

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 という可視化ライブラリを用いて可視化したサンプルは以下の様。

課題4

課題2では，それぞれの初期条件に関して別個のプログラムを作成したが，例えば，一つのプログラムではじめに初期値を数値で選べるようなプログラム（1 を選べば(a) waveeq2 の初期値となり 2 を選ぶと (b)waveeq3 の初期値となるような）を作成せよ．

差分近似について．