2.1. 差分

2.1.1. 差分公式

関数 f(x) の解析的な表式が与えられていれば微分 f'(x) は必ず計算できますが、 f(x) が数値的にしか得られない場合（複雑な積分を含むなど）には、微分を数値的に評価する必要がでてきます。 もし、f(x) が精度よく、滑らかに見える程度に十分なデータがあれば、スプライン補間をして微分を計算する手段が有効です（別ページ）。 固有値問題 や 偏微分方程式 の数値解法で示すように x を離散化する必要がある場合には、以下で説明する 差分公式 を使います。

例として、1階微分を考えます。 座標 x をある区間で等間隔に離散化し、間隔を h=x_{i+1}-x_{i} と表します。 関数 f(x) の各離散点 x_i 上における値 f(x_i) \equiv f_i を持っているとします。 すると、1階微分 f'(x) は

(1)f'(x_i) = \frac{f_{i+1} - f_{i}}{h} + \mathcal{O}(h^1)

と表せます。 これを 前進差分 (forward finite difference)と呼びます。 最後の \mathcal{O}(h^1) は、この近似式の誤差が h の1乗に比例することを表しています。 この公式は、点 x_{i} における微分を求めるために、x_{i} とひとつ先の点 x_{i+1} の情報を使っていますが、 x_{i} とひとつ前の点 x_{i-1} を使って微分を計算することもできます：

(2)f'(x_i) = \frac{f_{i} - f_{i-1}}{h} + \mathcal{O}(h^1)

これを 後退差分 (backward finite difference)と呼びます。 近似の精度は前進差分と同じ \mathcal{O}(h^1) です。

前進差分や後退差分の精度 \mathcal{O}(h^1) は、実用的には良いものではありません。どれだけ h を小さくしても必ず誤差が残ります。 また、対称性の悪い形になっている点も気になります（例えば、点 x_{i} における後退差分の式 (2) は 点 x_{i-1} における前進差分と同じ式です）。 そこで、x_{i+1} と x_{i-1} の2点を使うと、より対称性の良い公式が得られます：

(3)f'(x_i) = \frac{f_{i+1} - f_{i-1}}{2h} + \mathcal{O}(h^2)

これを 中心差分 (central finite difference)と呼びます。 この場合、誤差は \mathcal{O}(h^2) と次数が一つ下がります（つまり、精度が上がります）。 h に比例する項が消えることは、f_{i\pm1} \equiv f(x_i \pm h) を h についてテイラー展開することで確かめられます。 また、同じ計算から、誤差が -f'''(x_i) h^2 /6 となることも示せます。

さらに、x_{i\pm1} に加えて x_{i\pm2} の計4点を使ことで、より高精度の公式を導くこともできます：

f'(x_i) = \frac{-f_{i+2} + 8f_{i+1} - 8f_{i-1} + f_{i-2}}{12h} + \mathcal{O}(h^4)

同様にして、2階以上の微分も f_{i}, f_{i\pm1} などを適当な係数を掛けて足し合わせることで表現できます。 係数の一覧は、例えば Wikipedia - Finite difference coefficient にまとまっています。

2.1.2. 行列表示

座標 x を離散化して f(x) をベクトル表示 \bm{f} \equiv (f(x_0), f(x_1), \cdots, f(x_{N-1}))^\mathrm{T} すると（\mathrm{T} は転置の意味）、微分操作は行列・ベクトル積として表されます。 具体的に1階微分を考えると、 \bm{f}' \equiv (f'(x_0), f'(x_1), \cdots, f'(x_{N-1}))^\mathrm{T} は

\bm{f}' = D_x \bm{f}

と表されます。 ここで、D_x は N \times N 行列で、 中心差分公式 (3) を採用すると、

D_x = \frac{1}{2h} \begin{pmatrix} 0 & 1 & 0 & \cdots & \cdots & 0 & -1 \\ -1 & 0 & 1 & 0 & & & 0\\ 0 & -1 & 0 & 1 & \ddots & & \vdots \\ \vdots & 0 & -1 & 0 & \ddots & \ddots & \vdots \\ \vdots & & \ddots & \ddots & \ddots & \ddots & 0 \\ 0 & & & \ddots & \ddots & \ddots & 1 \\ 1 & 0 & \cdots & \cdots & 0 & -1 & 0 \end{pmatrix}

で与えられます。 左下の1と右上の-1は周期境界条件 f(x_{N}) = f(x_{0}) を仮定しているためです。

2.1.3. pythonでの実装方法

pythonで差分を計算する場合には、ベクトル表示を使って行列・ベクトル積を計算するとforループを使わずに実装できます。 特に、微分演算子は疎行列なので、疎行列モジュール scipy.sparse を使うと効率が良いです。 疎行列クラスとしていくつかありますが、scipy.sparse.csr_matrix (Compressed Sparse Row format)が適当であると思います。 右からベクトルを書ける演算が中心なので、行列を行ベクトルとして持っていた方が若干効率がよさそうという理由です（大差はないと思います）。