今回は１次元波動方程式(wave equation)を扱います。ここでは，１次元波動方程式をコンピュータを用いて解くために差分化を行い，コンピュータを用いて波動方程式を解くアルゴリズムを示します．
さて，天下り的ですが１次元波動方程式は弦の振動をあらわします．下記方程式(W)における u(x,t) は弦の（静止位置からの）変位をあらわします．演習では１次元波動方程式を数値計算を用いて解き，その解 u(x,t) を画面上にアニメーションとして表示する事を目標とします．つまり，コンピュータ内部に弦の振動を再現し，それをアニメーションとして表示するわけです．
さて，次の１次元波動方程式の初期値境界値問題、
(W)
但し、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 が計算できるという形をしています。したがって、初期条件として、 
が必要となります。まず 
は(W)の第三式より、
となります。 
を決める方法はいくつか考えられますが、例えば次のようにします。まず、テイラーの公式を用いると、
となります。(W)より
であり，さらに t=0 でも(W)の第一式が成り立つとすると、
となります。この式の右辺は２階差分商を用いて
と近似できます。つまり，
と求まります．
境界条件は（W）の第二式ですが，まずは境界条件について触れる必要があります．ここでは，境界の処理の為に，
という架空のセルを用意することにします．
様々な境界条件がありますが，今回の（W)における境界条件は Dirichlet 条件（この場合各時間 0 という値に固定されている．一般の場合，例えば a という値に固定されている場合には，次の式の右辺 0 が a となる．）と呼ばれるものです．空間離散化の方法として（１）を採用した事を思い出すと，
となればよいことが分かります．すなわち，
となります．
まとめると、以下のようになります．
(2) 
但し、 
が安定性のための必要十分条件（安定性条件とは，数値計算がうまく行く為の条件です）.
めでたく差分方程式が得られました。(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] が用意される．実際そのように用意すれば良い．
以下のアルゴリズムのイメージ．new_u という配列要素を u および old_u という配列要素から求める．
（１）初期パラメータ設定
N, T, λ, c, l を設定する (λ≦ 1)
h = l/N, τ=(λh)/c, M = T/τ（２） 初期値設定
j := 1,....,N の順に
を繰り返す
（境界条件）
j := 1,....,N の順に
を繰り返す
（境界条件）
（３）
k := 0,1,.....,M-1 の順に
u[1] 〜 u[N] を画面に表示 (g_move + g_plot, u[0], u[N+1]
は境界処理用なので表示しない)
j := 1,....,N の順に
を繰り返す
j := 1,....,N の順に
を繰り返す
j := 1,....,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) としてそれぞれ、
g(x) = 0,
g(x) = f'(x),
とした場合のシミュレーションをしてみよ．
課題３
弦の振動ではなく，膜の振動を考えたい．すなわち，１次元波動方程式ではなく２次元波動方程式を考える．どのような方程式になるであろうか？
二次元波動方程式を適当な初期値で解き，適当な方法で可視化した例は次のようになる（動画像）．
差分近似について．
関数 u およびその導関数が x の１価有限連続関数であるとき、Taylorの定理から、
(0.1)
また、
(0.2)
となります。これらの展開式を加えると、
(0.3)
が与えられます。ここで、 
は h の４次以上を含む項を示すものとします。これらは h のより低次のこうと比較して無視できるとすれば、次のようになります。
(0.4)
この式の右辺の主要誤差は h2 程度となります。また、(0.1)の h2 程度の項を無視すると、
(0.5)
が得られます。
