⚡ Pluto.jl ⚡

xxxxxxxxxx
 
begin 
    using Pkg
    Pkg.activate(mktempdir())
    Pkg.add("PyPlot")
    Pkg.add("PlutoUI")
    using PlutoUI
    using PyPlot
    using LinearAlgebra
end

24.2 s

Tridiagonal systems

In the previous lecture (nb08) we introudced the discretization matrix for the 1D heat conduction problem. In general form it can be written as a tridiagonal matrix $A = (\begin{matrix} b_{1} & c_{1} \\ a_{2} & b_{2} & c_{2} \\ a_{3} & b_{3} & ⋱ \\ ⋱ & ⋱ & c_{n - 1} \\ a_{N} & b_{N} \end{matrix})$

and stored in three arrays $a$ , $b$ , $c$ .

6.6 μs

Gaussian elimination using arrays $a, b, c$ as matrix storage ?

From what we have seen, this question arises in a quite natural way, and historically, the answer has been given several times and named differently
TDMA (tridiagonal matrix algorithm)
"Thomas algorithm" (Llewellyn H. Thomas, 1949 (?))
"Progonka method" (from Russian "прогонка": "run through"; Gelfand, Lokutsievski, 1952, published 1960)

6.9 μs

Прогонка: derivation

Write solution of $A u = f$ as

$a_{i} u_{i - 1} + b_{i} u_{i} + c_{i} u_{i + 1} = f_{i} (i = 1 \dots N)$

where we define $a_{1} = 0$ , $c_{N} = 0$ .

For $i = 1 \dots N - 1$ , assume there are coefficients $α_{i}, β_{i}$ such that

$u_{i} = α_{i + 1} u_{i + 1} + β_{i + 1}$

Re-arranging, we can express $u_{i - 1}$ and $u_{i}$ via $u_{i + 1}$ :

$(a_{i} α_{i} α_{i + 1} + b_{i} α_{i + 1} + c_{i}) u_{i + 1} + (a_{i} α_{i} β_{i + 1} + a_{i} β_{i} + b_{i} β_{i + 1} - f_{i}) = 0$

