Scientific Computing TU Berlin Winter 2021/22 © Jürgen Fuhrmann
Notebook 11

html"""

Scientific Computing TU Berlin Winter 2021/22 © Jürgen Fuhrmann

Notebook 11

"""

1.4 ms

begin

using LinearAlgebra

end

4.5 ms

TableOfContents(title="")

6.5 ms

Tridiagonal systems

In the previous lecture (nb10) 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$ .

md"""

# Tridiagonal systems

In the previous lecture (nb10) we introudced the discretization matrix for the 1D heat conduction problem. In general form it can be written as a tridiagonal matrix

$$A= \left(\begin{matrix}

{b_1} & {c_1} & { } & { } & { } \\

{a_2} & {b_2} & {c_2} & { } & { } \\

{ } & {a_3} & {b_3} & \ddots & { } \\

{ } & { } & \ddots & \ddots & {c_{n-1}}\\

{ } & { } & { } & {a_N} & {b_N}\\

\end{matrix}\right)$$

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

"""

7.1 μ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)

md"""

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)

"""

2.4 ms

Прогонка: 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}$

md"""

## Прогонка: derivation

Write solution of $A u=f$ as

$a_i u_{i-1} + b_i u_i + c_i u_{i+1} = f_i \quad (i=1\dots N)$

where we define $a_1=0$, $c_N=0$.

- For $i=1\dots N-1$, assume there are coefficients $\alpha_i, \beta_i$ such that

$u_i=\alpha_{i+1}u_{i+1}+\beta_{i+1}$.

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

$$(a_i\alpha_i\alpha_{i+1}+ b_i \alpha_{i+1} + c_i) u_{i+1} + (a_i\alpha_i\beta_{i+1}+ a_i\beta_i +b_i \beta_{i+1} -f_i)=0$$

- This is true for arbitrary $u$ if

$\begin{cases}

a_i\alpha_i\alpha_{i+1}+ b_i \alpha_{i+1} + c_i&=0\\

a_i\alpha_i\beta_{i+1}+ a_i\beta_i +b_i \beta_{i+1} -f_i &=0

\end{cases}$

- Re-arranging gives for ``i=1\dots N-1``:

$\begin{cases}

\alpha_{i+1}&= -\frac{c_i}{a_i\alpha_i+b_i}\\

\beta_{i+1} &= \frac{f_i - a_i\beta_i}{a_i\alpha_i+b_i}

\end{cases}$

"""

15.4 μ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}$

md"""

## Прогонка: realization

- Initialization of forward sweep:

$\begin{cases}

\alpha_{2}&= -\frac{c_1}{b_1}\\

\beta_{2} &= \frac{f_i}{b_1}

\end{cases}$

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

$\begin{cases}

\alpha_{i+1}&= -\frac{c_i}{a_i\alpha_i+b_i}\\

\beta_{i+1} &= \frac{f_i - a_i\beta_i}{a_i\alpha_i+b_i}

\end{cases}$

- Initialization of backward sweep: $u_N=\frac{f_N - a_N\beta_N}{a_N\alpha_N+b_N}$

- Backward sweep: for $i= N-1 \dots 1$:

$u_i=\alpha_{i+1}u_{i+1}+\beta_{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$ possible 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.

statement"""

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$ possible 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.

"""

165 ms

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.

md"""

## 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.

"""

3.9 μs

N=5

68 ns

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

md"""

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

"""

7.3 μs

5×5 Tridiagonal{Float64, Vector{Float64}}:
 0.422241  0.810134   ⋅         ⋅         ⋅ 
 0.695929  0.941789  0.690428   ⋅         ⋅ 
  ⋅        0.596871  0.190061  0.95558    ⋅ 
  ⋅         ⋅        0.290692  0.250277  0.263156
  ⋅         ⋅         ⋅        0.124175  0.220931

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

4.6 μs

LU{Float64, Tridiagonal{Float64, Vector{Float64}}}
L factor:
5×5 Matrix{Float64}:
 1.0       0.0        0.0       0.0       0.0
 0.0       1.0        0.0       0.0       0.0
 0.606731  0.399955   1.0       0.0       0.0
 0.0       0.0        0.0       1.0       0.0
 0.0       0.0       -0.587351  0.207751  1.0
U factor:
5×5 Matrix{Float64}:
 0.695929  0.941789   0.690428   0.0       0.0
 0.0       0.596871   0.190061   0.95558   0.0
 0.0       0.0       -0.49492   -0.382189  0.0
 0.0       0.0        0.0        0.124175  0.220931
 0.0       0.0        0.0        0.0       0.217258

lu(A)

3.5 μs

Int64

lu(A).p

2.7 μs

Float64

3.31506

-0.493445

-1.22

1.59735

3.62851

A\ones(N)

4.1 μs

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

md"""

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

"""

2.2 μ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)

md"""

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

"""

2.8 μs

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

md"""

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

"""

3.8 μs

5×5 Tridiagonal{Float64, Vector{Float64}}:
  2.99421    -0.713217    ⋅          ⋅          ⋅ 
 -0.0690588   2.51996   -0.316398    ⋅          ⋅ 
   ⋅         -0.635394   2.9221    -0.156421    ⋅ 
   ⋅           ⋅        -0.102088   2.8711    -0.592794
   ⋅           ⋅          ⋅        -0.107685   2.44815

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

39.6 ms

LU{Float64, Tridiagonal{Float64, Vector{Float64}}}
L factor:
5×5 Matrix{Float64}:
  1.0         0.0        0.0         0.0        0.0
 -0.0230641   1.0        0.0         0.0        0.0
  0.0        -0.253801   1.0         0.0        0.0
  0.0         0.0       -0.0359235   1.0        0.0
  0.0         0.0        0.0        -0.0375801  1.0
U factor:
5×5 Matrix{Float64}:
 2.99421  -0.713217   0.0        0.0        0.0
 0.0       2.50351   -0.316398   0.0        0.0
 0.0       0.0        2.8418    -0.156421   0.0
 0.0       0.0        0.0        2.86548   -0.592794
 0.0       0.0        0.0        0.0        2.42588

lu(A1)

2.4 μs

Int64

lu(A1).p

4.3 μs