Numerical Differentiation

Why Numerical Differentiation Matters

A simulator can evaluate a system's response to a control setting, but it does not expose a derivative. To estimate how sensitive the response is to a small change, evaluate the simulator at nearby settings and form a slope. Position samples used to estimate velocity, finite-difference gradient checks, and sensitivity analysis all use this idea.

This lesson estimates a first derivative of a real scalar function, $f'(x)$, using forward, backward, and central differences, then improves a sequence of central differences with Richardson extrapolation. You will trace every sample and table entry on an exact example, implement the methods in Python, C++, and Java, and distinguish a small change between estimates from evidence of actual accuracy.

Numerical differentiation is sensitive to noise: subtracting nearby values can discard significant digits, and dividing by a small step amplifies what remains. A smaller step reduces truncation error only until evaluation error or rounding becomes dominant. If a trustworthy analytic derivative or automatic differentiation is available, it often avoids that tradeoff.

Core Concepts and Terminology

For a positive step $h$, define

$$ D_+(h)=\frac{f(x+h)-f(x)}{h},\qquad D_-(h)=\frac{f(x)-f(x-h)}{h},\qquad D_c(h)=\frac{f(x+h)-f(x-h)}{2h}. $$

Forward and backward differences use a sample at the target and one on a chosen side. Central difference uses two samples straddling the target and does not evaluate $f(x)$. These are secant slopes approximating a tangent slope. They are also dot products between sample vectors and stencil weights; they do not require solving a linear system or pivoting a matrix.

For $f(x)=x^3$ at $x=1$ with $h=1/2$, the sample values are $1/8$, $1$, and $27/8$. The backward slope is $7/4$, the forward slope is $19/4$, and the central slope is $13/4$. The true derivative is $3$. Central difference has less error on this example, though a one-sided rule may be necessary near a domain boundary.

The computed point x + h may differ slightly from its mathematical value, or even equal x. The implementations use the textbook nominal denominators $h$ and $2h$, check that required sample locations are finite and strictly separated, and report when a stencil cannot be formed. Separation alone does not establish that the realized spacings remain accurately symmetric; this matters near floating-point resolution.

How Numerical Differentiation Works, Step by Step

1. Validate the target and step. Require finite $x$ and a positive finite $h$. Choose a side compatible with the function's domain. For a central stencil, both $x-h$ and $x+h$ must be valid, finite, distinct points on the appropriate sides, and $2h$ must be finite. The code can check representability, but only the caller knows the function's true domain and smoothness.

2. Evaluate the selected stencil and form its weighted difference. Each basic rule uses exactly two calls. A nonfinite value or checked arithmetic result raises an exception. A function can also throw its own domain exception. When several schemes are compared independently, the supplied programs do not share their evaluations; a production implementation can cache common values.

Taylor expansion explains the difference in accuracy. For sufficiently smooth $f$,

$$ \begin{aligned} D_+(h)&=f'(x)+\frac{h}{2}f''(x)+O(h^2),\\ D_-(h)&=f'(x)-\frac{h}{2}f''(x)+O(h^2),\\ D_c(h)&=f'(x)+\frac{h^2}{6}f'''(x)+\frac{h^4}{120}f^{(5)}(x)+O(h^6). \end{aligned} $$

The symmetric central stencil cancels the even terms in the two function expansions before division, leaving an error series in even powers of $h$. That structure is what the central Richardson algorithm uses.

3. Compute two central differences and cancel their leading error. If $D_c(h)=d+c_2h^2+c_4h^4+\cdots$, where $d=f'(x)$, then

$$ R=\frac{4D_c(h/2)-D_c(h)}{3} =D_c(h/2)+\frac{D_c(h/2)-D_c(h)}{3} =d+O(h^4). $$

The correction form on the middle expression is used in the code. The first factor is $4$ because the central rule's leading error scales as $h^2$. A forward-difference sequence instead starts with a factor of $2$ for its $O(h)$ term; using the central factor indiscriminately would cancel the wrong model.

4. Extend to a triangular Richardson table. Let $h_i=h_0/2^i$. Start row $i$ with $R_{i,0}=D_c(h_i)$, then compute

$$ R_{i,j}=R_{i,j-1} +\frac{R_{i,j-1}-R_{i-1,j-1}}{4^j-1}, \qquad 1\le j\le i. $$

