Numerical methods for linear inhomogeneous sets of equations







Numerical Methods

Numerical methods for linear inhomogeneous sets of equations

Let us consider the $N$ equations

\[\begin{equation}\label{lin1} \begin{array}{c} A_{11}x_{1}+A_{12}x_{2}+ \cdots +A_{1N}x_{N}=b_{1} \\ A_{21}x_{1}+A_{22}x_{2}+ \cdots +A_{2N}x_{N}=b_{2} \\ . \\ . \\ . \\ A_{N1}x_{1}+A_{N2}x_{2}+ \cdots +A_{NN}x_{N}=b_{N} \end{array} \end{equation}\]

where the $A_{ij}$ and the $b_{i}$ are real quantities. We further require that:

\[\begin{equation} \mid b_{1} \mid + \mid b_{2} \mid + \cdots + \mid b_{N} \mid \ne 0 \end{equation}\]

With these assumptions (\ref{lin1}) constitutes a real, linear, inhomogeneous set of equations (or linear system) of $N$-th order. The quantities $x_{1},\ldots,x_{N}$, which satisfy simultaneously all the above equations, are called the solutions of the system.

(\ref{lin1}) is usually expressed in matrix form:

\[\begin{equation}\label{lin2} A \cdot \vec{x} = \vec{b} \end{equation}\]

The solution of this kind of systems is a central problem in numerical mathematics, since a wide range of numerical methods, such as methods for the numerical interpolation, least-square methods, differential methods etc. can be reduced to the problem of solving a set of inhomogeneous linear equations.

From a theoretical point of view, the solution of (\ref{lin1}) does not present any difficulty, as long as the determinant of the matrix of coefficients $A$ does not vanish, i.e. as long as the problem is non-singular: \[\begin{equation} \mbox{det} (A) \ne 0 \end{equation}\]

In this case the problem can be solved with Cramer's rule. In practice, however, using this rule for $N \geq 4$ is very complicated; for several other reasons this procedure is not appropriate for use on a computer.

As in other fields of numerical mathematics, also for the present case there are two classes of methods:

Direct methods:
These do not include any methodological error and therefore always return the exact solution - except for rounding errors. Indeed, rounding errors can have severe effects in direct methods, which employ algorithms that are computationally expensive.

As an example in the following we will illustrate the elimination method of Gauss, in the formulation of LU decomposition.
Iterative methods:
are often characterized by a particularly simple algorithm, stable with respect to rounding errors, and economical in terms of memory. However, there are also cases where iterative methods do not converge. Furthermore, also when they converge, they do not give the exact solution to the problem, but only an approximate one, which is also affected by a truncation error.

As an example of iterative method for a linear set of equations, in this chapter we will discuss the Gauss-Seidel method.

Aim of the direct methods: transformation of the matrix of coefficients into a triangular matrix

The matrix $A$ can be written as the product $A=LU$ where $L$ is a lower triangular matrix and $U$ is an upper triangular matrix of the form,

\[\begin{equation} \left[ \begin{array}{cccc} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \\ A_{41} & A_{42} & A_{43} & A_{44} \end{array} \right]= \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ L_{21} & 1 & 0 & 0 \\ L_{31} & L_{32} & 1 & 0 \\ L_{41} & L_{42} & L_{43} & 1 \end{array} \right]\left[ \begin{array}{cccc} U_{11} & U_{12} & U_{13} & U_{14} \\ 0 & U_{22} & U_{23} & U_{24} \\ 0 & 0 & U_{33} & U_{34} \\ 0 & 0 & 0 & U_{44} \end{array} \right]. \end{equation}\]

This is a convenient form for solving equations such as $A \cdot \vec{x} = \vec{b}$.

\[\begin{equation} A \vec{x} = LU\vec{x}=\vec{b}. \end{equation}\]

To find the solution we define $\vec{y} =U\vec{x}$. Then the equation $L\vec{y}=\vec{b}$ and easily be solved by letting $L$ act on $\vec{y}$,

