Lesson 1: Numerical Differentiation

Lesson 1-1: Fundamentals of Numerical Differentiation: Approximation Using Finite Difference Methods

Before numerically evaluating the derivative of a function, let us review the mathematical definition of differentiation. The derivative of a function \(f(x)\) is defined as: \[f'(x)=\lim_{h\rightarrow 0} \frac{f(x+h)-f(x)}{h}\] Based on this definition, the derivative of a function \(f'(x)\) can be approximated using a sufficiently small \(h\) as follows: \[f'(x) \approx \frac{f(x+h) - f(x)}{h}. \tag{1} \] This method of approximating differentiation using a finite (but sufficiently small) width \(h\) is called the finite difference method. The error of the finite difference approximation can be evaluated using a Taylor expansion. For instance, the Taylor expansion of the function \(f(x)\) is as follows: \[f(x+h) = f(x) + f'(x)h +O(h^2). \tag{2}\] Here, \(O(h^2)\) represents terms of order \(h^2\) or higher and indicates the truncation error of the Taylor series at a finite order. Rearranging equation (\(2\)) gives the following: \[f'(x) = \frac{f(x+h)-f(x)}{h}+O(h). \] Therefore, the error of the finite difference approximation formula, equation (\(1\)), is proportional to \(h\), and this is referred to as an \(O(h)\) error. From this, it is clear that the derivative \(f'(x)\) can be accurately approximated by using a sufficiently small \(h\).

By using the Taylor expansion, it is also possible to derive more accurate approximation formulas. Here, we consider the second-order Taylor expansions of the function \(f(x+h)\) and \(f(x-h)\) as follows: \[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}\] By taking the difference between equations (\(3\)) and (\(4\)), we can derive the following approximation formula: \[f'(x) = \frac{f(x+h)-f(x-h)}{2h} +O(h^2). \tag{5}\] The error of this finite difference approximation formula, equation (\(5\)), is proportional to \(h^2\), and in the limit of small \(h\), it converges to the exact derivative more rapidly than the difference formula in equation (\(2\)). From this, it is evident that for sufficiently small \(h\), equation (\(5\)) is a more accurate approximation formula than equation (\(2\)). Difference formulas like equation (\(2\)) are referred to as forward difference formulas, while those like equation (\(5\)) are referred to as central difference formulas.

Numerical Differentiation Code Using C Language

In this section, we will learn how to write numerical differentiation in C language using the forward difference formula and the central difference formula. As a specific example, we consider the differentiation of the exponential function, \( f(x) = e^x \), at \( x = 0 \). The following code calculates the derivative of the function, \( f'(0) \), while varying the step size, \( h \), and evaluates the error from the exact derivative value. The error in numerical differentiation is defined as the deviation from the exact derivative value: \( 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;
}

In the code above, the derivative value at \(x=0\) is being calculated. For this purpose, in line 7 of the code, the variable x is initialized by assigning x=0. Additionally, in line 8, the exact derivative value of the function is computed using exact = exp(x).

In line 13 of the code, a file is opened to output the results. The printf statement is used to format and display the results.

In the code mentioned above, a for loop begins, where the numerical differentiation is repeatedly evaluated while varying the step size \(h\). The initial value of h is 1.0, and it is halved in each iteration.

The difference calculations are performed using the forward_diff function and the central_diff function. The former employs the forward difference formula, while the latter uses the central difference formula.