The divisors are $3,15,63,\ldots$. Each new column removes one further even-power error term, under the smooth expansion assumption. In exact arithmetic, for fixed column $j$, $R_{i,j}=f'(x)+O(h_i^{2j+2})$ as the step decreases and sufficient derivatives exist. Higher columns do not promise better floating-point answers when the underlying samples are noisy.

5. Stop with a reason. After each row beyond the first, compare consecutive diagonal entries using

$$ \delta_i=\lvert R_{i,i}-R_{i-1,i-1}\rvert,\qquad \delta_i\le\text{atol}+\text{rtol}\,\lvert R_{i,i}\rvert. $$

If this holds, the supplied solver returns status estimated. This explicitly names a diagnostic, rather than a certified error bound. If the row limit is reached first, return limit. If a finer stencil becomes unrepresentable after at least one usable row, return the latest diagonal with stalled; if the very first stencil is unusable, raise an exception. When only one row exists, no diagonal change can be measured, so its reported change is infinity.

flowchart TD
    A["Validate x, initial step, row limit, and tolerances"] --> B{"Current central stencil representable?"}
    B -->|No| C{"Any completed row?"}
    C -->|No| D["Raise invalid-stencil error"]
    C -->|Yes| E["Return latest row with stalled status"]
    B -->|Yes| F["Evaluate two samples and central difference"]
    F --> G["Build the row with factors 4, 16, 64, ..."]
    G --> H{"Previous diagonal available?"}
    H -->|Yes| I["Compute diagonal change and tolerance threshold"]
    I --> J{"Change within threshold?"}
    J -->|Yes| K["Return estimated status"]
    J -->|No| L{"Row budget remains?"}
    H -->|No| L
    L -->|No| M["Return limit status"]
    L -->|Yes| N["Halve the step"]
    N --> B

The diagram assumes finite evaluated values and arithmetic; otherwise an exception interrupts the run. The programs check the combined tolerance threshold for finiteness as well. Both tolerances may be zero, in which case the stopping diagnostic requires exact equality of consecutive computed diagonals.

Worked Example with a Full Trace

Differentiate $f(x)=x^3$ at $x=1$, where $f'(1)=3$. Choose $h_0=1/2$, maximum rows $6$, $\text{atol}=10^{-12}$, and $\text{rtol}=0$. These dyadic positions and cubic values are exactly representable, making the cancellation easy to inspect without incidental rounding.

First, the coarse basic rules give

$$ D_+(1/2)=\frac{27/8-1}{1/2}=\frac{19}{4},\qquad D_-(1/2)=\frac{1-1/8}{1/2}=\frac74,\qquad D_c(1/2)=\frac{27/8-1/8}{1}=\frac{13}{4}. $$

The Richardson run is separate from those three basic calls. Its complete evaluation trace contains six samples, in the following order within the three rows; $f(1)$ is not needed.

Now fill every table entry and test the stopping rule:

For a cubic, $D_c(h)=3+h^2$ at this target, so the first extrapolation column cancels its entire error. The algorithm nevertheless needs the third row to observe agreement between successive diagonals. The returned value is $3$, the change diagnostic is $0$, and the count is $6$ function evaluations over $3$ rows.

The actual errors of the coarse forward, backward, and central rules are $7/4$, $5/4$, and $1/4$. These relatively large errors are intentional: the example uses a large, readable step to show the error terms. They should not be hidden behind a single success label.

Implementation in Python, C++, and Java

Each program provides difference(f, x, h, scheme) for forward, backward, or central, plus richardson based specifically on the central rule. Richardson accepts a maximum row count between $2$ and $12$ and finite nonnegative absolute and relative tolerances. It returns value, change, evaluations, status, and the completed triangular table.

The table uses plain nested arrays, with row $i$ containing $i+1$ entries. The computation follows a recurrence rather than matrix elimination; no pivoting is part of this algorithm. Every central row evaluates two new points. The code preserves every completed row so the complete trace is reviewable, rather than storing only a final estimate.

Finite function values and checked arithmetic results are required. Invalid parameters, an unusable first stencil, unknown basic scheme, or a nonfinite result raise an exception. A later loss of representable points returns stalled without making further function calls. A nonfinite function evaluation is an exception even if earlier rows were usable. The latest diagonal on a limit or stalled exit is a reported estimate, not a claim that it is the most accurate row.

