Appendix D: Mathematics Refresher

A compact reminder of the mathematics this book leans on — the linear algebra, calculus, and numerical intuition that Part V puts to work. It is written for a reader who once learned this material and wants it back quickly, in the exact form the code needs. Everything here is in service of the computation, not of rigour: there are no proofs, the notation matches the chapters, and each idea points to where the book actually uses it. If a topic is already comfortable, skip it. If a formula in Chapter 21, 22, or 24 looks unfamiliar, this is the page to read first.

Notation follows the book: scalars in italics ($x$, $\alpha$), vectors bold-lower ($\mathbf{v}$, $\mathbf{x}$), matrices capital ($A$). In Fortran a vector is a rank-1 array and a matrix a rank-2 array — the mathematics below maps one-to-one onto the arrays of Chapter 5.

D.1 Linear Algebra

Vectors and matrices. A vector $\mathbf{v} = (v_1, v_2, \dots, v_n)$ is an ordered list of $n$ numbers — a point (or arrow) in $n$-dimensional space, and in Fortran the rank-1 array v(n). A matrix $A$ is a rectangular grid of numbers with entry $A_{ij}$ in row $i$, column $j$; an $m \times n$ matrix has $m$ rows and $n$ columns and is the rank-2 array a(m,n). Fortran stores that grid column-major — first index fastest — which is exactly the layout LAPACK expects (Chapter 21).

Matrix–vector product. Multiplying a matrix by a vector produces a vector, $\mathbf{y} = A\mathbf{x}$, whose entries are

$$ y_i = \sum_{j=1}^{n} A_{ij}\, x_j . $$

Read it two ways: each $y_i$ is the dot product of row $i$ of $A$ with $\mathbf{x}$; equivalently, $A\mathbf{x}$ is a weighted sum of the columns of $A$. It costs about $n^2$ multiply-adds — an $O(n^2)$ operation (Level-2 BLAS in Chapter 21).

Matrix–matrix product. Multiplying two matrices, $C = AB$, gives

$$ C_{ij} = \sum_{k} A_{ik}\, B_{kj}, $$

the dot product of row $i$ of $A$ with column $j$ of $B$. Matrix multiplication is associative ($A(BC) = (AB)C$) but not commutative ($AB \neq BA$ in general). In Fortran it is the intrinsic matmul(a, b); it costs $O(n^3)$ operations on $O(n^2)$ data, the ratio that makes it compute-bound and the reason a tuned dgemm beats a hand-written loop (Chapters 21 and 29).

The linear system $A\mathbf{x} = \mathbf{b}$. This is the central problem of numerical linear algebra: given a known $n \times n$ matrix $A$ and a known right-hand side $\mathbf{b}$, find the unknown vector $\mathbf{x}$ that satisfies $n$ linear equations at once. For example,

$$ \begin{aligned} x + y + z &= 6 \\ 2y + 5z &= -4 \\ 2x + 5y - z &= 27 \end{aligned} \qquad\Longleftrightarrow\qquad A = \begin{bmatrix} 1 & 1 & 1 \\ 0 & 2 & 5 \\ 2 & 5 & -1 \end{bmatrix},\quad \mathbf{b} = \begin{bmatrix} 6 \\ -4 \\ 27 \end{bmatrix}, $$

whose solution is $\mathbf{x} = (5, 3, -2)$. Chapter 21 solves exactly this with one call to LAPACK's dgesv.

Identity, transpose, inverse. The identity matrix $I$ has $1$s on its diagonal and $0$s elsewhere; it is the "do nothing" matrix, $I\mathbf{x} = \mathbf{x}$ and $AI = IA = A$. The transpose $A^{\mathsf{T}}$ swaps rows and columns, $(A^{\mathsf{T}})_{ij} = A_{ji}$; a matrix with $A = A^{\mathsf{T}}$ is symmetric. The inverse $A^{-1}$ is the matrix that undoes $A$, meaning $A A^{-1} = A^{-1} A = I$. Formally the solution of the system above is $\mathbf{x} = A^{-1}\mathbf{b}$ — and in code you should almost never compute it that way.

