数理解析研究所講究録1032 巻 1998 年 86-100

86Langevin 型拡散方程式 (KPZ 方程式) の数値解法
齊藤善弘 (聖徳学園女子短期大学)
新宮康平 ((株) 電通国際情報サービス)
三井斌友 (名古屋大学人間情報学研究科)

1

はじめに

拡散現象を記述するモデルの中に確率拡散方程式がある. これは拡散方程式に確率的なノイズを加えたモデルであって, 特に微粒子や量子の動き, 界面の現象など, ノイズに左右されやすい現象を記述することができる. そのような現象の具体例をあげれば, 障子紙の下端を黒インクにつけ, 立てておくと, インクが紙に染み込んでいき, 黒インクのぎざ
ぎざの境界線が少しずつ形を変えながら, 次第に上昇していくような現象をあげることができる. これは, 基盤が 1 次元のときの “荒れた成長界面” に関する現象とといわれる.
現在, このような現象のダイナミックスについての研究が非常に盛んである. その理由の–つとして, 臨界現象の研究において, その重要性が明らかになったスケーリング則(scaling rule) が成立し, 定量的な理解が可能であることが挙げられる $[2, 7]$ . すなわち時刻 とシステムの大きさ (length scale of the system)
における界面の乱れの平均 $(w_{L}(t))$
は, 次のようなスケーリング則に従うことが知られている [4].
$L$

$t$

(1.1)

$\langle w_{L}(t)\rangle\simeq L^{\alpha}f(t/L^{z}.)$

ここでスケーリング函数 $f(x)$ は$f(x)\simeq x^{\beta}$$f(x)\simeq$

const.

( $x\ll 1$ のとき)
( $x\gg 1$ のとき)を満たし, かつ $z=\alpha/\beta$ であるような函数である.

しかも, 粗さ指数$\alpha$

と動的指数

$z$

は,

空間次元や秩序変数の自由度など非常に少数のパラメータに依存していると考えられてと 2 は互いに独立ではなく, 基盤が 1 次元の場合に\alpha +z $=2$ という関係式を満たすことがわかっている. この事実から, 基盤の次元に関わりなく, $\alpha+z=2$ といる.

また$\alpha$いう関係式を満たすと予想されている.
当初このような荒れた成長界面の研究には, 荒れた界面をアルゴリズムによって記述する ballistic deposition のような “離散モデル” が用いられてきた. その後 Kardar , Parisi,
Zhang はこのような界面を記述する連続体モデルを提案し, この研究に飛躍的前進をもたらした [8]. それは界面の運動を Langevin 型の確率拡散方程式で表すもので, これに87より初めて界面の運動を解析的に扱うことが可能となった. この方程式は著者たちの頭文字を取って “KPZ 方程式” と呼ばれる. Kardar たちは動的くりこみ群の方法を適用して, 基盤の次元が 1 次元 $(d= 1)$ のときの指数の値を求めた. このときに得られた値$(\alpha=1/2, \beta=1/3)$ は離散モデルからの結果と–致し, KPZ 方程式がこのような界面の本質を捉えたものであることを明らかにした. しかし $d\geq 2$ のときには解析的には指数の値は求まらないので, $d\geq 2$ での KPZ 方程式に対しては, 離散的な数値解法が必要になってくる. なお動的くりこみ群の理論によると, $d=3$ のときには非線型項の係数 の値により相転移が起こり, slnootllと roughの二つの異なる相が存在すると予想されている.
本論文では, この KPZ 方程式を数値的に解く際の問題点を含め, 数値シミュレーショ
ンの結果について述べる. KPZ 方程式も含めて, 殆どの非線型な方程式では解析解は困難であり, 何らかの数値解法を使わざるを得ない. 従来用いられてきた数値解法の多くは,偏微分方程式の前進差分スキームに確率的ノイズを加えるものである. 方程式の非線型性のため, 陽的解法が適用しやすいが, 陽的解法は数値的不安定性を生じる可能性がある.
そこで, 陽的解法ではあるが安定性を考慮した数値スキームを提案し, 数値的不安定性が緩和されることを示す. まず, 2 節において KPZ 方程式について簡単に紹介する. 3 節
では, 加法的ノイズをもつ確率常微分方程式 (Stochastic Ordinary Differential Equation
以下 SODE と略す) の数値解法について述べる. 4 節では, KPZ 方程式に対する数値シ
ミュレーションの結果を与える. また, 5 節では数値的不安定性を生じた原因について,
線型安定性の観点から解析を行う. 最後に今後の課題について述べる (6 節).
$\lambda$$\mathrm{p}1_{1}\mathrm{a}\mathrm{s}\mathrm{e}(\alpha=0)$2$\mathrm{p}\mathrm{h}_{\mathrm{a}}\mathrm{S}\mathrm{e}(\alpha\neq 0)$

