Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Chapter 7

Numerical Integration

These exercises supplement those in the book. They range from short practice problems to more challenging comparisons and can be used for extra practice, self-study, or homework assignments. The physics projects use numerical integration as part of a larger physical investigation.

Chapter 7 codes Chapter 7 resources Chapter 7 updates & errata Chapter 7 overview

Additional exercises

  1. Convergence of Newton–Cotes rules

    Use the composite trapezoidal and Simpson rules for the three smooth test integrals below. Start with a small number of subintervals and repeatedly double it. Plot the absolute error versus the number of subintervals and see whether the expected convergence rate appears in your results.

    ∫0π sin(x) dx = 2
    ∫01 ex dx = e − 1
    ∫01 dx/(1 + x2) = π/4

    For a more challenging case, try ∫01 √x dx = 2/3. The integrand is continuous but has an infinite derivative at x = 0. Estimate the observed convergence rate from the slope of the error curve. How does it differ from the nominal rates for smooth integrands?

  2. Simpson versus Gauss–Legendre quadrature

    Compare composite Simpson and Gauss–Legendre quadrature on the integrals below. Give the two methods roughly comparable numbers of function evaluations and compare the resulting errors.

    ∫01 ex dx = e − 1
    ∫01 dx/(1 + x2) = π/4
    ∫01 e−x2 dx = (√π/2) erf(1)

    Which method gains accuracy more quickly for these smooth functions? Is the answer the same for all three integrals? As a hint, the complex singularities of 1/(1+x²) at x = ±i are relatively close to the real integration interval and can slow high-order convergence.

  3. Two peaks on very different scales

    Consider a function with one broad peak and one much narrower feature,

    f(x) = e−40(x − 0.30)2 + 0.05e−4000(x − 0.72)2,    0 ≤ x ≤ 1.

    Compare composite Simpson on a uniform grid with adaptive Simpson or adaptive Gauss–Kronrod integration. Use several grid sizes or tolerances and plot the integrand together with the points or subintervals used by each method.

    Can a uniform-grid calculation appear nearly converged while still missing part of the narrow peak? What changes if the narrow peak is made even narrower? For the parameters above, the narrow peak contributes only about 0.5% of the total integral, so a calculation can look convincing while still missing it.

  4. Requested tolerance versus actual error

    Run an adaptive integration routine with requested tolerances from about 10−3 down to 10−12. Use the two integrals

    ∫01 ex dx = e − 1
    ∫01 e−20x dx = (1 − e−20)/20.

    For each tolerance, record the numerical result, the routine's error estimate, the actual error, and the number of function evaluations if it is available. Does asking for a tolerance of 10−10 actually produce an error close to 10−10?

  5. An endpoint singularity that can be removed

    Evaluate

    ∫01 ln(1 + x)/√x dx = π + 2 ln(2) − 4.

    First apply an adaptive routine directly. Then use the substitution x = t2 and integrate the transformed function on 0 ≤ t ≤ 1. Compare accuracy and computational effort.

    Plot both integrands. Why is the transformed problem easier numerically even though the value of the integral is unchanged? The square-root behavior here is closely related to the turning-point behavior in the WKB project below.

  6. Two transformations for an infinite interval

    Evaluate

    ∫0∞ e−x sin(x) dx = 1/2

    by mapping the semi-infinite interval onto a finite interval in two different ways:

    x = t/(1 − t),    0 ≤ t < 1
    x = tan(πt/2),    0 ≤ t < 1.

    Plot the transformed integrands and compare convergence, function-evaluation counts, and sensitivity near t = 1. Which transformation works better for this particular integral? Explain the result by examining how rapidly each transformed integrand approaches zero or oscillates as t → 1.

  7. Cauchy principal value

    Consider

    PV ∫−11 cos(x)/(x − 0.2) dx.

    First remove a symmetric interval (0.2 − δ, 0.2 + δ) around the pole and study the result as δ becomes smaller. Then evaluate the same integral using analytical subtraction or a library routine designed for principal values.

    The reference value is approximately −0.5912784964342435, matching the principal-value benchmark on the Chapter 7 resources page.

    For additional practice, repeat the calculation for

    PV ∫−11 ex/(x − 0.3) dx.

    Its reference value is approximately 1.620314024361904.

    Why is it unsafe simply to give an ordinary quadrature routine an interval containing the pole?

  8. What happens as the oscillation frequency grows?

    Study

    I(ω) = ∫01 e−x cos(ωx) dx
    I(ω) = [1 + e−1 (ω sinω − cosω)]/(1 + ω2).

    for ω = 10, 25, 50, 100, and larger values. Compare composite Simpson with an oscillatory-specific method such as Filon quadrature or a weighted library routine.

    Track both the error and the number of function evaluations. How quickly does a general-purpose method become expensive as the oscillations get faster? The case ω = 50 is also used as a benchmark on the Chapter 7 resources page.

  9. Cancellation inside the integrand

    Evaluate

    I(δ) = ∫01 [e−x − e−(1+δ)x] dx,    δ = 10−8.

    A high-precision reference value is

    I(δ) ≈ 2.642411168541014 × 10−9.

    Compare direct evaluation of the integrand with the stable form −e−xexpm1(−δx). Also derive the leading-order approximation I(δ) ≈ δ(1 − 2/e). Use increased-precision arithmetic to evaluate the reference, since subtracting the two separately integrated exponentials in ordinary double precision suffers from the same cancellation as the direct integrand.

    This exercise is less about the quadrature formula itself and more about an important lesson: an integration routine cannot recover accuracy that has already been lost while evaluating the integrand.

  10. Two-dimensional product quadrature

    Apply an n × n Gauss–Legendre product rule to

    ∫01 ∫01 dy dx /(1 + x2 + y2).

    Use several values of n and plot the absolute error against the total number n2 of function evaluations. Compare your result with a two-dimensional library routine if one is available. A high-accuracy reference value is 0.6395103518703110.

    Notice how quickly the number of points grows when a one-dimensional rule is extended by taking products.

  11. How accurate should the inner integral be?

    Use the variable-boundary integral

    ∫01 ∫0sin(x) x2/(y2 + 2) dy dx.

    Hold the outer tolerance fixed while changing the tolerance used for the inner integration. Find the point at which making the inner calculation more accurate no longer improves the final answer, and then make the inner tolerance looser until it begins to affect the result noticeably.

    The goal is to see that nested integration has two accuracy controls that must work together rather than independently.

  12. Companion implementation versus a library routine

    Choose one of the Chapter 7 companion programs and compare it with a widely used routine in your preferred language: MATLAB integral or quadgk, Python/SciPy quad, C++/GSL, or Fortran/QUADPACK.

    Test both routines on at least one easy integral and one more demanding case from the exercises above. Compare accuracy, error estimates, function-evaluation counts when available, and ease of use.

    Which implementation helps you understand the algorithm better? Which would you choose for a larger scientific calculation, and why?