All three drivers use $f(x)=x^3$, $x=1$, $h=1/2$, a maximum of $6$ Richardson rows, absolute tolerance $10^{-12}$, and relative tolerance $0$. They check the three expected finite-step values, then independently compare Richardson's result with the analytic derivative $3$. This does not claim the coarse one-sided values are within $10^{-12}$ of the derivative.

import math
from dataclasses import dataclass


def finite(value):
    if not math.isfinite(value):
        raise ValueError("nonfinite function value or arithmetic result")
    return value


def validate(x, h):
    if not math.isfinite(x) or not math.isfinite(h) or h <= 0.0:
        raise ValueError("finite x and positive finite h required")


def stencil_ok(x, h, scheme):
    if not math.isfinite(h) or h <= 0.0:
        return False
    plus, minus = x+h, x-h
    if scheme in ("forward", "central") and (not math.isfinite(plus) or plus <= x):
        return False
    if scheme in ("backward", "central") and (not math.isfinite(minus) or minus >= x):
        return False
    return scheme != "central" or math.isfinite(2.0*h)


def difference(f, x, h, scheme):
    validate(x, h)
    if scheme not in ("forward", "backward", "central"):
        raise ValueError("unknown difference scheme")
    if not stencil_ok(x, h, scheme):
        raise ValueError("unrepresentable stencil or denominator")
    if scheme == "forward":
        fx, fp = finite(f(x)), finite(f(x+h))
        result = (fp-fx)/h
    elif scheme == "backward":
        fx, fm = finite(f(x)), finite(f(x-h))
        result = (fx-fm)/h
    else:
        fp, fm = finite(f(x+h)), finite(f(x-h))
        result = (fp-fm)/(2.0*h)
    return finite(result)


@dataclass
class Richardson:
    value: float
    change: float
    evaluations: int
    status: str
    table: list


def richardson(f, x, h, levels=8, atol=1e-12, rtol=0.0):
    validate(x, h)
    if (not isinstance(levels, int) or not 2 <= levels <= 12
            or not math.isfinite(atol) or not math.isfinite(rtol)
            or atol < 0.0 or rtol < 0.0):
        raise ValueError("levels in [2, 12] and finite nonnegative tolerances required")
    table, change, step = [], math.inf, h
    for i in range(levels):
        if not stencil_ok(x, step, "central"):
            if not table:
                raise ValueError("first stencil is unrepresentable")
            return Richardson(table[-1][-1], change, 2*len(table), "stalled", table)
        # Row i starts with the central difference at h / 2**i.
        row = [difference(f, x, step, "central")]
        factor = 4.0
        for j in range(1, i+1):
            correction = finite(row[j-1]-table[i-1][j-1])/(factor-1.0)
            row.append(finite(row[j-1]+correction))
            factor *= 4.0
        table.append(row)
        if i > 0:
            change = abs(finite(row[-1]-table[i-1][-1]))
            threshold = finite(atol+rtol*abs(row[-1]))
            if change <= threshold:
                return Richardson(row[-1], change, 2*len(table), "estimated", table)
        step *= 0.5
    return Richardson(table[-1][-1], change, 2*len(table), "limit", table)


def main():
    def f(x):
        return (x*x)*x

    x, h, exact = 1.0, 0.5, 3.0
    # Rows are [computed difference, expected finite-step value].
    # These are scalar stencils and an extrapolation table, with no matrix pivoting.
    checks = [[difference(f,x,h,"forward"),4.75],
              [difference(f,x,h,"backward"),1.75],
              [difference(f,x,h,"central"),3.25]]
    result = richardson(f,x,h,6,1e-12,0.0)
    ok = all(abs(row[0]-row[1]) <= 1e-12 for row in checks)
    ok = ok and result.status == "estimated" and abs(result.value-exact) <= 1e-12
    for name, row in zip(["forward","backward","central"],checks):
        print(f"{name}: derivative = {row[0]:.6f}, actual error = {abs(row[0]-exact):.6f}")
    print(f"richardson: derivative = {result.value:.6f}, actual error = {abs(result.value-exact):.6f}")
    print(f"richardson status = {result.status}")
    print(f"richardson rows = {len(result.table)}, evaluations = {result.evaluations}")
    print(f"diagonal change = {result.change:.6f}")
    print(f"check = {'PASS' if ok else 'FAIL'}")