Why we solve rather than invert. Forming $A^{-1}$ and then multiplying is both slower and less accurate than solving $A\mathbf{x} = \mathbf{b}$ directly. The direct route factors $A$ into a product of a lower- and an upper-triangular matrix, $A = LU$ (Gaussian elimination with pivoting), which costs about $\tfrac{2}{3}n^3$ operations; the two triangular solves that follow cost only $O(n^2)$. Building the full inverse costs several times more arithmetic and folds extra rounding into every entry. This is why dgesv returns $\mathbf{x}$, never $A^{-1}$, and why the professional habit (Chapter 21) is state the system, solve it, check the residual $A\mathbf{x} - \mathbf{b}$ — not invert.

Determinant (briefly). The determinant $\det A$ is a single number attached to a square matrix. For a $2 \times 2$ matrix it is

$$ \det \begin{bmatrix} a & b \\ c & d \end{bmatrix} = ad - bc, $$

and geometrically $|\det A|$ is the factor by which the map $A$ scales area (in 2D) or volume (in higher dimensions). The one fact that matters for computation: $A$ is invertible — the system $A\mathbf{x} = \mathbf{b}$ has a unique solution — if and only if $\det A \neq 0$. A matrix with $\det A = 0$ is singular. In practice you never compute a large determinant to test this; the factorization inside dgesv reports singularity directly through its info status.

Eigenvalues and eigenvectors. For a square matrix $A$, a nonzero vector $\mathbf{v}$ that $A$ merely stretches — without rotating it — is an eigenvector, and its stretch factor $\lambda$ is the eigenvalue:

$$ A\mathbf{v} = \lambda \mathbf{v}. $$

Eigenvalues expose the intrinsic behaviour of the operator $A$: the natural frequencies of a vibrating system, the principal axes of a dataset, the growth or decay rates of a dynamical process. Chapter 21 computes them for a real symmetric matrix with LAPACK's dsyev. There is a direct payoff in this book: the eigenvalues of the discrete Laplacian are what set the stability limit of the explicit heat scheme in Chapter 24 — the same $r \le 1/4$ you meet below.

D.2 Calculus

Derivatives. The derivative of a function measures its instantaneous rate of change — the slope of the tangent to its graph:

$$ f'(x) = \frac{df}{dx} = \lim_{h \to 0} \frac{f(x+h) - f(x)}{h}. $$

The notations $f'(x)$ and $\tfrac{df}{dx}$ are interchangeable; for a rate of change in time the dot form $\dot{u} = \tfrac{du}{dt}$ is common. The second derivative $f''(x) = \tfrac{d^2 f}{dx^2}$ is the rate of change of the slope — the curvature — and it is positive where the graph bends upward (a valley) and negative where it bends down (a hill). A computer cannot take the limit $h \to 0$, so it stops at a small finite $h$; that single act of stopping early is the finite-difference idea of §D.3 and Chapter 22.

The chain rule (briefly). To differentiate a composition — a function of a function — multiply the outer derivative by the inner:

$$ \frac{d}{dx} f\big(g(x)\big) = f'\big(g(x)\big)\, g'(x). $$

For instance $\tfrac{d}{dx}\sin(x^2) = \cos(x^2)\cdot 2x$. It is the workhorse behind most derivatives you meet, and the conceptual root of automatic differentiation.

Integrals as area. The definite integral $\int_a^b f(x)\,dx$ is the signed area between the curve $y = f(x)$ and the $x$-axis from $a$ to $b$ — area above the axis counts positive, below negative. The Fundamental Theorem of Calculus ties it to the derivative: if $F' = f$, then

$$ \int_a^b f(x)\,dx = F(b) - F(a). $$

When no closed-form antiderivative $F$ exists — the common case in science — you approximate the area by a weighted sum of samples, which is exactly the quadrature of Chapter 22 (the trapezoidal, Simpson, and Gauss rules).

