Solving PDEs in C++: Numerical Methods in a Unified Object-Oriented Approach

The finite-element discretization method produces a large system of algebraic equations called the stiffness system. The number of equations in this system, N, is as large as the number of nodes in the mesh. Actually, the system can be written algebraically as
where A is the coefficient (stiffness) matrix of order N, f is the given N-dimensional right-hand-side vector, and x is the N-dimensional vector of unknowns. The matrix A is sparse in the sense that most of the elements in it are zero. Once this system is solved numerically, the ith component in x, x i is the value of the numerical solution of the partial differential equation (PDE) at the node of index i in the list of nodes in the mesh. The approximate solution at points that are not in the grid can also be obtained by some interpolation.
The numerical solution of a linear system of the form Ax = f is, thus, an essential part of the numerical solution of PDEs. In the past, before computers became available, this task had to be done manually using Gaussian elimination (Chapter 2, Section 21). Of course, this approach is impractical unless the order of the system is very small.
Fortunately, the digital computer freed the human mind from the exhausting...