Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 5

Interpolation

These exercises supplement those in the book. They range from short practice problems to broader comparisons of interpolation methods and can be used for extra practice, self-study, or homework assignments. The physics projects use interpolation as it commonly appears in computational work: reconstructing physical quantities from values known only at discrete points.

Chapter 5 codes Chapter 5 resources Chapter 5 updates & errata Chapter 5 overview

Additional exercises

  1. Interpolation order and local stencil size

    Sample

    f(x) = e−x cos(4x),    0 ≤ x ≤ 3

    at 11 equally spaced points. Using the companion polynomial interpolation routine, compare polynomial degrees 1 through 5. For degree d, use the d + 1 data points nearest the evaluation point. Evaluate the interpolant on a fine grid of at least 500 points and record the maximum absolute error for each degree.

    Repeat with 21 base points. Does increasing the local polynomial degree always improve the result? Use direct evaluation of f(x) as the reference solution.

  2. Unequally spaced data

    Begin with 15 equally spaced points for

    f(x) = 1/(1 + x2),    0 ≤ x ≤ 4.

    Let Δx denote the original spacing. Move each interior point by a random amount chosen uniformly between −0.15Δx and 0.15Δx, leaving the two endpoints fixed. These perturbations are too small to reverse the ordering of neighboring points. Compare linear, cubic local polynomial, and cubic spline interpolation on the regular and irregular grids.

    Repeat for several random grids and compare the maximum interpolation errors with the exact function. Does irregular spacing affect all methods in the same way?

  3. Where is the interpolation error largest?

    Tabulate

    f(x) = tanh(3x),    −1 ≤ x ≤ 1

    at 13 equally spaced points. Compare local linear, quadratic, cubic, and cubic spline interpolation. Plot the pointwise error rather than only a single average error, and identify where the maximum error occurs.

    Relate the error distribution to the regions in which the function changes most rapidly. Use direct evaluation of tanh(3x) as the reference.

  4. Chebyshev nodes beyond the standard Runge example

    Compare equally spaced and Chebyshev nodes for

    f(x) = |x|,    −1 ≤ x ≤ 1.

    Construct global interpolating polynomials using 9, 17, and 33 nodes. Plot the interpolants and the pointwise errors on a fine grid. Use the Chebyshev extrema xj = cos[jπ/(n − 1)], j = 0, …, n − 1. When available, use barycentric evaluation to avoid instability in the polynomial evaluation itself. Chebyshev spacing reduces the endpoint oscillations associated with equally spaced nodes, but f(x) is not differentiable at x = 0.

    Determine which difficulty remains even after the node distribution has been improved, and compare how the maximum error changes as the number of nodes increases. Estimate the convergence rate from a log–log plot of maximum error versus the number of nodes.

  5. Spline boundary conditions

    Interpolate

    f(x) = ex/2,    0 ≤ x ≤ 4

    using 9 equally spaced data points. Compare natural, not-a-knot, clamped, and the Forsythe–Malcolm–Moler (FMM)-type endpoint conditions used by the companion spline program. For the clamped spline use the exact endpoint derivatives f′(x) = (1/2)ex/2.

    Plot the interpolation error near both ends of the interval and in the interior. Which boundary condition is most accurate for this particular function, and why?

  6. Shape preservation and spline overshoot

    Consider the monotonic tabulated data

    x = [0, 1, 2, 3, 4, 5],    y = [0, 0.02, 0.10, 0.85, 0.97, 0.99].

    Compare a standard cubic spline with a shape-preserving method such as PCHIP, Steffen interpolation, or another monotonic piecewise cubic interpolant available in your language.

    Check whether either interpolant introduces artificial extrema or values outside the data range 0 ≤ y ≤ 0.99. Plot both curves together with the original data points.

  7. Interpolation versus noisy data

    Generate 21 exact values from

    f(x) = e−x/2,    0 ≤ x ≤ 6,

    and add independent, normally distributed zero-mean random perturbations with standard deviations 0.01 and 0.05 to create two synthetic measurement sets. Use a fixed random seed so that your calculation can be reproduced.

    Interpolate each noisy data set exactly with a cubic spline and compare the result with the underlying exact function. The purpose is not to introduce curve fitting, which is treated later in the book, but to demonstrate why exact interpolation is inappropriate when the data themselves are uncertain.

    Note that f(6) = e−3 ≈ 0.05, so the larger noise standard deviation is about as large as the signal at the end of the interval. As an optional extension, differentiate the interpolating spline and compare its derivative with the exact derivative. This provides a direct link to the noise-amplification problem discussed in Chapter 6.

  8. Interpolation is not extrapolation

    Tabulate

    f(x) = ln(1 + x),    0 ≤ x ≤ 4

    at 9 equally spaced points and construct several interpolants. First test them only between known data points. Then evaluate them at x = 4.5, 5, and 6. Compare the extrapolated values with the exact values ln(5.5), ln(6), and ln(7).

    How quickly does the error grow outside the sampled interval? Does a method that interpolates well inside the data range necessarily extrapolate well?

  9. Companion implementation versus a library routine

    Choose either polynomial or cubic spline interpolation and compare the Chapter 5 companion implementation with a standard routine in your preferred language: MATLAB interp1 or spline, SciPy interpolation routines, GSL interpolation, or one of the specific Fortran options listed on the Chapter 5 resources page.

    Use one smooth test problem and one more demanding problem from the exercises above. Compare accuracy, behavior near the boundaries, extrapolation policy, ease of use, and any available endpoint options. Explain when a transparent implementation helps you understand the algorithm and when a library routine is the better practical choice.