Partial derivatives, the gradient, and the Laplacian. When a quantity depends on several variables — temperature $u(x, y, t)$ on a plate, say — its partial derivative with respect to one variable is its rate of change with the others held fixed. The curly $\partial$ signals this: $\tfrac{\partial u}{\partial x}$ varies $x$ alone, $\tfrac{\partial u}{\partial t}$ varies $t$ alone. Collect the first spatial partials into a vector and you have the gradient,

$$ \nabla u = \left( \frac{\partial u}{\partial x},\ \frac{\partial u}{\partial y} \right), $$

which points in the direction of steepest increase of $u$ and whose length is that steepest slope (heat flows down the gradient, from hot to cold). Take the divergence of the gradient — sum the pure second partials — and you get the Laplacian, the single most important operator in this book:

$$ \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}. $$

Intuitively $\nabla^2 u$ at a point compares that point to the average of its neighbours: it is positive when the surroundings are, on balance, higher (the point sits in a valley of the field), negative at a local peak, and zero when the point equals its neighbourhood average. That reading is the whole physical content of the heat equation in §D.4, and §D.3 shows how to compute it on a grid.

D.3 Taylor Series and the Finite-Difference Formulas

The Taylor series expresses the value of a smooth function a small step $h$ away from $x$ in terms of the function and its derivatives at $x$:

$$ f(x+h) = f(x) + h\,f'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(x) + \frac{h^4}{24}f^{(4)}(x) + \cdots $$

The general term is $\tfrac{h^n}{n!}f^{(n)}(x)$. Replacing $h$ with $-h$ flips the sign of every odd-power term:

$$ f(x-h) = f(x) - h\,f'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{6}f'''(x) + \frac{h^4}{24}f^{(4)}(x) - \cdots $$

These two expansions are the source of every finite-difference formula — and of the error in each one. The trick is always the same: combine them so the derivative you want survives and the rest is a known power of $h$.

Forward difference. Rearrange the first expansion by dropping $f(x)$ and dividing by $h$:

$$ \frac{f(x+h) - f(x)}{h} = f'(x) + \underbrace{\frac{h}{2}f''(x) + \cdots}_{\text{error}} . $$

The approximation $f'(x) \approx \big(f(x+h) - f(x)\big)/h$ carries a leftover error whose leading term is proportional to $h$: it is first-order accurate, $O(h)$. (Looking backward to $x-h$ instead gives the equally first-order backward difference.)

Central difference. Now subtract the second expansion from the first. The $f(x)$ terms cancel, the even-power terms cancel (they share a sign), and the odd terms reinforce:

$$ f(x+h) - f(x-h) = 2h\,f'(x) + \frac{h^3}{3}f'''(x) + \cdots $$

Divide by $2h$:

$$ \frac{f(x+h) - f(x-h)}{2h} = f'(x) + \underbrace{\frac{h^2}{6}f'''(x) + \cdots}_{\text{error}} . $$

The leading error is now proportional to $h^2$: the symmetric formula is second-order accurate, $O(h^2)$. The $f''$ term that dominated the forward difference's error has cancelled, which is why the central difference is markedly more accurate for the same $h$.

Second difference — the stencil. Add the two expansions instead. Now the odd terms cancel and the $f''$ terms survive:

$$ f(x+h) - 2f(x) + f(x-h) = h^2 f''(x) + \frac{h^4}{12}f^{(4)}(x) + \cdots $$

Divide by $h^2$ to get the three-point approximation to the second derivative:

$$ f''(x) \approx \frac{f(x+h) - 2f(x) + f(x-h)}{h^2}, \qquad \text{error } \frac{h^2}{12}f^{(4)}(x) = O(h^2). $$

This is the one-dimensional heart of the Laplacian, with the memorable coefficient pattern $(1, -2, 1)$. Apply it in $x$ and again in $y$ and add, and — on a square grid with $\Delta x = \Delta y = h$ — the two second differences combine into the five-point stencil of Chapter 24:

