Alexander Godunov Personal website

Companion Website

Computational Physics

From Equations to Simulation

Resources

This page provides practical material that complements the book and the companion source codes: links to maintained numerical libraries, benchmark problems for testing, method-selection guidance, and selected references.

The companion codes are intended primarily for learning, experimentation, and situations where a numerical method is useful as a transparent component of a larger program. For general-purpose scientific computing, well-tested library routines are usually the preferred choice because they provide mature error control, extensive testing, and robust implementations.

Appendix B of the book provides a broader survey of numerical routines and scientific software. The links on this page are meant to provide convenient access to current documentation.

Chapter 5

Interpolation

The Chapter 5 companion programs provide transparent implementations of local polynomial interpolation, the Runge/Chebyshev comparison, and cubic spline interpolation. The resources below connect those implementations with maintained numerical libraries and provide practical guidance for selecting an interpolation method.

Software and library resources

The table lists representative library options corresponding to the interpolation methods discussed in Chapter 5. High-level environments provide interpolation directly, while C++ and Fortran users commonly rely on established numerical libraries. Exact boundary conditions, extrapolation behavior, and calling conventions differ, so the linked documentation should be consulted for details.

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

Numerical task MATLAB Python C++ Fortran
One-dimensional local interpolation Linear and low-order interpolation for tabulated data. interp1 numpy.interp for linear interpolation;
scipy.interpolate for higher-order methods 		GSL interpolation NAG E01 interpolation routines
Global polynomial interpolation A single polynomial through all selected data points. Companion implementation for direct study;
low-degree polynomial construction can also be performed with standard polynomial tools 		BarycentricInterpolator gsl_interp_polynomial NAG E01 polynomial interpolation
Cubic spline interpolation Piecewise cubic interpolation with continuous first and second derivatives. spline
interp1(...,'spline') 		CubicSpline
make_interp_spline 		gsl_interp_cspline NAG E01 cubic spline interpolation
Shape-preserving piecewise cubic interpolation Useful when overshoot or loss of monotonicity is undesirable. pchip
makima 		PchipInterpolator
Akima1DInterpolator 		gsl_interp_steffen or gsl_interp_akima Netlib PCHIP
NAG E01 monotonic interpolation

Choosing an interpolation method

Character of the problem Method to consider Practical note
Only a few neighboring points are needed Local linear, quadratic, or cubic polynomial interpolation Low-order local interpolation is simple, inexpensive, and usually more stable than a high-degree global polynomial.
Many points describe a smooth function Cubic spline A good general-purpose choice when a smooth interpolant and smooth first derivative are important.
The data are monotonic or overshoot must be avoided PCHIP, Steffen, or Akima-type piecewise cubic interpolation Shape-preserving methods sacrifice second-derivative continuity in order to reduce artificial extrema and oscillations.
A global polynomial is required and node locations can be chosen Chebyshev nodes with a stable polynomial representation Clustering nodes near the interval endpoints greatly reduces the Runge oscillations seen with equally spaced points.
Many equally spaced points would produce a high-degree global polynomial Prefer local interpolation or a spline Adding points does not automatically improve a global polynomial; endpoint oscillations may become worse.
Smooth derivatives of the interpolant are needed Cubic spline The usual cubic spline is twice continuously differentiable; boundary conditions should be chosen with the application in mind.
Data contain measurement noise or random fluctuations Curve fitting or smoothing rather than exact interpolation Exact interpolation follows the noise. Use fitting or smoothing methods when the data are uncertain.
A value is required outside the data interval Extrapolate only with strong physical justification Extrapolation is much less reliable than interpolation and should be kept close to the known data whenever possible.

Selected references and documentation

  • G. E. Forsythe, M. A. Malcolm, and C. B. Moler, Computer Methods for Mathematical Computations, Prentice-Hall, 1977. The spline routine used as the basis for the companion implementation comes from this work.
  • C. de Boor, A Practical Guide to Splines, Springer-Verlag, 1978. A classic reference on spline construction and B-spline representations.
  • J.-P. Berrut and L. N. Trefethen, “Barycentric Lagrange Interpolation,” SIAM Review 46 (2004), 501–517.
  • L. N. Trefethen, Approximation Theory and Approximation Practice, SIAM, 2013. Includes a modern treatment of polynomial approximation and Chebyshev points.
  • F. N. Fritsch and R. E. Carlson, “Monotone Piecewise Cubic Interpolation,” SIAM Journal on Numerical Analysis 17 (1980), 238–246.
  • MATLAB interpolation documentation
  • SciPy interpolation documentation
  • GNU Scientific Library: Interpolation
  • NAG Library E01: Interpolation
  • Netlib SLATEC PCHIP routines

