diff_order.py 技術ドキュメント

プログラムの動作

diff_order.py は、指数関数 $f(x) = e^x$ の導関数を、異なる次数（2点前方差分、3点、5点、7点中心差分）の数値微分スキームで計算し、その精度を比較するPythonスクリプトです。

このプログラムは以下の機能を提供します。

特定の点 $x_0 = 1.0$ における関数 $e^x$ の導関数を計算します。
差分幅 $h$ を初期値 0.5 から順次減少させながら、各数値微分スキームを適用します。
計算された数値微分値と、解析解である $e^{x_0}$ との絶対誤差を評価します。
計算結果を標準出力に表示するとともに、diff_order.csv というCSVファイルに保存します。
matplotlib を用いて、差分幅 $h$ に対する各数値微分値、および絶対誤差の変化を示すグラフを生成し、表示します。これにより、差分幅の選択やスキームの次数が精度に与える影響を視覚的に確認できます。

原理

本プログラムは、関数 $f(x) = e^x$ の導関数 $f'(x) = e^x$ を、様々な数値微分スキームで近似し、その誤差を評価します。微分を評価する中心点 $x_0 = 1.0$ は固定されています。

差分幅 $h$ は、初期値 h = 0.5 から kh = 0.8 の比率で nh = 101 回減少させながら計算が行われます。各ステップ $ih$ における差分幅 $h_{ih}$ は以下の式で計算されます。

\[h_{ih} = h \cdot kh^{ih}\]

プログラムで実装されている主な数値微分スキームは以下の通りです。

2点前方差分 (Forward Difference): $$f'(x) \approx \frac{f(x+h) - f(x)}{h}$$ これは1次精度 ($O(h)$) の近似です。
3点中心差分 (Central Difference): $$f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$$ これは2次精度 ($O(h^2)$) の近似です。
5点中心差分: $$f'(x) \approx \frac{-f(x+2h) + 8f(x+h) - 8f(x-h) + f(x-2h)}{12h}$$ これは4次精度 ($O(h^4)$) の近似です。
7点中心差分: $$f'(x) \approx \frac{f(x+3h) - 9f(x+2h) + 45f(x+h) - 45f(x-h) + 9f(x-2h) - f(x-3h)}{60h}$$ これは6次精度 ($O(h^6)$) の近似です。

各スキームで計算された数値微分値と、解析解 $f'(x_0) = e^{x_0}$ との絶対誤差 $|f'_{numerical}(x_0) - f'_{analytical}(x_0)|$ を算出し、その変化をグラフで可視化します。両対数プロットでは、誤差が $O(h^n)$ で減少する場合、その傾きが $n$ となることが期待されます。

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

本プログラムの実行には、以下の非標準ライブラリが必要です。

numpy: 数値計算を効率的に行うためのライブラリです。
matplotlib: グラフ描画のためのライブラリです。

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

pip install numpy matplotlib

必要な入力ファイル

このプログラムは、外部からの入力ファイルを必要としません。微分対象の関数 $f(x) = e^x$ およびその解析解 $f'(x) = e^x$ はプログラム内に直接定義されています。

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

プログラムは以下のファイルを生成します。

diff_order.csv:
- 各行は、異なる差分幅 $h$ における数値微分結果と、その差分幅の情報を記録します。
- カンマ区切りのCSV (Comma Separated Values) 形式です。
- ヘッダー行は以下の通りです。 Ndiv,h,Df/Dx(2-point),Df/Dx(3-point),Df/Dx(5-point),Df/Dx(7-point)
- 各列の内容は以下の通りです。
  - Ndiv: 差分幅のインデックス（1から nh までの整数）。
  - h: 計算に使用された差分幅の値。
  - Df/Dx(2-point): 2点前方差分スキームで計算された数値微分値。
  - Df/Dx(3-point): 3点中心差分スキームで計算された数値微分値。
  - Df/Dx(5-point): 5点中心差分スキームで計算された数値微分値。
  - Df/Dx(7-point): 7点中心差分スキームで計算された数値微分値。
プロット画像 (インタラクティブ表示):
- プログラム実行中に matplotlib によって生成されるグラフウィンドウが表示されます。
- このウィンドウには2つのサブプロットが含まれます。
  - 上段: 横軸を $h$、縦軸を $\Delta f / \Delta x$ とするプロット。各スキームの数値微分値と解析解が示されます。
  - 下段: 横軸を $h$ (対数スケール)、縦軸を絶対誤差 $| \Delta f / \Delta x - df/dx |$ (対数スケール) とするプロット。各スキームの誤差の変化が示され、その収束の次数を視覚的に確認できます。
- この画像はファイルとして自動保存はされず、実行中にインタラクティブに表示されるだけです。画像をファイルとして保存したい場合は、プログラムコードを修正し、plt.savefig() 関数を呼び出す必要があります。

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

本プログラムはコマンドライン引数を受け付けません。実行する際は、Pythonインタープリターを使って直接スクリプトファイルを実行します。

python diff_order.py

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

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

python diff_order.py

上記のコマンドを実行すると、以下のような情報が標準出力に表示されます（出力の一部を抜粋）。

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.297442541400 : 2.762198084530 : 2.718765270560 : 2.718290352934 
 2    :  4.0000e-01  : 3.084196885360 : 2.736990230239 : 2.718370503022 : 2.718282305596 
 3    :  3.2000e-01  : 2.951508240579 : 2.724773229671 : 2.718297471208 : 2.718281903901 
 ...
101   :  1.6961e-09  : 2.718281828460 : 2.718281828459 : 2.718281828459 : 2.718281828459 
Press ENTER to exit>> 

同時に、プログラムが実行されたディレクトリに 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.2974425414,2.76219808453,2.71876527056,2.718290352934
2,0.4,3.08419688536,2.736990230239,2.718370503022,2.718282305596
3,0.32,2.95150824057929,2.724773229671448,2.7182974712081514,2.718281903901416
...
101,1.696144575815049e-09,2.718281828459046,2.718281828459045,2.718281828459045,2.718281828459045

さらに、数値微分結果と誤差を示す2つのグラフが新しいウィンドウに表示されます。ユーザーがそのグラフウィンドウを閉じるか、コンソールで Enter キーを押すとプログラムが終了します。