Classic physics applications

These projects use familiar physics problems to practice numerical integration in settings where the physical interpretation and limiting behavior are well understood. They are useful for building confidence with the numerical methods before moving to less standard applications.

1. Debye heat capacity of a solid

Objective: use numerical integration to calculate how the heat capacity of a solid changes with temperature in the Debye model.

CV = 9NkB(T/ΘD)3 ∫0ΘD/T x4ex/ (ex − 1)2 dx.

Here N is the number of atoms, kB is Boltzmann's constant, ΘD is the Debye temperature, and x is the dimensionless integration variable. Calculate CV over a broad range of temperature and plot CV/(NkB) as a function of T/ΘD.

At x = 0, define the integrand by its limiting value, zero. For small x, evaluate ex − 1 with expm1(x). At very low temperatures the upper limit can be large; use the equivalent form x4e−x/(1 − e−x)2 to avoid overflow.

Validation: at low temperature, verify CV ≈ (12π4/5)NkB(T/ΘD)3. At high temperature, verify the Dulong–Petit limit CV → 3NkB.

Possible extension: choose the Debye temperature of a real material and compare the model with experimental heat-capacity data.

2. Quantum tunneling in the WKB approximation

Objective: study how the tunneling probability depends on particle energy and on the shape of the potential barrier.

P(E) ≈ exp{−2 ∫x1x2 √[2m(V(x) − E)]/ℏ dx}.

Here m is the particle mass, E is its energy, V(x) is the barrier potential, x1 and x2 are the classical turning points, and ℏ is the reduced Planck constant. Compare a rectangular barrier of height V0 and width a with a smooth barrier such as V(x) = V0 exp(−x2/a2). A convenient first calculation is to work in dimensionless units with m = ℏ = V0 = 1, choose a wider barrier such as a = 5, and vary 0 < E < V0. The larger width gives an exponent well above unity over much of the energy range, where the WKB approximation is more reliable. For the Gaussian barrier, the turning points are known exactly: x = ±a√ln(V0/E).

