植物数理モデリング

最初に一次元空間の場合を考える。このとき拡散方程式は    

     (Eq 1.1)

で与えられる。D は拡散定数。数値計算を行なうために、Eq 1.1の左辺と右辺を、それぞれ時間の１階差分と空間の２階差分で置き換える。ここで前進差分（右辺に時刻 t での値）を用いると

     (Eq 1.2)

が得られる。ただし uit は時刻 t  位置（細胞）i の u の値、Δt は時間刻み幅、Δx は空間刻み幅。この式を uit+Δtについて解くと

     (Eq 1.3)

が得られる。この式により、時刻 t  の情報から、時間が Δt 進んだ時の値（uit+Δt）を直ちに得ることができる。この方法は陽解法と呼ばれ、直観的にわかりやすく容易に計算できる利点がある反面、解の安定性が低いという欠点がある。今回の場合、解が安定であるためには時間刻み幅に関しての制約条件

     (Eq 1.4)

が課せられることになる。この条件は、拡散は濃度の高い方から低い方へと流れる現象であり、その逆には移動しないことを反映している。この条件を満たしていない場合、数値計算は破綻することになるので、陽解法を用いる場合は常にこの制約条件（時間刻み幅の上限）に留意する必要がある。

　　一般に、時間刻み幅が小さくなるにつれて、計算精度は高くなるのに対し、計算量すなわち計算時間は増大することになる。従ってEq 1.4の制約条件は、計算時間を強く拘束する要因になり、実際に計算を行なう上で不都合を引き起こす可能性がある。この制約条件は、次項の陰解法を用いることにより回避することができる。

Eq 1.2 において、前進差分ではなく後退差分（右辺に時刻 t+Δt の値）を用いると

     (Eq 2.1)

となり、これを変形することにより

     (Eq 2.2)

が得られる。この式に従って計算される uit+Δt の値を用いるのが陰解法である。しかし、この値は直ちに得られる訳ではなく、一旦Eq 2.2を解く必要があり、陽解法に比べて手続きが煩雑になる。また、一般に陰解法においては解を求めること自体が困難な場合もあり、常にこの手法が使えるとは限らない。この様な欠点があるものの、陰解法は解の安定性が陽解法に比べて高く、今回においてもEq 1.4 の制約条件がなくなるという非常に大きな利点を有している。つまり、どのような時間刻み幅（Δt）をとろうが、この方法を用いている限り解が発散することはない。

　　しかしながら、無制限にΔt を大きくしても良いという訳ではない。上述のように、時間刻み幅を大きくしていくと計算量すなわち計算時間が減少していくが、それと同時に計算精度が悪くなることに注意しなければならない。従って、Δt を大きくする場合は、自分の目的に照らし合わせて、それが適切かどうかを確認する必要がある。より計算精度の高い方法を次項で紹介する予定です。

Eq 2.2 を解く際には、行列表記をした方が便利である。

     (Eq 2.3a)

     (Eq 2.3b)

     (Eq 2.3c; 周期境界条件)

     (Eq 2.3d; ノイマン境界条件)

ただし α = D Δt/Δx2、n は細胞数（空間分割数）。これは線形の連立方程式であるので、ut+Δt に関して解くことができるが、実用的には行列 A をLU分解することにより容易に求めることができる。（境界条件およびLU分解はいづれ説明するかも）

前項では陰解法を紹介したが、より一般化させて前進差分と後退差分を組み合わせてみる。

     (Eq 3.1)

右辺第１項は前進差分、第２項は後退差分に対応し、パラメータθ（0 ≤ θ ≤ 1）は後退差分の比率を表している。この式において、θ = 0 のときはEq 1.2 すなわち陽解法に対応しており、それ以外の場合（θ > 0）は広い意味での陰解法に対応する。特に θ = 1 のときはEq 2.1 と一致し、ここでは特に断らない限りこの場合（θ = 1）を単に陰解法と呼ぶことにする。これらに加え、θ = 1/2 の場合は特にCrank-Nicolson法と呼ばれており、以下のように書き換えられる（α = D Δt/Δx2）。

     (Eq 3.2)