if __name__ == "__main__":
    main()

#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
using Fn=std::function<double(double)>;

double finite_value(double value) {
    if (!std::isfinite(value)) throw std::invalid_argument("nonfinite function value or arithmetic result");
    return value;
}
void validate(double x,double h) {
    if (!std::isfinite(x) || !std::isfinite(h) || h<=0)
        throw std::invalid_argument("finite x and positive finite h required");
}
bool stencil_ok(double x,double h,const std::string& scheme) {
    if (!std::isfinite(h) || h<=0) return false;
    double plus=x+h,minus=x-h;
    if ((scheme=="forward" || scheme=="central") && (!std::isfinite(plus) || plus<=x)) return false;
    if ((scheme=="backward" || scheme=="central") && (!std::isfinite(minus) || minus>=x)) return false;
    return scheme!="central" || std::isfinite(2.0*h);
}
double difference(const Fn& f,double x,double h,const std::string& scheme) {
    validate(x,h);
    if (scheme!="forward" && scheme!="backward" && scheme!="central")
        throw std::invalid_argument("unknown difference scheme");
    if (!stencil_ok(x,h,scheme)) throw std::invalid_argument("unrepresentable stencil or denominator");
    double result;
    if (scheme=="forward") {
        double fx=finite_value(f(x)),fp=finite_value(f(x+h));result=(fp-fx)/h;
    } else if (scheme=="backward") {
        double fx=finite_value(f(x)),fm=finite_value(f(x-h));result=(fx-fm)/h;
    } else {
        double fp=finite_value(f(x+h)),fm=finite_value(f(x-h));result=(fp-fm)/(2.0*h);
    }
    return finite_value(result);
}
struct Richardson {
    double value,change;int evaluations;std::string status;
    std::vector<std::vector<double>> table;
};
Richardson richardson(const Fn& f,double x,double h,int levels=8,double atol=1e-12,double rtol=0.0) {
    validate(x,h);
    if (levels<2 || levels>12 || !std::isfinite(atol) || !std::isfinite(rtol) || atol<0 || rtol<0)
        throw std::invalid_argument("levels in [2, 12] and finite nonnegative tolerances required");
    std::vector<std::vector<double>> table;
    double change=std::numeric_limits<double>::infinity(),step=h;
    for (int i=0;i<levels;++i) {
        if (!stencil_ok(x,step,"central")) {
            if (table.empty()) throw std::invalid_argument("first stencil is unrepresentable");
            return {table.back().back(),change,2*static_cast<int>(table.size()),"stalled",table};
        }
        // Row i starts with the central difference at h / 2^i.
        std::vector<double> row={difference(f,x,step,"central")};double factor=4.0;
        for (int j=1;j<=i;++j) {
            double correction=finite_value(row[j-1]-table[i-1][j-1])/(factor-1.0);
            row.push_back(finite_value(row[j-1]+correction));factor*=4.0;
        }
        table.push_back(row);
        if (i>0) {
            change=std::abs(finite_value(row.back()-table[i-1].back()));
            double threshold=finite_value(atol+rtol*std::abs(row.back()));
            if (change<=threshold)
                return {row.back(),change,2*static_cast<int>(table.size()),"estimated",table};
        }
        step*=0.5;
    }
    return {table.back().back(),change,2*static_cast<int>(table.size()),"limit",table};
}
int main() {
    Fn f=[](double x){return (x*x)*x;};
    double x=1.0,h=0.5,exact=3.0;
    // Rows are [computed difference, expected finite-step value].
    // Scalar stencils and an extrapolation table require no matrix pivoting.
    std::vector<std::vector<double>> checks={ {difference(f,x,h,"forward"),4.75},
        {difference(f,x,h,"backward"),1.75}, {difference(f,x,h,"central"),3.25} };
    Richardson result=richardson(f,x,h,6,1e-12,0.0);bool ok=true;
    for (const auto& row:checks) ok=ok && std::abs(row[0]-row[1])<=1e-12;
    ok=ok && result.status=="estimated" && std::abs(result.value-exact)<=1e-12;
    std::vector<std::string> names={"forward","backward","central"};
    std::cout<<std::fixed<<std::setprecision(6);
    for (int i=0;i<3;++i)
        std::cout<<names[i]<<": derivative = "<<checks[i][0]<<", actual error = "<<std::abs(checks[i][0]-exact)<<'\n';
    std::cout<<"richardson: derivative = "<<result.value<<", actual error = "<<std::abs(result.value-exact)<<'\n';
    std::cout<<"richardson status = "<<result.status<<'\n';
    std::cout<<"richardson rows = "<<result.table.size()<<", evaluations = "<<result.evaluations<<'\n';
    std::cout<<"diagonal change = "<<result.change<<'\n';
    std::cout<<"check = "<<(ok?"PASS":"FAIL")<<'\n';
    return 0;
}

