Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 6

Numerical Differentiation

These exercises and projects supplement those in the book. They introduce different functions, data sets, and physical applications, with particular attention to irregular grids, noisy data, error estimation, boundary treatment, and derivatives in more than one dimension.

Chapter 6 codes Chapter 6 resources Chapter 6 updates & errata Chapter 6 overview

Additional exercises

  1. Truncation error and roundoff

    Calculate the first derivative of f(x) = exp(x) at x = 1 using the forward-difference formula

    f′(x) ≈ [f(x+h)−f(x)]/h

    and the centered-difference formula

    f′(x) ≈ [f(x+h)−f(x−h)]/(2h).

    Use h = 10−1, 10−2, ..., 10−16 and plot the absolute error versus h on logarithmic axes. Identify the ranges with the expected first- and second-order convergence, find the step size that gives the smallest error for each formula, and explain why further reductions in h eventually increase the error. In particular, explain why 1 + 10−16 rounds to 1 in double precision and examine the representable increment hactual = (x + h) − x.

  2. Differentiation on an irregular grid

    Generate 41 points on 0 ≤ x ≤ 3, then displace each interior point randomly by up to 20% of the uniform-grid spacing; keep the two endpoints fixed. Tabulate

    f(x) = ln(1 + x2).

    Calculate the first derivative using local Lagrange interpolation with the actual coordinates. Compare quadratic three-point and quartic five-point local interpolants, including suitable shifted stencils near the endpoints. Compare both with a calculation that incorrectly assumes uniform spacing. Repeat the experiment for several random grids and summarize the RMS and maximum errors.

  3. When differentiation amplifies noise

    Sample the signal

    y(t) = sin(2πt) + 0.20 sin(10πt),    0 ≤ t ≤ 2,

    using Δt = 0.005, and add normally distributed random noise with standard deviation σ = 0.001, 0.01, and 0.05. Although the high-frequency term has only 20% of the signal amplitude, its derivative amplitude is 0.20(10π) = 2π, the same as that of the main term. Estimate dy/dt using centered differences, then repeat after applying a moving polynomial or Savitzky–Golay filter. Vary the sampling interval, noise level, and smoothing-window width. For centered differences, estimate the derivative noise using σ/(√2Δt); it is about 7 when σ = 0.05 and Δt = 0.005, slightly larger than the true derivative amplitude 2π. Determine when smoothing improves the derivative and when it removes genuine high-frequency structure from the signal.

  4. A multilevel Richardson table

    For f(x) = xe−x at x = 0.5, begin with h = 0.4 and calculate centered first- and second-derivative approximations using h, h/2, h/4, and h/8. Construct a Richardson table with more than one extrapolation level. The exact derivatives are f′(x) = (1−x)e−x and f″(x) = (x−2)e−x.

    Compare each extrapolated result with the exact derivative. Estimate the observed order from successive differences and determine when further extrapolation stops improving the result. Compare the levels at which roundoff begins to limit the first- and second-derivative estimates; they need not be the same.

  5. Boundary formulas and periodic data

    Tabulate f(x) = cos(3x) on 0 ≤ x < 2π using N = 32 and N = 64 equally spaced points. Compare three treatments near the ends of the stored array: first-order one-sided differences, second-order one-sided differences, and centered differences with periodic wraparound.

    Plot the pointwise error over the entire grid. Explain why the periodic treatment is appropriate here and why it would be incorrect for a general nonperiodic data set.

  6. Gradient and Laplacian in two dimensions

    Use the manufactured field

    u(x,y) = ex+y + sin(πx) cos(2πy),    0 ≤ x,y ≤ 1.

    Compute both components of the gradient, the Laplacian, and the mixed derivative ∂2u/∂x∂y on a sequence of two-dimensional grids. Compare with the exact derivatives and determine the observed convergence rate. Repeat the calculation with unequal spacings in the two coordinate directions. First report errors on interior points only. Then include boundary points using second-order one-sided formulas and determine whether the boundary treatment changes the observed convergence rate.

  7. Finding stationary and inflection points from data

    Sample the function

    f(x) = e−(x−1)2 + 0.35e−8(x+0.5)2,    −2 ≤ x ≤ 3.

    Use numerical first and second derivatives to locate local maxima, minima, and inflection points. Refine each location by interpolating the derivative through zero. Study how the inferred locations change with the grid spacing. Begin with 101 equally spaced points and repeat with 201 and 401 points. Explain how you verified that no stationary or inflection points were missed.

    For validation, use a root finder on the analytic derivatives. The stationary points are maxima near x = −0.43164 and x = 1.00000 and a minimum near x = −0.11021. The inflection points are near x = −0.73227, −0.28090, 0.34641, and 1.70711. The maximum near x = 1 lies almost exactly on every specified grid, so the other stationary points provide the more demanding location tests.

  8. Companion implementation versus a library routine

    Choose one Chapter 6 example and compare a companion implementation with a numerical-differentiation routine in the language you normally use. Compare accuracy, step-size control, endpoint behavior, reported error estimates when available, and ease of use. Consult the Chapter 6 resources page for specific routines. If you use nonuniform data, also check the documented accuracy of the routine on a nonuniform grid; it may differ from its uniform-grid accuracy.

    Discuss when a transparent finite-difference implementation is preferable and when a tested library routine is the better scientific choice.

Physics projects