← Back to Chapter 5

Chapter 7

Numerical Integration

The Chapter 7 companion programs implement the principal algorithms explicitly in MATLAB, Python, C++, and Fortran. The resources below connect those implementations with general-purpose library routines and provide reference problems for verification.

Software and library resources

The table lists typical production-quality options corresponding to the numerical tasks treated in Chapter 7. The exact capabilities and calling conventions differ between libraries, so the linked documentation should be consulted for details.

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

Numerical task MATLAB Python C++ Fortran
General adaptive one-dimensional integration Automatic quadrature with user-controlled tolerances. integral
quadgk 		scipy.integrate.quad GSL qag / qags QUADPACK DQAG / DQAGS
Infinite intervals Improper integrals with one or both limits infinite. integral
quadgk 		quad with infinite bounds GSL qagi / qagiu / qagil QUADPACK DQAGI
Cauchy principal-value integrals Integrals containing a simple pole inside the interval. Analytical subtraction with integral or quadgk quad
weight='cauchy' 		GSL qawc QUADPACK DQAWC
Oscillatory integrals Rapid oscillations or explicit sine/cosine kernels. quadgk or specialized oscillatory methods quad
weight='sin' or 'cos' 		GSL qawo / qawf QUADPACK DQAWO / DQAWF
Fixed Gaussian quadrature Efficient high-order quadrature for smooth integrands. quadgk for adaptive Gaussian-type quadrature scipy.integrate.fixed_quad GSL fixed quadrature routines Library-specific Gaussian rules; custom Gauss–Legendre implementations are also common.
Multidimensional integration Low-dimensional nested quadrature or higher-dimensional integration. integral2
integral3 		nquad, dblquad, tplquad Nested one-dimensional quadrature;
GSL VEGAS / MISER for higher dimensions. 		Nested QUADPACK calls for low dimensions;
specialized or Monte Carlo methods for higher dimensions.

Benchmark integrals and reference values

These examples are useful for checking a new implementation, comparing algorithms, or testing changes to the companion programs. Whenever possible, compare both the integral and the reported error behavior.

Test problem Reference value What it tests
∫0π sin(x) dx 2 Basic composite and Gaussian quadrature.
∫01 dx / [1 + 400(x - 0.2)2] [atan(16) + atan(4)] / 20
≈ 0.1417097590		 Localized narrow peak and adaptive subdivision.
∫0∞ e-x sin(x) dx 1/2 Infinite interval and variable transformation.
PV ∫-11 cos(x)/(x - 0.2) dx ≈ -0.5912784964342436 Cauchy principal value and analytical subtraction.
∫01 e-x cos(50x) dx ≈ -1.6717740468665e-3 Rapid oscillation; comparison of general and oscillatory-specific methods.
∫01∫01 1/(1+x2+y2) dy dx ≈ 0.6395103518703110 Two-dimensional product Gauss–Legendre quadrature.
∫01∫0sin(x) x2/(y2+2) dy dx ≈ 0.1034464976429315 Nested integration with a variable inner limit.

Choosing an integration method

Character of the problem Method to consider Practical note
Smooth function on a finite interval Gauss–Legendre or a general adaptive library routine High-order Gaussian rules can be very efficient for smooth integrands.
Automatic error control is needed Adaptive Gauss–Kronrod Use absolute and relative tolerances appropriate to the physical problem.
Localized peak or rapidly varying region Adaptive Simpson or adaptive Gauss–Kronrod Plot the integrand first; adaptive subdivision should concentrate work where needed.
Infinite interval Variable transformation or a library routine supporting infinite limits Inspect the transformed integrand as well as the original one.
Simple pole inside the interval Analytical subtraction or a dedicated principal-value routine Do not treat an interior pole as an ordinary integrand.
Rapid oscillations Filon-type, weighted, or other oscillatory-specific quadrature General adaptive methods may require many function evaluations.
Smooth integral in two dimensions Product Gaussian quadrature or adaptive multidimensional routine Tensor-product rules are effective at low dimension but scale poorly as dimension increases.
Variable multidimensional boundary Nested integration Control the accuracy of inner integrations so they do not dominate the outer error.

Selected references and documentation

← Back to Chapter 7