$$ \nabla^2 u \big|_{i,j} \approx \frac{u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4\,u_{i,j}}{h^2}. $$

You can check it by hand on $u = x^2 + y^2$, whose exact Laplacian is $2 + 2 = 4$ everywhere. Each neighbour contributes a full field value — the east neighbour is $(x+h)^2 + y^2$, the north neighbour is $x^2 + (y+h)^2$, and so on — so the four neighbours sum to $4x^2 + 4y^2 + 4h^2$. Subtracting $4u_{i,j} = 4(x^2 + y^2)$ leaves exactly $4h^2$, which divided by $h^2$ gives $4$ for any spacing. The stencil is exact for quadratics because its error involves the fourth derivative, which vanishes here. That the whole stencil really converges at $O(h^2)$ is what the Chapter 22 Project Checkpoint verifies numerically.

D.4 The Heat / Diffusion Equation

The book's running project solves the heat equation (identical mathematics to the diffusion equation). In two dimensions:

$$ \frac{\partial u}{\partial t} = \alpha \nabla^2 u = \alpha \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right). $$

Read each piece physically. Here $u(x, y, t)$ is the temperature at a point of the plate at a moment in time. The left side $\tfrac{\partial u}{\partial t}$ is the rate at which that point heats up or cools down. The Laplacian $\nabla^2 u$ on the right is the spatial curvature — the neighbourhood comparison of §D.2 — and the constant $\alpha > 0$ is the thermal diffusivity, how readily the material conducts heat (large for copper, small for wood). The equation therefore says: a point's temperature changes at a rate proportional to how much hotter or colder its surroundings are than it is. A point cooler than its neighbours warms ($\nabla^2 u > 0$); a local hot spot cools ($\nabla^2 u < 0$). That is diffusion in one sentence — everything drifts toward the average of what surrounds it.

The one-dimensional version, a heated rod, keeps a single spatial term:

$$ \frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}. $$

Steady state. When the plate stops changing — $\tfrac{\partial u}{\partial t} = 0$ everywhere — the equation collapses to

$$ \nabla^2 u = 0, $$

Laplace's equation, whose solution is the smooth, settled temperature distribution the plate relaxes into. In one dimension $u'' = 0$ means a straight line, which is why a rod held at $0$ and $100$ degrees settles to a linear profile.

To turn this into a simulation, Chapter 24 discretizes: the Laplacian becomes the five-point stencil of §D.3, and the time derivative becomes a forward difference in time (the forward-Euler step you can also read as the method of lines from Chapter 23). The explicit update is

$$ u^{n+1}_{i,j} = u^n_{i,j} + r\left(u^n_{i+1,j} + u^n_{i-1,j} + u^n_{i,j+1} + u^n_{i,j-1} - 4u^n_{i,j}\right), \qquad r \equiv \frac{\alpha\,\Delta t}{h^2}, $$

where the superscript $n$ counts time steps. The single dimensionless number $r$ — physics, timestep, and grid rolled together — controls everything, and §D.5 explains the hard limit on how large it may be.

D.5 A Little Numerical-Analysis Intuition

Four ideas recur across Part V. None is difficult; together they are the difference between a number you can trust and a plausible lie.

Order of accuracy. A discretization has order $p$ when its error shrinks like a fixed power of the step size, $E(h) \approx C\,h^{p}$, for some constant $C$. The exponent $p$ — not the constant — is what matters, because it decides how fast refining the grid pays off: halving $h$ divides an $O(h^p)$ error by $2^p$. So a first-order method's error halves, a second-order method's error quarters, and a fourth-order method's error drops by sixteen. This gives a universal correctness check (Chapter 22): compute the error at $h, h/2, h/4, \dots$, take successive ratios, and confirm they approach $2^p$. If the ratio is wrong, your code has a bug. The forward difference is $p = 1$; the central difference, the trapezoidal rule, and the second-difference stencil are $p = 2$; Simpson's rule is $p = 4$.

