コンテンツ

第１０回

１次元拡散方程式の数値解法

これまでは変数が時間にのみ依存するものを扱っていましたが、今回からは、時間以外に空間にも依存する変数を扱います。つまり、変数としては $u(x,t)$ と書けます。１次元拡散方程式（１次元熱伝導方程式とも呼ばれる）の初期値境界値問題を数値的に解きます。
それは、次のようなものです。（教科書や文献によって記号など多少異なるかもしれません。）

問題１.１

$\displaystyle \frac{\partial u(x, t)}{\partial t} = D\frac{\partial^2u(x, t)}{\partial x^2} ~~ (0 < x < L, t> 0),$
$u(0, t) = 0,$
$u(L, t) = 0~~(t>0),$
$u(x, 0) = 2x(L-x)~~(0\le x \le L)$

1行目の式は拡散方程式、もしくは熱伝導方程式と呼ばれます。式中の $D$ は拡散係数または熱伝導率とよばれます。
拡散とは物が散らばっていく様子（例：コップにインクを垂らすと、時間とともに広がって最後には一様になる）であり、熱伝導とは熱が伝わる様子（例：スプーンを熱湯に浸けていると徐々に熱が伝わって持っている部分まで熱くなる）です。
これらの現象を表現しているのが上式なのです。
「物の散らばり」と「熱の伝わり」という異なる現象が同じ方程式で表されるというところは、数式の面白さといって良いでしょう。
熱伝導方程式と見ると $u(x,t)$ は時刻 $t$ での温度分布、拡散方程式とすると $u(x,t)$ は粒子の密度分布に対応し、何れにしても、分布が時間的な変化を考えることができます。

2行目と3行目の式は境界条件と呼ばれ、空間の端の部分（ $x=0$ および $x=L$ ）での拘束条件になります。
問題1では、 $x = 0$ と $x = L$ でいつも $u = 0$ となるディリクレ境界条件を境界に課しています。
温度の分布の言葉でいうと、両端を氷か何かでつねに冷やしている状況を考えていることになります。
つまり、熱が両端からどんどんと奪われているわけです。

4行目の式は、初期条件（ $t=0$ ）です。

以下では、拡散方程式・熱伝導方程式という呼称を拡散方程式に統一します。

空間と時間の分割

さて、計算機で扱えるよう離散化し、数値的に近似解を求めます。
これまでの常微分方程式の解法では、数値計算において、時間という本来連続な変数をそのまま扱うわけにはいかないため、 $t=t_i~(i=0,1,\cdots,N)$ のように離散化をしていたわけですが、もう一つ増えた空間を表す変数 $x$ についても同様に離散化をする必要があります。ただし、以下のように、空間方向の分割方法は時間方向の分割方法とは少し異なります。

上図の(1)は時間方向の分割方法で、これまでと変わりません。数直線で書き表すと、考える範囲を $I$ 等分した目盛り上に位置します。一方、空間方向の分割は、考える範囲を $I$ 等分するところまでは同じようにしますが、そこから目盛と目盛の中間の点に空間座標がきます。なぜこのような面倒なことをするのか、というと、境界条件について考えやすいためです（後述）。

空間の区切り方のイメージとしては、アパートのように同じ間取りの部屋で区切り、隣との壁が目盛り、物質の濃度や温度は部屋の中心で測る、といった感じでしょうか。アパートでは隣の部屋の物質が壁を通過して流れ込んだりすることはそうそうないですが、断熱材が少ないアパートは熱は多少なりとも伝わっていったりします。両隣の部屋がエアコンをガンガンにつけていると、真ん中の部屋の光熱費は少なくて済むかも？

これまで $h$ を時間刻み幅としてきましたが、今後は $h$ を空間刻み幅、 $\tau$ を時間刻み幅とします。すなわち、空間方向では $[0, L]$ 区間を $N$ 等分（ $h=L/N$ ）、時間方向では $[0, T]$ 区間を $M$ 等分（ $\tau=T/M$ ）して、 $x_j = (j - 1/2)h, t_k = k\tau$ おくことができ、

$U^k_j \approx u(x_j, t_k), \ \ (j=1, \cdots, N, k= 0, 1, \cdots , M)$

と書くことができます。

差分方程式

上記の拡散方程式と境界条件・初期条件を以下のように差分化します。

拡散方程式の差分化

拡散方程式は、上記の差分化を使えば、

$\frac{\partial u(x, t)}{\partial t} = D\frac{\partial^2u(x, t)}{\partial x^2}$

は次の差分方程式で近似されます。

$\frac{1}{\tau} (U_j^{k+1} - U_j^k) = \frac{D}{h^2}(U_{j-1}^k - 2U_j^k + U_{j+1}^k), \ \ (j=1, \cdots, N, k= 0, 1, \cdots ,M)$

