Using Numerical Mathematics Computing 7th Edition in Practice
The book itself is a reference work first and foremost. It covers iterative methods for linear systems, numerical linear algebra, optimization, and differential equations. The code examples are in MATLAB. That is the most important thing to understand upfront, because MATLAB's floating point behavior does not exactly match what you will see in production C or Fortran routines. The companion code is available from the SIAM website under the book's product page. You will need to register or log in. The download is a zip file, roughly 85 megabytes, containing M-files for every chapter. Extract it to a dedicated folder. Do not scatter the files across your MATLAB path or you will run into naming conflicts. The chapters use custom wrapper functions that shadow built-in MATLAB names in some cases, and that causes debugging headaches if you are not paying attention. Before running anything, check your MATLAB version. Some of the later examples use `fzero` with tighter tolerances and `ode15s` options that assume a version from 2014 onward. If you are on an older install, set `MaxOrder` to 5 in the `ode15s` options structure and you will avoid a common error on stiff test problems.
Iterative Methods and Why They Fail
Chapter 4 on iterative methods for linear systems is where most people hit their first wall. The book derives Richardson, Jacobi, Gauss-Seidel, and SOR from first principles. The derivations are clean. The implementation details are where things get messy in practice. The SOR implementation assumes you have the diagonal of the matrix available as a separate vector. If your matrix comes from a sparse factorization or an implicit discretization, the diagonal might be scattered or partially zero. I encountered this when working with a discretized reaction-diffusion system where the coupling term produced near-zero diagonal entries at boundary nodes. The textbook code stalled immediately because the relaxation parameter omega required a non-zero pivot. My workaround was to add a diagonal shift of 1e-12 before passing the matrix to the solver. This is not elegant, but it keeps the iteration going long enough to converge. You should monitor the residual after adding the shift because it changes the condition number slightly. A counter-intuitive point that the book does not emphasize enough: SOR convergence is not monotonic in the way introductory courses suggest. The residual norm can oscillate wildly before settling, especially when omega is near the optimal value. I spent two days debugging a code that appeared to diverge until I stopped plotting the residual every iteration and instead logged only every twentieth step. The behavior was normal under-relaxation dynamics, not a bug.
Numerical Linear Algebra Section
Chapter 3 covers Gaussian elimination with partial pivoting and the LU factorization. The theory section is standard. The practical warnings matter more. MATLAB's backslash operator uses different algorithms depending on matrix structure. When the book shows a custom LU implementation, do not confuse it with `mldivide`. The custom implementation includes explicit pivoting arrays and column norms for scaling. Running the textbook code alongside `A\b` on the same matrix gives slightly different results because of the order of operations and threshold settings. That is expected. The difference is on the order of machine epsilon for moderate sized systems, but can grow for ill-conditioned matrices. The condition number estimation routine in the companion code uses a one-norm estimate. For matrices larger than 5000 by 5000, this can become expensive. I found that switching to `condest` with a relaxed tolerance gave sufficient accuracy for my purposes and cut the estimation time from roughly 40 seconds down to under 6 seconds on a standard workstation. The book does not mention this trade-off, but it is worth knowing if you are working with large sparse systems.
Get the Full Details
Ordinary Differential Equations
Chapters 8 and 9 treat ODEs. The explicit and implicit Runge-Kutta methods are explained thoroughly. The test problem involving the Robertson chemical kinetics system is a standard benchmark. The issue here is stiffness detection. The textbook assumes the reader understands when a method is stiff without explicitly mapping out every case. A variable-step explicit method will take thousands of tiny steps on the Robertson system. Using `ode15s` with the default settings resolves it in under a hundred steps. The code provided in the book uses a fixed-step BDF method, which is accurate but slower for this particular problem. If you are following along, replace the fixed-step integrator with `ode15s` for the Robertson example and you will see the speed difference immediately. Another subtlety: the Jacobian sparsity pattern matters more than its values for stiff solvers. I wasted an afternoon running a custom Jacobian routine that computed every entry as nonzero when the true sparsity pattern had over 90 percent zeros. Switching to a pattern-correct sparse Jacobian reduced runtime from twelve minutes to forty seconds for a single trajectory. The book's examples compute dense Jacobians for simplicity, which is fine for small systems but becomes a liability quickly.
Common Pitfalls
Three things to watch for when working through this book: First, the indexing. Several examples use one-based indexing consistent with MATLAB, but a few derivations in the text assume zero-based indexing from a C or Fortran perspective. If you port the code elsewhere, recheck loop bounds. Second, tolerance settings. The default stopping criteria in the companion code are often tighter than what is needed for practical work. Setting a loose tolerance of 1e-6 instead of 1e-12 for a preconditioned iterative solve typically makes no meaningful difference to the final answer while cutting runtime significantly.
Third, the book does not cover parallel implementations extensively. Most of the methods are serial. If you are working with large-scale problems, you will need to look beyond the companion code for MPI or GPU variants. There are third-party toolboxes that extend these methods, but they are not covered here.
When This Book Falls Short
For modern production code, the MATLAB examples are useful for understanding the algorithms, but they are not ready for deployment. The code lacks error handling for edge cases like singular matrices or convergence failure. You will need to add your own guards if you plan to use any of these routines in a larger system. For classroom learning or self-study, the book remains solid. For anything beyond that, supplement it with library documentation from packages like SuperLU or PETSc.