この方法は、陰解法と同様に計算手続きは煩雑であるものの、解は常に安定である。しかしこれだけではこの手法を使用する意味はあまりないが、Crank-Nicolson法は陰解法に比べて計算精度が高い（つまり計算誤差が小さい）という大きな利点があり、それによってしばしば用いられている。次項では、その点を含めてこれらの手法間の比較をしようと思う。Eq 3.2 の計算は、陰解法と同様にして行なうことができる。

Table 1

陽解法、陰解法、およびCrank-Nicolson法の基本的特徴をTable 1にまとめておく。陽解法は数値計算の手続きが非常に簡潔であるのに対し、陰解法およびCrank-Nicolson法はそれに比べて遥かに煩雑であり、それだけ計算量（計算時間）も余計にかかる。その反面、陽解法は時間刻み幅 Δt に強い制約（Eq 1.4）があるのに対し、他の２手法ではそれが回避できる利点を有する。これは状況にもよるが、計算手続きが煩雑であるという欠点を上回る利点となりうる。

Table 2

ここで、陽解法での制約条件Eq 1.4について考えてみる。今この制約を満たすほぼ上限の Δt 値を用いているとして、そこから空間分割数を a (>1) 倍にした時、それが計算量（つまり計算時間）にどのような影響をもたらすかを考えてみる（Table 2参照）。分割数が増えることにより、計算量は１次元空間の場合は a 倍に、２次元空間の場合は a2 倍に増加することになる。一方で、空間分割数が a 倍になることは、空間刻み幅 Δx がΔx/a になることを意味するので、制約条件Eq 1.4を満たすためには時間刻み幅 Δt を Δt/a2 以下にする必要がある。これは計算量を a2 倍に増加させる効果を及ぼすことになる。これらの効果を合わせると、例えば空間分割数を a = 10 倍にした場合、計算量は１次元空間の場合は a3 = 1,000 倍に、２次元空間の場合には a4 = 10,000 倍に膨れ上がることになる。この制約条件が、数値計算にいかに重大な影響を与えうるかが想像できると思う。特に２次元空間を扱う際にはこの効果は顕著になるので、可能ならば陰解法を用いることを検討した方が良い。

数値計算の安定性と計算精度（計算誤差）は似て非なるものです。計算誤差（Table 1）において、O(Δtn) は Δtn のオーダーで計算値に誤差が生じることを意味しており、nが大きいほど計算誤差は小さくなる、つまり計算精度が高くなることを意味している。従って、陽解法と陰解法は同程度の計算精度で、それらに比べてCrank-Nicolson法は計算精度が高いことになる。例えば、時間刻み幅Δt を a (= 10) 分の１（つまり Δt → Δt/a）に変更した場合、これにより計算量は a (= 10) 倍に増えることになる。一方で計算誤差は、陽解法と陰解法では 1/a (= 1/10) になるのに対し、Crank-Nicolson法では 1/a2 (= 1/100) になり、遥かに計算精度が向上することが分かると思う。

　　プログラムの作成を考えたとき、陽解法から陰解法への変更に比べて、陰解法からCrank-Nicolson法への変更ははるかに簡単に行なうことができる。したがって、もし陰解法のプログラムが既にできているのならば、折角なのでCrank-Nicolson法への変更も考えても良いと思う。

どの手法や条件を選択するかは、目的に依存する問題であり、一概に決めることはできない。多くの場合の判断基準としては、その目的に必要な計算精度を満たしつつ、なるべく計算時間が短くなる様な条件を選ぶということである。その過程で、Table 1およびこの項の説明を参考にしていただければ幸いです。

２次元空間の場合、拡散方程式は

     (Eq 5.1)

と書ける。この式を前進差分で離散化すると

     (Eq 5.2)

となり、さらに変形して

     (Eq 5.3)

