Lesson 1: 数値微分

Lesson 1-1: 数値微分の基礎: 差分法を用いた近似

数値的に関数の微分を評価する前に、数学的な微分の定義を復習してみましょう。関数 \(f(x)\)の微分は次の式によって定義されています: \[f'(x)=\lim_{h\rightarrow 0} \frac{f(x+h)-f(x)}{h}\] この式に基づき、十分小さな\(h\)を用いることで関数の微分\(f'(x) \)を次のように近似することが出来ます: \[f'(x) \approx \frac{f(x+h) - f(x)}{h}. \tag{1} \] このような有限の大きさの(ただし十分小さい)幅\(h\)を使って微分を近似する方法は有限差分法と呼ばれます。有限差分法による近似の誤差は、Taylor展開を用いることで次のように評価することが出来ます。例えば、関数\(f(x\))のTaylor展開は次のようになります。 \[f(x+h) = f(x) + f'(x)h +O(h^2). \tag{2}\] ここで、\(O(h^2)\)は、\(h\)の２次以上の項からなり、Taylor展開を有限の次数で打ち切った場合の誤差を表しています。式(\(2\))を変形することで次のような式を得ます。 \[f'(x) = \frac{f(x+h)-f(x)}{h}+O(h). \] したがって、有限差分近似公式, 式(\(1\)), の誤差は\(h\)に比例する事が分かります。これを\(O(h)\)[読み:オーターh]の誤差と呼びます。このことから、十分小さな\(h\)を用いることで関数の微分値\(f'(x)\)を正確に近似できることが分かります。

Talyor展開を用いることで、さらに高精度な近似公式を作ることも可能です。ここでは、関数\(f(x+h)\)と\(f(x-h)\)の2次までのTaylor展開を次のように考えてみます: \[f(x+h) = f(x) + f'(x)h + \frac{1}{2}f''(x)h^2 + O(h^3), \tag{3}\] \[f(x-h) = f(x) - f'(x)h + \frac{1}{2}f''(x)h^2 + O(h^3). \tag{4}\] ここで式(\(3\))と式(\(4\))の差を取ることによって次のような近似公式を導くことが出来ます: \[f'(x) = \frac{f(x+h)-f(x-h)}{2h} +O(h^2). \tag{5}\] この有限差分近似公式,式(\(5\)),の誤差は\(h^2\)に比例しており、\(h\)が小さい極限で式(\(2\))の差分公式よりも速く厳密な微分値に収束することが分かります。このことから、十分小さな\(h\)に対して、式(\(5\))の方が式(\(2\))よりも高精度な近似公式となっていることが分かります。式(\(2\))のような差分公式は前進差分公式と呼ばれ、式(\(5\))のような差分公式は中心差分公式と呼ばれています。

C言語を使った数値微分の計算コード

この節では、前進差分公式と中心差分公式を用いた数値微分のC言語の書き方について学びます。具体的な例題として、指数関数,\(f(x)=e^x\) の\( x=0 \)における微分を考えます。次のコードでは、関数の微分値, \( f'(0) \),を差分幅, \(h\),を変えながら計算し、厳密な微分値からの誤差を評価します。数値微分の誤差を厳密な微分値からのずれとして定義しています: \( Error = \left|f'_{numerical}(0)-f'_{exact}(0) \right| \).

(lesson1_1.c)

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

int main() {
    double x = 0.0;
    double dfdx_e = exp(x);
    double h = 1.0;
    int n = 16;
    FILE *fp;

    // Open file for output
    fp = fopen("fd_error.out", "w");

    for (int i = 1; i <= n; i++) {
        double ff_0 = exp(x);
        double ff_p = exp(x + h);
        double ff_m = exp(x - h);

        double dfdx_f = (ff_p - ff_0) / h;
        double dfdx_m = (ff_p - ff_m) / (2 * h);

        fprintf(fp, "%16.6e %16.6e %16.6e %16.6e %16.6e\n", h, dfdx_f, fabs(dfdx_f - dfdx_e), dfdx_m, fabs(dfdx_m - dfdx_e));
        
        h = h / 2.0;
    }

    // Close file
    fclose(fp);

    return 0;
}

上のコードでは、\(x=0\)における微分の値を計算しています。このため、コードの7行目で変数xにはx=0を代入して初期化しています。また、8行目では、関数の厳密な微分値をexact = exp(x)によって計算しています。

コードの13行目では、結果を出力するためのファイルを開いています。printf文を使って結果をフォーマット化して表示しています。

上記のコードではforループが始まり、刻み幅\(h\)を変えながら繰り返し数値微分の評価が行われています。hの初期値は1.0で、各反復で半分にされています。

差分計算はforward_diff関数とcentral_diff関数で行われています。前者は前進差分公式を、後者は中心差分公式を使用しています。