Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 6

Numerical Differentiation

The Chapter 6 companion programs illustrate differentiation from tabulated data and the improvement of centered finite differences through Richardson extrapolation. The resources below connect these examples with widely used numerical tools and provide benchmark functions and practical guidance for choosing a method and step size.

Chapter 6 codes Chapter 6 exercises & projects Chapter 6 updates & errata Chapter 6 overview

Software and library resources

The table below lists representative tools for numerical differentiation. Their formulas, endpoint treatments, and error estimates are not necessarily identical, so check the documentation before comparing results or using a routine in a larger calculation.

On a narrow screen, scroll horizontally to view all language columns.

Numerical task MATLAB Python C++ Fortran
Derivatives of tabulated data First derivatives on uniform or nonuniform grids. gradient
Centered differences in the interior and first-order one-sided differences at the endpoints. On a nonuniform grid, do not assume second-order accuracy; use the companion code or nonuniform finite-difference weights when needed 		numpy.gradient
Use edge_order=2 when second-order endpoint formulas are desired 		Companion implementation or a custom local interpolation formula Companion implementation or a custom local interpolation formula
Derivative of a supplied function Finite differences with method-specific step selection and error estimates. Evaluate the function on a grid and use gradient, or use the companion implementation scipy.differentiate.derivative
Adaptive step size and an error estimate; requires SciPy 1.15 or later 		Boost finite_difference_derivative
Automatically selected step and optional error estimate
Boost complex_step_derivative for analytic functions;
GSL gsl_deriv_* (C library, callable from C++) 		Companion Richardson implementation, an implementation of Ridders' extrapolation method, or GSL through the FGSL Fortran bindings
Multidimensional gradient Partial derivatives of data sampled on a regular or specified grid. gradient numpy.gradient Apply one-dimensional finite-difference stencils along each coordinate direction Apply one-dimensional finite-difference stencils along each coordinate direction
Discrete Laplacian Approximation of the sum of second partial derivatives. del2
Multiply by 4 in two dimensions to obtain ∇2u; in N ≥ 2 dimensions multiply by 2N. At boundaries, del2 linearly extrapolates second differences from the interior 		scipy.ndimage.laplace
Uses unit index spacing and defaults to mode='reflect'; for uniform physical spacing divide the result by h2 		Custom finite-difference stencil with explicitly chosen boundary conditions Custom finite-difference stencil with explicitly chosen boundary conditions
Derivatives of noisy data Smoothing or local polynomial fitting before differentiation. sgolay for differentiation filters;
sgolayfilt for smoothing; equally spaced samples required
csaps for smoothing splines
Toolboxes required 		scipy.signal.savgol_filter
Set deriv and the constant sample spacing delta;
make_smoothing_spline 		Boost discrete_lanczos_derivative for equally spaced data Smoothing or local least-squares fitting followed by differentiation

For functions implemented as compatible computer code, automatic differentiation avoids finite-difference truncation error and does not require a step size. Examples include jax.grad, MATLAB dlgradient, and Boost.Math autodiff. The differentiated operations must be supported by the selected tool.

Benchmark functions and reference derivatives

These examples can be used to test a new implementation, compare finite-difference formulas, or study how the error changes with step size. Compare each numerical result with the exact value shown in the table.

Function and point Exact derivatives Why it is useful
f(x) = exp(x) at x = 1 f'(1) = e ≈ 2.718281828459045
f''(1) = e ≈ 2.718281828459045		 All derivatives equal e, so the exact value and leading error constants are especially simple.
f(x) = sin(x) at x = π/4 f'(π/4) = √2/2 ≈ 0.7071067811865476
f''(π/4) = -√2/2 ≈ -0.7071067811865476		 Checks signs and trigonometric derivatives.
f(x) = ln(x) at x = 2 f'(2) = 1/2
f''(2) = -1/4		 Matches the Richardson-extrapolation companion program.
f(x) = 1/(1+x2) at x = 1 f'(1) = -1/2
f''(1) = 1/2		 Tests a rational function with opposite derivative signs.
f(x) = x4-2x2+1 at x = 1/2 f'(1/2) = -3/2
f''(1/2) = -1		 The centered first- and second-derivative errors are both exactly 2h2. A five-point stencil or one Richardson step is exact apart from roundoff.
f(x) = tanh(20x) at x = 0.05 f'(0.05) = 20 sech2(1) ≈ 8.399486832280521
f''(0.05) = -800 sech2(1)tanh(1) ≈ -255.8800033796898		 A steep-gradient test in which both derivatives are nonzero. At x = 0, the centered second derivative would vanish identically by symmetry and would not test numerical accuracy.

When an exact derivative is zero, report absolute error; relative error is undefined or misleading.

Choosing a differentiation method

Character of the problem Method to consider Practical note
Smooth function available at arbitrary points Centered finite differences Usually more accurate than a forward or backward formula with the same step size and number of nearby evaluations.
Higher accuracy is needed Higher-order stencil or Richardson extrapolation Verify the expected convergence by repeating the calculation with smaller step sizes.
Point is at or near a boundary One-sided finite difference or a shifted interpolation stencil Do not evaluate the function outside its physical or tabulated domain.
Data are tabulated on a nonuniform grid Differentiate a local interpolation polynomial Use the actual coordinate values rather than formulas derived for constant spacing.
Data contain measurement noise Local polynomial fit, smoothing spline, or Savitzky–Golay method Differentiation amplifies high-frequency noise; decreasing the spacing usually makes the derivative estimate noisier.
Choosing the step size Plot error against h on a log–log scale Truncation error decreases with h, but subtractive cancellation and roundoff eventually increase the total error.
Function of several variables One-dimensional stencils applied along each coordinate Handle grid spacing and boundary conditions separately in each coordinate direction.
Mixed partial derivative Tensor-product stencil or successive one-dimensional derivatives Use consistent stencils, coordinate spacing, and boundary treatments when applying successive one-dimensional derivatives.
Analytic function that accepts complex arguments Complex-step first derivative Avoids subtractive cancellation, but every operation must preserve the complex perturbation. Operations such as abs, max, or real-valued conditional branches can invalidate the result.
Function implemented with operations supported by an automatic-differentiation tool Automatic differentiation No finite-difference step is required. The result differentiates the executed computer program, so numerical conditioning and unsupported operations still matter.

Starting estimates for the step size

For double-precision calculations, let ε denote machine epsilon and let L be a characteristic scale for x, often max(1,|x|). Useful initial estimates are:

  • forward difference: h ≈ ε1/2L, typically about 10-8L;
  • centered first derivative: h ≈ ε1/3L, typically about 10-5L;
  • centered second derivative: h ≈ ε1/4L, typically about 10-4L;
  • fourth-order first derivative: h ≈ ε1/5L, typically about 10-3L;
  • fourth-order second derivative: h ≈ ε1/6L, typically a few times 10-3L.

These are starting values rather than universal optima; the best step also depends on the function and its derivative scales. In floating-point arithmetic, use the representable increment hactual = (x + hrequested) - x.

A complex-step first derivative permits a much smaller step because it avoids subtractive cancellation. The practical lower limit depends on underflow and on whether every operation preserves the imaginary perturbation, so an extremely small universal value should not be prescribed.

Selected references and documentation

← Back to Chapter 6