diff_richardson.py の技術ドキュメント

プログラムの動作

diff_richardson.py は、リチャードソン補外法を用いて関数の数値微分を行うPythonスクリプトです。具体的には、指数関数 \(f(x) = \exp(x)\) の \(x=1.0\) における導関数を数値的に計算し、その結果を解析解と比較します。

主な機能:

関数定義: 微分対象の関数 \(f(x) = \exp(x)\) とその解析的な導関数 \(f'(x) = \exp(x)\) を内部で定義しています。
初期近似: 初期ステップサイズ \(h_0\) を用いて、中央差分公式に基づいた数値微分の初期近似値を計算します。
ステップサイズ調整: 計算の繰り返しごとにステップサイズ \(h\) を半分にしていきます。
リチャードソン補外: 異なるステップサイズで得られた数値微分値から、リチャードソン補外法によってより高精度な近似値（D行列）を計算し、行列を構築します。
収束判定: D行列の対角要素（最も精度の高い近似値）が収束条件（隣接する近似値の差が許容誤差 eps 以下）を満たした場合に計算を終了します。
結果出力: 解析解と計算された数値微分値、および計算過程でのD行列の各要素と収束判定の差分を標準出力に表示します。

解決する課題: 単純な数値微分（例：中央差分）は、ステップサイズ \(h\) を小さくするほど真の値に近づきますが、ある程度以上小さくすると浮動小数点数の丸め誤差が支配的になり、かえって精度が悪化するという問題があります。リチャードソン補外法は、この問題に対し、複数のステップサイズで得られた近似値から誤差項をキャンセルすることで、丸め誤差の影響を抑えつつ、より高精度な近似値を効率的に導き出すことを可能にします。

原理

このプログラムは、数値微分をより高精度に行うために「リチャードソン補外法」を利用しています。

1. 中央差分公式

まず、数値微分の最も基本的な手法の一つである中央差分公式を使用します。関数 \(f(x)\) の \(x_0\) における導関数 \(f'(x_0)\) は、ステップサイズ \(h\) を用いて次のように近似されます。

\[ D(h) = \frac{f(x_0+h) - f(x_0-h)}{2h} \]

この近似の誤差は \(O(h^2)\) となります。つまり、真の値との差は \(h^2\) に比例します。

2. リチャードソン補外法

中央差分公式の誤差は、\(h\) の偶数次の項で展開できることが知られています。

\[ f'(x_0) = D(h) + C_1 h^2 + C_2 h^4 + C_3 h^6 + \dots \]

リチャードソン補外法は、この誤差展開を利用して、異なるステップサイズで計算された複数の近似値から、より高精度な近似値を外挿する手法です。

ステップサイズ \(h\) と \(h/2\) で中央差分公式を適用して得られた近似値をそれぞれ \(D(h)\) と \(D(h/2)\) とします。これらの値から、\(h^2\) の誤差項をキャンセルした新しい、より高精度な近似値 \(D'\) を次のように計算できます。

\[ D' = \frac{4D(h/2) - D(h)}{3} \]

この \(D'\) は \(O(h^4)\) の誤差を持ちます。

このプロセスを一般化したものがD行列を構築するリチャードソン補外法の再帰式です。プログラムでは、このD行列の各要素 \(D_{m,k}\) を計算するために以下の漸化式が用いられています。

\[ D_{m,k} = \frac{4^m D_{m-1,k+1} - D_{m-1,k}}{4^m - 1} \]

ここで、

\(D_{0,k}\) は、ステップサイズ \(h = h_0 / 2^k\) で計算された中央差分による初期近似値です。
\(m\) は補外の「次数」を表し、\(m\) が大きいほど高次の誤差項をキャンセルしています。
\(k\) はD行列の行（または列）インデックスに対応し、異なるステップサイズでの初期近似値を示します。

プログラムのコードでは、D[m1][k] が \(D_{m1, k}\) に対応し、m1 が \(m\) に対応します。D[0][m] は \(D_{0,m}\) に対応します。この漸化式を繰り返し適用することで、行列の対角方向 (\(D_{m,0}\)) に高精度の近似値が得られます。

必要な非標準ライブラリとインストール方法

このプログラムは、標準ライブラリ math の exp 関数のみを使用しており、外部の非標準ライブラリは特に必要ありません。 import numpy as np が記述されていますが、コード内で numpy の機能は一切使用されていません。

必要な入力ファイル

diff_richardson.py は、外部からの入力ファイルを必要としません。微分対象の関数 func(x)、その解析導関数 dfdx(x)、および計算パラメータ (h0, nhmax, eps, x0) はすべてスクリプトのソースコード内にハードコードされています。

生成される出力ファイル

diff_richardson.py は、いかなるファイルも生成しません。すべての計算結果と途中経過は、実行時に標準出力（コンソール）に表示されます。

コマンドラインでの使用例 (Usage)

diff_richardson.py は、コマンドライン引数を必要としません。スクリプトを直接実行するだけで動作します。

python diff_richardson.py

コマンドラインでの具体的な使用例

以下のコマンドを実行します。

python diff_richardson.py

このコマンドを実行すると、リチャードソン補外法による数値微分の計算過程と結果が標準出力に表示されます。初期パラメータは以下の通りです。

h0 = 1.0 (初期ステップサイズ)
nhmax = 15 (最大計算回数またはD行列の最大サイズ)
eps = 1.0e-6 (収束判定の許容誤差)
x0 = 1.0 (微分を行う点)

\(f(x) = \exp(x)\) および \(f'(x) = \exp(x)\) なので、x0 = 1.0 における解析値は \(f(1.0) = \exp(1.0) \approx 2.718281828\) です。

以下は実行結果の抜粋と説明です。実際の出力はさらに多くのD行列の要素を含みますが、ここでは主要な部分を示します。

Numerical differentiation using by Richardson extrapolation method

eps for convergnece: 1e-06
Analytical values:
  f(1.0)=2.718281828459045
  df/dx(1.0)=2.718281828459045

D00: 3.1945280494548487
D01: 2.822467367375211
D10: 2.7231464403487396
D02: 2.738127533088012
D11: 2.718300262423377
D20: 2.7182818318451163
  diff = D[2][0] - D[1][0] = 2.7182818318451163 - 2.7231464403487396 = -0.004864608503623276
D03: 2.722510344463372
D12: 2.718282305018659
D21: 2.718281828494884
D30: 2.7182818284590635
  diff = D[3][0] - D[2][0] = 2.7182818284590635 - 2.7182818318451163 = -3.38605282216503e-09

出力結果の説明:

最初の数行で、解析値と収束条件 eps が表示されます。f(1.0) と df/dx(1.0) は解析的に \(\exp(1.0) \approx 2.718281828459045\) です。
D00 は、初期ステップサイズ \(h_0 = 1.0\) での中央差分による近似値です。真の値から大きく外れています。
D01, D02, ... は、ステップサイズを半分にしながら計算された中央差分による近似値です。ステップサイズが小さくなるにつれて真の値に近づきます。
D10, D11, D12, ... および D20, D21, ... は、リチャードソン補外法によって計算された高精度の近似値です。行インデックスが大きくなるほど（例: D10 -> D20 -> D30）、誤差項がより多くキャンセルされ、精度が向上します。
diff は、現在の最も高精度な近似値 (D[m][0]) と、一つ前の (D[m-1][0]) との差分です。この差分の絶対値が eps (1.0e-6) より小さくなると、計算が終了します。
上記の抜粋例では、D[3][0] が計算された時点で diff = -3.38605282216503e-09 となり、これは eps (\(10^{-6}\)) よりも小さいため、計算が停止します。
最終的に得られた近似値 D[3][0] は 2.7182818284590635 であり、解析解 2.718281828459045 と非常に近い値であることがわかります。これはリチャードソン補外法によって高い精度が達成されたことを示しています。