Classic physics applications

These projects use familiar physics examples in which a physical quantity is known analytically but is treated as if only a sparse table were available. That makes it possible to evaluate the interpolation error directly.

1. Optical dispersion from tabulated refractive-index data

Objective: reconstruct the wavelength dependence of a refractive index from a sparse optical data table.

For fused silica, generate reference values from the Sellmeier relation

n2(λ) = 1 + B1λ2/(λ2 − C12) + B2λ2/(λ2 − C22) + B3λ2/(λ2 − C32),

where λ is measured in micrometers and B1 = 0.6961663, B2 = 0.4079426, B3 = 0.8974794, C1 = 0.0684043, C2 = 0.1162414, and C3 = 9.896161. Use the seven wavelengths 0.4, 0.6, …, 1.6 µm as the available data, then interpolate n(λ) at intermediate wavelengths using local polynomial interpolation and cubic splines. These coefficients are from I. H. Malitson, “Interspecimen Comparison of the Refractive Index of Fused Silica,” Journal of the Optical Society of America 55 (1965), 1205–1208, doi:10.1364/JOSA.55.001205.

Validation: compare the interpolated values with the Sellmeier relation evaluated directly on a fine wavelength grid.

Possible extension: estimate the group index ng = n − λ dn/dλ from the interpolated curve and examine how interpolation accuracy affects this derived quantity.

2. Resolving a resonance in a scattering cross section

Objective: investigate how data spacing affects reconstruction of a narrow physical resonance.

Use a smooth background plus a Breit–Wigner resonance,

σ(E) = σ0 + A(Γ/2)2 / [(E − Er)2 + (Γ/2)2].

Here E is energy, Er is the resonance energy, Γ is the full width at half maximum, σ0 is the background cross section, and A is the resonance amplitude. In arbitrary units, begin with σ0 = 1, A = 4, Er = 5.23, and Γ = 0.30. Generate one table with spacing ΔE = 0.5 and another with ΔE = 0.1 over 3 ≤ E ≤ 7.

Compare linear interpolation, cubic local interpolation, and cubic splines, paying particular attention to the reconstructed resonance position, height, and width.

Validation: compare each interpolated curve with the analytic cross section used to generate the table. The exact maximum occurs at E = Er and has value σ0 + A. Because the analytic cross section never falls below σ0, also check whether an interpolant dips below the background or becomes negative.

Possible extension: use nonuniform sampling with more points near the resonance and determine how much the error is reduced for the same total number of tabulated points.

3. Reconstructing a quantum-mechanical wavefunction

Objective: examine interpolation of a quantity that changes sign and contains a physically important node.

Use the normalized first excited harmonic-oscillator wavefunction in the dimensionless coordinate ξ,

