**By R. Belmans**

## Linear system equation solvers

**There are various methods known** to solve a system of linear equations. In this section, we will give an overview of the most important methods suitable to solve equation systems obtained from finite element models of electromagnetic field problems with its special properties. We will consider a system of the form

The matrix A iws derived from a linear elliptic partial differential equation (PDE) by using a finite difference or finite element discretisation. Most of the numerical methods that will be discussed in this section are applicable to a wide range of problems.Here, we will focus on the standard example of Poisson’s equation

**Discretised by the** finite element method, matrix A will usually be extremely sparse. For the present problem, the matrix is also

**• symmetric**

**• positive definite**

**• structured.**

**These are properties that are not always satisfied for more general problems. The solution of linear systems derived from a PDE can be approached from two different points of view:**

**• algebraic:**

**One forgets about the nature of the** original problem and deals with (8.1) as an abstract matrix equation. Algebraic methods are usually very robust. They tend to be sub-optimal, however, but easy to apply, and widely available in standard software packages.

**• geometric:**

One tries to employ all available knowledge about the mesh and the nature of the PDE and boundary conditions. The system is seen as a collection of discrete finite element equations, whose molecules (stencils) reflect the connectivity of the unknowns and their relative strengths. Geometric methods are more specialised towards one particular class of problems. They are usually of more limited applicability, often hard to implement but possibly extremely effective.

**To study the theoretical aspects of the** different equation solver, refer to the literature (Briggs 23, George & Liu Heath 55, Saad 97). In this topic it is intended to give a rough overview of possible methods and with the properties of the particular algorithms discussed here, the reader is enabled to choose an appropriate equation solver for his class of problem.

The various methods to solve a linear system of equations are traditionally distinguished by direct and iterative methods.

## Direct methods

**Direct methods can be characterised as follows:**

**• They produce a solution to** the system of equations in a finite number of operations.

**• They are mainly of algebraic nature,** and widely applicable under extremely weak conditions (e.g. non-singularity). Some direct methods are geometric, e.g. most fast direct solvers. Certain methods use the geometric information reflected in the matrix adjacency graph, e.g., to determine an optimal renumbering of the unknowns.

**• The accuracy of the solution depends** on the conditioning of the problem, the numerical stability of the solver and the precision of the computer arithmetic. The accuracy can sometimes be improved by an iterative refinement technique.

**• They usually require large** amounts of memory, and may be very time-consuming.

**• They do not require an initial** estimate for the solution. On the other hand, they cannot take advantage of good initial approximations, which are sometimes available.

**• The user has little control over the** accuracy of the solution. In some cases, a low accuracy solution suffices. This is especially the case for PDE problems, where the discretisation error is often in the range of 10"2 (or even worse). Direct methods cannot take advantage of this fact.

## Iterative methods

**Iterative methods have the following qualities and drawbacks:** • Usually, they require little memory: storage for the nonzero matrix elements, the right-hand side, the solution, and perhaps a few additional vectors for temporary storage. In the geometric approach no matrix has to be stored explicitly. The storage for the matrix elements can be integrated into the mesh data structure.

**• Often such methods can be implemented** as soon as a procedure is available that computes the product of the matrix with an arbitraiy vector. Hence, such methods can be applied in cases where one does not want to construct the matrix. This occurs, for example, when solving non-linear problems with the Newton procedure (construction of Jacobian matrix can be avoided).

**• When optimally tuned,** they usually require much less work than direct methods. For good performance they may require accurate estimation of various problem-dependent parameters, such as the extreme-eigenvalues of the discretisation matrix. They may lose much of their potential performance, or even stop to converge, if the parameters are not chosen carefully.

**• They are less robust than direct methods.** They often require different conditions on the matrix to be satisfied (e.g., symmetry, positive definiteness,…).

**• A small perturbation to the system matrix** or to its structure may have a fatal effect on the convergence. In certain cases, one may want to resort to two-level iterative schemes, with the classical iteration used within a block-level iteration.