import java.util.ArrayList;
import java.util.Locale;
import java.util.function.DoubleUnaryOperator;

public class Main {
    static double finite(double value) {
        if (!Double.isFinite(value)) throw new IllegalArgumentException("nonfinite function value or arithmetic result");
        return value;
    }
    static void validate(double x,double h) {
        if (!Double.isFinite(x) || !Double.isFinite(h) || h<=0)
            throw new IllegalArgumentException("finite x and positive finite h required");
    }
    static boolean stencilOk(double x,double h,String scheme) {
        if (!Double.isFinite(h) || h<=0) return false;
        double plus=x+h,minus=x-h;
        if ((scheme.equals("forward") || scheme.equals("central")) && (!Double.isFinite(plus) || plus<=x)) return false;
        if ((scheme.equals("backward") || scheme.equals("central")) && (!Double.isFinite(minus) || minus>=x)) return false;
        return !scheme.equals("central") || Double.isFinite(2.0*h);
    }
    static double difference(DoubleUnaryOperator f,double x,double h,String scheme) {
        validate(x,h);
        if (!scheme.equals("forward") && !scheme.equals("backward") && !scheme.equals("central"))
            throw new IllegalArgumentException("unknown difference scheme");
        if (!stencilOk(x,h,scheme)) throw new IllegalArgumentException("unrepresentable stencil or denominator");
        double result;
        if (scheme.equals("forward")) {
            double fx=finite(f.applyAsDouble(x)),fp=finite(f.applyAsDouble(x+h));result=(fp-fx)/h;
        } else if (scheme.equals("backward")) {
            double fx=finite(f.applyAsDouble(x)),fm=finite(f.applyAsDouble(x-h));result=(fx-fm)/h;
        } else {
            double fp=finite(f.applyAsDouble(x+h)),fm=finite(f.applyAsDouble(x-h));result=(fp-fm)/(2.0*h);
        }
        return finite(result);
    }
    static class Richardson {
        double value,change;int evaluations;String status;ArrayList<double[]> table;
        Richardson(double value,double change,int evaluations,String status,ArrayList<double[]> table) {
            this.value=value;this.change=change;this.evaluations=evaluations;this.status=status;this.table=table;
        }
    }
    static Richardson richardson(DoubleUnaryOperator f,double x,double h,int levels,double atol,double rtol) {
        validate(x,h);
        if (levels<2 || levels>12 || !Double.isFinite(atol) || !Double.isFinite(rtol) || atol<0 || rtol<0)
            throw new IllegalArgumentException("levels in [2, 12] and finite nonnegative tolerances required");
        ArrayList<double[]> table=new ArrayList<>();double change=Double.POSITIVE_INFINITY,step=h;
        for (int i=0;i<levels;++i) {
            if (!stencilOk(x,step,"central")) {
                if (table.isEmpty()) throw new IllegalArgumentException("first stencil is unrepresentable");
                double[] last=table.get(table.size()-1);
                return new Richardson(last[last.length-1],change,2*table.size(),"stalled",table);
            }
            // Row i starts with the central difference at h / 2^i.
            double[] row=new double[i+1];row[0]=difference(f,x,step,"central");double factor=4.0;
            for (int j=1;j<=i;++j) {
                double correction=finite(row[j-1]-table.get(i-1)[j-1])/(factor-1.0);
                row[j]=finite(row[j-1]+correction);factor*=4.0;
            }
            table.add(row);
            if (i>0) {
                double[] previous=table.get(i-1);
                change=Math.abs(finite(row[i]-previous[previous.length-1]));
                double threshold=finite(atol+rtol*Math.abs(row[i]));
                if (change<=threshold) return new Richardson(row[i],change,2*table.size(),"estimated",table);
            }
            step*=0.5;
        }
        double[] last=table.get(table.size()-1);
        return new Richardson(last[last.length-1],change,2*table.size(),"limit",table);
    }
    public static void main(String[] args) {
        Locale.setDefault(Locale.ROOT);DoubleUnaryOperator f=x -> (x*x)*x;
        double x=1.0,h=0.5,exact=3.0;
        // Rows are [computed difference, expected finite-step value].
        // Scalar stencils and an extrapolation table require no matrix pivoting.
        double[][] checks={ {difference(f,x,h,"forward"),4.75},
            {difference(f,x,h,"backward"),1.75}, {difference(f,x,h,"central"),3.25} };
        Richardson result=richardson(f,x,h,6,1e-12,0.0);boolean ok=true;
        for (double[] row:checks) ok=ok && Math.abs(row[0]-row[1])<=1e-12;
        ok=ok && result.status.equals("estimated") && Math.abs(result.value-exact)<=1e-12;
        String[] names={"forward","backward","central"};
        for (int i=0;i<3;++i)
            System.out.printf("%s: derivative = %.6f, actual error = %.6f%n",names[i],checks[i][0],Math.abs(checks[i][0]-exact));
        System.out.printf("richardson: derivative = %.6f, actual error = %.6f%n",result.value,Math.abs(result.value-exact));
        System.out.println("richardson status = "+result.status);
        System.out.printf("richardson rows = %d, evaluations = %d%n",result.table.size(),result.evaluations);
        System.out.printf("diagonal change = %.6f%n",result.change);
        System.out.println("check = "+(ok?"PASS":"FAIL"));
    }
}

