Numerical Differentiation & Finite Differences

1
#include <iostream>
2
#include <cmath>
3
#include <functional>
4
#include <limits>
5
#include <iomanip>
6
7
enum class Scheme { Forward, Backward, Central };
8
9
// Choose a reasonable step size based on machine epsilon and the x-scale
10
double optimal_step(double x, Scheme s) {
11
    const double eps = std::numeric_limits<double>::epsilon();
12
    // Scale with |x| to keep relative spacing; ensure at least 1.0
13
    const double scale = std::max(1.0, std::abs(x));
14
    if (s == Scheme::Central) {
15
        return std::cbrt(eps) * scale; // ~ eps^{1/3}
16
    } else {
17
        return std::sqrt(eps) * scale; // ~ sqrt(eps)
18
    }
19
}
20
21
// Generic first-derivative finite difference
22
// f: function to differentiate; x: point; h: step; s: scheme
23
// If h <= 0, an automatic step is chosen.
24
template <class F>
25
double derivative(F f, double x, double h, Scheme s) {
26
    double hh = (h > 0) ? h : optimal_step(x, s);
27
28
    switch (s) {
29
        case Scheme::Forward:
30
            return (f(x + hh) - f(x)) / hh; // O(h)
31
        case Scheme::Backward:
32
            return (f(x) - f(x - hh)) / hh; // O(h)
33
        case Scheme::Central:
34
        default:
35
            return (f(x + hh) - f(x - hh)) / (2.0 * hh); // O(h^2)
36
    }
37
}
38
39
int main() {
40
    // Test function: f(x) = sin(x), exact derivative f'(x) = cos(x)
41
    auto f = [](double x) { return std::sin(x); };
42
    auto fprime_exact = [](double x) { return std::cos(x); };
43
44
    double x = 1.0; // point of differentiation
45
46
    for (Scheme s : {Scheme::Forward, Scheme::Backward, Scheme::Central}) {
47
        double h = 0.0; // let the function choose a default h
48
        double d = derivative(f, x, h, s);
49
        double exact = fprime_exact(x);
50
        double abs_err = std::abs(d - exact);
51
52
        std::string name = (s == Scheme::Forward) ? "Forward" : (s == Scheme::Backward ? "Backward" : "Central");
53
        std::cout << std::fixed << std::setprecision(12);
54
        std::cout << name << " difference at x=" << x << ": " << d
55
                  << ", exact=" << exact << ", |error|=" << abs_err << "\n";
56
    }
57
58
    // Show sensitivity to h for central difference
59
    std::cout << "\nCentral difference error vs h:" << std::endl;
60
    for (double h = 1e-1; h >= 1e-12; h *= 0.1) {
61
        double d = derivative(f, x, h, Scheme::Central);
62
        double err = std::abs(d - fprime_exact(x));
63
        std::cout << "h=" << std::setw(10) << std::scientific << h
64
                  << ", error=" << err << std::fixed << std::endl;
65
    }
66
67
    return 0;
68
}
69

1
#include <iostream>
2
#include <vector>
3
#include <cmath>
4
#include <iomanip>
5
6
// Compute first derivative from sampled data y[i] = f(x[i]) on a uniform grid
7
// Uses forward difference at i=0, backward difference at i=n-1, and central otherwise
8
std::vector<double> differentiate_uniform(const std::vector<double>& x, const std::vector<double>& y) {
9
    size_t n = x.size();
10
    std::vector<double> dydx(n, 0.0);
11
    if (n < 2) return dydx;
12
13
    // Assume uniform grid and compute h from endpoints
14
    double h = (x.back() - x.front()) / (n - 1);
15
16
    // Forward at left boundary
17
    dydx[0] = (y[1] - y[0]) / h;
18
19
    // Central in the interior
20
    for (size_t i = 1; i + 1 < n; ++i) {
21
        dydx[i] = (y[i + 1] - y[i - 1]) / (2.0 * h);
22
    }
23
24
    // Backward at right boundary
25
    dydx[n - 1] = (y[n - 1] - y[n - 2]) / h;
26
27
    return dydx;
28
}
29
30
int main() {
31
    // Sample f(x) = exp(x) on [0,1]
32
    size_t n = 11; // 11 points
33
    std::vector<double> x(n), y(n);
34
    double a = 0.0, b = 1.0;
35
    double h = (b - a) / (n - 1);
36
    for (size_t i = 0; i < n; ++i) {
37
        x[i] = a + i * h;
38
        y[i] = std::exp(x[i]);
39
    }
40
41
    std::vector<double> dydx = differentiate_uniform(x, y);
42
43
    std::cout << std::fixed << std::setprecision(6);
44
    std::cout << "i\tx\tdy/dx (FD)\texact\t|error|\n";
45
    for (size_t i = 0; i < n; ++i) {
46
        double exact = std::exp(x[i]); // derivative of exp(x) is exp(x)
47
        double err = std::abs(dydx[i] - exact);
48
        std::cout << i << '\t' << x[i] << '\t' << dydx[i] << '\t' << exact << '\t' << err << '\n';
49
    }
50
51
    return 0;
52
}
53

1
#include <iostream>
2
#include <cmath>
3
#include <limits>
4
#include <iomanip>
5
6
// Central difference O(h^2)
7
double central_diff(const std::function<double(double)>& f, double x, double h) {
8
    return (f(x + h) - f(x - h)) / (2.0 * h);
9
}
10
11
// Richardson extrapolation: for p=2 (central diff), D4 = (4 D(h/2) - D(h)) / 3
12
// Result is O(h^4) if f is sufficiently smooth
13
double richardson_central(const std::function<double(double)>& f, double x, double h) {
14
    double D1 = central_diff(f, x, h);
15
    double D2 = central_diff(f, x, h / 2.0);
16
    return (4.0 * D2 - D1) / 3.0;
17
}
18
19
int main() {
20
    auto f = [](double x) { return std::sin(x); };
21
    auto fprime_exact = [](double x) { return std::cos(x); };
22
23
    double x = 1.0;
24
25
    std::cout << std::scientific << std::setprecision(3);
26
    std::cout << "h\t\tcentral(O(h^2)) error\tRichardson(O(h^4)) error\n";
27
    for (double h = 1e-1; h >= 1e-6; h *= 0.1) {
28
        double d2 = central_diff(f, x, h);
29
        double d4 = richardson_central(f, x, h);
30
        double e2 = std::abs(d2 - fprime_exact(x));
31
        double e4 = std::abs(d4 - fprime_exact(x));
32
        std::cout << h << "\t" << e2 << "\t\t\t" << e4 << "\n";
33
    }
34
35
    return 0;
36
}
37