これを書き換えれば、

$U_j^{k+1} = \alpha U_{j-1}^k + (1 - 2\alpha) U_j^k + \alpha U_{j+1}^k, \ \ (j=1, \cdots, N, k= 0, 1, \cdots ,M)$

が得られます。ただし、 $\alpha = D\tau / h^2$ としました。
この式は見てのとおり、時刻 $t_k$ での $u$ から次の時刻 $t_{k+1}$ での $u$ の値が計算できる形をしています。

境界条件の差分化

今回は境界条件はディリクレ境界条件ですが、他にもよく使う境界条件であるノイマン境界条件・周期境界条件も導出しておきましょう。

ディリクレ境界条件

ディリクレ境界条件の一般的な形式は次の通りです。

$u(0, t) = c,u(L, t) = c~(t>0)$

$x=0$ について考えると、 $x=0$ の目盛り上で時間に依らず常にある一定値 $c$ （今回は $c=0$ ）となるためには、境界の両隣の2点の中間値が $c$ になれば良いので、

$\displaystyle \frac{U_0^k+U_1^k}{2} = c$

より、

$U_0^k=2c-U_1^k~ (k=0, 1, 2, \cdots ,M)$

とすれば良いです。反対側の $x=L$ の境界も同様に

$U_{N+1}^k = 2c-U_N^k~ (k=0, 1, 2, \cdots ,M)$

となります。

ノイマン境界条件

ノイマン境界条件は以下の通りです。

$\frac{\partial u(0, t)}{\partial x} = 0, \frac{\partial u(L, t)}{\partial x} = 0~~(t>0)$

$x=0$ について考えると、 $x=0$ の目盛り上で時間に依らず常に空間方向の微分が0とならなければならないのですが、空間方向は離散的に区切られています。そこで、

$\frac{\partial u(0, t)}{\partial x}\approx\frac{U_1^k-U_0^k}{h} = 0$

より、

$U_0^k=U_1^k~ (k=0, 1, 2, \cdots ,M)$

とすれば良いです。反対側の $x=L$ の境界も同様に

$U_{N+1}^k = U_N^k~ (k=0, 1, 2, \cdots ,M)$

となります。

周期境界条件

周期境界条件は、周期的になっていれば良いのですが、周期関数の定義

$u(x)=u(x+L)$

を使って考えます。境界を超えて、 $x<0$ や $x>L$ も同じように関数が繰り返されると考えれば良いので、

$U_0^k=U_N^k$
$U_{N+1}^k=U_1^k$

となります。

初期条件

初期条件は、 $x$ の離散化のみ考えればよく、

$U_j^0 = 2x_j(1-x_j) ~(j=1, 2, \cdots ,N)$

になります。

これで、準備はできました。

次に拡散方程式を計算機にて解くためのアルゴリズムを考えましょう。
解を求める時刻の上限を $T$ とし、時間刻み幅 $\tau$ に対応する時間刻み数を自然数 $M = T/\tau$ で定めます。空間方向も同様に、空間幅 $L$ および空間刻み幅 $h$ に対して、空間刻み数を自然数 $N = L/h$ で定めます。
問題１.１を解くアルゴリズムは次のようになります（ただし、 $D=1.0$ としています）。

アルゴリズム

ライブラリのインポート

（１）変数と初期パラメータを設定

T（時間上限）, τ（時間刻み幅）, L（空間幅）, h（空間刻み幅）, D（拡散係数）を設定する

M = int(T/τ), N = int(L/h), α = Dτ/(h^2)　(プログラム中では, τは「tau」, αは「alpha」等とする.)

x[N+2], u[N+2], new_u[N+2]（偏微分方程式の変数をリストとして準備する。）

（２）初期値設定

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

x[j] = (j - 0.5)*h

u[j] = 2*x[j]*(1-x[j])

を繰り返す

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

描画する準備

figとaxをつくる。

（３）時間発展を計算

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

t = k*τ

（３.１）数ステップ毎に(x[1], u[1]) ～ (x[N], u[N]) を画面に表示

(ただし x[0], u[0], x[N+1], u[N+1] は境界処理用なので表示しない)

更新されたuを描画

（３.２）計算

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

new_u[j] = αu[j-1] + (1 - 2α)u[j] + αu[j+1] # t=(k+1)*τのときのu

を繰り返す

u[1:-1] = new_u[1:-1] （リストの更新）

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

を繰り返す

