数値微分における刻み幅 \(h\) の影響解析

プログラムの動作

このプログラム diff_h.py は、数値微分における刻み幅 \(h\) が結果の精度に与える影響を解析するPythonスクリプトです。

目的

数値微分の基本的なアルゴリズムである2点前方差分法を用い、定義された関数 \(f(x)\) の導関数を計算します。複数の異なる刻み幅 \(h\) を適用し、各 \(h\) における数値微分結果と厳密な導関数（解析解）を比較することで、刻み幅 \(h\) の選択が数値微分の精度にどのように影響するかを視覚化することを目的としています。

主な機能

関数の定義: 微分対象となる関数 \(f(x) = x^3\) を定義します。
厳密な導関数の定義: \(f(x)\) の厳密な導関数 \(f'(x) = 3x^2\) を定義します。これは数値微分結果との比較対象となります。
数値微分の実行: 2点前方差分法 \(f'(x) \approx \frac{f(x+h) - f(x)}{h}\) を用いて、複数の異なる刻み幅 \(h\) で数値微分を実行します。プログラム内部では、h = [1.0, 0.1, 0.01, 0.001, 1.0e-6] の5種類の刻み幅が使用されます。
結果の出力: 計算された \(x\) の値、関数値 \(f(x)\)、厳密な導関数 \(f'(x)\)、および各刻み幅 \(h\) での数値微分結果を、標準出力とCSVファイルの両方に出力します。

解決する課題

数値微分の精度は、使用する刻み幅 \(h\) に大きく依存します。 \(h\) が大きすぎると離散化誤差が大きくなり、小さすぎると浮動小数点演算の丸め誤差の影響が顕著になります。このプログラムは、これらの誤差の挙動を具体的に示し、適切な \(h\) の選択の重要性を理解するのに役立ちます。

原理

このプログラムは、関数 \(f(x) = x^3\) の導関数を、その厳密な解析解 \(f'(x) = 3x^2\) と、数値計算によって得られる近似解を比較することで、刻み幅 \(h\) の影響を評価します。

対象関数と厳密な導関数

微分対象の関数は以下の通りです。

\[ f(x) = x^3 \]

この関数の厳密な導関数（解析解）は以下の通りです。

\[ f'(x) = \frac{d}{dx}(x^3) = 3x^2 \]

数値微分アルゴリズム

プログラムでは、2点前方差分法 (two-point forward difference method) と呼ばれる最も基本的な数値微分アルゴリズムを使用しています。この方法は、ある点 \(x\) における関数の導関数を、その点とわずかに前方（\(x+h\)）の関数値を用いて近似します。

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

この近似は、テイラー展開に基づくと、一次の誤差 \(O(h)\) を持つことが知られています。すなわち、刻み幅 \(h\) が小さいほど、近似の精度は向上する傾向がありますが、特定のしきい値を超えて \(h\) を小さくしすぎると、浮動小数点演算の丸め誤差により、かえって精度が低下する可能性があります。

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

このプログラムは、数値計算を効率的に行うために numpy ライブラリを使用しています。csv ライブラリはPythonの標準ライブラリであるため、別途インストールする必要はありません。

`numpy`

numpy は、Pythonで数値計算を効率的に行うための基本的なパッケージです。多次元配列オブジェクトと、それらを操作するための様々な関数を提供します。

インストール方法

numpy がシステムにインストールされていない場合、以下の pip コマンドを使用してインストールできます。

pip install numpy

必要な入力ファイル

diff_h.py は、外部の入力ファイルを必要としません。微分対象の関数 \(f(x)\) とその厳密な導関数 \(f'(x)\) は、プログラムのソースコード内で直接定義されています。

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

プログラムを実行すると、以下のCSVファイルが生成されます。

ファイル名: diff_h.csv

ファイルの内容

diff_h.csv は、計算された \(x\) の各点における関数値、厳密な導関数、および複数の刻み幅 \(h\) での数値微分結果をカンマ区切り形式で保存します。

ヘッダー行: x, f(x), df/dx, Df/Dx(h=1.0), Df/Dx(h=0.1), Df/Dx(h=0.01), Df/Dx(h=0.001), Df/Dx(h=1.0e-6)

x: 計算点の \(x\) 座標。
f(x): \(x\) における関数 \(f(x)=x^3\) の値。
df/dx: \(x\) における厳密な導関数 \(f'(x)=3x^2\) の値。
Df/Dx(h=...): 指定された刻み幅 \(h\) を用いた2点前方差分法による数値微分結果。括弧内には使用された刻み幅 \(h\) の値が示されます。

データ行: ヘッダー行に続く各行は、x の各点に対応する上記の値を含みます。 \(x\) は x0 = 0.0 から nx = 21 点まで、xstep = 0.1 刻みで計算されます（\(x = 0.0, 0.1, ..., 2.0\)）。

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

diff_h.py は、コマンドライン引数を必要としません。プログラムは、内部で定義されたパラメータ（刻み幅 \(h\) のリスト、 \(x\) の範囲など）に基づいて処理を実行します。

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

python diff_h.py

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

以下のコマンドを実行することで、プログラムを実行し、結果を標準出力で確認できます。同時に、diff_h.csv というファイルがプログラム実行ディレクトリに生成されます。

python diff_h.py

実行結果（標準出力の例）

プログラムを実行すると、以下のような形式で計算結果が標準出力に表示されます。これは、diff_h.csv に書き込まれる内容とほぼ同じですが、整形された表形式で表示されます。

Numerical differentiation using two-point forward difference method
Write to [diff_h.csv]

  x  :   f(x)    df/dx       h=1.0       h=0.1      h=0.01      h=0.001     h=1.0e-6   
0.000:  0.000000  0.000000   1.00000000  0.01000000  0.00010000  0.00000100  0.00000000  
0.100:  0.001000  0.030000   1.33100000  0.06310000  0.03310000  0.03030100  0.03000000  
0.200:  0.008000  0.120000   1.72800000  0.18760000  0.12760000  0.12060100  0.12000000  
...

（以降、同様の形式で \(x=2.0\) まで21行のデータが続きます。）

生成される出力ファイル `diff_h.csv` の内容

上記の実行により、diff_h.csv という名前のCSVファイルがプログラム実行ディレクトリに生成されます。このファイルは、スプレッドシートソフトウェアなどで開いて詳細な解析やグラフ化を行うことができます。

diff_h.csv の内容は以下のようになります。

x,f(x),df/dx,Df/Dx(h=1.0),Df/Dx(h=0.1),Df/Dx(h=0.01),Df/Dx(h=0.001),Df/Dx(h=1.0e-6)
0.0,0.0,0.0,1.0,0.01,0.0001,1.0000000000000001e-06,0.0
0.1,0.001,0.03,1.3310000000000002,0.0631,0.0331,0.030301,0.030000000000000006
0.2,0.008,0.12,1.728,0.1876,0.1276,0.120601,0.12000000000000001
...