数値計算 ~差分法~

解析解(数学解)が求められない場合に、数値解析で数値解を計算するが、計算機で扱うためには離散化が必要である。ここでは、微分方程式を解く時に必要な離散化を考える。微分方程式を解くための離散化には、時間の離散化と、空間の離散化がある。離散化で最もよく使われるのが差分法であり、前進差分、後退差分、中央差分がある。空間の離散化には、最も直感的な中央差分が使われることが多いように思う。もっとも、扱う系によって異なるし、格子状に離散化しない場合、例えば粒子法による流体計算は最も狭い意味での中央差分ではない(概念的には中央での離散化を行っているのは同じ)。時間の離散化には、陽解法では前進差分、陰解法では後退差分、クランク・ニコルソン法というものでは中央差分を使うものとなっていて、一長一短がある。陽解法はシンプルに現在の状態を元に漸化式を解くイメージで次の時刻の状態を計算するもので、実装が簡単だが、必ずしも数値的に安定で収束するわけではないため、時間ステップを大きく取れない。陰解法は現在の状態がどのようなものであれば一つ前の状態と適合するかを解くもので、連立方程式を解く必要があり、実装が複雑だが、時間ステップを大きく取れるため、シミュレーション時間が長い時には有利。クランク・ニコルソン法は、これも実装が複雑だが、離散化誤差は最も小さい。陰解法、クランク・ニコルソン法はともに数値的に安定な(情報が伝播する速さが、実際の物理現象の伝播する速さよりも早くある必要があるという、CLF条件による制約がない)ので、時間ステップは離散化誤差が許容できる範囲で大きく取れることになる。クランク・ニコルソン法では、ステップの中央に注目して微分方程式を解くが、離散化した時にそこには格子点が存在しないため、両脇の格子点の平均を近似として用いる。なお、CLF条件は必要条件であるが、それ以外に安定性の解析方法にフォン・ノイマンの安定性解析がある。これは、差分法でどの精度の差分式を使うかによって異なる打ち切り誤差が、計算の過程で収束するか発散するかを解析するものである。これは解きたい微分方程式によって異なる結果を出すが、時間ステップ幅と空間ステップ幅の関係がどうであれば打ち切り誤差が収束するかがわかる。ここでは書かないが、連立方程式を解くときの行列の固有値が関係してくる。関数の微分値を差分法で表すには、関数のテイラー展開を元にして計算を行えばよい。何次の項まで使うかによって精度が異なる。これによって、差分法には、(前進, 後退, 中央)差分×(N回微分)×(M次精度)というような組み合わせが言われうる。まず、前進と後退のテイラー展開から、1階微分の差分式を導出してみる。$f(x)$のn階微分を$f^{(n)}(x)$と書くことにする。前進差分: $$ f(x+h)=\sum^{\infty}_{n=0}\frac{f^{(n)}(x)}{n!}h^n=f(x)+f^{(1)}(x)\cdot h+\frac{f^{(2)}(x)}{2!}h^2+\frac{f^{(3)}(x)}{3!}h^3+… \tag{1}$$ ここから、1階1次精度の前進差分を求めることができる。 $$f^{(1)}(x)=\frac{f(x+h)-f(x)}{h}+\mathcal{o}(h) \tag{2}$$ 後退差分: $$ f(x-h)=\sum^{\infty}_{n=0}\frac{f^{(n)}(x)}{n!}(-h)^n=f(x)-f^{(1)}(x)h+\frac{f^{(2)}(x)}{2!}h^2-\frac{f^{(3)}(x)}{3!}h^3+… \tag{3}$$ ここから、1階1次精度の後退差分を求めることができる。 $$f^{(1)}(x)=\frac{f(x-h)-f(x)}{-h}+\mathcal{o}(h) \tag{4}$$ 中央差分: 前進差分と後退差分のテイラー展開の差をとると、右辺の奇数番目の項が消える。そこから1階2次精度の中央差分が計算できる $$f(x+h)-f(x-h)=2(f^{(1)}(x) \cdot h+\frac{f^{(3)}(x)}{3!}h^3+\frac{f^{(5)}(x)}{5!}h^5+…) \tag{5}$$ $$f^{(1)}(x)=\frac{f(x+h)-f(x-h)}{2h}+\mathcal{o}(h^2) \tag{6}$$