Convergence. A method converges if its answer approaches the true one as the grid is refined — $E(h) \to 0$ as $h \to 0$. Convergence is the minimum any usable method must have; the order of accuracy then says how fast it converges. A method that fails to converge, or "converges" only until round-off takes over, is telling you something is wrong with the method, the code, or the problem.

Stability. A computation is stable when errors already present — from rounding, from the input data — stay bounded as it proceeds, instead of growing without limit. Stability wears two faces in this book. As a property of an algorithm (Chapter 20), a stable method does not manufacture extra error through its own arithmetic (for example, by internally subtracting nearly equal numbers). As a property of a time-stepping scheme (Chapter 24), stability means errors do not amplify from one step to the next — and for the explicit heat scheme it is a hard threshold on the timestep, the CFL condition:

$$ r = \frac{\alpha\,\Delta t}{h^2} \le \frac{1}{4} \quad (\text{2D}). $$

The limit is $\tfrac12$ in 1D and $\tfrac16$ in 3D — the pattern is $\tfrac{1}{2d}$ in $d$ dimensions. Stay under it and the simulation marches sensibly toward steady state; cross it and errors double every step until the numbers overflow. Because $r$ depends on $h^2$, halving the grid spacing forces a four-times-smaller timestep — the reason explicit diffusion codes get slow on fine grids.

Conditioning. Where stability is a property of the method, conditioning is a property of the problem: how much the answer changes when the input is perturbed a little. A well-conditioned problem amplifies small input errors only a little; an ill-conditioned one amplifies them enormously. The amplification factor is the condition number $\kappa$; a rough rule (Chapter 20) is that $\kappa \approx 10^{k}$ costs you about $k$ of your roughly $16$ significant decimal digits in double precision. The slogan to remember: conditioning is the problem's fault; stability is yours. If a linear system's matrix is ill-conditioned, no solver — not even LAPACK — can recover the digits the problem itself destroys; a small residual $A\mathbf{x} - \mathbf{b}$ can then still hide a wrong $\mathbf{x}$, which is why Chapter 21 estimates the condition number rather than trusting info = 0 alone.

D.6 Big-O and Orders of Magnitude

When you ask how expensive an algorithm is, you rarely want an exact instruction count; you want to know how the cost grows as the problem gets bigger. Big-O notation captures that growth by keeping only the dominant term and dropping constants: an algorithm is $O(n^3)$ if, for large $n$, its work grows in proportion to $n^3$. The everyday shapes:

The exponent dominates once $n$ is large: double the size of an $O(n^2)$ job and it costs four times as much; double an $O(n^3)$ job and it costs eight times as much. Memory scales the same way — a dense $n \times n$ matrix needs $O(n^2)$ storage, which is why the million-unknown heat system of Chapter 21 would need terabytes if you stored every zero, and why sparse formats exist.

This cost model is the backbone of Part VII (Performance). It is why you profile first (Chapter 28) to find the term that actually dominates the runtime — an untouched $O(n^2)$ loop can own ninety percent of the time while you tune everything else in vain. It is why the arithmetic intensity of a kernel — how many operations it does per byte of memory it moves — decides whether it can run near the processor's peak (Chapters 27 and 29). And it is why two algorithms with the same Big-O can still differ by a large constant factor: a hand-rolled $O(n^3)$ matrix multiply and a tuned $O(n^3)$ dgemm do the same number of operations, yet the library wins by an order of magnitude because it moves memory more cleverly (Chapter 29). Big-O tells you which battles are worth fighting; the constant factor tells you how to fight the one that is.

D.7 Quick Reference: where each idea is used

For the machine-arithmetic side of these ideas — machine epsilon $\varepsilon_{\text{mach}} \approx 2.2\times10^{-16}$, catastrophic cancellation, and why a smaller step is not always better — Chapter 20 is the companion to this refresher, and Appendix C covers the compiler flags that turn these methods into fast, checkable programs.