偏微分方程式とセルオートマトン
(1)偏微分方程式の典型
未知数(従属変数)ではなく独立変数が2つ以上の場合は偏微分方程式となる。理工学でよく現れる線形偏微分方程式の典型には以下のようなものがある。
偏微分記号∂が印刷のとき文字化けを起こすようなので,ここでは単に d で書くことにする。
u が (x, y, z) の関数であるとして
ラプラス方程式:Δu = 0 ここで Δ = d2/dx2+d2/dy2+d2/dz2
ポアソン方程式:Δu = f(x,y,z) (線形・非同次方程式)
ヘルムホルツ方程式:(Δ+k2)u = 0
拡散方程式:du/dt = DΔu
波動方程式:d2u/dt2 = c2Δu
いずれも,これに初期条件,境界条件が与えられて初めて解を求めることができる。また,未知数がベクトル(3変数)でしかも非線形方程式である流体力学方程式
ナビエ-ストークス方程式:
du/dt + [u・▽]u = R-1Δu - ▽p + F ただし ▽ = (d/dx, d/dy, d/dz)
[u は流体の速度場,p は圧力,F は外力。R はレノルズ数と呼ばれるパラメータ]や,線形ではあるが電場・磁場ベクトルについて連立方程式であるマクスウェル方程式(省略)のようなものもある。あるいはナビエ-ストークス方程式を1次元に簡単化したモデル方程式
バーガース方程式:du/dt + u du/dx = D d2u/dx2
もある。同類のもので,ソリトン(孤立波)解を持つことで有名な
KdV方程式:du/dt + u du/dx = m d3u/dx3
もあり,こちらは3階微分を含んでいる。
いずれの場合も数値解法の基本は,各変数について差分化することである。差分化する際に必ずしも均一な間隔(メッシュ)を採用する必要はなく,連続変数に関する情報を有限のデータからいかにして効率よく得ることができるかについて様々な工夫が行われており,有限要素法とか境界要素法とか呼ばれている。また,線形方程式に関してはフーリエ変換・ラプラス変換の方法もある。
(2)拡散方程式
簡単のため1次元の拡散方程式
du/dt = D d2u/dx2
だけを扱うことにする。時間 t についての差分は前章のとおり u(t+h) - u(t) にとる。空間 x についての2階微分は,方程式の x について正負対称である性質を保持する意味で
[(u(x+a)-u(x))/a - (u(x)-u(x-a))/a] /a = [u(x-a) - 2u(x) + u(x+a)]/a2
とするのが妥当である。このとき,x = na とし,時間についての積分はオイラー法を用いるとして,漸化式は
u(n,t+h) = u(n,t) + (Dh/a2) [u(n-1,t) - 2u(n,t) + u(n+1,t)]
となる。もちろん、空間変数を差分化した以上,{u(n,t)}についての常微分連立方程式であるから,時間についてはルンゲ-クッタ法を用いてもよい。
これを与えられた初期条件と境界条件(両端での値)のもとで繰り返せばよい。ただし,対称性を保持するためには,この式を逐次的に,たとえばそのまま左から(n の小さい方から)順番に更新するのではなく,すべての要素を同時に「イッセイノーで」更新(並列処理と同等の処理)しなければならないことはいうまでもない。
この漸化式では係数 Dh/a2 ( = b とする)を十分小さくしなければならず,したがって a を無制限に小さくできないことが予想されよう。この漸化式を
u(n,t+h) = b u(n-1,t) + (1-2b) u(n,t) + b u(n+1,t)
と書けば,ベクトルの変換式となり,変換の行列は
( ............. )
(....0 0 0 b 1-2b b 0 0 0 0 ..............)
(....0 0 0 0 b 1-2b b 0 0 0 ..............)
(....0 0 0 0 0 b 1-2b b 0 0 ..............)
( .......... )
である。この行列を対角化したときの対角成分,すなわち行列の固有値の絶対値が全て1より小であることが収束の条件である。この条件は
b<1/2 すなわち a2>Dh/2
である。
(フーリエ級数法) この対角化を実現する方法がフーリエ級数法である。x の範囲を [0, L] とし,簡単のため境界条件は u(0) = u(L)= 0 とする。このとき,関数 u(x) は以下のように展開することができる。
∞ L
u(x) = u[k] sin(πkx/L) , u[k] = 2∫ u(x)sin(πkx/L)dx
k=1 0
u(x,t) に関する差分漸化式を用いれば,u[k,t] に対する漸化式が得られ
u[k,t+h] = [1 - 4b sin2(πka/2L)] u[k,t]
となる。つまり,上のように u(x) の一次結合で u[k] を定義することにより対角化される(連立方程式でなくなる)。この係数すなわち行列の固有値が
1 - 4b sin2(πka/2L)>-1
を満たしていなければならないことから,収束条件 「b<1/2」 が得られる。
この方法は他の線形方程式にも適用でき,実用的にもたいへん有用である。
(1次元拡散方程式の一般解)
点 x0 にだけ重み1で集中した分布(デルタ関数という)から出発したときの特解は
u(x,t) = (4πDt)-1/2 exp[-(x-x0)2/4Dt]
である。したがって,t=0 での分布が与えられれば,その後の解は
u(x,t) = (4πDt)-1/2∫u(x0,0) exp[-(x-x0)2/4Dt] dx0
と書くことができる。
(バーガース方程式について)
従属変数の変換
x
v(x,t) = exp[ -(1/2D)∫ u(s,t) ds ] or u = -2D (d/dx)log v
により,線形の拡散方程式
dv/dt = D d2v/ dx2
に帰着する。拡散方程式の解を利用すれば,
x0
U(x0) = ∫u(s, 0) ds
として
v(x,t) = (4πDt)-1/2∫exp{-(1/2D)[(x-x0)2/2t + U(x0)]} dx0
となる。
(3)差分方程式−−−どっちが先か?
ノートNo.3において,ランダムウオークは時間と空間について少しならした見方をすれば,酔っぱらいの滞在確率が拡散方程式に従うことを説明した。つまり,ランダムウオークは拡散現象の差分化モデルになっているのである。
[No.3で示した式は上の漸化式で b=1/2 の安定限界に相当するため,振動してしまうじゃないかと思うかもしれない。この振動は波長が格子点間隔の a に相当するフーリエ成分の振る舞いであって,離散的な跳躍モデルでは避けられない現象である。]
しかしながら,実際のブラウン運動(花粉から飛び出した微粒子が水分子に弾かれることにより激しく右往左往する運動)や,煙の粒子が空気分子に弾かれながら拡散していく様子を思い浮かべれば,離散化モデルであるランダムウオークの方がより忠実な基本方程式であり,これを連続近似したのが偏微分方程式の拡散方程式であると考える方が妥当かもしれない。
[注:ただし,数学(確率過程)でブラウン運動というときは,跳び移る確率がガウス分布に従い,確率分布が拡散方程式に従う運動のことを指す。]
(4)セル・オートマトン
独立変数である時間や空間座標を差分化するだけでなく,従属変数(関数)の取りうる値も離散化 (ディジタル化) したモデルで法則を表現することも可能である。有限系しか扱えない点を除けば,丸め誤差のことを心配せずに正確な計算ができ,計算機をもっとも効率的に駆使することができる。実際,この方法はビット演算を基本とするディジタル回路網の発想そのものと関わりを持っており,計算機が開発された当時から研究が行われてきた。その後「ライフゲーム」などを経て,今日ではセルオートマトン法と呼ばれ様々な分野で応用されている。
(1次元モデル)
(0, 1) のビットを1列に並べておき,次の瞬間の各ビットの値を「現在の自分と両隣 (場合によってはさらにその隣) の値のセットから決める規則」を与え,ビット列を同時的に更新する。
(規則の例)
111 110 101 100 011 010 001 000
| | | | | | | |
V V V V V V V V
1 0 1 1 1 0 0 0 <---- 184の2進数表現
(10111000)と読む
このような規則は,「どんな場合でも0に変わる」「どんな場合でも1に変わる」というトリビアルなものも含めて全部で256種類作ることができる。これに通し番号をつけて,上の例なら「ルール184」のように名付ける。
得られるパターンは,・一様なパターン,・周期的パターン,・進行パターン,・繰り返しパターン,などに分類することができる。この「ルール184」は0(1)の数を保存する進行パターンになり,先にあげたバーガース方程式の離散化モデルとも解釈される。
(多次元モデル)
多次元への拡張は明白であろう。しかしパターンの変化を図示し目で見て確認できるのは2次元に限られる。この場合は,2次元の正方格子網の各格子点にビットを配置しておき,自分と左右上下(斜め上下)の5つ(9つ)の値のセットから次の値を決めるルールを与える。
「ライフゲーム」の例
・各格子点は「生」「死」の2値をとる。
・左右上下斜め上下の8つの近隣の点のうち
「生」が1以下または4以上の場合には,「生」は「死」へ
「生」がちょうど3個の場合だけ,「死」は「生」へ
変わる。
プログラムは簡単である。
(例) 最初に中心に1個だけ 「生」# があり,ルールが 「8つの近隣のうち1個が『生』のときだけ新たに『生』に転じ,それ以外は変わらない」 の場合:
# # # # # # # #
# # # # # # # # # # # # # # # # # #
# # # # # # # #
# # # # # # # # # # # # # # # #
# # # # # # # # # # # # # #
# # # # # # # # # # # # # # # #
# # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # #
# # #
# # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# # # #
# # # # # # # # # # # # # # # #
# # # # # # # # # # # # # #
# # # # # # # # # # # # # # # #
# # # # # # # #
# # # # # # # # # # # # # # # # # #
# # # # # # # #
INTEGER :: b(-40:40, -40:40) = 0, a(-40:40, -40:40) = 0, i, j, t
b(0,0) = 1
DO t = 1, 38
DO i = -39, 39
DO j = -39, 39
a(i,j) = b(i+1,j) + b(i+1,j+1) + b(i,j+1) + b(i-1,j+1) &
+ b(i-1,j) + b(i-1,j-1) + b(i,j-1) + b(i+1,j-1)
END DO
END DO
!
DO i= -40, 40
DO j= -40, 40
IF( a(i,j)==1 ) b(i,j)=1
END DO
END DO
END DO
.................
Fortranでは配列演算式を用いれば,もっと簡単に書けるが,Cプログラムでも解釈可能な形で書いた。作業を2段階に分けることにより,配列 b の全要素は(結果的には)同時的に更新されるようになっていることに注意しよう。
<課題5−1>
(交通渋滞) 上で1次元の例にあげた 「ルール184」 は,「高速道の自然渋滞モデル」 と呼ばれるもので,(0,1) は各点に 「車がいない」 「いる」 を表し,「すぐ前 (右) が空いているときには車は1つずつ前 (右) へ進み,他の車がいるときには止まったまま」 という規則に相当する。
端は周期的に繋ぐ(=サーキットだ!)とし,適当な長さの格子点上に最初に(0,1)をランダムにばらまいておいて,その後のパターン変化の様子を(一部分でよいから)表示してみよ。
プログラムはすべてのルール番号(≦255)について適用可能な形で書くことにせよ。[ルールは関数ではなく,予め要素数8の配列で表すようにする方が効率的である。]
<課題5−2>
(森林火災) 森林にまばらに木が生えている。それぞれの木は,隣接する木のうち1本でも燃えだしたら自分も発火する。どこかで火災が発生したとき,どれくらいの密度で木が生えていたら森林のほとんどが類焼してしまうであろうか? デモ→ 森林火災(=サイト問題) 浸水(=ボンド問題)
(1) まず,L×Lの正方形の格子点に,pの割合(つまり全部でpL2本)でランダムに木をばらまく。[ Lは大きいほどよいが,L=200 くらいが限度であろう。]
(2) 最初にある木を任意に一本選んで発火させる(=値を1にする)。
(3) 上下左右斜め上下の隣接する8箇所のうち,1本でも発火している木があったら,次のステップで自分も発火するルールを与える。もちろん全ての木を同時的に更新。
(4) 燃えている木の数が一定になるまでステップを繰り返す。[ヒント:燃えている木の総数は減ることはなく,また同じ値が2回続けばそれ以上増えることはない。]
(5) 定常に達したときの燃えている木の占める率(最初にばらまいた本数に対する率)とpの関係を求め,これが0から1へと急激に変わるpのおよその値(自然鎮火の限界値)を見いだせ。特に 0.35 と 0.45 の間はくわしく調べること。各pの値に対して木の配置を変えて10回くらい繰り返して平均を取ること。
隣接点を上下左右のみにした例:限界値は上の場合と異なり,0.59 程度。
(ヒント) 上下辺上および右端の処理だけプログラム上は特別にしなければならないが,用意する配列を実際の森の大きさより1つずつ大きめにとり,ここは木が生えていないようにしておけば簡単化できる。なお,木が生えていることを表す配列と,燃えていることを表す配列は別にしておく方がプログラムしやすいであろう。
<課題5−3>
(相分離モデル) 上下左右斜め上下と自分をあわせて9つの格子点のうち
・5つ以上が「↑」であれば「↑」は生き続け,「↓」は「↑」に転じ
・4以下の場合は「↓」は「↓」のままで,「↑」は「↓」に転じる。
100×100程度の格子点に最初に半数ずつ「↑」と「↓」をランダムに配置しておき,このルールを適用してしばらくたった後の時刻におけるパターンをいくつか図示せよ。(定常パターンにはならない。)[画面や印刷にできない場合は,表示可能な範囲に切り取って示すこと。可能ならグラフィックスで上下左右が同じ尺度,すなわち正方形になるようにし,(0,1) の代わりに塗りつぶしパターンで描く方が望ましい。]