In-class: Second-Derivative

Let’s investigate a second-derivative, \(d^2f/dx^2\). We can start with the same Taylor expansions:

\[f(x_0 + \Delta x) = f(x_0) + \left . \frac{df}{dx} \right |_{x_0} \Delta x + \frac{1}{2} \left . \frac{d^2f}{dx^2} \right |_{x_0} \Delta x^2 + \mathcal{O}(\Delta x^3)\]

\[f(x_0 - \Delta x) = f(x_0) - \left . \frac{df}{dx} \right |_{x_0} \Delta x + \frac{1}{2} \left . \frac{d^2f}{dx^2} \right |_{x_0} \Delta x^2 - \mathcal{O}(\Delta x^3)\]

Now, instead of subtracting them, we add them:

\[f(x_0 + \Delta x) + f(x_0 - \Delta x) = 2 f(x_0) + \left . \frac{d^2f}{dx^2} \right |_{x_0} \Delta x^2 + \mathcal{O}(\Delta x^4)\]

Notice that the \(\Delta x\) and \(\Delta x^3\) terms cancel out. Now solving for the second derivative, we have:

\[\left . \frac{d^2f}{dx^2} \right |_{x_0} = \frac{f(x_0 + \Delta x) - 2 f(x_0) + f(x_0 - \Delta x)}{\Delta x^2} + \mathcal{O}(\Delta x^2)\]

We see that this is second-order accurate.

try it…

Let’s code up the second derivative and test it on

\[f(x) = x^4 + x^2 + 1\]

This has the second derivative:

\[f^{\prime\prime}(x) = 12 x^2 + 2\]

solution

Listing 68 second_derivative.cpp

#include <iostream>
#include <format>
#include <cmath>

// the function we are differentiating
double f(double x) {
    return std::pow(x, 4.0) + std::pow(x, 2.0) + 1.0;
}

// the analytic second-derivative of f(x)
double df2dx2(double x) {
    return 12.0 * x * x + 2.0;
}

int main() {

    // the x-value where we want to compute the second-derivative
    double x0{1.5};

    // the analytic second derivative at x0
    double df2dx2_true = df2dx2(x0);

    // output a header -- note the use of "^" to center these headings
    // in the column
    std::cout << std::format("{:^12} {:^12} {:^12}\n",
                             "dx", "d2f/dx2", "error");

    for (double dx : {0.1, 0.05, 0.025, 0.0125}) {

        // our 2nd-order accurate second-derivative finite-difference
        double df2dx2_approx = (f(x0 + dx) - 2.0 * f(x0) + f(x0 - dx)) /
            (dx * dx);

        // the absolute error
        double err = std::abs(df2dx2_true - df2dx2_approx);

        std::cout << std::format("{:12.6g} {:12.6g} {:12.6g}\n",
                                 dx, df2dx2_approx, err);
    }

}