6.2. Solving systems of linear equations using LU decomposition#

Systems of linear equations often appear in the topics of numerical analysis and numerical solutions to differential equations. Examples include the solution of the stage values of an implicit Runge-Kutta method and the solution to a boundary value problem using the finite-difference method. The methods that are applied to solve systems of linear equations fall into one of two categories: direct methods that use an algebraic approach and indirect methods that use an iterative approach. Here we will look at some common direct methods.

6.2.1. Systems of linear equations#

A system of linear equations with \(m\) equations and \(n\) unknowns is expressed as

\[\begin{split} \begin{align*} 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 ,\\ &\vdots \\ a_{m1} x_1 +a_{m2} x_2 +\cdots +a_{mn} x_n &=b_n , \end{align*} \end{split}\]

where \(x_i\) are the unknown variables, \(a_{ij}\) are coefficients and \(b_i\) are constant terms. It is often more convenient to express a system of linear equations as a matrix equation

\[ \begin{align*} A \vec{x} = \vec{b}, \end{align*} \]

where \(A\) is the coefficient matrix, \(\vec{x}\) is the variable vector and \(b\) is the constant vector

\[\begin{split} \begin{align*} A &= \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{pmatrix},& \vec{x} &= \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_m \end{pmatrix}, & \vec{b} &= \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{pmatrix}. \end{align*} \end{split}\]

The solution (if it exists) is the vector \(\vec{x}\) for which the equation \(A\vec{x} = \vec{b}\) is satisfied. The solution to a linear system may be one of the following

  • The system has infinitely many solutions. This usually occurs when the number of unknown variables exceeds the number of equations in the system.

  • The system has a single unique solution.

  • The system has no solution. This usually occurs where one equation in the system contradicts another such that not value of the variables could satisfy both.

6.2.2. Solving systems of linear equations using LU decomposition#

Given a system of linear equations of the form \(A \vec{x} = \vec{b}\) where LU decomposition has been applied to the coefficient matrix then since \(A = LU\) we have

\[ LU \vec{x} = \vec{b}.\]

Let \(\vec{y} = U \vec{x}\) then

\[ L \vec{y} = \vec{b}. \]

\(L\) is a lower triangular matrix so the solution of \(L \vec{y} = \vec{b}\) is easily calculated using forward substitution. Once \(\vec{y}\) has been calculated the solution to \(U\vec{x} = \vec{y}\) is calculated using back substitution.

The advantage of using LU decomposition to solve systems of linear equations is that once the LU decomposition of the coefficient matrix has been calculated it can be used for any values of the constant vector \(\vec{b}\). This is unlike Gaussian elimination where row reduction will need to be repeated if the values of \(\vec{b}\) change.

Algorithm 6.2 (Solving systems of linear equations using LU decomposition)

Inputs: A system of linear equations expressed using the coefficient matrix \(A\) and variable vector \(\vec{b}\) such that \(A \vec{x} = \vec{b}\).

Inputs: The variable vector \(\vec{x}\).

  • Calculate the LU decomposition of \(A\) to determine \(L\) and \(U\) such that \(A = LU\).

  • For \(i = 1, \ldots, n\) do

    • \(y_i \gets b_i - \displaystyle\sum_{j=1}^i \ell_{ij} y_j\)

  • \(x_n \gets \dfrac{y_i}{u_{nn}}\)

  • For \(i = n - 1, \ldots, 1\) do

    • \(x_i \gets \dfrac{1}{u_{ii}} \left( y_i - \displaystyle \sum_{j=i+1}^{n} u_{ij} x_j \right)\)

  • Return \(\vec{x} = (x_1, x_2, \ldots, x_n)\).

Example 6.2

Solve the following system using LU decomposition

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 3 & 0 \\ 2 & -4 & -1 \\ -3 & 1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}= \begin{pmatrix} -7 \\ 11 \\ 1 \end{pmatrix}. \end{align*} \end{split}\]
Solution (click to show)

We saw in the Example 6.1 above that the LU decomposition of the coefficient matrix is

\[\begin{split} \begin{align*} L & =\begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -3 & -1 & 1 \end{pmatrix}, & U &= \begin{pmatrix} 1 & 3 & 0 \\ 0 & -10 & -1 \\ 0 & 0 & 1 \end{pmatrix}. \end{align*} \end{split}\]

Solving \(L \vec{y} = \vec{b}\) using forward substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -3 & -1 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix} = \begin{pmatrix} -7\\ 11\\ 1 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} y_1 &=-7,\\ y_2 &=11-2y_1 =-2(-7)=25,\\ y_3 &=-1+3y_1 +y_2 =-1+3(-7)+1(25)=5. \end{align*} \end{split}\]

Solving \(U \vec{x} = \vec{y}\) using back substitution

\[\begin{split} \begin{align*} \begin{pmatrix} 1 & 3 & 0\\ 0 & -10 & -1\\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} -7 \\ 25 \\ 5 \end{pmatrix}, \end{align*} \end{split}\]

gives

\[\begin{split} \begin{align*} x_3 &= 5,\\ x_2 &= \frac{1}{-10}(25 + x_3) = -\frac{1}{10}(25 + 5) = -3,\\ x_1 &= -7 - 0x_3 - 3x_2 = -7 + 9 = 2. \end{align*} \end{split}\]

So the solution is \(\vec{x}=(2,-3,5)\).

6.2.3. Code#

The code below defines two functions called forward_substitution() and back_substitution() that calculate perform forward and back substitution. These are used to solve the system of linear equations from Example 6.2.

def forward_substitution(L, b):
    n = L.shape[0]
    x = np.zeros(n)
    for i in range(n):
        sum_ = 0
        for j in range(i):
            sum_ += L[i,j] * x[j]
            
        x[i] = (b[i] - sum_) / L[i,i]
    
    return x


def back_substitution(U, b):
    n = U.shape[0]
    x = np.zeros(n)
    x[-1] = b[-1] / U[-1,-1]
    for i in range(n - 2, -1, -1):
        sum_ = 0
        for j in range(i + 1, n):
            sum_ += U[i,j] * x[j]
            
        x[i] = (b[i] - sum_) / U[i,i]
        
    return x


# Define system
A = np.array([[1, 3, 0], 
              [2, -4, -1], 
              [-3, 1, 2]])
b = np.array([-7, 11, 1])

# Calculate LU decomposition
L, U = lu(A)

# Solve system
y = forward_substitution(L, b)
x = back_substitution(U, y)
print(f"x = {x}")

% Define system
A = [1, 3, 0 ;
     2, -4, -1 ;
     -3, 1, 2 ];
b = [-7 ; 11 ; 1];

% Calculate LUP decomposition
P = partial_pivot(A);
[L, U] = lu(P * A);

% Solve system
y = forward_substitution(L, P * b);
x = back_substitution(U, y)

% --------------------------------------------------------------
function x = forward_substitution(L, b)

n = size(L, 1);
x = zeros(n, 1);
for i = 1 : n
    sum_ = 0;
    for j = 1 : i - 1
        sum_ = sum_ + L(i,j) * x(j);
    end
    x(i) = (b(i) - sum_) / L(i,i);
end

end

% --------------------------------------------------------------
function x = back_substitution(U, b)

n = size(U, 1);
x = zeros(n, 1);
x(n) = b(n) / U(n,n);
for i = n - 1 : -1 : 1
    sum_ = 0
    for j = i + 1 : n
        sum_ = sum_ + U(i,j) * x(j);
    end
    x(i) = (b(i) - sum_) / U(i,i);
end

end