This book is meant to teach an advanced undergraduate physics major how to solve problems, that they can easily pose and understand, using numerical methods. It is self-contained (within reason) in its presentation of physical problems, so that a reader need not be previously familiar with all of the physical topics covered. The numerical techniques that are presented are complete, without being exhaustive. I have chosen a set of tools that can be motivated easily – some of these, like Runge–Kutta methods for solving ordinary differential equations, are themselves the “end of the story” for certain problems. In other cases, I choose to present a simplified method that will work for at least some problems, even if it is not what is actually implemented in practice. A homogenous and appropriate level of presentation is the goal. Sometimes that allows for a complete and formal discussion, in other cases, I motivate and inform while ensuring that a successful (if not competitive) method is carefully described.

The chapters are defined by computational technique, with chapters covering “Partial differential equations,” “Integration,” “Fourier transform,” “Matrix inversion,” etc. Each chapter introduces a disparate set of physical problems, most of which are “unsolvable” (meaning that there is no closed-form solution expressible in terms of simple functions). I have attempted to draw problems from as wide an array of physical areas as possible, so that in a single chapter, there may be examples from quantum mechanics, E&M, special relativity, etc.

We will start off our investigation of numerical techniques by solving ordinary differential equations (ODEs). Our first method is a relatively direct discretization of Newton's second law. This is called the Verlet method, and is an efficient way to find solutions (i.e. a vector function x(t)) given a force that depends on time and position (and appropriate initial data). The method has generalizations, but we will think about it only in the context of equations of motion. As a first method, it has the advantage of being straightforward, while relying in its derivation on ideas that can be extended to develop more involved numerical differential equation solvers.

In the more general context that follows the Verlet method, we note that any ODE can be re-written in the form: f′(x) = G(x, f(x)) for vector-valued function G and vector-valued target function f(x). The punch-line will be that if we are given initial conditions f(0), we can use Runge–Kutta methods to find values approximating f(x) on a grid. First, we'll set up a few ODEs of interest, then introduce the process of discretization and the (special) Runge–Kutta family of methods that can be used to solve the discretized problem.

Many problems that we have seen already require that we “minimize” a certain function, or find the “minimum” of a function of many variables. Formally, the problem is: given F(x), find a value x0 such that F(x) > F(x0) for all x, for a global minimum, or F(x) > F(x0) where x ϵ [a, b] (some constrained range) for a local minimum. We will review these familiar problems, and introduce some new ones. Our first job will then be to write a function that can find the minimum of a function of one variable, a sort of minimization analogue to the bisection method from Section 3.2.1 for finding roots of F(x). Then we will study methods that can be used to minimize a function u(x) where x ϵ ℝn.

Physical motivation

There are many places in physics where finding the minimum of a function can be useful. We have already done linear minimization in solving the least squares problem for data fitting from Section 10.2 – there we turned the minimization of the Pythagorean length of the residual vector into a matrix inverse problem with a unique solution. There are nonlinear extensions of this process that rely on iterative techniques directly (as opposed to the matrix inversion that solves the least squares problem, which can be performed iteratively or not).

Recall the harmonic approximation from Chapter 8: Given the potential energy U(x), some complicated function, we know that in the vicinity of a local minimum, the motion of a particle will be harmonic with effective spring constant set by the second derivative of the potential evaluated at the minimum.

A programming language useful to this book must provide a minimal set of components that can be used to combine numbers, compare quantities and act on the result of that comparison, repeat operations until a condition is met, and contain functions that we can use to input data, and output results. Almost any language will suffice, but I have chosen to use Mathematica's programming environment as the vehicle. The reasoning is that 1. The input/output functions of Mathematica are easy to use, and require little additional preparation 2. We will be focused on the ideas, numerical and otherwise, associated with the methods we study, and I want to draw a clear distinction between those ideas and issues of implementation. This book is not meant to teach you everything you need to know about programming – we will discuss only the bare essentials needed to implement the methods. Instead, we will focus on the physical motivation and tools of analysis for a variety of techniques. My hope is that the use of Mathematica allows us to discuss implementation in a homogeneous way, and our restriction to the basic programming structure of Mathematica (as opposed to the higher-level functionality) allows for easy porting to the language of your choice.

The Fourier transform plays a major role in the solution of ordinary and partial differential equations, and so is of theoretical utility in physics. Just as important is its applications to data analysis, and signal processing, making the Fourier transform an indispensible tool in experimental physics as well. The reason for its centrality is the connection the transform makes between temporal data and frequency information. A time-varying sinusoidal signal may have complicated patterns of repetition, and yet be made up of only a few frequencies. A signal with ten thousand data points may be describable in terms of two or three frequency values. We know how the frequencies combine, so we can reconstruct the signal.

Where do physical oscillations come from? The harmonic potential shows up far more than its springy roots would imply. Any local minimum in a complicated potential has an effective spring constant, and that effective spring constant determines a time scale of oscillation, i.e. a period (whose inverse defines the frequency). If we know the frequencies associated with a signal, then, we know something about the local minimum defined by some complicated potential.

There is also a human element to the Fourier transform. Much of the information we process can be described in terms of frequencies – both light and sound come to mind. So our ability to isolate and manipulate the individual frequencies of physiological input can lead to interesting changes in perception.

We have seen some direct methods for inverting, and calculating the eigenvalues and eigenvectors of matrices. In many cases, however, these direct methods will take too long to solve a particular problem. In addition, we may only have interest in approximate solutions. When solving the Poisson problem, for example, we may already be on a coarse grid, and require only a qualitative description of the solution – high precision and “exact” matrix inverses are unnecessary. On the eigenvalue side, we may only need part of the spectrum of a matrix, maybe we just want a few bound states for a potential in quantum mechanics, for example. In these cases, what we would like is a process that ultimately would give the full inverse, or the complete set of eigenvalues, but a process that can be truncated “along the way” and still provide partial information.

There are two broad schemes for these approximate methods, and we'll describe and see examples of both. The first approach is relevant to matrix inversion, and involves decomposing a matrix into a simple (to invert) part, and a (hopefully small) “other” part. We proceed to invert the simple part and use that inversion to drive an iteration that will converge to the exact numerical solution (computed using QR factorization, for example). The second approach involves constructing a particular subspace of ℝn, called a Krylov subspace, and we invert matrices and find eigenvalues within that subspace.

Coupled, nonlinear sets of ODEs, of the sort that describe most physical processes, do not always have well-behaved numerical solutions (as we saw in the last chapter). In addition, there is certain physical behavior that cannot be (or at any rate has not been) rendered into well-defined, deterministic, ODEs. Whatever the source, there are some problems where we have a large amount of data, but no good rule for generating output from input. Suppose we have an electromagnetic signal that is supposed to propagate through space to an antenna on the ground – as the signal enters the Earth's atmosphere … things happen. The signal is scattered by the atmosphere, for example, and the properties of the air that govern that scattering change as the light makes its way to an antenna on the Earth. Now there is nothing physically obscure about this process, and yet one can imagine that the particular model of the atmosphere plays a large role in taking us from an observed signal to properties of its source. If we had associated pairs of known source signals together with ground-based measurements of them, we could side-step the model-dependent portion of the problem, and simply estimate the signal given the measurement by comparing with the previously observed data.