The identical printed output is:

• forward: derivative = 4.750000, actual error = 1.750000
• backward: derivative = 1.750000, actual error = 1.250000
• central: derivative = 3.250000, actual error = 0.250000
• richardson: derivative = 3.000000, actual error = 0.000000
• richardson status = estimated
• richardson rows = 3, evaluations = 6
• diagonal change = 0.000000
• check = PASS

The three separate basic computations use six calls in total; Richardson uses a further six. Its reported evaluation count covers only its own run. In applications where no analytic derivative is known, actual error cannot be printed this way; inspect a step-size sweep, compare independent methods, and account for the function evaluator's error instead.

Complexity, Correctness, and Invariants

Let $C_f$ be the cost of evaluating the scalar function and $m$ the number of completed Richardson rows, with a configured limit $L$. The stated orders concern a smooth function, exact arithmetic, and suitably symmetric nodes. A finite step error bound additionally requires derivative bounds on the stencil's neighborhood.

For a valid basic stencil the best, average, and worst work are identical: always two evaluations. Richardson's earliest estimated exit uses two rows and four evaluations. A later representability failure can return after just one row; an invalid first stencil fails before evaluation. The worst completed run uses $2L$ evaluations, $O(LC_f+L^2)$ time, and $\Theta(L^2)$ table space. Its average row count depends on the function, step, tolerances, and evaluator error; there is no distribution-free average. Keeping only two rows reduces working space to $O(L)$ if the full table is not required.

Why each stencil approximates the derivative. Taylor expansion around $x$ identifies the derivative term and the first omitted remainder. Forward and backward formulas retain an $O(h)$ contribution; symmetric subtraction cancels the $h^2f''$ terms in the function values, yielding an $O(h^2)$ derivative error. A central derivative stencil can even be exact for a quadratic despite using only two off-center values.

Why Richardson cancels the error model. Suppose a column has leading error $c h^{2j}$. The value at half the step has leading error $c h^{2j}/4^j$. Combining the two as $[4^jA(h/2)-A(h)]/(4^j-1)$ removes that term while preserving the common derivative. Induction across columns proves the formal order invariant when sufficiently many even-power terms exist. This is an algebraic cancellation argument; the stopping test itself does not prove a derivative exists or that rounding is small.

