Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Exercises & Projects

The material on this page extends the exercises and projects in the book. The additional exercises emphasize numerical experimentation, testing, and comparison of methods. The projects use numerical integration to investigate physical problems.

Most exercises and projects can be completed in any programming language. Companion codes may be used as starting points when appropriate, but they should be modified, tested, and validated rather than treated as black boxes.

Chapter 7

Numerical Integration

The exercises below focus on the numerical behavior of integration methods. The physics projects use the same ideas in mechanics, thermal physics, quantum mechanics, electromagnetism, and optics.

Chapter 7 codes Chapter 7 resources Chapter 7 overview

Additional exercises

  1. Convergence of Newton–Cotes rules

    Use the composite trapezoidal and Simpson rules for several smooth functions. Double the number of subintervals successively and determine the observed convergence rate from the numerical error.

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

    As an extension, try ∫01 √x dx = 2/3 and examine how reduced smoothness changes the convergence.

  2. Simpson versus Gauss–Legendre quadrature

    Compare composite Simpson and Gauss–Legendre quadrature using approximately comparable numbers of function evaluations. Suggested test integrals are

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

    Examine how rapidly the error decreases as the computational effort increases.

  3. Move and sharpen a narrow peak

    Starting from the peaked integrand

    I = ∫01 dx/[1 + A(x − x0)2],

    vary both the peak position x0 and the parameter A. Compare a uniform-grid method with adaptive Simpson or adaptive Gauss–Kronrod integration. The companion example corresponds to A = 400 and x0 = 0.2.

  4. Requested tolerance versus actual error

    Run an adaptive method with requested tolerances ranging from approximately 10−3 to 10−12. Use one smooth and one sharply peaked test integral:

    ∫01 ex dx = e − 1
    ∫01 dx/[1 + 400(x − 0.2)2] = [atan(16) + atan(4)]/20

    Compare the requested tolerance, estimated error, actual error, and number of function evaluations. As an extension, record the accepted subintervals and show where the adaptive method concentrates its work.

  5. Transformations for an improper integral

    Evaluate

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

    using two mappings of the semi-infinite interval to a finite interval:

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

    Plot the transformed integrands and compare convergence and numerical effort.

  6. Cauchy principal value

    Consider

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

    First exclude a symmetric interval around the pole and study convergence as its width approaches zero. Then compare the result with analytical subtraction. A second test problem is

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

  7. Increasing oscillation frequency

    Study

    I(ω) = ∫01 e−x cos(ωx) dx

    for ω = 10, 25, 50, 100, and larger values. Compare composite Simpson and Filon integration and determine how the computational effort changes with oscillation frequency.

  8. Two-dimensional product quadrature

    Apply an n × n Gauss–Legendre product rule to

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

    Plot the absolute error versus the total number n2 of function evaluations.

  9. Tolerance allocation in nested integration

    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. Determine when the inner calculation becomes unnecessarily accurate and when insufficient inner accuracy begins to affect the final result.

  10. Companion implementation versus a library routine

    Work in the programming language you normally use. Choose one Chapter 7 example and compare the companion implementation with a corresponding general-purpose routine: MATLAB integral or quadgk, Python/SciPy quad, C++/GSL, or Fortran/QUADPACK.

    Compare numerical accuracy, the reported error estimate, function-evaluation count when available, and ease of use. Discuss when a transparent implementation is useful and when a well-tested library routine is the better choice.

Physics projects

These projects use numerical integration as a tool for exploring physical behavior. They are intended to involve both computation and interpretation rather than only evaluation of an integral.

1. Nonlinear pendulum

Objective: determine how the period of a pendulum changes when the small-angle approximation is no longer valid.

T(θ0) = 4√(L/g) ∫0π/2 dφ / √[1 − sin2(θ0/2) sin2φ].

Calculate T/T0 as a function of the initial angle, where T0 = 2π√(L/g). Compare with the small-angle approximation and investigate the behavior as θ0 approaches π.

Validation: recover T/T0 → 1 for small amplitudes.

Possible extension: compare the numerical result with one or more analytical approximations for the amplitude-dependent period.

2. Blackbody radiation

Objective: investigate the integral underlying the Stefan–Boltzmann law and the distribution of radiated energy over frequency.

∫0∞ x3/(ex − 1) dx = π4/15.

Evaluate the improper integral numerically using an interval transformation or controlled truncation. Then calculate the fraction of the total radiated power lying below or above selected values of the dimensionless frequency x.

Validation: compare the total integral with π4/15.

Possible extension: convert the dimensionless calculation to a physical wavelength or frequency range for a selected temperature.

3. Debye heat capacity of a solid

Objective: use numerical integration to calculate the temperature dependence of the heat capacity in the Debye model.

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

Calculate CV over a broad temperature range and plot CV/(NkB) versus T/ΘD.

Validation: demonstrate the T3 behavior at low temperature and recover the Dulong–Petit limit CV → 3NkB at high temperature.

Possible extension: use the Debye temperature of a real solid and compare with experimental heat-capacity data.

4. Quantum tunneling in the WKB approximation

Objective: study how the transmission probability depends on particle energy and barrier shape.

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

Compare a rectangular barrier with a smooth barrier such as V(x) = V0 exp(−x2/a2). Determine the classical turning points and calculate the transmission probability as a function of E.

Validation: for the rectangular barrier, compare the numerical WKB exponent with the analytic constant-barrier result.

Possible extension: investigate how the result changes when the barrier width or height is varied.

5. Electrostatic potential of a finite charged disk

Objective: calculate a genuinely two-dimensional electrostatic integral and compare an off-axis numerical result with an on-axis analytic benchmark.

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

First calculate the potential on the symmetry axis (ρ = 0), then move the observation point off axis and evaluate the full two-dimensional integral.

Validation: for z > 0, 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.

6. Fresnel diffraction by a single slit

Objective: calculate a diffraction pattern directly from an oscillatory wave integral and follow the transition toward the far-field limit.

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

Evaluate the real and imaginary parts of the amplitude numerically and calculate the intensity for several propagation distances.

Validation: at sufficiently large z, compare with the Fraunhofer single-slit result I(θ) ∝ [sin(β)/β]2, where β = (πa/λ) sinθ.

Possible extension: vary the slit width or wavelength and determine how the diffraction pattern changes.

← Back to Chapter 7