This is true for arbitrary $u$ if ${\begin{cases} a_{i} α_{i} α_{i + 1} + b_{i} α_{i + 1} + c_{i} & = 0 \\ a_{i} α_{i} β_{i + 1} + a_{i} β_{i} + b_{i} β_{i + 1} - f_{i} & = 0 \end{cases}$
Re-arranging gives for $i = 1 \dots N - 1$ :

${\begin{cases} α_{i + 1} & = - \frac{c_{i}}{a_{i} α_{i} + b_{i}} \\ β_{i + 1} & = \frac{f_{i} - a_{i} β_{i}}{a_{i} α_{i} + b_{i}} \end{cases}$

15.7 μs

Прогонка: realization

Initialization of forward sweep:

${\begin{cases} α_{2} & = - \frac{c_{1}}{b_{1}} \\ β_{2} & = \frac{f_{i}}{b_{1}} \end{cases}$

Forward sweep: for $i = 2 \dots N - 1$ :

${\begin{cases} α_{i + 1} & = - \frac{c_{i}}{a_{i} α_{i} + b_{i}} \\ β_{i + 1} & = \frac{f_{i} - a_{i} β_{i}}{a_{i} α_{i} + b_{i}} \end{cases}$

Initialization of backward sweep: $u_{N} = \frac{f_{N} - a_{N} β_{N}}{a_{N} α_{N} + b_{N}}$
Backward sweep: for $i = N - 1 \dots 1$ :

$u_{i} = α_{i + 1} u_{i + 1} + β_{i + 1}$

12.6 μs

Прогонка: properties

$N$ unknowns, one forward sweep, one backward sweep $\Rightarrow$ $O (N)$ operations vs. $O (N^{2.75})$ for algorithm using full matrix
No pivoting $\Rightarrow$ stability issues
Stability for diagonally dominant matrices where $| b_{i} | > | a_{i} | + | c_{i} |$
Stability for symmetric positive definite matrices
In fact, this is a realization of Gaussian elimination on a particular data structure.

10.3 μs

Tridiagonal matrices in Julia

In Julia, solution of a tridiagonal system is based on the LU factorization in the LAPACK routine dgtsv which also does pivoting.

4.0 μs

xxxxxxxxxx
 
N=5

994 ns

LU Factorization in the case of a tridiagonal matrix with random diagonal entries

3.9 μs

5×5 Tridiagonal{Float64,Array{Float64,1}}:
 0.186793  0.985412   ⋅         ⋅         ⋅ 
 0.201945  0.336733  0.671338   ⋅         ⋅ 
  ⋅        0.518801  0.191529  0.416203   ⋅ 
  ⋅         ⋅        0.889957  0.147218  0.769763
  ⋅         ⋅         ⋅        0.726546  0.470719

xxxxxxxxxx
 
A=Tridiagonal(rand(N-1),rand(N),rand(N-1))

138 ms

LU{Float64,Tridiagonal{Float64,Array{Float64,1}}}
L factor:
5×5 Array{Float64,2}:
 1.0       0.0       0.0       0.0       0.0
 0.924968  1.0       0.0       0.0       0.0
 0.0       0.0       1.0       0.0       0.0
 0.0       0.0       0.0       1.0       0.0
 0.0       0.769797  0.752337  0.420408  1.0
U factor:
5×5 Array{Float64,2}:
 0.201945  0.336733   0.671338  0.0        0.0
 0.0       0.673945  -0.620967  0.0        0.0
 0.0       0.0        0.889957  0.147218   0.769763
 0.0       0.0        0.0       0.726546   0.470719
 0.0       0.0        0.0       0.0       -0.777014

xxxxxxxxxx
 
lu(A)

50.1 ms

Int641

xxxxxxxxxx
 
lu(A).p

46.2 μs

Float641

1.5038

0.729746

0.671173

1.18418

0.296653

xxxxxxxxxx
 
A\ones(N)

329 ms

Solving this system with a positive right hand side can yield negative solution components.

2.6 μs

We see that the in order to maintain stability, pivoting is performed: the LU factorization is performed as $P A = L U$ where $P$ is a permutation matrix. The underlying permutation can be obtained as lu(A).p)

2.6 μs

Define a diagonally dominant matrix with random entries with positive main diagonal and nonpositive off-diagonal elements:

3.4 μs

5×5 Tridiagonal{Float64,Array{Float64,1}}:
  2.04343   -0.907038    ⋅          ⋅          ⋅ 
 -0.265936   2.25515   -0.597263    ⋅          ⋅ 
   ⋅        -0.739934   2.24558   -0.790491    ⋅ 
   ⋅          ⋅        -0.701657   2.47384   -0.233202
   ⋅          ⋅          ⋅        -0.662966   2.22106

xxxxxxxxxx
 
A1=Tridiagonal(-rand(N-1),rand(N).+2,-rand(N-1))

69.5 ms

LU{Float64,Tridiagonal{Float64,Array{Float64,1}}}
L factor:
5×5 Array{Float64,2}:
  1.0        0.0        0.0        0.0       0.0
 -0.130142   1.0        0.0        0.0       0.0
  0.0       -0.346232   1.0        0.0       0.0
  0.0        0.0       -0.344154   1.0       0.0
  0.0        0.0        0.0       -0.301103  1.0
U factor:
5×5 Array{Float64,2}:
 2.04343  -0.907038   0.0        0.0        0.0
 0.0       2.1371    -0.597263   0.0        0.0
 0.0       0.0        2.03879   -0.790491   0.0
 0.0       0.0        0.0        2.20179   -0.233202
 0.0       0.0        0.0        0.0        2.15084

xxxxxxxxxx
 
lu(A1)

3.6 μs

Int641

xxxxxxxxxx
 
lu(A1).p

4.7 μs

Here we see, that no permutation is needed to maintain stability, confirming the statement made. In this case, the underlying algorithm is equivalent to Progonka, and the resulting LU factorization can be stored in three diagonals.

2.4 μs

Float641

0.844489

0.800026

0.97042

0.742815

0.671959

xxxxxxxxxx
 
A1\ones(N)

5.6 μs

Here we get only nonnegative solution values, though the matrix off-diagonal elements are nonpositive. Later we will see that this is a theorem for this type of matrices.

2.4 μs

5×5 Array{Float64,2}:
 0.519495    0.231457   0.0686099  0.0225583  0.00236853
 0.0678612   0.52144    0.154568   0.0508208  0.00533598
 0.024921    0.191491   0.55307    0.181845   0.019093
 0.00727301  0.0558852  0.16141    0.469004   0.0492435
 0.00217093  0.0166812  0.0481793  0.139993   0.464934

xxxxxxxxxx
 
inv(A1)

41.9 ms

The inverse is a nonnegative full matrix! This is a theorem as well.

2.4 μs