These projects use numerical differentiation as part of a larger physical investigation. Each project includes an analytical or limiting result that can be used for validation.

1. Electric field and charge density from a potential map

Objective: reconstruct an electric field and charge density from potential values sampled on a two-dimensional grid.

V(x,y) = V0 sin(πx/L) sin(πy/L),    0 ≤ x,y ≤ L.

Calculate E = −∇V and use Poisson's equation, ρ = −ε0∇2V, to infer the corresponding charge density. Plot the potential, electric-field vectors, and charge-density map. Examine how the errors change with grid resolution and near the boundaries. Use L = V0 = ε0 = 1 and grids with 33, 65, and 129 points in each direction as an initial dimensionless test. For the main error comparison, compute the charge density on interior nodes. As a separate boundary test, use second-order one-sided formulas to evaluate it on the edges.

Validation: derive the exact field components and verify that ρ = 2π2ε0V0 sin(πx/L)sin(πy/L)/L2. Because the exact charge density is zero on the boundary, use absolute error there rather than relative error.

Possible extension: add small random errors to the potential values and compare how strongly the field and charge density are affected.

2. Group velocity and effective mass in an energy band

Objective: extract transport quantities from a numerically tabulated dispersion relation.

E(k) = E0 − 2J cos(ka),    vg = (1/ℏ) dE/dk,    m* = ℏ2/(d2E/dk2).

Use dimensionless values ℏ = a = J = 1 and E0 = 0. Tabulate E(k) at 64 equally spaced points on the endpoint-exclusive grid −π/a ≤ k < π/a across the first Brillouin zone and calculate the group velocity and effective mass numerically. Locate the extrema of the band, the maximum magnitude of the group velocity, and the points where the effective mass changes sign. Because m* diverges when the band curvature passes through zero, plot 1/m* as well as, or instead of, plotting m* directly. Compare second- and fourth-order formulas and use periodic wraparound at the Brillouin-zone boundaries. Do not include both −π/a and π/a, because they represent the same physical state and duplicate an endpoint in a wrapped stencil. This is the same convention used in the periodic cos(3x) exercise above.

Validation: differentiate the cosine band analytically and verify the numerical velocity and curvature. Show that |vg| is maximal at ka = ±π/2, precisely where the curvature is zero and the effective mass changes sign through infinity.

Possible extension: perturb the tabulated energies to imitate calculated or measured band-structure uncertainty and test a local polynomial fit before differentiating.

3. Galaxy mass density from a rotation curve

Objective: infer a spherical mass distribution from a circular-speed profile and study the sensitivity to noisy data.

v(r) = v0r/√(r2 + a2),    M(r) = rv2(r)/G,    ρ(r) = [4πr2]−1dM/dr.

Generate a synthetic rotation curve on both uniform and irregular radial grids. Calculate the enclosed mass and then differentiate it to obtain the density. Compare direct finite differences, local polynomial differentiation, and a smoothed calculation after adding 1% and 5% normally distributed relative errors to the speeds. For an initial dimensionless calculation, use G = v0 = a = 1 and 201 points over 0 ≤ r ≤ 10a.

Validation: verify that ρ(r) = v02(r2+3a2)/ [4πG(r2+a2)2]. It approaches 3v02/(4πGa2) near the origin and decreases as r−2 at large radius. Treat the origin by its limiting value because direct division by 4πr2 is numerically fragile. Also compare direct differentiation of M(r) with dM/dr = [v2+2rv(dv/dr)]/G. For the noisy cases, plot relative density error versus radius and explain why noise is especially strongly amplified near the origin.

Possible extension: use a rotation curve from the public SPARC database and discuss the uncertainty of the inferred density. Since SPARC galaxies are disks whereas the formula above assumes spherical symmetry, interpret the result as a spherical-equivalent density, not the actual three-dimensional density of the disk. See F. Lelli, S. S. McGaugh, and J. M. Schombert, “SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves,” The Astronomical Journal 152 (2016), 157, doi:10.3847/0004-6256/152/6/157.

4. Divergence and vorticity of a two-dimensional flow

Objective: identify local expansion, compression, and rotation from velocity-field data.

∇·u = ∂u/∂x + ∂v/∂y,    ωz = ∂v/∂x − ∂u/∂y.

Use the regularized radial source

ur(r) = [Q/(2πr)] [1 − e−r2/σ2],    ∇·usource = [Q/(πσ2)] e−r2/σ2.

Construct a velocity field by combining solid-body rotation, u = −Ωy, v = Ωx, with the smooth regularized radial source above. Convert its radial velocity to Cartesian components and add them to the rotational field. When constructing the source on the grid, set its Cartesian components to their limiting values u = v = 0 at r = 0 to avoid division by zero. Use Ω = Q = 1, σ = 0.30, and 51- and 101-point grids on −1 ≤ x,y ≤ 1 as initial values. Calculate the divergence and vorticity numerically, and use contour maps to distinguish the rotating and expanding regions.

Validation: verify that the combined field has the Gaussian divergence shown above and uniform vorticity ωz = 2Ω. Examine the convergence of both quantities and the treatment of the removable r = 0 limit.

Possible extension: add a Lamb–Oseen vortex with uθ = [Γ/(2πr)] [1−e−r2/σv2] and exact vorticity ωz = [Γ/(πσv2)] e−r2/σv2, or use a measured particle-image-velocimetry field and investigate the influence of spatial resolution and noise.

← Back to Chapter 6