さて、ここからは、1階前進差分の高次精度の差分式をみたいと思う(後退差分については$h$を$-h$に読み換えれば良い)。前進差分で$2h$先をみたときのテイラー展開は $$ f(x+2h)=\sum^{\infty}_{n=0}\frac{f^{(n)}(x)}{n!}(2h)^n=f(x)+f^{(1)}(x)\cdot 2h+\frac{f^{(2)}(x)}{2!}(2h)^2+\frac{f^{(3)}(x)}{3!}(2h)^3+… \tag{7}$$ であるが、(1)式の4倍から(7)式を引いて$f^{(1)}(x)$について整理すると$h^2$の項が消えて、2次精度の差分式ができる。 $$f^{(1)}(x)=\frac{-f(x+2h)+4f(x+h)-3f(x)}{2h}+\mathcal{o}(h^2) \tag{8}$$ 同様に、前進差分で$3h$先をみたときのテイラー展開は $$ f(x+3h)=\sum^{\infty}_{n=0}\frac{f^{(n)}(x)}{n!}(3h)^n=f(x)+f^{(1)}(x)\cdot 3h+\frac{f^{(2)}(x)}{2!}(3h)^2+\frac{f^{(3)}(x)}{3!}(3h)^3+… \tag{9}$$ であるから、($h$以外の)$h^2$と$h^3$の項を消すために、(1)式-(7)式*(-1/2)+(9)式*(1/9)を計算し、$f^{(1)}(x)$について整理すると $$f^{(1)}(x)=\frac{2f(x+3h)-9f(x+2h)+18f(x+h)-11f(x)}{6h}+\mathcal{o}(h^3) \tag{10}$$ なお、$-1/2, 1/9$という数字をどうやって出すかというと、繰り返しになるが$h^2$と$h^3$の項の係数がゼロとなるように、係数についての連立方程式を解いて求める。具体的には次の連立方程式を解いた。 $$\left\{\begin{array}{l}\frac{f^{(2)}(x)}{2!} \cdot(1+4B+9C)=0\\\frac{f^{(3)}(x)}{3!} \cdot(1+8B+27C)=0\end{array}\right. \tag{11}$$ これを解くと、$B=-1/2, C=1/9$となる。ここまでやると、N階M次精度微分の前進差分式を(同じように後退差分も)一般的に決定できることがわかる。つまり、$h, 2h, 3h, …Mh$まで先をみたテイラー展開を行い、$h^N$以外のh冪乗の項を消去できるように係数についての連立方程式を立てて解く。このとき、$M \geq N$でないと連立方程式が解けない。これは直感的にも、N階微分の計算にN+1個以上の点が必要であるということと一致する。

さて、2階微分の中央差分式も考えたい。まずは$h$まで先をみた前進と後退のテイラー展開を行った和を考える(和にしたのは$h^2$の項を消さないためである)。これは、$h^2$の項まで捉えているので、2次精度である。 $$f^{(2)}(x)=\frac{f(x+h)-2f(x)+f(x-h)}{h^2} +\mathcal{o}(h^2) \tag{12}$$ 続いて、1階微分の4次精度の中央差分式を考える。$2h$まで先をみた前進と後退のテイラー展開を行った差を考えると(5)式と同様に $$f(x+2h)-f(x-2h)=2(f^{(1)}(x) \cdot 2h+\frac{f^{(3)}(x)}{3!}(2h)^3+\frac{f^{(5)}(x)}{5!}(2h)^5+…) \tag{13}$$ (5)式と(13)式を足し引きして、右辺$h^3$の項の係数がゼロになるようにしたい。(右辺では求めたい$f^{(1)}(x)$に$h$がかかっていることに注意。4次精度なので$h^5$の項は消さないし、消そうにも式の本数が足りない、消すなら$3h$まで先をみたテイラー展開の式との連立が必要である。そのとき6次精度となる。)結局のところ、(13)式では$h^3$の項の係数が(5)式の8倍なので、(5)式*8-(13)式を計算して$f^{(1)}(x)$について整理すると $$f^{(1)}(x)=\frac{f(x+2h)-8f(x+h)+8f(x-h)-f(x-2h)}{12h}+\mathcal{o}(h^4) \tag{14}$$ ここでも、前進差分の時と同様に、N階の2M-1または2M次精度の中央差分式が計算できるとわかる。手順としては、前後に$h, 2h, 3h, … Mh$までみたテイラー展開を行い、Nが奇数なら前進と後退のテイラー展開の差を、Nが偶数なら和を取るところから始める。そして、求めたい$h^N$の項以外の係数を消去できるように、係数についての連立方程式を立てて解けばよい。

さて、前進差分、後退差分、中央差分の違いは何であるかというと、文字通りに読めば、どこを差分の中心とするかということである。問題によって使い分けるのが良く、必要に応じてこれらの差分式を混合することで、差分の中心を任意の場所に置くことができる。