diff_order.py 技術ドキュメント

プログラムの動作

diff_order.py は、指数関数 \(f(x) = \exp(x)\) の導関数を、様々なステップサイズ \(h\) と異なる数値微分近似法を用いて計算し、その精度を評価するPythonスクリプトです。

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

  1. 関数の定義: 微分対象の関数 \(f(x) = \exp(x)\) と、その解析的な導関数 \(f'(x) = \exp(x)\) を定義します。

  2. 数値微分法の実装: 以下の5種類の数値微分近似法を関数として実装しています。

    • 2点前方差分

    • 2点後方差分

    • 3点中心差分

    • 5点中心差分

    • 7点中心差分

  3. 誤差評価: 初期ステップサイズ h から、kh 倍ずつ減少させながら nh 回にわたってステップサイズを変化させます。各ステップサイズにおいて、定義された各数値微分法で導関数を計算し、その結果と解析解との絶対誤差を求めます。

  4. 結果の出力:

    • 計算された数値微分値と誤差をコンソールに表示します。

    • すべての計算結果を diff_order.csv というCSVファイルに保存します。

    • matplotlib を用いて、ステップサイズ h に対する各数値微分値と、各数値微分法の誤差を視覚的に表示するグラフを生成します。

このプログラムは、数値微分におけるステップサイズと近似法の選択が、結果の精度にどのように影響するかを理解し、可視化することを目的としています。

原理

このプログラムは、関数のテイラー展開に基づいた数値微分近似法を利用しています。関数 \(f(x)\) の点 \(x\) における導関数 \(f'(x)\) は、以下のテイラー展開から導き出されます。

\[ f(x+h) = f(x) + hf'(x) + \frac{h^2}{2!}f''(x) + \frac{h^3}{3!}f'''(x) + \dots \]

これを変形することで、様々な差分近似式が得られます。

  1. 2点前方差分 (2-point forward difference): \(f(x+h) = f(x) + hf'(x) + O(h^2)\) から \(f'(x)\) について解きます。

    \[ f'(x) \approx \frac{f(x+h) - f(x)}{h} \]

    この近似は \(O(h)\) のオーダーで、ステップサイズ \(h\) が減少すると線形に誤差が減少します。

  2. 2点後方差分 (2-point backward difference): \(f(x-h) = f(x) - hf'(x) + O(h^2)\) から \(f'(x)\) について解きます。

    \[ f'(x) \approx \frac{f(x) - f(x-h)}{h} \]

    この近似も \(O(h)\) のオーダーです。

  3. 3点中心差分 (3-point central difference): \(f(x+h) = f(x) + hf'(x) + \frac{h^2}{2!}f''(x) + \frac{h^3}{3!}f'''(x) + O(h^4)\)\(f(x-h) = f(x) - hf'(x) + \frac{h^2}{2!}f''(x) - \frac{h^3}{3!}f'''(x) + O(h^4)\) を用います。 これらを引くと、偶数次の項が消去され、高次の精度が得られます。

    \[ f'(x) \approx \frac{f(x+h) - f(x-h)}{2h} \]

    この近似は \(O(h^2)\) のオーダーで、2点差分法よりも高精度です。

  4. 5点中心差分 (5-point central difference): さらに多くの点を使い、テイラー展開のより高次項を考慮して導出される、より高精度な近似です。

    \[ f'(x) \approx \frac{-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)}{12h} \]

    この近似は \(O(h^4)\) のオーダーです。

  5. 7点中心差分 (7-point central difference): 5点中心差分よりもさらに多くの点を使用し、非常に高精度な近似を提供します。

    \[ f'(x) \approx \frac{f(x+3h) - 9f(x+2h) + 45f(x+h) - 45f(x-h) + 9f(x-2h) - f(x-3h)}{60h} \]

    この近似は \(O(h^6)\) のオーダーです。

プログラムでは、関数 func(x) として \(e^x\) を、解析解 dfdx(x) としても \(e^x\) を用いてこれらの数値微分法を評価し、それぞれの誤差を解析解との絶対差として計算します。

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

diff_order.py の実行には、以下の非標準ライブラリが必要です。

  • numpy: 数値計算を効率的に行うためのライブラリです。

  • matplotlib: グラフ描画のためのライブラリです。

これらのライブラリは、Pythonのパッケージインストーラ pip を使ってインストールできます。コマンドプロンプトやターミナルで以下のコマンドを実行してください。

pip install numpy matplotlib

必要な入力ファイル

diff_order.py は、実行に必要な外部入力ファイルはありません。 すべての計算パラメータ(初期ステップサイズ h、ステップサイズ減少率 kh、計算回数 nh、微分中心点 x0)はスクリプトのソースコード内部で定義されています。

もしこれらのパラメータを変更したい場合は、diff_order.py のソースコードを直接編集する必要があります。

生成される出力ファイル

diff_order.py は、以下のファイルと視覚出力を生成します。

  1. diff_order.csv:

    • このファイルは、各ステップサイズ h における数値微分結果をカンマ区切り形式で保存します。

    • ファイルはプログラム実行ディレクトリに生成されます。

    • 内容は以下の列で構成されます。

      • Ndiv: ステップサイズの計算回数(1から nh まで)。

      • h: その回のステップサイズ。

      • Df/Dx(2-point): 2点前方差分法による数値微分値。

      • Df/Dx(3-point): 3点中心差分法による数値微分値。

      • Df/Dx(5-point): 5点中心差分法による数値微分値。

      • Df/Dx(7-point): 7点中心差分法による数値微分値。

  2. Matplotlibグラフ:

    • プログラムの実行が完了すると、matplotlib によって生成されたグラフウィンドウが表示されます。

    • このウィンドウには2つのサブプロットがあります。

      • 上段: ステップサイズ h に対する各数値微分法の計算結果と、解析解 df/dx (exact) がプロットされます。

      • 下段: ステップサイズ h に対する各数値微分法の誤差(解析解との絶対差)が両対数スケールでプロットされます。これにより、誤差が h に対してどのように減少するかが視覚的に確認できます。

    • このウィンドウは、ユーザーがターミナルで Enter キーを押すまで開いたままになります。

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

diff_order.py はコマンドライン引数を受け付けません。そのため、特別なオプションなしで直接実行します。

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

python diff_order.py

プログラムは、ソースコード内で定義されているデフォルトのパラメータ(x0 = 1.0, h = 0.5, kh = 0.8, nh = 101)を使用して計算を実行します。

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

diff_order.py を実行すると、まずコンソールに出力が行われ、その後グラフウィンドウが表示されます。

  1. 実行コマンド:

    python diff_order.py
    
  2. 実行結果の説明:

    • コンソール出力: プログラムはまず、解析値と、各ステップサイズ h に対する各数値微分法の計算結果を標準出力に表示します。 出力例(一部抜粋):

      Numerical differentiation using differnet approximations
      Write to [diff_order.csv]
      
      Analytical values:
        f(1.0)=2.718281828459045
        df/dx(1.0)=2.718281828459045
      Ndiv:       h           Df/Dx(2-point)  Df/Dx(3-point)  Df/Dx(5-point)  Df/Dx(7-point)  
       1 :  5.0000e-01      3.376839353998  2.839849503463  2.718362635445  2.718282710323  
       2 :  4.0000e-01      3.136890356616  2.771986422204  2.718296316270  2.718281830493  
      ...
      101 : 1.6385e-08      2.718281804797  2.718281828461  2.718281828459  2.718281828459  
      
    • diff_order.csv ファイルの生成: プログラムが実行されたディレクトリに diff_order.csv という名前のファイルが生成されます。このファイルには、上記のコンソール出力と対応するデータがCSV形式で格納されます。 内容例(一部抜粋):

      Ndiv,h,Df/Dx(2-point),Df/Dx(3-point),Df/Dx(5-point),Df/Dx(7-point)
      1,0.5,3.376839353998135,2.839849503463375,2.718362635444983,2.718282710322521
      2,0.4,3.136890356616147,2.771986422203649,2.718296316270388,2.718281830492801
      ...
      101,1.6384597395029463e-08,2.7182818047968594,2.718281828461159,2.718281828459045,2.718281828459045
      
    • Matplotlibグラフの表示: コンソール出力に続いて、matplotlib によって生成されたグラフウィンドウが表示されます。 このグラフは、上段に各差分法の数値微分値(\(y\)軸)と解析解がステップサイズ h\(x\)軸)に対してプロットされ、下段には各差分法の誤差(\(y\)軸)がステップサイズ h\(x\)軸)に対して両対数プロットされます。 グラフの例(テキストでの説明):

      • 上段プロット: \(h\) が大きいときは各差分法の結果が大きく異なり、解析解から離れています。\(h\) が小さくなるにつれて、より高次の差分法(5点、7点)の結果が解析解に急速に近づくことが示されます。

      • 下段プロット: 誤差が両対数スケールでプロットされており、各差分法のオーダー(\(O(h^1)\), \(O(h^2)\), \(O(h^4)\), \(O(h^6)\))に対応する傾きで誤差が減少していく様子が確認できます。ただし、非常に小さい \(h\) では浮動小数点演算の丸め誤差により、誤差が再び増加する傾向が見られることがあります。

    グラフウィンドウは、コンソールで Press ENTER to exit>> と表示された後、ユーザーが Enter キーを押すまで表示され続けます。Enter を押すと、プログラムは終了します。