計算数理A演習第13回・第14回

今回の演習

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

波動方程式の数値解法

今回は１次元波動方程式(wave equation)を扱います。
ここでは，１次元波動方程式をコンピュータを用いて解くために差分化を行い，コンピュータを用いて波動方程式を解くアルゴ リズムを示します。
今日の内容に関する詳しい情報は講義で説明があります（ました）。
演習では細かい事は気にせずに一気にアルゴリズムの紹介まで行きたい と思います。
（なんとか，今週の演習，講義および来週の演習で１次元波動方程式は解けるようになってほしい。）

さ て，天下り的ですが１次元波動方程式は弦の振動をあらわします。
（導出等は講義で説明されます。また前回の「連結バネ系」からも導出できます。時間があれば紹介します。）
下記方程式 (W) におけるu(x, t)は弦の変位をあらわします．
演習では１次元波動方程式を数値計算を用いて解き，その解u(x, t)を画面上にアニメーションとして表示する事を目標とします。
つまり，コンピュータ内部に弦の振動を再現し，それをアニメーションとして表示するわけです。

時間ー空間２階偏微分方程式の差分化

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

(W)   \begin{equation*}\left\{ \begin{array}{lc}\frac{\partial^2 u(x, t)}{\partial t^2} = c^2\frac{\partial^2u(x, t)}{\partial x^2} & (0 < x < \ell, t> 0) \\u(0, t) = u(\ell, t) = 0 & (t>0) \\u(x, 0) = f(x), \frac{\partial u(x, 0)}{\partial t} = g(x) & (0<x<\ell)\end{array}\right.\end{equation*}

但し、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, \tau>0とし，x_j = (j-1/2)h, t_k=k\tauとして，

    \[ U_j^k \approx u(x_j, t_k) \ \ (j=1, \cdots, N, k=0, 1, \cdots) \]

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

    \[ \frac{1}{\tau^2} (U_j^{k-1} - 2U_j^k + U_j^{k+1}) = \frac{c^2}{h^2} (U_{j-1}^{k} - 2U_j^k + U_{j+1}^{k})\ \ (j=1, \cdots, N, k=1, 2, \cdots) \]

\lambda = c\tau/hとして、まとめると

    \[U_j^{k+1} = 2U_j^k - U_j^{k-1} + \lambda^2(U_{j-1}^{k} - 2U_j^k + U_{j+1}^{k})\]

という式が得られます。
２階微分の差分化については，このページ最後のスライドを参照してください。
（右辺最後の\lambda^2(\cdots)という項は, どこかで見覚えがあると思います。（連結バネ系と同じです。））

初期条件

常微分方程式同様, 初期条件を与える必要があります。
これは空間全域に対して与える必要があります。
また上記の式はこの式は、時刻t_{k-1}, t_kでのUから次の時刻t_{k+1}でのUが計算できるという形をしています。
したがって、初期条件として、 が必要となります。
まず は(W)の第三式より、

    \[ U_j^0 = f(x_i) \ \ (i=1, 2, \cdots, N) \]

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

    \[ u(x, \tau) = u(x, 0) + \tau\frac{\partial u(x, 0)}{\partial t} + \frac{\tau^2}{2} \frac{\partial^2 u(x, 0)}{\partial t^2} + O(\tau^3) \]

となります。(W)より

    \[ \frac{\partial u(x, 0)}{\partial t} = g(x) \]

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

    \[ \frac{\partial^2 u(x, 0)}{\partial t^2} = c^2 \frac{\partial^2 u(x, 0)}{\partial x^2} \]

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

    \[\frac{\partial^2 u(x, 0)}{\partial x^2} \approx \frac{c^2}{h^2} \{ u(x-h, 0) - 2u(x,0) + u(x+h, 0) \} \]

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

    \[ U_j^1 = U_j^0 + \tau g(x_j) + \frac{\lambda^2}{2} (U_{j-1}^0 - 2U_j^0 + U_{j+1}^0) \]

と求まります．

境界条件

空間に関しても端を適切に扱う必要があり，そのための境界条件が与えられます。
境界条件は（W）の第二式ですが，まずは境界条件について触れる必要があります。
詳しくは，講義で説明がありますが，ここでは，境界の処理の為に，U_0^k, U_{N+1}^kという架空のセルを用意することにします．（以下講義のプリント）

様々な境界条件がありますが，今回の（W)における境界条件は Dirichlet条件と呼ばれるものです．
空間離散化の方法として（２）を採用した事を思い出すと，

    \[ \frac{U_0^k + U_1^k}{2} = 0,\frac{U_N^k + U_{N+1}^k}{2} = 0 \ \ (k=1, 2, \cdots)\]

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

    \[ U_0^k = - U_1^k, U_{N+1}^k = - U_N^k \ \ (k=1, 2, \cdots) \]

となります．

解くべき（差分）方程式

