186

数理解析研究所講究録
第2076巻 2018年 186-196

3次元スクロール波フイラメントの負張力状態の数理モデル
筑波大学・名誉教授1, 統計数理研究所2

金野秀敏1,2, 田村義保2

Hidetoshi Konno and Yoshiyasu Tamura
University of Tsukuba, Tsukuba
Institute of Statistical Mathematics, Tachikawa

1

はじめに

近年,位相特異点のダイナミクスの解析が益々注目されるようになってきている.心室細動時の回転らせん波が生成死滅を繰り返す状態では,位相特異点も生成死滅過程が現れることはよく知られているが,数多くの非平衡系

で位相特異点の生成死滅過程が知られている: (1) 心室細動 [1, 2, 3, 4, 5, 6], (2) 傾斜のある系での RB 対流や回転するシリンダー中の \mathrm{R}\mathrm{B} 対流系の欠陥 [7], (3) Hepta‐Penta 欠陥カオス [8], (4) 液晶における欠陥乱流 [9], (5) 生
物系での大変形や移動のダイナミクス [10], (6) 量子乱流 [11], (7) プラチナ表面上での一酸化炭素の酸化反応 [12]など枚挙にいとまがない.

位相特異点のダイナミクスを記憶のない(Markov)マスター方程式を使って記述する方法の開発は Gil, Lega,
Menuer [13] によってなされた.彼らは空間2次元の複素ギンスブルグ.ランダウ方程式 (Complex Ginzburg‐Landau
方程式(CGLE) のシミュレーションデータ (回転らせん波乱流状態) を用いて,位相特異点の生成率や死滅率を特定した$\lambda$_{n}= $\nu$+ $\lambda$ n, $\mu$_{n}= $\mu$ n+ $\mu$ n^{2}(1)相互作用を書き下すと次のように書ける:
X\rightarrow 0 (death); 0\rightarrow X (birth); 2X\rightarrow 0 (pair extinction); X\rightarrow 2X (wave break)(2)

しかし,非線形の生成死滅率 (対消滅) が内在するため,マスター方程式の平衡解のみが解析的に求まっているが,時間に依存した解は求められていない.従って,寿命分布の解析的表現なども得られていないし,長期記憶のある場合の解も求められていない.

空間1次元の複素ギンスブルグ.ランダウ方程式の振幅 (欠陥) 乱流領域 (Uchiyama and Konno) [14, 15] では
変調波 (Modulated Amphtude Wave, MAW), ホール (Hole,H), 欠陥 (Defect,D) の混在した複雑な非線形素励起が内在しており,各粒子の生成死滅則はボアソン分布に従うが,各粒子の寿命分布や出現の待ち時間には幕状の厚い裾があることが明らかにされた.各粒子の生成死滅則は次のように同定され,時間に依存する解や寿命分布さ
らには,特異点の速度の揺らぎの解析も実行できた.
$\lambda$_{n}= $\nu$

及び$\mu$_{n}= $\mu$ n .

(3)

単一の素励起 X(MAW, H, D) に関する生成死滅の確率法則 (周辺分布 P(\mathrm{X}) に関する法則) は次のように簡単なものとなるため,様々な時間依存解析が可能となったわけである.
X\rightarrow 0

及び 0\rightarrow \mathrm{X}(4)

これは空間1次元の特徴として位相特異点に 「右巻き」 と 「左巻き」 の区別は無いので同種の素励起の衝突による対消滅は存在しないからである.記憶効果は各粒子の寿命分布に顕著に現れるが,この記憶効果を導入するためには2つの方法があり,(a) $\nu$(t) 及び $\mu$(t) に畳み込み積分の無い時間依存係数を用いて記憶を導入する方法; 及び187(b)畳み込み積分で記憶を導入する方法がある.寿命分布の長期記憶の効果のみを記述するならどちらの定式化でも良いが,寿命分布の共鳴構造なども再現するためには (a) の定式化を使うしかない.
一方,空間3次元では空間1次元や2次元の系とは異なり位相特異点が (i) 点ではなく フィラメントの構造を形成している; (ii) リング状のフイラメント,曲がったフイラメント,ならびに端点が表面に固定されているフィラメ
ントが存在する.CGLE [16], Barkley model [3], 3\mathrm{V}‐SIM[17] などの計算では次のような3種類の異なる生成死滅率が報告されている:(A) 1種(hnear line tension, Reid, 2011 [2])(5)$\lambda$_{n}= $\lambda$ n+ $\nu,\ \mu$_{n}= $\mu$ n+ $\epsilon$(B) 2種(linear line tension , Clayton, 2009 [5])(6)$\lambda$_{n}= $\lambda$ n_{\rangle} $\mu$_{n}= $\mu$ n(C) 3種(scroll wave turbulence, Reid , 2011 [2])$\lambda$_{n}= $\nu,\ \mu$_{n}= $\mu$ n^{2}(7)(C) の空間3次元の乱流状態では主として対生成過程と対消滅過程がバランスしていると考えられる.
0\rightarrow 2X (pair creation); 2X\rightarrow 0 (pair extinction)

(B)

(8)

のClayton の計算したモデルでは次のような相互作用が支配的であると考えられる:
X\rightarrow 0 (death and loss); X\rightarrow 2X (wave break)(9)

おそらく,異なる形状の特異点の紐の間には相互転換が存在すると推定されるが詳しい解析の報告はなされていない.

空間2次元の相互作用で対消滅を無視して書き下すと,次のように書ける:
X\rightarrow 0 (death); 0\rightarrow X (birth);

X\rightarrow 2X (wave break)(10)

すなわち,吸い込み項が無い.一方,DavidsenやReidらの空間3次元の数値実験では 「吸い込み項は正値でなく,負値の吸い込み項」 としてパラメータが同定されている.

昨年度には (5) 式のような生成死滅率を持つ生成死滅過程の場合の確率母関数の解析はフラクショナル微分の形
で長期記憶が導入された場合でもホインの微分方程式の方法を使えば解析できることを報告した [18]. しかし,記憶効果がフラクショナル微分で表されるモデルの場合のマスター方程式の場合数学的な特異点の存在のために,確率質量,絶滅確率や待ち時間分布 (寿命分布) が求められないことに気がついた.そこで本報では通常のマルコフ系に話しを限定し,ポアソン変換とフォッカー プランク近似を組み合わせることにより近似解析解を得ることができ,また,絶滅確率や待ち時間 (寿命) 分布なども求められたので報告する.

2

マスター方程式

2.1

線形生成死滅率を持つモデル

時刻 t に n 粒子が存在する確率質量 P(n, t) に関するマスター方程式を考える:\displaystyle \frac{d}{dt}P(n, t)=[$\lambda$_{n-1}P(n-1, t)-$\lambda$_{n}P(n, t)]+[$\mu$_{n+1}P(n+1, t)-$\mu$_{n}P(n, t)] ,(11)

ただし,生成率 $\lambda$_{n} 及び死滅率 $\mu$_{n} は次式で与えられるとする:$\lambda$_{n}= $\lambda$ n+ $\nu$ 及び $\mu$_{n}= $\mu$ n+ $\epsilon$ .

ここで, $\lambda$, $\mu$, $\nu$,$\epsilon$(12)は定数である.初期条件は次式で与えられるとする:

P(n, 0)=$\delta$_{n.n_{\mathrm{O}}} .

(13)

188

(12) 式で $\epsilon$=0 の場合の応用例は多くの書籍や文献に見られる [19], Í20], [21], [22], [23]. 直交多項式を基礎とするマスター方程式 (11) の解の定式化はKerlin及びMcGregor[24] に報告されている.例えば煽 = $\lambda$ 及び $\mu$_{n}= $\mu$
なら,初期条件 P(n, t_{0})=$\delta$_{n,n_{0}} 下の確率質量 P(n, t|n0,t_{0}) は Chevishev 多項式; $\lambda$_{n}= $\nu$ 及び $\mu$_{n}= $\mu$ n な
ら Charier 多項式, $\lambda$_{n}= $\lambda$ n+ $\nu$ 及び $\mu$_{n}= $\mu$ n[25] なら Meixner 多項式で与えられることが知られている.その他の非線形生成死滅率の場合の解はIsmail[26], Sasaki [27] やSchoutens [28] などに与えられている.

