Richardson 外挿法による数値微分

プログラムの動作

このPythonプログラム diff_richardson.py は、リチャードソン外挿法（Richardson Extrapolation Method）を用いて関数の数値微分を行うことを目的としています。具体的には、定義済みの関数 \(f(x) = e^x\) の \(x=1.0\) における1階微分を計算します。

主な機能は以下の通りです。

  • 関数と解析解の定義: 微分対象の関数 func(x) とその解析的な1階微分 dfdx(x) をそれぞれ exp(x) として定義しています。

  • 中心差分近似: 逐次的にステップサイズ \(h\) を半分にしながら、中心差分近似 \(D(h) = \frac{f(x_0+h) - f(x_0-h)}{2h}\) を計算します。

  • リチャードソン外挿: 計算された複数の中心差分近似値から、より高次の精度を持つ微分値を外挿によって算出します。

  • 収束判定: 計算された微分値の差分が事前に定義された許容誤差 eps (デフォルト \(1.0 \times 10^{-6}\))を下回った場合に、計算を停止します。

  • 結果の出力: 計算の途中経過（各ステップでの近似値）と、解析解、そして最終的な収束結果を標準出力に表示します。

このプログラムは、数値微分において単純な差分近似では \(h\) を小さくしすぎると浮動小数点数の丸め誤差が支配的になるという課題を、リチャードソン外挿法を用いることで高精度な結果を効率的に得ることを目指しています。

原理

このプログラムは、リチャードソン外挿法を用いて数値微分を行います。基本的なアイデアは、異なるステップサイズ \(h\) で計算された粗い近似値を組み合わせて、より高い精度の近似値を生成することです。

  1. 中心差分近似: 関数の1階微分 \(f'(x_0)\) の中心差分近似は、以下の式で与えられます。

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

    この近似の誤差は、ステップサイズ \(h\) の偶数乗の項で展開できます。

    \[f'(x_0) = D(h) + K_1 h^2 + K_2 h^4 + K_3 h^6 + \dots\]

    ここで \(K_1, K_2, \dots\) は定数です。

  2. リチャードソン外挿: リチャードソン外挿は、上記の誤差展開を利用して、より高次の精度を持つ近似値を生成します。ステップサイズ \(h\) を半分にした \(h/2\) を用いた中心差分近似 \(D(h/2)\) も同様に表現できます。

    \[f'(x_0) = D(h/2) + K_1 (h/2)^2 + K_2 (h/2)^4 + \dots\]
    \[f'(x_0) = D(h/2) + K_1 \frac{h^2}{4} + K_2 \frac{h^4}{16} + \dots\]

    これらの2つの式から \(h^2\) の項を消去することで、より高次の近似値を得ることができます。

    \[D^{(1)}(h) = \frac{4 D(h/2) - D(h)}{4 - 1}\]

    この \(D^{(1)}(h)\) は \(h^4\) のオーダーの誤差を持つため、単純な中心差分近似よりも高精度です。このプロセスを繰り返すことで、さらに高次の近似値を生成できます。

    一般化された漸化式は以下の通りです。

    • 初期値として、ステップサイズ \(h_m = h_0 / 2^m\) を用いた中心差分近似を \(D_{0,m}\) とします。

      \[D_{0,m} = \frac{f(x_0+h_m) - f(x_0-h_m)}{2h_m}\]

    • リチャードソン外挿による更新式は、以下のように定義されます。

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

      ここで、\(k\) は外挿の次数（\(k=1\) で \(h^2\) の項を消去、\(k=2\) で \(h^4\) の項を消去、といった具合）、\(m\) はステップサイズ \(h_m\) に対応するインデックスです。コードでは、D[m1][k] が \(D_{m1, k}\) に対応しており、m1 が外挿の次数、\(k\) がベースとなる中心差分近似のインデックスを表します。

プログラムは、nhmax 回までこの外挿プロセスを繰り返し、隣接する近似値 D[m][0] と D[m-1][0] の差 diff が eps より小さくなった時点で収束したとみなし、計算を終了します。

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

このプログラムは、以下の非標準ライブラリをインポートしていますが、そのうち実際に数値計算に使用されているのは math.exp のみであり、numpy はインポートされているものの、明示的にその機能が利用されている箇所はありません。

  • numpy:

    • 目的: Pythonで数値計算を効率的に行うためのライブラリです。多次元配列オブジェクトや、それらを操作するための関数群を提供します。

    • インストール方法: コマンドプロンプトやターミナルで以下のコマンドを実行してインストールします。

      pip install numpy

必要な入力ファイル

このプログラムは、実行時に外部の入力ファイルを必要としません。 微分対象の関数 func(x)、その解析解 dfdx(x)、および計算パラメータ（初期ステップサイズ h0、最大反復回数 nhmax、収束判定の許容誤差 eps、微分点 x0）はすべてプログラムのソースコード内にハードコードされています。

生成される出力ファイル

このプログラムは、実行時にファイル形式の出力を生成しません。 すべての計算結果と途中経過は、標準出力（コンソール）に直接表示されます。これには、解析的な関数値と導関数値、各ステップで計算されるリチャードソン外挿の近似値 D[i][j]、および収束判定に使われる差分値などが含まれます。

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

diff_richardson.py は、特別なコマンドライン引数を必要としません。 プログラムを単体で実行するだけで、内部に定義されたパラメータに基づいて数値微分が実行され、結果が標準出力に表示されます。

基本的な実行コマンドは以下の通りです。

python diff_richardson.py

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

以下のコマンドを実行することで、プログラムを起動します。

python diff_richardson.py

実行すると、以下のような出力がコンソールに表示されます（出力は一部を抜粋）。

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: 2.859296856455551
D01: 2.753517452654637
D10: 2.717657651354336
D02: 2.7262441951912976
D11: 2.717665393246857
D20: 2.718282729904268
  diff = D[2][0] - D[1][0] = 2.718282729904268 - 2.717657651354336 = 0.0006250785499319088
D03: 2.719602073981882
D12: 2.717666293961805
D21: 2.718281804246875
D30: 2.718281829051877
  diff = D[3][0] - D[2][0] = 2.718281829051877 - 2.718282729904268 = -9.008527878297405e-07

この出力は、リチャードソン外挿法の進行状況を示しています。

  • D00, D01, D02, ...: これらはステップサイズ \(h_0, h_0/2, h_0/4, \dots\) を使用した最も基本的な中心差分近似の結果です。

  • D10, D11, D12, ...: これらは1回目のリチャードソン外挿を適用して、誤差の \(h^2\) の項を消去した近似値です。

  • D20, D21, ...: これらは2回目のリチャードソン外挿を適用して、誤差の \(h^4\) の項を消去した近似値です。

  • D[m][0]: 各ステップで得られる、最も高精度な近似値の候補です。

  • diff: 現在の D[m][0] と前のステップの D[m-1][0] の差分です。この値が eps (この例では \(1.0 \times 10^{-6}\))を下回ると、計算が終了します。

上記の出力例では、D[3][0] と D[2][0] の差分が -9.008...e-07 となり、eps (\(1.0 \times 10^{-6}\))よりも小さくなったため、プログラムはここで収束して停止します。最終的に得られた D[3][0] の値 \(2.718281829051877\) は、解析解である \(e^1 \approx 2.718281828459045\) に非常に近い値を示していることがわかります。