Why the step cannot shrink indefinitely. A useful simplified total-error model is $E(h)\approx Ah^p+B uF/h$, where $p$ is the truncation order, $u$ characterizes relative floating-point error, and $F$ represents a local function-value scale. Balancing the terms gives $h$ proportional to $(uF/A)^{1/(p+1)}$, with constants depending on the stencil. Under well-scaled smooth behavior, this suggests square-root precision scaling for a first-order rule and cube-root precision scaling for a central rule. The constants, units, derivative magnitudes, and noise level matter; these are starting heuristics rather than universal step formulas.

For bounded additive errors of magnitude at most $\eta$ in each function value, the induced worst-case error is at most $2\eta/h$ for a basic one-sided difference and $\eta/h$ for a central difference. Extrapolation combines additional noisy values and can amplify their errors. The programs make no claim that their change tolerance controls these contributions.

Variants, Pitfalls, and Optimizations

Check representable positions before evaluating. In binary64 arithmetic, $1+10^{-20}$ and $1-10^{-20}$ both round to $1$. Blindly evaluating a difference there could produce zero and suggest a flat function. The code rejects an initial collapsed stencil or reports stalled if refinement reaches that point after usable rows. Even before collapse, realized offsets may differ from the nominal $h$ or from each other; exact-grid order statements then cease to describe all rounding effects.

Separated nodes still do not guarantee useful values. For $f(x)=10^{16}+x$ at $x=0$, steps $h=1/2$ and $h=1/4$ produce distinct input points, but the computed function values all round to $10^{16}$. Central differences and Richardson return zero with an estimated status, although the mathematical derivative is $1$. Rewriting the evaluator to avoid losing the small varying component, using more precision, or using a larger step can help. Subtracting a constant after the values have already rounded cannot recover the lost component.

Smoothness is essential. For $\lvert x\rvert$ at $0$, symmetric estimates agree at zero even though the derivative is undefined. Compare one-sided behavior around suspected kinks and respect domain boundaries. At a left boundary, the second-order forward formula $[-3f(x)+4f(x+h)-f(x+2h)]/(2h)$ offers higher order if all three forward points are valid and the function is smooth. That variant is not implemented in the supplied programs.

Use the right Richardson factor. The general cancellation for leading order $p$ uses $2^p$ when the step halves. Central differences have an even-power series, so successive factors are $4,16,64,\ldots$. Forward differences generally have all integer powers, so their corresponding sequence starts $2,4,8,\ldots$. Do not substitute $4^j$ into an arbitrary black-box approximation without deriving its error expansion.

Measured data needs a noise model. If positions are noisy, their finite-difference velocities can be much noisier. Local polynomial fits, smoothing splines, or Savitzky–Golay differentiation can trade bias against noise amplification. For uneven samples, derive weights for the actual coordinates rather than pretending they form the grid $x-h,x,x+h$. Choosing a larger stencil can improve formal order but does not automatically improve noisy-data accuracy.

Consider alternatives and extensions. Automatic differentiation propagates derivatives through supported program operations and avoids a finite-difference step; branch behavior and differentiability still need attention. Complex-step differentiation, $\operatorname{Im}f(x+ih)/h$, avoids subtractive cancellation for a compatible analytic complex extension, but real-only code, absolute values, and nonanalytic branches can invalidate it. Higher derivatives and Jacobians require their own stencils, scaling, and error analysis; for $d$ input coordinates, a basic central gradient needs $2d$ scalar function calls.

Optimize only after validating behavior. Cache repeated center values when comparing one-sided rules, batch independent stencil evaluations when calls are expensive, and use scaled coordinates when input units differ greatly. A step sweep that reveals an initial convergence regime followed by worsening estimates is more informative than increasing Richardson's row limit indefinitely. The supplied solver caps rows and exposes its table so that behavior remains inspectable.

Practice Problems and Real-World Use

Start with functions whose derivatives are known independently, then introduce boundaries, rounding, and noise one at a time. Record the step, actual stencil points, function values, estimate, and analytic error when available. For Richardson, inspect the entire triangle rather than only its last diagonal.

For further reading, NIST's numerical differentiation chapter develops derivative formulas from interpolation at equally spaced nodes. SciPy's derivative documentation describes finite-difference derivative estimation, termination diagnostics, and tolerance considerations. Its implementation and defaults differ from the explicit central Richardson table used here; compare the documented contracts rather than assuming the methods are interchangeable.