以上の考察から解くべき式をまとめると、以下のようになります．

(2)   \begin{equation*}\left\{ \begin{array}{l}U_j^{k+1} = 2U_j^k - U_j^{k-1} + \lambda^2(U_{j-1}^{k} - 2U_j^k + U_{j+1}^{k}) \ \ (j=1, \cdots, N, k=1, 2, \codts) \\\frac{U_0^k + U_1^k}{2} = 0,\frac{U_N^k + U_{N+1}^k}{2} = 0 \ \ (k=1, 2, \cdots) \\U_j^0 = f(x_j) \ \ (j=1, \cdots, N) \\U_j^1 = U_j^0 + \tau g(x_j) + \frac{\lambda^2}{2} (U_{j-1}^0 - 2U_j^0 + U_{j+1}^0)\end{array}\right.\end{equation*}

ただし、\lambda \le 1 \ (\lambda=c\tau/h)が安定性のための必要十分条件です。
拡散方程式では1/2が上限でしたが，波動方程式では違う値になります。
（安定性条件とは，数値計算がうまく行く為の条件ですなぜにこういった条件が必要であるかは，講義で説明があります）.

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

アルゴリズム

あくまで，参考用です．始めに自分で解き方を考えてください。

１次元波動方程式を解く為のアルゴリズムは次のようになります。
安定性条件を満たすため\lambda \le 1となるようにh, \tauを与えます。

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

float u[N+2];

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

アルゴリズム

(１)初期パラメータ設定
T（時間上限）, M（時間刻み数）, L（空間幅）, N（空間刻み数）, c を設定する
h = L/N, τ=T/M, λ = τ*c/h (λ ≦ 1 となるように. M, N を調整, プログラム中では, λは「lamb」, τは「tau」などとするとよい.)

（２） 初期値設定
j := 1,....,N の順に
    xj = j*h - 0.5*h
    old_u[j] = f(xj)
を繰り返す ( f(x) は課題１,２を参照 )

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

j := 1,....,N の順に
    u[j] = old_u[j] + tau*g(xj) + (λ*λ/2) * (old_u[j-1] - 2 * old_u[j] + old_u[j+1])
を繰り返す ( g(x) は課題１,２を参照 )

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

（３）計算, 動画描画
k := 0,1,.....,M-1 の順に
    (x[1], u[1]) ～ (x[N], u[N]) を画面に表示 (x[0], u[0], x[N+1], u[N+1] は境界処理用なので表示しない)
    j := 1,....,N の順に
        new_u[j] = 2*u[j] - old_u[j] +λ*λ(u[j-1] - 2 * u[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] （境界条件）
を繰り返す

課題１

L=1, c=1, f(x) = 2x(1-x), g(x) = 0, u(0, t)=u(1, t) = 0として (2) を解くプログラムを作成してください。
(Nは100程度，T \sim 100で\lambda \ll 1となるようにMを適度に選び, 描写頻度も適度に選ぶ事.)

課題２

L=1, c=1, u(0,t)=u(1,t)=0(N = 100程度) はそのままで、f(x)およびf(x)を他のもにかえてみよ。例えば関数g(x)=0、関数f(x)を

f(x) = \exp(-100(x-0.4)^2)

としたり

f(x) = \exp(-100(x-0.4)^2) - \exp(-50(x-0.7)^2) + \exp(-10(x-0.5)^2)

とする等。
（局在化した初期値からはじめると、波の伝播が観察される。）

課題３

他の境界条件，Neumann境界条件，周期境界条件について復習し（講義を受けていない人は図書館で調べる），それぞれの境界条件を適応したプログラムを作成し，実際に数値計算してみよ．

周期境界条件u(x,t)=u(x+1,t)のもと、進行波解\sin 2\pi (x-t)を実現するためには、どのようなパラメータで、どのような初期値から計算する必要が有るか

課題４

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

課題5

２次元波動方程式を差分化し、プログラムを作成し、数値計算せよ

 

差分近似について．

導関数の差分近似について

関数 u(x) を x のまわりでTaylor展開すると

(S1)   \begin{equation*}u(x+d) = u(x) + du'(x) + \frac{1}{2}d^2 u''(x) + \frac{1}{6}d^3 u'''(x) + O(d^4)\end{equation*}

(S2)   \begin{equation*}u(x-d) = u(x) - du'(x) + \frac{1}{2}d^2 u''(x) - \frac{1}{6}d^3 u'''(x) + O(d^4)\end{equation*}

となり、両辺を加えると、

(S3)   \begin{equation*}u(x+d) + u(x-d) = 2u(x) + d^2 u''(x) + O(d^4)\end{equation*}

となります。O(d4)を無視すると、

(S4)   \begin{equation*}u'' \approx \frac{1}{d^2} \{ u(x+d) - 2u(x) + u(x-d) \}\end{equation*}

この式の右辺の誤差はd^2程度となります。