コンテンツ

第１３回

波動方程式の数値解法

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

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

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

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

(W) $\begin{equation*}\left\{ \begin{array}{lc}\displaystyle\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), \displaystyle\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)を採用してきました。今回も同様です。

$[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})$

という漸化式が得られます。
２階微分の差分化については、このページ最後のスライドを参照してください。

初期条件

初期条件も差分化する必要がありますが、拡散方程式のときと異なる点があります。
先ほどの離散化した漸化式において、次の時刻 $t_{k+1}$ での $U$ を計算するためには2つの時刻 $t_{k-1}$ 及び $t_k$ での $U$ が必要です。
したがって、漸化式を順次計算する上で、まずは $U_j^0$ と $U_j^1$ が必要となります。
まずは(W)の第三式より、

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

となることは良いでしょう。では $U_J^1$ はどのように決めれば良いでしょうか？
$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)における境界条件は Dirichlet条件ですので、

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

となります。
他の境界条件の場合については、第10回を参考にして下さい。

解くべき（差分）方程式

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

(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$ を与えます。

拡散方程式のときには「u_j 」、「 new_u_j 」のみ定義していましたが、今回は「old_u_j 」も必要になります。また、それぞれのリストの個数は境界条件の計算用に2つ加えたN+2個となります。（new_u_jについては、実際には境界処理用の要素は用いないですが…）

(１)初期パラメータ設定
T（時間上限）, M（時間刻み数）, L（空間幅）, N（空間刻み数）, c を設定する
λ = τ*c/h (λ ≦ 1 となるように. h, τ 決め、N = L/h, M = T/τ（いずれも自然数）とする。 プログラム中では, λは「lamb」, τは「tau」などとするとよい.)
 
（２） 初期値設定
j := 1,.....,N の順に
    x_j = (j - 0.5)*h
    old_u[j] = f(x_j) # t=0のときのu
を繰り返す ( f(x_j) は課題１,２を参照 )
old_u[0] = -old_u[1], old_u[N+1] = -old_u[N] （境界条件）
j := 1,.....,N の順に
    u[j] = old_u[j] + tau*g(x_j) + (λ*λ/2) * (old_u[j-1] - 2 * old_u[j] + old_u[j+1]) # t=hのときのu
を繰り返す ( g(x) は課題１,２を参照 )
u[0] = -u[1], u[N+1] = -u[N] （境界条件）

（３）計算, 動画描画
k := 0,....,M の順に
    t = k*τ
    （i）描画
    一定間隔毎に(例えば k%100==0 の時に, 等)
        (x[1], old_u[1]) ～ (x[N], old_u[N]) を画面に表示
        (x[0], old_u[0], x[N+1], old_u[N+1] は境界処理用なので表示しない)

    （ii）計算
    j := 1,.....,N の順に
        new_u[j] = 2*u[j] - old_u[j] +λ*λ(u[j-1] - 2 * u[j] + u[j+1]) # t=(k+2)*τのときのu
    を繰り返す
    new_u[0] = -new_u[1]
    new_u[N+1] = -new_u[N]
 
    （iii）更新（順番に注意）
    old_u[:] = u[:]
    u[:] = new_u[:]
    

を繰り返す

課題1

$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 \le 1$ となるように $M$ を適度に選び, 描写頻度も適度に選ぶ事。)

課題2

$L=1, c=1, u(0,t)=u(1,t)=0$ ( $N = 100$ 程度) はそのままで、 $f(x)$ および $g(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)$

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

課題3

境界条件を、Neumann境界条件や周期境界条件を適応したプログラムを作成し、実際に数値計算してみよ。

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

発展課題1

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

発展課題2

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

Pythonでの3次元プロットの方法は各自調べてみてください。

補足：差分近似について

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

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

(S1) $\begin{equation*}u(x+d) = u(x) + d\dot{u}(x) + \frac{1}{2}d^2 \ddot{u}(x) + \frac{1}{6}d^3 \dddot{u}(x) + O(d^4)\end{equation*}$

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

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

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

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

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

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

Masashi FUJII

計算数理A演習第13回