特異点のダイナミクスを解析するマスター方程式を基礎としたアプローチは Gil ら [13] によって開発された.非線形の生成死滅率の場合の解析が行われているが,平衡分布鳥 (n) のみ解析解が得られているが時間依存解はほとんど得られていない.本論文では (11) 式のマスター方程式で線形の生成死滅率 (12) で平衡解鳥 (n) 存在する
( $\mu$> $\lambda$) の場合に焦点を絞った解析結果を紹介する.詳細釣り合いが成立しているので,
P_{s}(n+1)$\mu$_{n+1}=P_{s}(n)$\lambda$_{n} ,(14)平衡分布は次式のように求まる:

P_{8}(n)=\displayst le\frac{(\frac{$\nu$}{ \lambda$})_{n} (\frac{$\epsilon$}{ \mu$}+1)_{n}\frac{1} _{2}F_{1}( ,\frac{$\nu$}{ \lambda$},\frac{$\epsilon$}{ \mu$}+1;\frac{$\lambda$}{ \mu$})(\frac{$\lambda$}{ \mu$})^{n}(15)

ここで, {}_{2}F_{1} ([a, b], [\mathrm{c}];x) は超幾何関数, (x)_{n}(\displaystyle \equiv\frac{ $\Gamma$(x+n)}{ $\Gamma$(x)} はPoc—er の記号である.定常状態での平均,分散及び3次のキュムラントは次式で与えられる:\displaystyle\langlen(\infty)\rangle=\frac{$\nu$-$\epsilon$}{$\mu$-$\lambda$}, $\sigma$_{n}^{2}(\displayst le\infty)=\frac{$\nu\mu$- \epsilon\lambda$}{($\mu$- \lambda$)^{2} 及び $\kap a$_{3}(\displaystyle\infty)=\frac{($\lambda$+$\mu$)($\nu\mu$- \epsilon\lambda$)}{($\mu$- \lambda$)^{3} .
2.2

(16)

特性関数法による確率母関数の厳密解マスター方程式 (11) の母関数 (GF) [19, 20] は次式で与えられる :\displaystyle \frac{d}{dt}g(z,t)=( $\lambda$ z- $\mu$)(z-1)\frac{\partial}{\partial z}g(z,t)+( $\nu$-\frac{ $\epsilon$}{z})(z-1)g(z, t) ,(17)

ここで,母関数は次式で定義される :g(z, t)=\displaystyle \sum_{n=0}^{\infty}z^{n}P(n, t)

.

(18)

(13) 式の初期条件より, g(z,0)=z^{n_{\mathrm{O}}} となる.特性関数法を用いて確率母関数の解析解を得ることができる.特性方程式は次式で与えられる:\displaystyle \frac{dt}{1}=-\frac{dz}{(z-1)( $\lambda$ z- $\mu$)}=\frac{dg}{(z-1)( $\nu$- $\epsilon$/z)g} .

(19)

この解 g(z, t) は次式で与えられる:巽g(z, t)=z^{-\frac{ $\epsilon$}{ $\mu$} (1+\displaystyle \frac{ $\lambda$}{ $\lambda$- $\mu$}(1-z)[e^{( $\lambda$- $\mu$)t}-1])^{-} (1-\displaystyle\frac{(1-z)e^{($\lambda$-$\mu$)t} {1+\frac{$\lambda$}{$\lambda$-$\mu$}(1-z)[e^{($\lambda$-$\mu$)t}-1]})^{n_{0}\times\displaystyle \times(1+\frac{ $\lambda$}{ $\lambda$- $\mu$}(1-z)[e^{( $\lambda$- $\mu$)t}-1]-(1-z)e^{( $\lambda$- $\mu$)t})^{\frac{e}{ $\mu$}(20)

この解は Kerlin と McGregor [24, 25] ならびに Goel と Dyn [20] などの有名な確率過程の文献や書籍にも記述が無い.この表式が正しいことはモーメントを導出してみれば明らかであり, \langle n(t)\rangle(= g'(z, t)|_{z=1}) 及び分散$\sigma$_{n}(t)^{2}(=9''(z, t)|_{z=1}+g^{r}(z,t)|_{z=1}-(_{g}\prime(z, t)|_{z=1})^{2}) ( g' 及び g'' は z に関する1階及び2階微分を表す) は次式のように求められる :\displaystyle \langle n(t)\rangle=n_{0}e^{( $\lambda$- $\mu$)t}+\frac{ $\nu$- $\epsilon$}{ $\mu$- $\lambda$}(1-e^{( $\lambda$- $\mu$)t})

(21)

189

及び

$\sigma$_{n}(t)^{2}=\displaystyle\frac{$\mu$+$\lambda$}{$\lambda$-$\mu$}n_{0}e^{($\mu$-$\lambda$)t}(1-e^{($\lambda$-$\mu$)t})+\frac{$\nu\lambda$-$\epsilon\mu$}{($\lambda$-$\mu$)^{2} (1-e^{($\lambda$-$\mu$)t})^{2}+\frac{$\epsilon$+$\nu$}{$\mu$-$\lambda$}(1-e^{($\lambda$-$\mu$)t}) .

(22)

これらの時間依存モーメントは時間無限大で平衡値 (16) に収束する.
しかし,ここで問題が起こる.母関数が求まっているにもかかわらず確率質量 (PM),1 \partial^{n}P(n,t)=\overline{n!}\overline{\partial z^{n} g(z, t)|_{z=0} ,(23)

絶滅確率, P(0,t) , 及び待ち時間分布, f( $\tau$)=-\displaystyle \frac{d}{d $\tau$}P(0, $\tau$) などが (20) 式の特異点 z=0 の存在のために求まらない.そこで,近似解を導出することは意味がある.

3

近似解

3.1

ボアソン変換マスター方程式の近似解を得る方法は様々なものが存在する [19]. van Kampen[29] によるシステムサイズ展開が有名である.本報では,1変数に注目したモデルの解析に関する限りシステムサイズ展開よりも簡単なボアソン変換を採用した解析を行う.ボアソン変換(Poisson transform, PT) は次式で定義される [19]:P(n,t)=\displaystyle \frac{\langle n\rangle^{n} {n!}\int_{0}^{\infty}\exp(-\{n\rangle y)y^{n} $\psi$(y, t)dy ,(24)

このPT にはパラメータ \langle n\rangle が含まれており,ポアソン分布の平均値 (分散値) に対応する.(11) 式に逆ポアソン変換$\psi$(y,t)=\displaystyle \frac{\{n\rangle}{2 $\pi$ q}e^{\langle n\rangle y}\int_{-\infty}^{+\infty}\sum_{n=0}^{\infty}(\frac{ix}{q})P(n,t)e^{-ix_{q} dx4n

,

(25)

を適用することにより次のようなFokker‐Plan& 方程式 (FPE) が得られる [21],\displaystyle \frac{\partial}{\partial t} $\psi$(y, t)=-\frac{\partial}{\partial y}[K(y) $\psi$(y, t)]+\frac{\partial^{2} {\partial y^{2} [D(y) $\psi$(y, t)] ,(26)

ここで,回帰項と拡散項は次式で与えられる:

K(y)= $\alpha$- $\beta$ y 及び D(y)= $\gamma$+ $\delta$ y ,ここで,$\alpha$, $\beta$) $\gamma$ 及び $\delta$(27)は次式で与えられる:$\alpha$=\displayst le\frac{$\nu$- \epsilon$}{\n\rangle}, $\beta$= $\mu$- $\lambda$, $\gam a$=\displayst le\frac{$\epsilon$}{\n\rangle^{2} 及び $\delta$=\displaystle\frac{$\lambda$}{\langen\rangle} .

(28)

注意すべき点は,マスター方程式 (11) と FP 方程式 (26) は一対一ではなく,FP 方程式の3次より大きな高階微分の項を無視していることである :\displayst le\frac{$\epsilon$}{\langlen\rangle^{3}\frac{\partial^{3}{\partialy^{3} $\psi$(y,t)+\frac{$\epsilon$}{\langlen\rangle^{4}\frac{\partial^{4}{$\vartheta$y^{4} $\psi$(y,t)+\cdots(29)

展開式には $\epsilon$ 及び \langle n } が含まれているから,近似の精度に $\epsilon$ 及び \langle n\rangle が関係する.また,(ii) \langle n } はシステムサイズ
展開のシステムサイズと同様な役割を果たす [29]. 興味深いことに,FP 方程式に対応するランジュバン方程式は相乗性雑音を有するシステムに対応する:

dy=( $\alpha$- $\beta$ y)dt+\sqrt{ $\delta$ y+ $\gamma$}dW .

(30)

190

3.2

定常解

FPE を簡単な型にするために,線形変換を行う:

z=\displaystyle \frac{ $\gamma$}{ $\delta$}+y\equiv \mathrm{c}+y .

(31)

その結果,回帰項及び拡散項はそれぞれ K(z)= $\alpha$+ 穿 - $\beta$ z 及び D(z)= $\delta$ z , となる.この 座標での FPE は第z

1種のガンマ過程 [19] に帰着する.数理経済学の分野では,これは CIR モデル [22] と呼ばれている.
パラメータ $\alpha$, $\beta$, $\gamma$ 及び $\delta$ が正値のとき, t\rightarrow\infty で定常状態が出現するのは簡単に確かめられる.
常分布はガンマ分布となる$\psi$_{s}(z)\displaystyle \ap rox\frac{b^{a}z^{a-1} { $\Gamma$(a,bc)}\exp(-bz) ,

z

座標での定

(32)

ここで,

a=\displayst le\frac{$\alpha$}{$\delta$}+\frac{$\beta\gam a$}{ \delta$^{2} 及び b=\displayst le\frac{$\beta$}{ \delta$} ,(33)及び$\Gamma$(z,v)=\displaystyle \int_{v}^{\infty}u^{z-1}e^{-u}du(34)

は第2種の不完全ガンマ関数である. y\in[0, \infty] であるので z\in [c, \infty] となる.なぜなら, z は (31) 式で定義されているからである.従って,規格化因子は $\Gamma$ ( a , bc) であり,r(a) でない. c>0 に関する限り, $\Gamma$(a, bc)\approx $\Gamma$(a) と近似できる.ガンマ分布のパラメータ ( $\alpha$, $\beta,\ \gamma$, $\delta$) は生成死滅率 ( $\lambda,\ \mu$, $\nu$, $\epsilon$) 及びポアソン変換のパラメータく n\rangle で書けている.
z

座標を y 座標に戻すと$\psi$_{s}(y)\displaystyle \approx\frac{b^{a} { $\Gamma$(a)}[y+c]^{a-1}\exp(-b[y+c]) ,ここで,

y\in [0, \infty] , また,(35)a_{)}b 及び c は ( $\lambda$, $\mu$, $\nu$, $\epsilon$) を使って書く と

a=\displaystyle\frac{$\nu$-$\epsilon$}{$\lambda$}+\frac{$\epsilon$($\mu$-$\lambda$)}{$\lambda$^{2}, b=\displaystyle\frac{$\mu$- \lambda$}{$\lambda$}\langlen\rangle 及び c=\displayst le\frac{$\epsilon$}{ \lambda$\langlen\rangle} .

(36)

パラメータ ( $\lambda$, $\mu$, $\nu$, $\epsilon$) が正の定数であるとすると, $\mu$> $\lambda$ 及び $\nu$> $\epsilon$ なら定常状態 (分布) が存在する. $\psi$_{s}(y) に対応する定常確率質量は

P_{FP,s}(n)=\displaystyle \frac{\langle n\rangle^{n} {n!}\frac{b^{a} { $\Gamma$(a)}e^{-b\mathrm{c} \int_{0}^{\infty}e^{-(\langle n\rangle+b)y}y^{n}(y+c)^{a-1}dy,=\displaystyle \frac{e^{-bc} {n! $\Gamma$(a)}$\theta$^{a}(1- $\theta$)^{n} $\Gamma$(n+1,1-a;(\langle n\rangle+b)c) ,(37)

ここで, $\Gamma$(z, u;v) は一般化ガンマ関数であり次式で与えられる ;$\Gamma$(z,u;v)\displaystyle \equiv\int_{0}^{\infty}e^{-y}y^{z-1}(y+v)^{-u}dy ,(38)

ここで, $\theta$=\displaystyle \frac{b}{\langle n\rangle+b} と置いた.確率質量を ( $\lambda$, $\mu$, $\nu$, $\epsilon$) を用いて書く と

P_{F,s}(n)=\displaystle\frac{\exp[-\frac{$\epsilon$( \mu$- \lambda$)}{\lambda$^{2}]$\Gam a$(n+1,-\frac{$\nu$- \epsilon$}{ \lambda$}-\frac{$\epsilon$( \mu$- \lambda$)}{\lambda$^{2},\notleq_{$\lambda$}^{ \epsilon$}){n!$\Gam a$(\frac{$\nu$- \epsilon$}{ \lambda$}+\frac{$\epsilon$( \mu$- \lambda$)}{\lambda$^{2})(1-\frac{$\lambda$}{ \mu$})^{\underlin{$\nu$}- \epsilon$}\mathrm{R}^+\frac{$\epsilon$( \mu$- \lambda$)}{\lambda$^{2}(\frac{$\lambda$}{ \mu$})^{n(39)

P_{\mathrm{s}(n)=\displaystle\frac{$\Gam a$(n+\frac{$\nu$}{ \lambda$}){$\Gam a$(\frac{$\nu$}{ \lambda$})n!(1-\frac{$\lambda$}{ \mu$})^{\frac{$\nu$}{ \lambda$}(\frac{$\lambda$}{ \mu$})^{n}(40)

これは修正パスカル分布 (または修正負の2項分布) である.パラメータ $\epsilon$ は表式に含まれている (cf. 厳密な平衡
分布 (15) と比較せよ). $\epsilon$=0 のとき ( a\displaystyle \rightar ow\frac{ $\nu$}{ $\lambda$}, b\displaystyle \rightar ow\frac{( $\mu$- $\lambda$)}{ $\lambda$}\langle n\rangle, c\rightarrow 0 及び $\Gamma$(n+1,1-a;(\langle n\rangle+b)c)\rightarrow $\Gamma$(n+aこの表式は, P_{FP,s}(n) , パスカル分布 (または,負の2項分布) に帰着する :

191

3.3

時間に依存した解

さて,FPE の時間依存解を初期条件(41)$\psi$(y,0)= $\delta$(y-y_{0})のもとで解くことを考えよう.この FPE の初期条件はマスター方程式 (11) では次の初期条件に対応する:

P_{FP}(n,0)=\displaystyle \frac{(\langle n\rangle y_{0})^{n} {n!}\exp(-\langle n\rangle y_{0})$\delta$_{n,n_{\mathrm{O} .

(42)

極限, $\xi$\equiv\langle n\rangle y_{0}\rightarrow 0 では,これはデルタ関数に一致する $\delta$_{n,0} (i.e., n_{0}=0 ).
られており [19, 22] 次式のようになる :z

座標での FP 方程式の解は与え$\psi$(z, t|z_{0},0)=b^{a}z^{a-1}e^{-bz}\displaystyle \times\sum_{m=0}^{\infty}\frac{m!e^{- $\beta$ mt} { $\Gamma$(m+a)}L_{m}^{a-1} (bz)L_{m}^{a-1}(bz_{0}) ,(43)

ここで, L_{n}^{ $\alpha$}(x) はLaguerre 多項式であり次式で定義される:

L_{n}^{ $\alpha$}(x)=\displaystyle \frac{( $\alpha$+1)_{\mathfrak{n} {n!}{ _{1}F_{1}(-n, $\alpha$+1, x) .

(44)

ラゲール多項式の和公式 [30] を用いると次式を得る:$\psi$(z, t|z_{0},0)=(\displaystyle\frac{$\beta$}{$\delta$})(\frac{z}{z_{0}e^{-$\beta$t} )^{\frac{a-1}{2} \frac{e^{-\not\inz} {1-e^{-$\beta$t} \exp(-\frac{e^{-$\beta$t} {1-e^{-$\beta$t} \frac{$\beta$}{$\delta$}(z+ _{0}) \displaystyle\timesI_{a-1}(2\frac{$\beta$}{$\delta$}\frac{\sqrt{z _{0} {1-e^{-$\beta$t}e^{-$\tau$^{t} $\beta$), (45)

ここで I_{q}(x) は第1種の修正ベッセル関数であり次式で定義される :I_{q}(x)=\displaystyle\sum_{k=0}^{\infty}\frac{(\frac{x}{2})^{q+2k} {k!$\Gam a$(q+k+1)}

.

(46)

y 座標での FP 近似の結果,時間に依存した確率質量は次式で与えられる :P_{FP}(n,t)=\displaystyle\frac{\langlen\rangle^{n} {n!}\frac{(_{$\delta$}^{$\var ho$})e^{\frac{$\beta$(a-1)}{2}t {1-e^{-$\beta$t} \int_{0}^{\infty}dye^{-\langlen\rangley} ^{n}(\frac{y+c}{y_{0}+c})^{2(a-1)_{e^{-} $\beta$7(y+c)\perp\displaystyle \times\exp(-\frac{e^{- $\beta$ t} {1-e^{- $\beta$ t} \frac{ $\beta$}{ $\delta$}(y+y_{0}+2c) \times I_{a-1}(2\frac{ $\beta$}{ $\delta$}\frac{\sqrt{(y+c)(y_{0}+c)} {1-e^{- $\beta$ t} e^{-g_{t} )

.

(47)

上式で n=0 と置けば,絶滅確率は次式で与えれられる:

P_{FP}(0,t)=\displaystyle\frac{(_{$\delta$}^{4})e^{\frac{(a-1)$\beta$}{2}t {1-e^{-$\beta$t} \int_{0}^{\infty}dye^{-\langlen\rangley}(\frac{y+c}{y_{0}+c})^{1}\mathrm{z}(a-1)_{$\beta$}e^{-$\tau$(y+\mathrm{c})\displaystyle \times\exp(-\frac{e^{- $\beta$ t} {1-e^{- $\beta$ t} \frac{ $\beta$}{ $\delta$}(y+y_{0}+2c) \times I_{a-1}(2\frac{ $\beta$}{ $\delta$}\frac{\sqrt{(y+c)(y_{0}+c)} {1-e^{- $\beta$ t} e^{-E_{t} 2)

.

(48)

直接数値計算によりこれらは評価でき , その時間微分 f_{FP}( $\tau$)=-dP_{FP}(0, $\tau$)/d $\tau$ , や待ち時間分布 f_{FP}( $\tau$) も評価できる.

$\epsilon$=0(c=0) の場合,ボアソン変換の積分の評価が解析的にに実行できる.すなわち, $\epsilon$=0 のときの時間に依存した確率質量は次式となる:

P_{FP}(n,t)=\displaystyle\frac{\langlen\rangle^{n}{n!}\frac{(\frac{$\beta$}{$\delta$})e^{\frac{$\beta$(a-1)}{2}t {1-e^{-$\beta$t}\int_{0}^{\infty}dye^{-\langlen\rangley} ^{n}(\frac{y}{y_{0})^{\frac{1}{2}(a-1)}e^{-\ovalbox{\t\smal REJ CT}_{(y)}\displaystyle\times\exp(-\frac{e^{-$\beta$t} {1-e^{-$\beta$t} \frac{$\beta$}{$\delta$}(y+y_{0}) \timesI_{a-1}(2\frac{$\beta$}{$\delta$}\frac{\sqrt{y _{0} {1-e^{-$\beta$t} e^{-\frac{$\beta$}{2}t )

.

(49)

(49) の積分の結果,次式を得る:

P_{FP}(n,t)=\displaystle\frac{$\Gam a$(\frac{$\nu$}{ \lambda$}+n){ !$\Gam a$(\frac{$\nu$}{ \lambda$})(\frac{1-\frac{$\lambda$}{ \mu$}{1-\frac{$\lambda$}{ \mu$}e^{-$\beta$ }1^{\frac{\mathrm{v} $\lambda$}(\frac{$\lambda$}{ \mu$}\mathrm{K}(t) ノ e’甲 (-\displaystyle\frac{e^{-$\beta$t}{1-e^{-$\beta$t}\frac{$\beta$}{$\delta$}y_{0}) {}_1}F_{1}(n+\displayst le\frac{$\nu$}{$\lambda$}),\frac{$\nu$}{$\lambda$};\frac{C^{2}(t)}{B(t)} (50)
n

192

ここで,

B(t)=\displaystle\frac{\langen\rangle(1-\frac{$\lambda$}{ \mu$}e^{-$\beta$ }){\frac{$\lambda$}{ \mu$}(1-e^{$\beta$ }) ,c(t)=\frac{(n\rangle(1-\frac{$\lambda$}{ \mu$})\sqrt{y_0}e^{-$\beta$ }{(\frac{$\lambda$}{ \mu$})(1-e^{$\beta$ }) K(t)=\frac{1-e^{$\beta$ }{1-\frac{$\lambda$}{ \mu$}e-$\beta$ }

であり,

1 F1

(a, b;z)=はKmmner 関数である.

\displaystyle\frac{$\Gam a$(b)}{$\Gam a$(b-a)$\Gam a$(a)}\int_{0}^{1}\exp(zu) u^{a-1}(1 -u)^{b-a-1}

FP 近似下での絶滅確率は初期条件 (41) の下で (49) 式で n

=

du ,

0 と置き {}_{1}F\mathrm{i}(a, a;z) =\exp(z) に注意すれば,次式を得る:

PFP (0, t)=(1-\displaystyle\frac{$\lambda$}{$\mu$})^{$\nu$}\mathrm{X}(1-\frac{$\lambda$}{$\mu$}e^{-$\beta$t})^{-$\tau$}\exp$\nu$\{-\frac{e^{-$\beta$t}{1-\frac{$\lambda$}{$\mu$}e-$\beta$t}(1-\frac{$\lambda$}{$\mu$})\langlen\rangley0\}y_{0} = 0 のとき,この表式は1Fi (a, b;0) = 1に注意すれば $\epsilon$

=

.

(51)

0 で初期条件 p(n, 0) = $\delta$_{n,0} のときの表式に帰着す

る.また,待ち時間分布は (51) の P_{FP}(0, t) を微分して ( f_{FP}( $\tau$) = -dP_{FP}(0_{T})/d_{T} ) , 次のように求まる :f_{FP}($\tau$)=(1-\displaystyle\frac{$\lambda$}{$\mu$})^{\frac{$\nu$}{$\lambda$}(1-\frac{$\lambda$}{$\mu$}e^{-$\beta\tau$})^{-\frac{$\nu$}{$\lambda$}-2 $\beta$e^{-$\beta\tau$}\exp\{-\frac{e^{-$\beta\tau$}{1-\frac{$\lambda$}{$\mu$}e-$\beta\tau$}(1-\frac{$\lambda$}{$\mu$})\{n\rangley_{0}\\displaystyle\times[\frac{$\nu$}{$\mu$}(1-\frac{$\lambda$}{$\mu$}e^{-$\beta\tau$})-(1-\frac{$\lambda$}{$\mu$})\langlen\rangley_{0}](52)$\epsilon$=0 のとき,初期条件 P(n, 0) =$\delta$_{n,n0} の下での絶滅確率は (20) 式より次式を得る

P(0,t|n0,0)=:(1-\displaystyle\frac{$\lambda$}{$\mu$})^{$\nu$}\mathrm{x}(1-\frac{$\lambda$}{$\mu$}e^{-($\mu$-$\lambda$)t} ^{-\mathrm{x}-n0}$\nu$(1-e^{-($\mu$-$\lambda$)t} ^{n0}(53)}n $\xi$)} でパラメータ $\xi$= \langle n\rangle y_{0} であるとして,我々は (51) をこの分布
ここで, n_{0} がボアソン分布恥 (n_{0}) \simeq。^{n0!^{\exp(の和をとることにより得ることができる : P_{FP}(0,t)=\displaystyle \sum_{n0=0}^{\infty} p0 (no)P(0, t|n_{0},0) .
=4

数値例と議論本解析で最も重要な結果は,時間依存解の近似解が求められ,それを使えば待ち時間分布が求められることである.以下ではこの結果を中心に述べる.

4.1

待ち時間分布

確率質量の時間依存解は (47) で与えられている.絶滅確率 (48) を使って,待ち時間分布は

f_{FP}( $\tau$)=-\displaystyle \frac{d}{d $\tau$}P_{FP}(0, $\tau$)(54)で評価することができる.

図1はパラメータの組み ( $\lambda$, $\mu$, $\nu$, $\epsilon$) に対する待ち時間分布 f( $\tau$) を示す: (i) (1, 2, 3,0) (実線),(ii) ( 1, 2, 12, 0)(破線) 及び (iii) ( 1, 2, 20, 0) (点線).パラメータ a の値が増加すると,単純な指数減衰が図に示すように2重指数減衰に変化する.(i) の場合の直接数値計算の待ち時間分布の曲線は (52) 式の待ち時間分布に完全に一致する.
図2はパラメータの組 ( $\lambda$, $\mu$, $\nu$, $\epsilon$)=(0.018,0.08,1.0,0.0) では待ち時間の分布f( $\tau$) は2つの異なるモードからできていることを示している: (i) 実線は y_{0}=0 のときの2つのモードが存在する場合; (五) 点線は $\tau$ が大きな場合の遅いモードを示す (f_{L}\exp(- $\beta \tau$)) そして (iii) 破線は $\tau$ の小さな場合の早いモードを示す f_{S}(1-\displaystyle \frac{ $\lambda$}{ $\mu$}\exp(- $\beta \tau$))^{-} 矢 -1.
図3はパラメータの組 ( $\lambda,\ \mu$, $\nu$, $\epsilon$) に対する待ち時間分布 \mathrm{f}_{FP}( $\tau$) を示す: (a) (1, 2, 20,0.7), (b) (1,2, 12, 0.7) 及び(c) (1, 2, 3, 0.7) 但し, $\epsilon$\neq 0 である.各図の実線 [(\mathrm{a}),(\mathrm{b}),(\mathrm{c})](48) 式の直接数値計算の結果を示す.各場合において,図のスロープはほとんど同じ値 $\beta$ をとる.大きな $\tau$ でのこのスロープは $\epsilon$ の値によらない普遍則であることを意味する.これが正しければ,後に示すようにパラメータの組 ( $\lambda$, $\mu$, $\nu$, $\epsilon$) を確率質量と待ち時間分布から推定可能
になる.

193

1

m①

10

\mathrm{i}\mathrm{o}1\mathrm{O}^{-}

10^{-}

ttme,

図1: パラメータの組み ( $\lambda$, $\mu$, $\nu$, $\epsilon$) に対する待ち時間分布 f( $\tau$) を示す: (i) ( 1, 2, 3, 0) (実践),(ii) ( 1, 2, 12, 0) (破線 ) 及び (iii) ( 1, 2, 20, 0) (点線).パラメータ a の値が増加すると,単純指数減衰が図に示すように2重指数減衰に変化する.(i) の場合の直接数値計算の待ち時間分布の曲線は (52) 式の待ち時間分布に完全に一致する.

$\xi leftarow
trme \mathrm{t}

図2: パラメータの組 ( $\lambda$, $\mu$, $\nu$, $\epsilon$)(0.018, 0.08, 1.0, 0.0) では待ち時間の分布 f( $\tau$) は2つの異なるモードか
0 のときの2つのモードが存在する場合; (ii) 点線は時間 $\tau$
らできていることを示している: (i) 実線は y_{0}が大きな場合の遅いモードを示す (f_{L}\exp(- $\beta \tau$))_{\rangle}. そして (iii) 破線は時間 $\tau$ の小さな場合の早いモードを示す==

f_{S}(1-\displaystyle \frac{ $\lambda$}{ $\mu$}\exp(- $\beta \tau$))^{-\mathrm{x}-1} $\nu$.

1\displaytemhr{g}fcD\Sapoxir=\\vteildwed{gfruancrl\im{Dtshmrb{Bet}\fa\cmathhrro}{{Bd}ms\geqdot{}subeq\mathr{}infy\geq\ap rox

F\mathrm{Q}\in

11$\iota_{frc9}X2dsplyemh\a{TiotrQ} 1\underlin{\ $Phi$}\chek{gqaproxundlie\mth{}ar$\vpi 1

.

1

tmot

1(

tmet

\mathrm{i}\mathrm{i}\mathrm{m}$\varepsilon$\mathrm{I}

図3: パラメータの組 ( $\lambda$, $\mu$, $\nu$, $\epsilon$) に対する待ち時間分布 f_{FP}( $\tau$) を示す: (a) (1,2, 20, 0.7), (b) (1, 2, 12, 0.7) 及び(c) (1, 2, 3, 0.7) 但し, $\epsilon$\neq 0 である.各図の実線 [(\mathrm{a}),(\mathrm{b}),(\mathrm{c})](48) 式の直接数値計算の結果を示す.各場合において,図のスロープはほとんど同じ値 $\beta$ をとる.数値計算はMaple©
時の値である.

で実行した.待ち時間分布の破線 $\epsilon$=0 の

194

4.2

パラメータ推定実験 (または数値実験) データ解析の観点からすると,パラメータの推定は重要である.パラメータの推定には様々な方法がある: (i) 最尤推定; (ii) モーメント法; (iii) 最小二乗法.(15) 式の超幾何分布に関しては,超幾何関数の微分が推定式に含まれるので最尤推定法は困難となる.ここではモーメント法による簡便法を紹介する.

(15) 式の確率質量は3つのパラメータを含んでいる:$\thea$_{1}=\displayst le\frac{$\nu$}{ \lambda$} , $\thea$_{2}=\displayst le\frac{$\lambda$}{ \mu$} 及び $\thea$_{3}=\displayst le\frac{$\epsilon$}{ \mu$}
これらのパラメータは平均(55)

M , 分散 V 及び3次のキュムラント K を使って推定できる

何故なら, $\theta$_{1}, $\theta$_{2} 及び $\theta$_{3} は次式で関係付けられているからである:

M=\displaystyle\frac{$\theta$_{1}$\theta$_{2}-$\theta$_{3}{1-$\theta$_{2}, V=\displayst le\frac{$\theta$_{2}($\theta$_{1}-$\theta$_{3}){(1-$\theta$_{2})^{2} 及び K=\displaystyle\frac{$\theta$_{2}(1+$\theta$_{2})($\theta$_{1}-$\theta$_{3}) {(1-$\theta$_{2})^{3} .
従って,3つのパラメータ $\theta$_{1}, $\theta$_{2} 及び(56)$\theta$_{3} は M, V 及び K を用いて次式で推定できる:$\theta$_{1}=\displaystyle \frac{2V^{2} {(K-V)}-M, $\theta$_{2}=\displaystyle \frac{K-V}{(K+V)} 及び $\theta$_{3}=\displaystyle \frac{2V^{2} {(K+V)}-M .

(57)

もし,パラメータ $\beta$ が$\beta$\equiv $\mu$- $\lambda$ ,T(58)の大きなところでのスケール指数として得られると,4つのパラメータ ( $\lambda$, $\mu$, $\nu$, $\epsilon$) が次式で推定される:$\lambda$=\displayst le\frac{$\beta\theta$_{2}{1-$\theta$_{2}, $\mu$=\displaystyle\frac{$\beta$}{1-$\theta$_{2}, $\nu$=\displaystle\frac{$\beta\heta$_{1}$\thea$_{2} 1-$\thea$_{2} 及び $\epsilon$=\displayst le\frac{$\beta\thea$_{3} 1-$\thea$_{2} ,4. 3

3\mathrm{D}(59)

スクロール波の数値実験との関連性

Davidsen ら [1], Reid ら [2] 及び Clayton [4, 5] は

3\mathrm{D} スクロール波の負張力状態及び乱流状態での数値計算

を3次元の Barkley モデル [3], 3次元の複素 Ginzburg‐Landau 方程式 [16] 及び3次元の 3\mathrm{V}‐SIMモデル [17] で行っている.負張力状態では,彼らは生成死滅率は次式で表現できることを報告している:$\lambda$_{n}= $\lambda$ n+ $\nu$

及び$\mu$_{n}= $\mu$ n+ $\epsilon$(60)

ここで ( $\lambda$, $\mu$, $\nu$) は正値で $\epsilon$ は負値となっている.Clayton の数値実験では ( $\lambda$>0, $\mu$>0, $\nu$=0 及び $\epsilon$=0 ) として解析が可能である.しかし,Reid et al [2] の典型的なパラメータは 負張力状態では $\lambda$=0.018, $\mu$=0.08, $\nu$=1
及び $\epsilon$=-2.8 となっている.母関数に基づいた解析的なアプローチに関する限り,(i) 平均と分散は (20) 式の母関数から計算できる; (ii) 定常確率質量 (15) は数値実験の分布を再現する.しかし,(iii) 時間依存解 P(n, t) , 絶滅確
率 P(0, t) 及び待ち時間分布 f( $\tau$) を求めることができない,何故なら 母関数 (20) の z=0 の数学的な特異点が存在するからである.もし $\epsilon$<0 の値が十分に小さく、パラメータ a>0 であれば z の値も正値に留まるから,われわれのこの近似も妥当であろう.しかし,我々の方法は大きな負値の $\epsilon$ に対する時間依存解を得ることには失敗する.定常状態に関する限り,(i) 測定された確率質量(Reid et al. [2]) は (15) で捉えることができ, (\d ot{\mathrm{u} )(55) -(59) のパラメータ推定は $\epsilon$ の絶対値が大きくても ( $\lambda$=0.018, $\mu$=0.08, $\nu$=1 及び $\epsilon$=-2.8 ) うまく働く.

5

まとめと諸言(2) 式の線形生成死滅率を有するマスター方程式 (1) を研究した.このマスター方程式に対応する母関数には
厳密解 (20) があるが, z=0 に存在する数学的な特異点のために,確率質量 P(n, t) , 絶滅確率 P(0, t) 及び待ち時間
分布 f( $\tau$)=-\displaystyle \frac{dP(0. $\tau$)}{d $\tau$} が求められていなかった.ボアソン変換を用いることにより,確率質量の近似解, P_{FP}(n, t)

195

を(47) 式に,近似絶滅確率 P_{FP}(0,t) を(48) 式にそして近似待ち時間分布か P(T)=-\displaystyle \frac{P_{FP}(0, $\tau$)}{d $\tau$} を求めた.確率質量の表式 PFP (n, t は n_{0} がポアソン則に従っているときのものである.我々が知る限り,初めて得られたものである.これらの結果は4つのパラメータ ( $\lambda,\ \mu$, $\nu$, $\epsilon$) が関係した,生成死滅率 $\lambda$_{n}= $\lambda$ n+ $\nu$ 及び $\mu$_{n}= $\mu$ n+ $\epsilon$ を有する (cf. [l9,20, 21, 23, 2]) 確率過程の解析に有用である.この確率質量と待ち時間の情報から,4つのパラメータの推定も可能である.

謝辞

本研究は科学研究費補助金 (JSPS, 挑戦的萌芽研究,課題番号 15\mathrm{K}11993 ) の援助を受けている.

参考文献[1] J. Davidsen, M. Zhan and R. Kapral, Phys. Rev. Lett., 101, 208302 (2008).
[2] J. C. Reid, H. Chate and J. Davidsen, European Phys. Lett., 94, 68003 (2011).

[3] M. Dowle, R. M. Mantel and D. Barkley, Int. J. Bif. Chaos 7, 2529 (1997).
[4] R. H. Clayton, Chaos, 18, 043127 (2008).
[5] R. H. Clayton, Physica \mathrm{D} , 238, 951 (2009).

[6] S. Alonso, R. Kahler, A. S. Mikhailov and S. Saugues, Phys. Rev. \mathrm{E} , 70, 056201 (2004).

[7] K. E. Daniels, C. Beck and E. Bodenschatz, Physica \mathrm{D} , 193, 208 (2004).
[8] Y.‐N.Young and H. Recke, Phys. Rev. Lett., 90, 134502 (2003).
[9] S. Kai and W. Zimmermann, Progress of Theoretical Physics, Supp. 99 , 458 (1989).
[10] M. M. Genkin, A. Sokolov, O. D. Laventovich and I. S. Aranson, Phys. Rev. X, 7, 011029 (2017).
[11] A. Vilhois and D. Proment and G. Krstulovic, Phys. Rev. \mathrm{E} , 93, 061103 (R) (2016).
[i2] C. Beta, A. S. Mikhailov, H. H. Rotermund and G. Ertl, Europhys. Lett., 75, 868 (2006).
[13] L. Gil, J. Lega, and J. Meunier,Phys. Rev. \mathrm{A} , 41, 1138 (1990).
[14] Y. Uchiyama and H. Konno, Phys. Lett. \mathrm{A} , 378, 1350 (2014).
[15] Y. Uchiyama and H. Konno, Phys. Rev. \mathrm{E} , 91, 022127 (2015).
[16] I. S. Aranson and L. Kramer, The world of the complex Ginzburg‐Landau equation, Rev. Mod. Phys. 74,
99 (2002).

[17] F. H. Fenton, E. M. Cherry, H. S. Hastings and S. J. Evans, Chaos, 12, 852 (2002).

[18] H. Konno, J. Phys. 8\mathrm{o}\mathrm{c} . Jpn, 86, 064001 (2017).
[19] C. W. Gardiner, Handbook of Stochastic Methods, (Springer, Berlin, 1983).
[20] N. S. Goel and N. R. Dyn, Stochastic Models in Biology, (Academic Press, New York, 1971).

[21] W. Feller, An Introduction to Probability Theory and Its Applications, (Wiley, New York, 1968).
[22] J. C. Cox, J. E. Ingersol and S. A. Ross, Econometrica 53, 385‐409 (1985).

[23] S. Pigolotti, et al. PNAS, 102, 15747‐15751 (2005).

196[24] S. Kerlin and J. McGregor, Trans. Amer. Math. Soc. 86, 366‐401 (1957).
[25] S. Kerlin and J. McGregor, J. Math. and Mech. 7, 643\mapsto 662 (1958).

[26] M. E. H. Ismail, Classical and Quantum Orthogonal Polynomials in One Variable, (Cambridge Umiv. Press,
Cambridge, 2009)[27] R. Sasaki, J. Math. Phys. 50, 103509 (2009).
[28] W. Schoutens, J. Math. Anal. Apph. 253, 515‐531 (2001).
[29] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, (North‐Holland, Amsterdam, 1992).

[30] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products (sixth edition). (Academic Press,
San Diego, 2000).