が得られる。Eq 5.1–5.3 は１次元空間の場合のEq 1.1–1.3 にそれぞれ相当し、これが２次元空間の場合の陽解法に対応する。１次元空間の場合と同様に、この方法はとても簡便であるものの、 Eq 1.4 と同様の時間刻み幅に関する制約条件が課せられることになる。前項で述べたように、この制約は２次元空間の場合により厳しく影響するので、実際に計算を行なう上でしばしば問題になることがある。その場合には、次項以降に述べる陰解的な手法を用いることをお勧めする。

Eq 5.1で後退差分を用いると

     (Eq 6.1)

となり、変形して

     (Eq 6.2)

が得られる。Eq 6.1–6.2 は１次元空間の場合のEq 2.1–2.2 にそれぞれ対応する。この連立方程式はとても繁雑なので、 ui,jt+Δt を解くのは通常は容易ではない。しかし周期境界条件を仮定できれば、フーリエ変換を利用して比較的簡単に解くことができる。

Eq 6.2 は、畳み込み（convolution）を用いることにより形式的に

     (Eq 6.3)

と書き換えることができる。ただし、２次元空間を M × N に分割した場合、（巡回）畳み込みは

     (Eq 6.4)

で与えられ（i = 0, ⋯ , M–1;  j = 0, ⋯ , N–1）、関数は周期的であるとする。 ut は  

     (Eq 6.5)

を満たす離散関数、k は拡散の効果を表す離散関数で、行列表記では

     (Eq 6.6)

に対応する（α = D Δt/Δx2）。畳み込みの重要な性質として、畳み込みのフーリエ変換はそれぞれのフーリエ変換の積に等しい（畳み込み定理）が知られている。

     (Eq 6.7)

そこで、Eq 6.3を（離散）フーリエ変換することにより

     (Eq 6.8)

となる。両辺をで割り、（離散）逆フーリエ変換することにより数値的に解くことができる。

     (Eq 6.9)

実装する場合は、高速フーリエ変換で計算することになるので、そのコードが使用できるのならば、この方法はとても有用だと思う。

（※）この方法は「発生の数理（三浦岳 著、京都大学学術出版会）」でMathematicaのコードとともに詳しく紹介されています。

ここで紹介した手法は、画像処理の分野では頻繁に用いられており、デコンボリューション（deconvolution: 逆畳み込み）と呼ばれている。その場合、ut は観察画像、ut+Δt は復元画像、k は画像劣化の作用に対応することになる。

次項では、別の陰的解法を紹介する。

ここでは境界条件に関係なく適用できるADI法を紹介する。この方法では、２次元の拡散を x 方向の拡散と y 方向の拡散に分解して交互に計算を行なう。これにより各ステップの計算は１次元の拡散のみを考慮すれば良いことになる。

     (Eq 7.1)

各ステップの拡散は１次元的であるので、それぞれを１次元の陰解法を用いることにより、Eq 5.1は

     (Eq 7.2a)

     (Eq 7.2b)

と離散化できる。これら２式の両辺を足すと

     (Eq 7.3)

となり、陽解法（Eq 5.2）や陰解法（Eq 6.1）との違いが認識できる。Eq 7.2を変形すると

     (Eq 7.4a)

     (Eq 7.4b)

が得られ、Eq 7.4a により時刻 t の値から (t + Δt/2)の値が、さらにEq 7.4b により (t + Δt/2)の値から (t + Δt)の値が計算できる。これらの式は、１次元空間の陰解法（Eq 2.2）での手法を繰り返し適用することにより計算できる。このADI法の特徴を簡単にまとめると

(１) 陽解法で見られる制約条件がないので、時間刻み幅によらず数値計算は常に安定である。

(２) 境界条件に関係なく適用できる。

(３) １次元空間の陰解法の計算手法をそのまま適用できるので、比較的容易にプログラムを構築できる。

(４) Crank-Nicolson法を容易に組み込むことができ、それにより計算精度を上げることができる。

次項では、これまで紹介した手法間の関係を、プログラム作成の観点から簡単にまとめようと思う。

（つづく）