KPZ 方程式
本節では, 不可逆的に発展しながら形成される界面の運動を記述するモデル方程式と,その空間離散近似を述べる. I く ardar, Parisi, Zhang は 次元のある基盤からの高さ
動きを表現するために, 次のような Langevin $\text{型の拡散方程}x$ を提案した [8].
.
$d$

(2.1)

$h$

の

$\frac{\partial h(r,t)}{\partial t}=\nu\nabla^{2}/l(r, t)+\frac{1}{2}\lambda(\nabla h(r, t))^{2}+\eta(r, t)$

ここで, $h(r, t)$ は界面の高さを表す. すなわち,次元基盤上の位置 $r(\in[0, L]^{d})$ と時刻における基盤からの高さを意味する.
は正の定数である. KPZ 方程式の右辺において, 第–項は表面張力によって界面が滑らかになろうとする傾向 (空間方向の拡散効果)を, 第二項は局所的な成長の度合いを表している. そして, 最後の項\eta \sim , ) はランダムな力による効果を表し, 次の性質をみたすガウス型白色ノイズである.

$t$

$d$

$\lambda$$\nu,$$t$(2.2)$\langle\eta(r, \mathrm{f})\rangle$

$=$

$0$

$\langle 7’(r, t)7’(r, t’J)\rangle$$=$$2D\delta(\{(r-r)J\delta(t-t’)$

$\{$

88

は標本空間上での平均値を表す.
界面の粗さ (surface width) は次のような尺度ではかる.

ここで,

(2.3)

$D$

はノイズの強さを,$\langle$$\cdot)$$w_{L}(t\mathrm{I}=(\overline{l^{2}\iota(r,t)}-(\overline{l\iota(r,t)})2)1/2$.

ここで ; は時間$t$を固定したときの空間平均を表す. つまり,
$t\backslash \overline{)}=\frac{1}{L^{d}}\int_{[0,L]}d^{/_{l}t})(r,dr$$l_{l(r},$は式 (2.3) が示すように, 界面の高さ $h(r, t)$ の平均自乗偏差である.
個の 1次元の空間をさて, KPZ 方程式 (2.1) に対して数値解法を導こう. まず,の等軸格子上のセルに分割する. 次に (2.1) 式の右辺の空間微分を, 格子定数辺である. また$w_{L}$$L^{d}$$d$$\Delta x$$\triangle x$.

での中心差分によって離散化すると次式を得る.
(2.4)

.

$\frac{\mathrm{d}.l_{l}n(t)}{(1t}$$=. \frac{1}{(\triangle\backslash l\mathrm{i})2}.‘\sum_{i=1}^{d}\{\nu[\dot{l}_{l}n+e.\cdot(t)-2h_{n}(t)+hn-e.\cdot(t)]$$+ \frac{1}{\mathrm{S}}\lambda[/ln+e,(t)-hn-^{e}.\cdot(t)]^{2}\}+\eta’(t)\cross\frac{1}{\sqrt{(\triangle x)^{d}}}$ここで, 格子点の位置を整数ベクト) $n=(r\iota 1, tl2, \cdots, rld)$ , 基本ベクトルを$e_{1},$ $\cdots,$ $e_{d}$で$\langle\eta’(t)\gamma|’(tJ)\rangle=2D\delta(t-t’)$ を満たすガウス型表すことにする. また, $7l’(t)$ は方程式, つまり加法的ノイズをもつ SODE になっ
白色ノイズを表す. 式 (2.4) はている. 論文 [9] では, 式 (2.4) の時間方向の積分に対し Euler スキームを使用して解いているが, その安定性は限界的である. 実際われわれの行った計算によっても, 次元数 ,定数\mu , の採り方により数値的不安定性を示すことがわかった (Fig. 1). そこで, Euler
スキームより高次で, より安定な数値スキームを提案する. 次客では SODE の数値解法とその基礎概念を述べ, 4 節でこれらの数値スキームを KPZ 方程式に適用した結果を示
そう.
$\langle\eta’(t)\rangle=0,$$\mathrm{L}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{g}\mathrm{e}\mathrm{v}\mathrm{i}\dot{\mathrm{n}}$$d$$\lambda$

SODE の数値解法

3

Langevin 方程式 (2.4) を簡単のため,(3.1)$\frac{\mathrm{d}_{\lrcorner}\lambda’(t)}{\mathrm{d}t}=f(X)+g_{7}l(t)$は定数で,と書くことにしよう. ここで $X\in R^{d’},$ $f$ は d’ 次元ベクトル函数,
は白色ノイズである. 白色ノイズ\eta (t) は, 次式のように Wiener 過程 $W(t)$ の形式的な導函数で$g$表現できる.
$\uparrow l(t)=\frac{\mathrm{d}\nu V(t)}{\mathrm{d}t}$$\eta(t)$

89

よって,

Langevin 方程式 (3.1) は, 次のような加法的ノイズをもつ SODE$\mathrm{d}X(t)=f(X)\mathrm{d}t+g\mathrm{d}W(t),$(3.2)$t\in[0, T]$

,

$\{$

$X(0)=x$ ,によって表される.

ここで,$W(t)$ は d’ 次元標準

Wiener 過程である.

SODE(3.2) に対して, 様々な時間離散近似解法が提案されているが, これらの時間離散近似解法は大きく分けて二つの種類に, すなわち強い意味の近似 (strong approximation)と弱い意味の近似 (weak approximation) に分けられる. 強い意味の近似は, 解 $X(t)$ の軌
道 (sample path) を最適に近似する解法である. 弱い意味の近似は, 解 $X(t)$ の統計量, すなわち平均や分散を求めるのに効率的な解法である.
KPZ 方程式 (21) を解くとき, 物理的に有意義な量は高さ の平均自乗偏差の平均値であるから, 弱い意味の数値スキームを使うべきである. 以降本論文では弱い意味の数値ス
キームに限定する. また, 偏微分方程式の空間離散化を細かくするに伴い, すなわちを大きくするのに伴い, (3.2) の次元
を小さくするかが増加し, 計算量が増えるので, 函数計算が複雑な数値スキームは好ましくない. また, KPZ 方程式は加法的ノイズ$h$$\Delta x$

$d’$

$L$

をもつので, 数値スキームの次数に関する制約条件は,で, 低次で函数計算の手間をとらない,

通常より幾分緩和される. そこ

Runge-Kutta 型の数値スキームを用いよう.

まず, 収束次数 1 の Euler スキームは次のように与えられる. 以降上付き添え字は, そのベク 1\) の成分を表す.

(3.3)$A\overline{\lambda^{\prime i}}_{7l+}1=\lambda_{n}^{\overline{r}\iota’}+f^{i}n+g\triangle t_{n}\triangle\nu V^{i}n$.

次に Runge-Kutta スキームの–つである, Heun スキームをあげる. これは 2 段 Runge法) の SODE 版である $[3, 5]$ .
Kutta 法 ($\mathrm{H}\mathrm{e}\mathrm{U}11$(3.4)$X \text{晶}=_{d}\overline{\lambda}_{7l}^{-i}+\frac{1}{2}[F_{1}^{i}+Fi|2+gW^{i}\triangle t_{n}\triangle n$’

ここで,$l_{2}^{i}\overline{r}^{1}=f^{i}(.\overline{\lambda}_{n}r+F_{1}\triangle t_{7l}+g\triangle Wn)$

$F_{1}^{i}=]^{i}.(\lrcorner\overline{\mathrm{x}}_{7}^{\vee})l$’.

Heun スキーム (3.4) は, 加法的ノイズの方程式に対して弱い意味で 2 次を達成することができる.

Heun スキームは,$t=t_{n+1}$ でEuler スキーム (3.3) で予測値 X-41 を求め, さらに台形

公式$\lambda_{n+1}^{i}\vee=\overline{\lambda}_{7l^{+}}^{ri}\wedge\overline{2}\perp[].i(\lambda_{?\iota}\overline{\prime})+fi(\lambda_{n+1}^{\overline{\prime}})]\triangle t_{n}+g\Delta W_{n}i$によって修正して,を求めるスキームとみることができる. 台形公式は SODEs に対して安定性に優れているとみられているので [11], この点に着目し,スキームで予$\overline{\lambda_{n+1}’}$$\mathrm{E}\mathrm{u}\mathrm{l}\dot{\mathrm{e}}\mathrm{r}$90

測し, 台形公式で 2 回修正するスキームを提案しよう.

(3.5)$\overline{\lambda}^{ri}=n+\iota\cdot 7l\overline{\mathrm{x}^{r}}i\frac{1}{2}+[F_{1}^{i}+Fi]3g\Delta\triangle t_{n}+\nu V^{i}n$.

ここで,$F_{2}^{i}.=f^{i}(X_{l}^{\overline{\prime}},+F_{1}\triangle t_{?l}\cdot+g\triangle\nu V.\iota)$$F_{1}^{i}=f^{i}(\overline{\lambda}_{?}^{-})1$’,$F_{3}^{i}=f^{i}(-n+ \frac{1}{2}[F_{1}+F_{2}]\Delta t_{n}+g\Delta W_{n})$.

上述の理由から, スキーム (3.5) をスキームと呼ぼう.
スキームは弱い意味で収束
次数が 2 であるが, Heun スキームより数値的安定性が優れていると期待できる.
$\mathrm{P}\mathrm{C}$次に$\mathrm{P}\mathrm{C}$$\mathrm{P}\mathrm{C}$スキームと同じ 3 段の Runge-Kutta スキームをあげる $[12, 13]$ .

(3.6)$\lambda_{n+1}^{\overline{r}i}=-\overline{\lambda}_{r1}’i+\frac{1}{4}[F_{1}^{i}+3\Gamma_{3}^{i}\triangle t_{n}+g\triangle W_{n}^{i}$ここで$F_{1}^{i}=f^{i}.\cdot(\lrcorner\lambda_{n}^{\overline{\prime}})$,$\Gamma_{2}^{i}’=f^{\dot{\iota}}(\Delta\lambda^{\overline{r}_{7\mathrm{t}}}+.\frac{1}{3}\int\tau.,\triangle\downarrow t,+‘\frac{1}{3}\mathit{9}\triangle \mathrm{t}.W_{7}1.)$,.

$l_{s^{i}}^{\^{\urcorner}=}.f^{\dot{l}}(-n+ \frac{2}{3}F_{2}\Delta tn+\frac{2}{3}g\Delta W_{n})$スキーム (3.6) は 3 段 Heun スキームと呼ばれ, 収束次数は Heun スキームと同じで, 弱い意味で 2 次である.
スキーム (3.5), そして 3 段 Heun スキー
Euler スキーム (3.3), Heun スキ一ム (3.4),
ム (3.6) を適用するとき, Wiener 過程の増分\triangle I/\iota \mbox{\boldmath $\gamma$}7i\iota は, 平均 , 分散 1 の正規乱数\xi S を使っ$\mathrm{P}\mathrm{C}$$0$て次のように実現すべきである.
$\triangle W_{n}^{i}=\xi_{n}^{i}\sqrt{\triangle t_{n}}$しかしながら, 平均や分散などを求める弱い意味の近似に対しては,

正規乱数\xi n を次ののいずれかで置き換えても, 精度が充分保証されることが示されている $[5]$ .
以下 $U$ は $[0,1)$ 上の–様乱数を示す.
弱い意味で 1 次の数値スキームに対しては$X\text{うな乱数_{}\hat{\xi}_{n}}$(3.7)$\hat{\xi}_{n}=\sqrt{12}(U-0.5)$

または

(3.8)
$(3.7),(3.8)$ に現れる

1

$U< \frac{1}{2}$$-1$$\frac{1}{2}\leq U$$\hat{\xi}_{71}=\{$$\hat{\xi}_{n}$は, 共通して次の性質を満たす.
$\langle\hat{\xi}_{n}\rangle=\langle\hat{\xi}_{n}^{3}\rangle=0,$$\langle\hat{\xi}_{n}^{2})=1$

91

弱い意味で 2 次の数値スキームに対しては$-\sqrt{3}$(3.9)$\hat{\xi}_{n}=$$\backslash$$U< \frac{1}{6}$$0$$\frac{1}{6}\leq U<\frac{5}{6}$-I$\frac{5}{6}\leq U$が使われる. 今度は, \xi 7’ は次の性質を満たす.
$\langle\hat{\xi}_{n}\rangle=\langle\hat{\xi}7\iota 3\rangle=\langle\hat{\xi}_{n}^{5}\rangle=0,$4$\langle\hat{\xi}_{n}^{2}\rangle=1,$$\langle\hat{\xi}_{n}^{4}\rangle--3$

KPZ 方程式の数値シミュレーション

本節では, KPZ 方程式に前節で述べた数値スキームを適用した結果を述べる. まず,
KPZ 方程式 (2.1) を空間離散化した式 (2.4) に Euler スキームを適用する. 1 次元の場合に, (2.4) に現れるパラメータを次のように設定する.
$d=1$ ,

$L=40$ ,$\Delta x=1$,$\nu=0.5$,

D- $=0.005$

後はパラメータ の設定を変えることで, 方程式を特徴づけることができる. 時間ステッ
プ幅は等間隔に $\triangle t_{n}=0.02$ と沖時刻 $t=10\dot{0}000$ まで計算した.
, 以降の計算に用いた計算機
Sun 製 SPARC Station 20 (Sun) である.
最初に\mbox{\boldmath $\lambda$}2 $=1500$ の場合, Euler スキ一ムを適用した結果を Fig. 1 に示すと, 10 本中 8本の軌道が発散する. 次に, 同じ設定で Heun スキームを適用したときの結果が Fig. 2 である. 今度は解が爆発せず, 数値解が安定に求まっていることがわかる.
次に 2 次元の場合に対して行ってみよう. パラメータは$\lambda$$\text{また}$$\mathrm{o}\mathrm{s}4.1.4$$|\mathrm{h},$,$d=2$ ,

$L=40$ ,$\triangle x=1$,$\nu=0.5$ ,$D=$ 0.005とする. ステップ幅は等間隔に $\triangle t_{n}=0.02$ とし, 時刻 $t=$ 100000 まで Euler スキームによって計算した.

$\lambda^{2}=1500$の場合, 1 次元のときに現れた数値的不安定性がみられない(Fig. 3).
今度は 1 次元で\mbox{\boldmath $\lambda$}の値を変えて, Heun スキームの数値的安定性の限界を調べよう.
1750 及び\mbox{\boldmath $\lambda$}2 $=2000$ のとき, Heun スキームを適用した結果を Fig. 4 及び Fig. 5 に示す. これよりの標本 10 本のうち 4 本あるいは 8 本が爆発した. よって, \mbox{\boldmath $\lambda$}が増えるに従い, Heun
スキームの不安定性が増すことがわかる. 同様に,スキーム及び 3 段 Runge-Kutta ス
キームで行った結果を表にまとめたのが Table 1 である. これより, Heun スキームより$\lambda^{2}=$$w$$\mathrm{P}\mathrm{C}$92

Table. 1. The nulnbel of explosive sample paths for KPZ equation with several values of$\cdot$$\lambda$

スキームが安定性において優れていることがわかる. 残念ながら, この表からは 3 段
スキームより良い結果を与えていない. .,.
Runge-Kutta スキームは以上の結果をまとめると, Euler スキームより Heun スキーム, Heun スキームより PC
スキームが数値的安定性において優れていることがわかる. また, 次元によって数値的不安定性の発生する状況が違う. そこで, 次節では数値解の不安定性の原因を, KPZ 方程式を半離散化した確率常微分方程式に線型化を施し, 固有値分布を調べることにより追求
する.
$\mathrm{P}\mathrm{C}$$\mathrm{P}\mathrm{C}$5

KPZ 方程式に対する線型安定性解析

加法的ノイズをもつ確率常微分方程式の安定性解析は文献[6] に詳しい. 結論からいう.
と, 加法的ノイズの場合は方程式のトリブト項のみ調べればよく, 常微分方程式の場合の線型安定性解析と同様に扱うことが出来る.
最初に 1 次元の場合を考察する.
方程式を $L=40,$ $\triangle x=1$ のとき半離散化した方程式 (2.4) のドリフト項に対して, 周期境界条件を考慮して, ヤコビ行列を計算すると$1\backslash \mathrm{P}\mathrm{Z}$

(5.1)

を得る.

ここで,

$\alpha_{i}---2\nu$,$\beta_{i}=\nu+\lambda[hi+1(t)-/li-1(t)]/4$,$\gamma_{i}=\nu-\lambda[hi+1(t)-h_{i}-1(t)]/4$

93

であって,
$2p_{i}$.

$w_{L}(t)$$l\iota_{0}=h_{40},$$h_{41}=/\iota_{1}$と仮定しよう.

と解釈する (周期境界条件). さて,$p_{i}=\{$

1
$-1$

とする.

$h_{:+1}(t)-h_{i-1}(t)=$ここで$( \text{確率}\frac{1}{2})$$( \text{確^{}\mathrm{R}}4\frac{1}{2})$

この仮定の根拠は, 次式が成り立つことによる.
$\overline{|l\iota_{i}(t)-/l_{j}(t)|}\leq\sqrt{\overline{(l\iota_{i}(t)-l\iota j(t))^{2}}}\leq 2w_{L}(t)$そこで, 線型化行列は行列 (5.1) において$\alpha_{i}=-2\nu$,$\beta_{i}=\nu+\lambda piwL(t)/2$,$\gamma_{i}=\nu-\lambda p:w_{L}(t)/2$と置き換えた行列となる.
1500, 1750, 2000 についてみるが, 不安定になる時の の値は, いずれの場合においても $w=0.‘ \mathit{2}0$ 程度である. よって, これらの値を線型化行列 (5.1) に代入し, 固有値の分布を複素平面上に函示したのが, Fig. 6 である. これより, \mbox{\boldmath $\lambda$}が増えれば, 固有値の分
布は $(-1,0)$ を中心に縦横に伸びることがわかる. Euler スキームと Heun スキームの絶対安定領域の有界性によって, \mbox{\boldmath $\lambda$}が大きくなれば, ステップ幅を相当に小さくしなければならないことがわかる. さらに, 固有値の分布は本来各ステップ毎に様相が変化することを考慮にいれると, 考察している方程式系では数値的不安定性が生じやすいことがわかる.
次に, 2 次元の場合を見ることにしよう. ヤコビ行列を計算するのに 1 次元のときと同じ仮定, すなわち$\lambda^{2}=$$w$$/_{\iota_{(i,j)}(t})-/_{l}(k,l)(t\mathrm{I}=2pi,j ..\omega_{L}(t)$をとる. ただし p 痘は 1 次元のときの$p_{i}$と同じ性質をもつ確率変数である. 線型化行列はブロックで表示すると

$A=$

94

となる.

ここで行列

$A_{i},$$B_{i}$, Ci は次のような行列である.

$A_{i}=$
diag$B_{i}=$$[\beta_{i,1}’, \ldots,\beta_{i,40}’]$,,$Ci=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[\gamma_{i,1}’, \ldots, \gamma_{i,40}’]$,

また各定数は$\alpha_{i,j}=-2_{\mathcal{U}}$,$\beta_{i,j}=\nu+\lambda Pi,jwL(t)/2$$\beta_{i,j}’=\nu+\lambda l^{J_{i}’},jw_{L}(t)/2$,,$\gamma_{i,j}=\nu-\lambda p_{i},jw_{L}(t)/2$$\gamma_{i,j}=\nu-\lambda p’i,jL(\prime tw)/2$である. さて, 1 次元の場合と同じ $\lambda^{2}=15.00$ ,1750, 2000 に対して, Euler スキームは不安定現象を生じなかった. そこで, いずれの場合でも $w=0.10$ 程度であった定常時の $w$の値を線型化行列に代入し, 各々の の値に対し, 固有値分布を計算し, 図示すると Fig.
7, 8, 9 になる. これらの図は, 1 次元とは対照的に, 円に近い分布を示している. また$\lambda$$w$の値も 1 次元の場合の半分の 0.10 程度であることも加えて, 2 次元の場合は数値的に安定に解けると考えられる. 本節での行列の固有値計算には, NUMPAC の HEQRVD を使用
した.

6

今後の課題

本論文では確率拡散方程式, 特に KPZ 方程式に対して陽的な Runge-Kutta 型数値ス
キームによる数値解法について考察した. 1 次元の場合, Euler スキームに比べ 2 段 Heun
スキームの方が数値解を安定に与える. しかし, Heun スキームも陽的スキームであるゆえ,の値が大きくなると数値解が不安定になる割合は高くなる. 同じ理由から,ス
キーム及び 3 段 Runge-Kutta スキームも, 安定性の面で Heun スキームより優るなどの特筆すべき結果は得られなかった. また, 第 5 節で数値解法の線型安定性解析を試みた.
確率的な要素が加わるため, 大胆な仮定をしたが, 1 次元と 2 次元の場合について, 不安
定性の–つの要因に触れることができた. さらに Heun スキームの使用により, Euler ス
キームでは困難だった\mbox{\boldmath $\lambda$}の高い値の KPZ 方程式について計算し,を求めることが可能になり, 統計力学の観点から興味深い結果を得ることができた.
今後は, 確率偏微分方程式の数値計算するにあたり, 確率常微分方程式においても同様であるが, 計算コストの軽減に注目せざるを得ない. 具体的には, 確率偏微分方程式は空$\lambda$$\mathrm{P}\mathrm{C}$$\langle$$w_{L})$

95

間離散化を行うことにより, 高次元の確率常微分方程式系に置き換わる. より精度よく求めようとするならば, 計算効率を上げるために工夫するしかない. よって, 標本 (軌道)を求めるアルゴリズムの並列化と時間方向の積分の効率化による計算コストの軽減を推し進める必要がある. そのためには並列計算用の乱数発生器の開発や, 自動ステップ幅調節機能をもつ数値スキームの開発が今後の課題といえよう. また, KPZ 方程式のような非線型方程式に対し, 陰的スキームは陽的スキームと比べ, 適用が難しいが, 実装化を考慮し, 安定性を保持する工夫も考えていかなければならない. 常微分方程式の数値解法で行われている研究に比べ, 残された課題は多いといえる.

謝辞 KPZ 方程式およびその物理的意義について教示いただき, また有益な討論をいた
だいた本田勝也氏 (信州大学・理学部) に深謝する.

参考文献[1] Alllal , J. and Family, F.,solution of a continium equation for interface$41$ (1990),
growth in 2+1 dimensions, Phys. Rev.
3399-3402.
$\cdot$$\mathrm{N}_{\mathrm{l}\mathrm{I}11\mathrm{e}}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}1$$\mathrm{A},$[2] Barab\’asi, A-L. and Stanley, 1-1.E., Fractal Concepts in Surface Growth, Cambridge
University Press, 1995.
[3] Chang, C. C., Nunlerical solution of stochastic differential equations with constantdiffusion coefficients, Math. Comp., 49(1987), 523-542.
[4] Family, F. and Vicsek, T. Scaling ofactive zone in the Eden process on percolation
networks and the ballistic deposition nlodel, J. Phys.
18(1985),.
$\mathrm{t}1_{1}\mathrm{e}$$A_{f}$[5]$\mathrm{G}1^{\cdot}\mathrm{e}\mathrm{i}_{1}1\mathrm{e}\mathrm{r}$, A.,$\mathrm{S}\mathrm{t}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{t}111\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{c}^{\mathrm{J}}1^{\cdot},$$\backslash ’\backslash ;.$tic differential equations., .
$i$

, and
$6’t$

$\mathrm{L}75-\mathrm{L}81$, J., Numerical integration of stochasatist. Physics, 51(1987), 95-108.
$\mathrm{I}^{-}1_{0}11\mathrm{e}\mathrm{r}\mathrm{k}\mathrm{a}1111^{)}$[6] Hernandez,. and Spigler, R., A-stability of Runge-Kutta Methods for systems
with additive noise, BIT,
, 620-633.
$\mathrm{D}.\mathrm{B}$$32(1^{\iota\underline{\cdot)}}\mathrm{J}9)$[7] 本田勝也, 荒れた界面のダイナミックス, 日本物理学会誌, $49(1994)$ , 819-826.
[8] Kardar, M., Parisi, . and, Y.-C., Dyllanlic scaling of growing interfaces,
Phys. Rev. Lett., 56 (1986), 342-345.
$\mathrm{G}^{1}$$r\Delta 1_{1\mathrm{a}}1\mathrm{l}\mathrm{g}$[ Wolf, D. E., Nulnerical solutions of the[9] Moser, K., Kert\’esz, J.
Kardar-ParisiZhang equation in one, two and three dinlensions, Physica $A,$ $178$ (1991), 215-226.
$\mathrm{a}\iota\iota \mathrm{C}$96[10] Moser, K. and Wolf, D.E., Vectorized and parallel simulations of the Kardar-ParisiZhang equation in 3+1, J. Phys. : Math. Gen., 27(1994), 4049-4054.
$A$$\mathrm{d}\mathrm{i}\mathrm{l}\mathrm{l}\mathrm{l}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{i}\backslash ,\mathrm{j}o1\mathrm{l}\mathrm{S}$[11] Pardoux, E., and Talay, D., Discretization and simulation of stochastic differential
equations, Acta Ap..pl. Math., $3(1985).$ ’ 23-47.
[12] 齊藤善弘三井斌友, 確率微分方程式の離散近似, 日本応用数理学会論文誌, 2(1992),

1-16.

Fig. 1. Euler$\mathrm{s}\mathrm{c}\mathrm{h}\mathrm{e}\mathrm{l}\mathrm{l}\iota \mathrm{e},$$d=1,$$\lambda^{2}=$1500, no. of samples $=10$ .

97Fig. 2. Heull scllellle, $d=1,$

化g$\lambda^{2}=$1500, no. of samples $=10$ .

$\lambda^{2}=$1500, no. of samples $=10$ .

$\mathrm{w}$Fig. 3. Euler scllellle, $d=2,$

98Fig. 4. Heun$\mathrm{S}\mathrm{C}1_{1\mathrm{e}\mathrm{n}}1\mathrm{e},$$d=1,$

Fig. 5 Heun scheme, $d=1,$$\lambda^{2}=$$\lambda^{2}=$1750, no. of samples $=10$ .

2000, no. of samples $=10$ .

99

1$\mathrm{m}$$\lambda^{2}=1500$$\lambda^{2}=\downarrow 750$$\lambda^{2}=20\infty$

Fig. 6. Eigenvalue distribution, $d=1$

1$\mathrm{m}$Fig. 7. Eigenvalue distribution, $d=2,$$\lambda^{2}=1500$

100

.

’

Fig. 8. Eigellvalue distfibution, $d=2,$ $\lambda^{2}=1750$

Fig. 9.

$\mathrm{E}\mathrm{i}\mathrm{g}\mathrm{e}\iota 1\mathrm{v}\mathrm{a}1_{\mathrm{U}}\mathrm{e}$distribution, $d=2,$ $\lambda^{2}=2000$