Solving Linear Equations

Serial Algorithm

To demonstrate the algorithm, we will consider a simple system of linear equations $Ax = b$:

A = [1.0 4.0 5.0 8.0 1.0; 
    2.0 -1.0 4.0 3.0 0.0; 
    7.0 6.0 3.0 -4.0 5.0; 
    -3.0 4.0 2.0 2.0 2.0; 
    0.0 -4.0 2.0 1.0 2.0]

b = [61.0; 24.0; 37.0; 29.0; 12.0];

The code in the following cell converts the general problem $Ax=b$ to the upper triangular equation system $Ux=y$. Note that this function assumes that the pivots are all nonzero. This function will be erroneos if any of the diagonal entries are zero!

function convert_to_upper_triangular!(A,b)
    n = size(A,1)
    # Upper Triangularization: convert Ax=b to Ux=y
    for k in 1:n
        for j in k+1:n
            # Divide by pivot
            A[k,j] = A[k,j] / A[k,k]
        end
        b[k] = b[k] / A[k,k]
        A[k,k] = 1
        # Substract lower rows
        for i in k+1:n 
            for j in k+1:n
                A[i,j]=A[i,j] - A[i,k] * A[k,j]
            end
            b[i] = b[i] - A[i,k] * b[k]
            A[i,k] = 0
        end
    end
    return A, b #U,y
end

convert_to_upper_triangular! (generic function with 1 method)

The function in the following cell solves the upper triangular equation system using backwards substitution. Note that the function alters the input values.

function solve_upper_triangular!(U,y)
    n = size(U,1)
    for step in reverse(1:n)
        if U[step,step] == 0
            if y[step] != 0
                return "No solution"
            else
                return "Infinity solutions"
            end
        else
        # Backwards substitution
            y[step] = y[step] / U[step,step]
        end
        for row in reverse(1:step-1)
            y[row] -= U[row,step] * y[step]
        end
    end
    return y 
end

solve_upper_triangular! (generic function with 1 method)

U,y = convert_to_upper_triangular!(A,b)
sol = solve_upper_triangular!(U,y)

5-element Vector{Float64}:
 1.0000000000000009
 1.999999999999999
 2.9999999999999964
 4.000000000000002
 5.000000000000005

We can test if the obtained solution is correct using @test:

using Test
@test sol ≈ [1.0; 2.0; 3.0; 4.0; 5.0]

Test Passed