Fundamentals of Numerical Computation: Julia Edition

The front matter includes the Title page, Copyright page, Dedication, Table of Contents, Preface to Julia Edition, Preface to the Original Edition, and Glossary.

Our first step is to discretize the real numbers—specifically, to replace them with a finite surrogate set of numbers. This step keeps the time and storage requirements for operating with each number at constant levels, but virtually every data set and arithmetic operation is perturbed slightly away from its idealized mathematical value. We can easily keep the individual roundoff errors very small, so small that simple random accumulation is unlikely to bother us. However, some problems are extremely sensitive to these perturbations, a trait we quantify using a condition number. Problems with large condition numbers are difficult to solve accurately using finite precision. Furthermore, even when the condition number of a problem is not large, some algorithms for solving it allow errors to grow enormously. We call these algorithms unstable. In this chapter we discuss these ideas in simple settings before moving on to the more realistic problems in the rest of the book.

One of the most frequently encountered tasks in scientific computation is the solution of the linear system of equations $\mathbf{Ax} = \mathbf{b}$ for a given square matrix $\mathbf{A}$ and vector $\mathbf{b}$. This problem can be solved in a finite number of steps, using an algorithm equivalent to Gaussian elimination. Describing the algorithm is mostly an exercise in organizing some linear algebra.

So far we have considered $\mathbf{Ax}=\mathbf{b}$ only when $\mathbf{A}$ is a square matrix. In this chapter we consider how to interpret and solve the problem for an $m\times n$ matrix where $m\gt n$—and in practice, $m$ is often much larger than $n$. This is called an overdetermined linear system because, in general, the system has more equations to satisfy than the variables allow. The complementary underdetermined case $m\lt n$ turns up less frequently and will not be considered in this book.

In this chapter we extend from linear algebra to deal with nonlinear algebraic problems. This kind of problem arises when there is a parameter or variable that can be changed in order to satisfy a constraint or achieve some goal. We start with scalar functions of a single variable, then generalize to $n$ variables and $n$ nonlinear equations. Finally, we generalize the problem of linear least squares to situations with more nonlinear constraints to satisfy than there are variables. In every case the strategy used is one of the cornerstones of numerical computing: replace a problem you can’t solve with an approximate one that you can. In the context of nonlinear algebraic problems, the particular tactic is to set up and solve a sequence of linear problems of the types covered in the two previous chapters.

In many scientific problems the solution is a function. Accordingly, our next task is to represent functions numerically. This task is more difficult and complicated than the one we faced in representing real numbers. With numbers it’s intuitively clear how one real value can stand for a small interval around it. But designating representatives for sets of functions is less straightforward—in fact, it’s one of the core topics in computing. The process of converting functions into numerical representations of finite length is known as discretization.

Quantities that change continuously in time or space are often modeled by differential equations. When everything depends on just one independent variable, we call the model an ordinary differential equation (ODE). Differential equations need supplemental conditions to define both the modeling situation and the theoretical solutions uniquely. The initial-value problem (IVP), in which all of the conditions are given at a single value of the independent variable, is the simplest situation. Often the independent variable in this case represents time.

In previous chapters, we have seen how matrices that represent square or overdetermined linear systems of equations can be manipulated into LU and QR factorizations. But matrices have other factorizations that are more intrinsic to their nature as mathematical linear transformations. The most fundamental of these are the eigenvalue and singular value decompositions.

What are the implications of the $O(n^3)$ work requirements for solving linear systems? Suppose tomorrow your computer became a thousand times faster. (Historically this has taken about 15 years in the real world.) Assuming you are willing to wait just as long today as you were yesterday, the size of the linear system you can solve has gone up only by a factor of 10. Nice, but not nearly the jump that you got in hardware power. In fact, there is an odd paradox: faster computers make faster algorithms more important, not less, because they demand that you work at larger values of $n$, where asymptotic differences are large.

In Chapter 5 we considered a few ways to map data values to functions via interpolation. The methods we deemed successful were piecewise low-degree polynomials. In this chapter we deal with approximations that are globally defined over the entire interval, not piecewise.

In Chapter 6 we examined how to solve initial-value problems (IVP) for ordinary differential equations (ODEs). In an IVP the supplemental conditions give complete information about the state of the system at one value of the independent variable. However, not all ODE problems come in this form. Instead we might have only partial information about the solution at two different points. Such ODE problems are called boundary-value problems (BVP).

To this point we have considered only ordinary differential equations, those having only one independent variable. In the rest of this book, we introduce the huge topic of solving partial differential equations (PDEs). We begin by pairing time with space.

Now that we have seen PDEs with both time and space variables, we have a new wrinkle to add. Some of these equations behave like those of the previous chapter, creating diffusive effects. Others, though, are about propagation or advection—generally, the behavior of waves.

We have graduated from ODEs, usually having either time or space as the independent variable, to PDEs in which both space and time are represented simultaneously. The final innovation in this book is to consider problems with two space dimensions. We will confine ourselves to the simplest possible class of two-dimensional regions, leaving aside the major issue of the geometry of the domain. Even so there are many important and valuable mathematical models within this constraint.

An ordinary number in $\mathbb{R}$ or $\mathbb{C}$ may be called a scalar. An $m\times n$ matrix $\mathbf{A}$ is a rectangular $m$-by-$n$ array of numbers called elements or entries. The numbers $m$ and $n$ are called the row dimension and the column dimension, respectively; collectively they describe the size or shape of $\mathbf{A}$. We say $\mathbf{A}$ belongs to the set $\mathbb{R}^{m\times n}$ if its entries are real, or $\mathbb{C}^{m\times n}$ if they are complex-valued. A square matrix has equal row and column dimensions. A row vector has dimension $1\times n$, while a column vector has dimension $m \times 1$.

This is the default way to start and run Julia. • Press ] to enter package management mode; backspace to leave it. • Press Ctrl-C to interrupt execution. • Press Ctrl-D to exit. • There are other interactive utilities.

The back matter includes bibliography, index, and back cover.