Numerical Methods  

Numerical methods for linear inhomogeneous sets of equationsLet 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, leastsquare 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 nonsingular: \[\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:
Aim of the direct methods: transformation of the matrix of coefficients into a triangular matrixThe 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 backsubstitution, \[\begin{equation} x_i = \frac{1}{U_{ii}}\left(y_i\sum \limits_{j=i+1}^{N}U_{ij}x_j\right)\qquad i=N1, N2,\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{LUMatrix}: \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 LUmatrix and a column vector $\vec{b}$ and returns the column vector $\vec{x}$. 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}$. 