ライブラリのインポート（１）変数と初期パラメータを設定 T（時間上限）, τ（時間刻み幅）, L（空間幅）, h（空間刻み幅）, D（拡散係数）を設定する M = int(T/τ), N = int(L/h), α = Dτ/(h^2)　(プログラム中では, τは「tau」, αは「alpha」等とする.) x[N+2], u[N+2], new_u[N+2]（偏微分方程式の変数をリストとして準備する。）　（２）初期値設定 j := 1,....,N の順に x[j] = (j - 0.5)*h u[j] = 2*x[j]*(1-x[j]) を繰り返す u[0] = -u[1], u[N+1] = -u[N] 　　（境界条件）描画する準備 figとaxをつくる。（３）時間発展を計算 k := 0,1,.....,M の順に t = k*τ （３.１）数ステップ毎に(x[1], u[1]) ～ (x[N], u[N]) を画面に表示 (ただし x[0], u[0], x[N+1], u[N+1] は境界処理用なので表示しない) 更新されたuを描画（３.２）計算 j := 1,2,....,N の順に new_u[j] = αu[j-1] + (1 - 2α)u[j] + αu[j+1] # t=(k+1)*τのときのu を繰り返す u[1:-1] = new_u[1:-1] （リストの更新） u[0] = -u[1], u[N+1] = -u[N] 　　（境界条件）を繰り返す

ライブラリのインポート
（１） 変数と初期パラメータを設定
T（時間上限）, τ（時間刻み幅）, L（空間幅）, h（空間刻み幅）, D（拡散係数） を設定する
M = int(T/τ), N = int(L/h), α = Dτ/(h^2)　(プログラム中では, τは「tau」, αは「alpha」等とする.)
x[N+2], u[N+2], new_u[N+2]（偏微分方程式の変数をリストとして準備する。）
　
（２） 初期値設定
j := 1,....,N の順に
    x[j] = (j - 0.5)*h
    u[j] = 2*x[j]*(1-x[j])
を繰り返す
u[0] = -u[1], u[N+1] = -u[N] 　　（境界条件）
描画する準備
figとaxをつくる。
（３）時間発展を計算
k := 0,1,.....,M の順に
    t = k*τ
    （３.１）数ステップ毎に(x[1], u[1]) ～ (x[N], u[N]) を画面に表示
        (ただし x[0], u[0], x[N+1], u[N+1] は境界処理用なので表示しない)
        更新されたuを描画
    （３.２）計算
    j := 1,2,....,N の順に
        new_u[j] = αu[j-1] + (1 - 2α)u[j] + αu[j+1] # t=(k+1)*τのときのu
    を繰り返す
    u[1:-1] = new_u[1:-1] （リストの更新）
    u[0] = -u[1], u[N+1] = -u[N] 　　（境界条件）
    
を繰り返す

このアルゴリズムを実現するプログラムを作成せよというのが当面の課題となります。

リストの更新や要素を限定したplotの部分については、第5回の小ネタを参考にしてください。

課題1

上のアルゴリズムに従ってプログラムを完成させよ。
ただし、 $T=0.5, \tau=0.0001, L=1.0, h=0.02, D=1.0$ 程度、
初期条件は $u(x,0) = 2x(1-x)$ 、境界条件はディリクレ条件とする。そのままだと描画に時間がかかるので、 $t$ が0.005経過するごとに描画せよ。

実際の計算部分を正しくプログラムにすることができれば、次のようなアニメーションを見ることができます。
両端から熱（物質）が奪われ、全体の温度（密度）が下がることがわかります。

課題1′

実は「要素を部分的に抜き出す」ことを使えば、アルゴリズムの18,19行目はfor文を使わずに書くことができる。どうすれば良いか？

ヒント：

import numpy as np

A = np.array([1., 2., 3., 4.])

D = A[:3] + A[1:]

import numpy as np A = np.array([1., 2., 3., 4.]) D = A[:3] + A[1:]

import numpy as np
A = np.array([1., 2., 3., 4.])
D = A[:3] + A[1:]

と

import numpy as np

A = np.array([1., 2., 3., 4.])

D = np.zeros(3)

for i in range(3):

D[i]= A[i]+A[i+1]

import numpy as np A = np.array([1., 2., 3., 4.]) D = np.zeros(3) for i in range(3): D[i]= A[i]+A[i+1]

import numpy as np
A = np.array([1., 2., 3., 4.])
D = np.zeros(3)
for i in range(3):
    D[i]= A[i]+A[i+1]

は同じ意味を持ちます。

課題1”

Plotされたグラフは、[0, x[1]]の範囲と[x[N], L]の範囲が描画されていない。どうすれば解決されるか？

Masashi FUJII

計算数理A演習第10回

第１０回

１次元拡散方程式の数値解法

問題１.１

空間と時間の分割

差分方程式

拡散方程式の差分化

境界条件の差分化

ディリクレ境界条件

ノイマン境界条件

周期境界条件

初期条件

アルゴリズム

課題1

課題1′

課題1”

Department of Mathematical and Life Sciences, Graduate School of Integrated Sciences for Life, Hiroshima University