Validation: for a rectangular barrier of width a, recover the WKB exponent 2a√[2m(V0 − E)]/ℏ. Compare the WKB result with the exact rectangular-barrier transmission coefficient T = {1 + V02sinh2(κa)/ [4E(V0−E)]}−1, where κ = √[2m(V0−E)]/ℏ. Relate the square-root behavior at smooth turning points to the endpoint-singularity exercise above.

Possible extension: construct two barriers with the same height and approximately the same width but different shapes, and compare their tunneling probabilities.

3. Electrostatic potential of a finite charged disk

Objective: move from a one-dimensional analytic benchmark to a genuinely two-dimensional electrostatic integral.

V(ρ,z) = (σ/4πε0) ∫0R ∫02π r′ dφ dr′ / √[ρ2 + r′2 − 2ρr′cosφ + z2].

Here σ is the surface charge density, ε0 is the vacuum permittivity, R is the disk radius, (ρ,z) specifies the observation point in cylindrical coordinates, and (r′,φ) specifies a source point on the disk. First calculate the potential on the symmetry axis (ρ = 0). Then move off axis and evaluate the full two-dimensional integral. A useful set of dimensionless observation points is ρ/R = 0, 0.5, 1, 1.5 and z/R = 0.25, 0.5, 1. Keep z > 0 so that the observation point does not lie on the charged surface.

Validation: compare the on-axis result with V(0,z) = (σ/2ε0) [√(R2 + z2) − z].

Possible extension: calculate equipotential curves in a plane containing the disk axis, or obtain an electric-field component by numerical differentiation of the potential.

4. Fresnel diffraction by a single slit

Objective: evaluate an oscillatory wave integral and follow the transition from near-field to far-field diffraction.

U(x,z) ∝ ∫−a/2a/2 exp[i k(x − x′)2/(2z)] dx′,    I(x,z) ∝ |U(x,z)|2.

Here a is the slit width, λ is the wavelength, k = 2π/λ is the wave number, x′ is a position across the slit, x is the observation coordinate, z is the propagation distance, U is the complex wave amplitude, and I is the intensity. Evaluate the real and imaginary parts of the amplitude separately and calculate the intensity profile for several values of z. As a guide, choose distances corresponding roughly to F = a2/(λz) = 10, 1, and 0.1 to sample the near-field, transition, and far-field regimes.

Validation: at sufficiently large z, compare the numerical pattern with the Fraunhofer single-slit result I(θ) ∝ [sin(β)/β]2, where θ is the diffraction angle and β = (πa/λ) sinθ. For every propagation distance, also compare the numerical complex amplitude with the exact expression in terms of the Fresnel integrals C and S, available in standard scientific libraries. The case F = 0.1 is only moderately far from the slit, so small visible differences from the Fraunhofer pattern are expected.

Possible extension: vary the slit width or wavelength and identify how the distance needed to reach the far-field regime changes.

Additional physics projects

These projects use the same numerical ideas in less routine settings. They are somewhat more open-ended and are intended for students who want a more challenging problem or for instructors looking for a project that goes beyond the most common textbook examples.

5. Dispersion from an absorption resonance

Objective: use a principal-value integral to connect an absorption line with the associated dispersive response.

Let ω be angular frequency, ω0 the resonance frequency, and γ the linewidth parameter. Begin with the Lorentzian absorption line

A(ω) = γ/[(ω − ω0)2 + γ2].

Form the corresponding dispersive function D(ω) through the Hilbert-transform form of the Kramers–Kronig relation,

D(ω) = (1/π) PV ∫−∞∞ A(ω′)/(ω − ω′) dω′.

A convenient first calculation uses dimensionless units with ω0 = 0 and γ = 1. Evaluate D(ω) over a range spanning the resonance. Compare direct exclusion of a small symmetric interval around the pole with analytical subtraction or a principal-value library routine. This project combines the two numerical difficulties treated separately in the principal-value and infinite-interval exercises above.

Validation: for the Lorentzian above and the sign convention used here, the exact result is D(ω) = (ω − ω0)/ [(ω − ω0)2 + γ2]. The result is antisymmetric about the resonance and changes sign at ω = ω0.

Possible extension: replace the single resonance by two nearby absorption lines and examine how the dispersive response changes.