ψ1(ξ) = (√2/π1/4) ξ e−ξ2/2.

Sample the wavefunction on −3.9 ≤ ξ ≤ 4.1 with spacing Δξ = 0.5 and reconstruct ψ1(ξ) between grid points. Compare local polynomial and cubic spline interpolation for both the wavefunction and the probability density |ψ1|2.

Validation: compare with the analytic wavefunction, locate the interpolated zero crossing numerically, and compare it with ξ = 0. On the common interval −3.9 ≤ ξ ≤ 3.9, check the expected odd symmetry ψ1(−ξ) = −ψ1(ξ).

Possible extension: repeat for a higher excited state, where additional nodes make the spatial sampling requirement more demanding.

Additional physics projects

These projects place interpolation inside a larger physical calculation or use data whose structure suggests a particular interpolation strategy.

4. Gamma-ray detector efficiency calibration

Objective: compare ordinary and transformed interpolation for a positive quantity that changes substantially over its energy range.

Use the model detector efficiency

ε(E) = 0.30(E/100 keV)−0.8 exp[−E/(3000 keV)],

where E is the photon energy and ε is the full-energy detection efficiency. This simplified model represents the decreasing part of an efficiency curve above 100 keV. Generate a calibration table at E = 100, 150, 250, 400, 700, 1000, and 1500 keV, and then treat only the tabulated values as known.

Interpolate the efficiency at intermediate energies using linear interpolation, a cubic spline, and linear interpolation of log ε versus log E. Compare the methods over the full energy range.

Validation: evaluate the model function directly at the test energies and compare with the interpolated values. The efficiency should remain positive and decrease monotonically with energy.

Possible extension: perturb the tabulated efficiencies by small measurement uncertainties and examine whether a higher-order interpolant introduces artificial structure.

5. Tabulated drag coefficient and aerodynamic drag

Objective: examine how interpolation of a tabulated aerodynamic coefficient affects the calculated drag force.

Use the model drag coefficient

CD(v) = 0.25 + 0.20 exp{−[(v − 340 m/s)/(70 m/s)]2},

where v is the speed. Sample CD(v) every 50 m/s between 0 and 600 m/s and then treat only the tabulated values as known. Compare linear interpolation and a cubic spline for CD(v) over the full speed range.

Use the interpolated drag coefficient to calculate the drag-force magnitude

FD(v) = (1/2)ρ CD(v) A v2,

where ρ is the air density and A is the cross-sectional area. Use ρ = 1.225 kg/m3 and A = 0.010 m2. Plot both CD(v) and FD(v) using the two interpolation methods.

Validation: evaluate the analytic CD(v) model directly on a fine speed grid and use it to compute a reference drag force. Compare the interpolated and reference values, and identify where the largest absolute and relative errors occur. Because FD(0) = 0, use absolute error near zero and report relative error only for v > 0 (for example, from 1 m/s upward).

Possible extension: repeat the calculation with table spacings of 25, 75, and 100 m/s and determine how the spacing affects the error. After studying numerical solutions of differential equations in Chapter 9, the same tabulated CD(v) can be used in a projectile trajectory calculation.

6. Two-dimensional equation-of-state table

Objective: extend one-dimensional interpolation to a physical quantity depending on two independent variables.

Use the reduced van der Waals equation

Pr = 8Tr/(3vr − 1) − 3/vr2,

where Pr = P/Pc, Tr = T/Tc, and vr = v/vc are pressure, temperature, and molar volume scaled by their critical values. Generate a 6 × 10 rectangular table using six uniformly spaced Tr values over 1.05 ≤ Tr ≤ 2 and ten uniformly spaced vr values over 0.6 ≤ vr ≤ 4. Treat only the tabulated pressures as known.

Reconstruct Pr at intermediate (Tr,vr) points using bilinear interpolation and, if available, a bicubic method.

Suitable library tools include MATLAB interp2 and SciPy RegularGridInterpolator.

Validation: compare interpolated pressures with direct evaluation of the reduced van der Waals equation at randomly selected points inside the table. Report both typical and maximum relative errors.

Possible extension: compare a uniform table with a nonuniform grid that places more points where the pressure changes most rapidly.

← Back to Chapter 5