\[\begin{equation} \left[ \begin{array}{y} y_1 \\ L_{21}y_1 + y_2 \\ L_{31}y_1 + L_{32}y_2 + y_3 \\ L_{41}y_1 + L_{42}y_2 + L_{43}y_3 +y_4 \end{array} \right]= \left[ \begin{array}{b} b_1 \\ b_2 \\ b_3 \\ b_4 \end{array} \right]. \end{equation}\]

The first component of the vector $\vec{y}$ is $y_1=b_1$ and the equations can then be solved from top to bottom by forward substitution,

\[\begin{equation} y_i = b_i-\sum \limits_{j=1}^{j < i}L_{ij}y_j. \end{equation}\]

Next let $U$ act on $\vec{x}$,

\[\begin{equation} \left[ \begin{array}{y} U_{11}x_1 + U_{12}x_2 + U_{13}x_3 + U_{14}x_4\\ U_{22}x_2 + U_{23}x_3 + U_{24}x_4\\ U_{33}x_3 + U_{34}x_4 \\ U_{44}x_4 \end{array} \right]= \left[ \begin{array}{b} y_1 \\ y_2 \\ y_3 \\ y_4 \end{array} \right]. \end{equation}\]

The last component of the vector $\vec{x}$ is $x_N=y_N/U_{NN}$ and the equations can be solved from bottom to top by back-substitution,

\[\begin{equation} x_i = \frac{1}{U_{ii}}\left(y_i-\sum \limits_{j=i+1}^{N}U_{ij}x_j\right)\qquad i=N-1, N-2,\cdots,1. \end{equation}\]

It is convenient to define an $LU$ matrix where the elements of $L_{ij}$ occupy the positions $j < i$ and the elements of $U_{ij}$ occupy the other positions.

\[\begin{equation} \mbox{LU-Matrix}: \left( \begin{array}{ccccc} U_{11} & U_{12} & U_{13} &....... & U_{1N} \\ L_{21} & U_{22} & U_{23} &....... & U_{2N} \\ L_{31} & L_{32} & U_{33} &....... & U_{3N} \\ . & & & & . \\ . & & & & . \\ . & & & & . \\ L_{N1} & L_{N2} & L_{N3} & ....... & U_{NN} \end{array} \right) \end{equation}\]

Given the $LU$ matrix for $A$, the code below solves $A\vec{x}=\vec{b}$ for $\vec{x}$ using forward- and back- substitution. A function LU_solve(LU,b) is defined that is passed an LU-matrix and a column vector $\vec{b}$ and returns the column vector $\vec{x}$.

<script> function LU_solve(LU,b) { //LU is a square NxN LU-matrix, the LU decomposition of matrix A //b is NxM matrix of M column vectors //LU_solve(LU,b) returns x, a NxM matrix of M column vectors //each column of x is a solution of Ax=b //initialize x and y to the same size as b x = new Array(b.length); y = new Array(b.length); for (i=0; i<b.length; i++) { x[i] = new Array(b[0].length); y[i] = new Array(b[0].length); } N = LU[0].length; Nb = b[0].length; //number of column vectors for (q=0; q<Nb; q++) { for (i=0; i<N; i++) { //forward substitution y[i][q]=b[i][q]; j = 0; while (j<i) { y[i][q] += - LU[i][j]y[j][q]; j++; } } for (i=N-1; i>-1; i--) { //back substitution x[i][q]=y[i][q]; j = i+1; while (j<N) { x[i][q] += - LU[i][j]x[j][q]; j++; } x[i][q]=x[i][q]/LU[i][i]; } } return x; } LU = [[2,7,6,2],[9,5,1,3],[4,3,8,4],[5,6,7,8]]; N = LU[0].length b = [[3],[7],[9],[3]]; x = LU_solve(LU,b); document.write("x = <br>"); for (i=0; i<N; i++) { for (j=0; j<x[0].length; j++) { document.write(x[i][j]+" "); } document.write("<\/br> "); } </script>

If you need to solve $A\vec{x}=\vec{b}$ for several vectors $\vec{b}$, these vectors can be passed to the function LU_solve(LU,b) where b is as an $N\times M$ matrix consisting of $M$ column vectors. Each column represents a different vector $\vec{b}$. Substitute b = [[3,4],[7,5],[9,1],[3,3]]; in the code and it will return two column vectors for $\vec{x}$.