6. The Gamow peak in thermonuclear reactions

Objective: see how two competing exponential factors create a narrow energy window that dominates a thermal nuclear reaction rate.

Let E be the relative kinetic energy of the reacting nuclei, T the temperature, kB Boltzmann's constant, and EG the Gamow-energy parameter that characterizes Coulomb-barrier penetration. In a simplified model with a slowly varying astrophysical S-factor, the important energy dependence can be written as

g(E) = exp[−E/(kBT) − √(EG/E)].

Plot g(E) and numerically evaluate its integral from E = 0 to ∞ for several temperatures. For a dimensionless first calculation, set EG = 1 and use kBT/EG = 0.01, 0.03, and 0.1. Compare a uniform-grid calculation with adaptive integration and examine where most of the contribution comes from. Define g(0) = 0 by its limiting value rather than evaluating √(EG/E) at the endpoint.

Validation: compare the numerical location of the maximum with E0 = [EG(kBT)2/4]1/3. Also compare the numerical integral with the saddle-point approximation (√π/2)Δ exp[−3E0/(kBT)], where Δ = 4√[E0kBT/3]. Examine how the approximation improves as kBT/EG decreases.

Possible extension: allow the S-factor to vary slowly with energy and determine when that variation noticeably shifts the effective reaction window.

7. Voigt spectral line profile

Objective: calculate a spectral line shape produced by the convolution of Doppler and collisional broadening.

Let x measure displacement from line center, let σ be the standard deviation of the Gaussian Doppler profile, and let γ be the half-width at half maximum of the Lorentzian collisional profile. Use the normalized forms

G(x) = exp[−x2/(2σ2)]/ (σ√(2π)),    L(x) = γ/[π(x2 + γ2)].

The Voigt profile is their convolution,

V(x) = ∫−∞∞ G(x′)L(x − x′) dx′.

For a simple starting calculation, set σ = 1 and use γ/σ = 0.1, 1, and 10. Calculate V(x) and compare how the line center and the wings change as one broadening mechanism becomes dominant. Remember that γ is a half-width, whereas σ is a standard deviation; γ/σ = 1 does not mean that the two component profiles have equal half-widths.

Validation: verify that the profile approaches a Gaussian when γ/σ is very small and a Lorentzian when γ/σ is very large. Check numerically that the area under the normalized profile remains close to one. Because Lorentzian wings decrease only as 1/x2, the computed area can depend strongly on the chosen finite integration range. For example, when γ/σ = 10, restricting the area calculation to |x| ≤ 50σ misses more than 10% of the total. Explain any observed shortfall rather than interpreting it immediately as a quadrature failure.

Possible extension: compare your direct numerical convolution with a Voigt-function routine available in a scientific library.

8. Exoplanet transit with stellar limb darkening

Objective: calculate a transit light curve by integrating the brightness over the visible part of a stellar disk.

Let the stellar radius be R, let r be the projected distance from the center of the stellar disk, let I0 be the intensity at the disk center, and let u be the linear limb-darkening coefficient. Define μ as the cosine of the angle between the local surface normal and the line of sight. For a spherical star,

I(r) = I0[1 − u(1 − μ)],    μ = √(1 − r2/R2).

In this simple model, 0 ≤ u ≤ 1: u = 0 gives a uniformly bright disk, while larger values produce stronger limb darkening. A useful initial choice is u = 0.6 and Rp/R = 0.1, where Rp is the planet radius. Place the dark circular planet at different positions across the stellar disk. Numerically integrate I(r) over the part of the star that remains visible and plot the normalized stellar flux as a function of the planet's position.

Validation: first set u = 0. For a planet completely inside the stellar disk, the fractional loss of light should equal the fractional area blocked, (Rp/R)2. With no planet and arbitrary linear limb darkening, verify the normalization F0 = πR2I0(1 − u/3). The stellar-limb integrand has square-root behavior, and a planet crossing the limb creates variable integration boundaries.

Possible extension: compare central and off-center transits, or compare several limb-darkening values such as u = 0, 0.3, 0.6, and 0.9. Study how the planet radius and the limb-darkening coefficient change the depth and shape of the light curve. Compare with the analytic transit models of K. Mandel and E. Agol, “Analytic Light Curves for Planetary Transit Searches,” The Astrophysical Journal 580 (2002), L171–L175, doi:10.1086/345520.

← Back to Chapter 7