Finite Differences

Numerical differentiation via finite differences for C++20.

Overview¶

The derivative is defined as a limit:

We cannot take h = 0 on a computer, but we can choose small h. The art of numerical differentiation lies in choosing h wisely---small enough that the approximation is good, but not so small that floating-point errors dominate.

Quick Start¶

#include <finite-diff/finite_diff.hpp>
#include <cmath>
#include <iostream>

int main() {
    using namespace finite_diff;

    // Derivative of sin(x) at x = 1
    auto f = [](double x) { return std::sin(x); };

    double x = 1.0;
    double df = central_difference(f, x);  // Uses optimal h

    std::cout << "f'(1) = " << df << "\n";         // ~0.5403
    std::cout << "cos(1) = " << std::cos(1.0) << "\n";  // 0.5403...
}

The Two Errors¶

Every finite difference approximation has two sources of error:

1. Truncation error: From the Taylor series approximation. Decreases as \(h \to 0\).
2. Round-off error: From floating-point arithmetic. Increases as \(h \to 0\).

where \(\varepsilon \approx 2.2 \times 10^{-16}\) (machine epsilon for double).

The optimal \(h\) minimizes total error:

API Reference¶

First Derivatives¶

// Forward difference: (f(x+h) - f(x)) / h
// O(h) accuracy
double df = forward_difference(f, x);

// Backward difference: (f(x) - f(x-h)) / h
// O(h) accuracy
double df = backward_difference(f, x);

// Central difference: (f(x+h) - f(x-h)) / (2h)
// O(h^2) accuracy - the workhorse
double df = central_difference(f, x);

// Five-point stencil: higher accuracy formula
// O(h^4) accuracy
double df = five_point_stencil(f, x);

Second Derivatives¶

// Second derivative: (f(x+h) - 2f(x) + f(x-h)) / h^2
double d2f = second_derivative(f, x);

// Five-point second derivative: O(h^4)
double d2f = second_derivative_five_point(f, x);

Multivariate Functions¶

// Gradient: grad(f)(x) for f: R^n -> R
auto grad = gradient(f, x);

// Directional derivative: D_v f(x)
double Dvf = directional_derivative(f, x, v);

// Hessian matrix: H[i][j] = d^2f/dx_i dx_j
auto H = hessian(f, x);

// Jacobian matrix for f: R^n -> R^m
auto J = jacobian(f, x);

Deriving the Formulas¶

Forward Difference¶

Taylor expand \(f(x+h)\):

Solve for \(f'(x)\):

Central Difference¶

Taylor expand \(f(x+h)\) and \(f(x-h)\):

Subtract:

The even-order error terms cancel, doubling our accuracy!

Five-Point Stencil¶

By taking linear combinations of more points, we can cancel more error terms:

This achieves \(O(h^4)\) by canceling the \(h^2\), \(h^3\), and \(h^4\) error terms.

Optimal Step Size Derivation¶

For central difference with round-off:

Minimize by taking derivative with respect to \(h\) and setting to zero:

For double precision: \(h \approx 6 \times 10^{-6}\).

Comparison with Automatic Differentiation¶

Use finite differences when: - Function is a black box (no source code) - Quick approximation needed - Validating other methods

Use automatic differentiation when: - Accuracy is critical - Computing many derivatives - Source code is available