User:HerrHartmuth/sandbox

Similarity to symmetric tridiagonal matrix

Given a given real tridiagonal, unsymmetic matrix

$T={\begin{pmatrix}a_{1}&b_{1}\\c_{1}&a_{2}&b_{2}\\&c_{2}&\ddots &\ddots \\&&\ddots &\ddots &b_{n-1}\\&&&c_{n-1}&a_{n}\end{pmatrix}}$

where $b_{i}\neq c_{i}$.

Assume that the product of off-diagonal entries is strictly positive $b_{i}c_{i}>0$ and define a transformation matrix $D$ by

$D:=\operatorname {diag} (\delta _{1},\dots ,\delta _{n})\quad {\text{for}}\quad \delta _{i}:={\begin{cases}1&,\,i=1\\{\sqrt {\frac {c_{i-1}\dots c_{1}}{b_{i-1}\dots b_{1}}}}&,\,i=2,\dots ,n\,.\end{cases}}$

The similarity transformation $J:=D^{-1}TD$ yields a symmetric tridiagonal matrix $J$ by

$J:=D^{-1}TD={\begin{pmatrix}a_{1}&{\sqrt {b_{1}c_{1}}}\\{\sqrt {b_{1}c_{1}}}&a_{2}&{\sqrt {b_{2}c_{2}}}\\&c_{2}&\ddots &\ddots \\&&\ddots &\ddots &{\sqrt {b_{n-1}c_{n-1}}}\\&&&{\sqrt {b_{n-1}c_{n-1}}}&a_{n}\end{pmatrix}}\,.$

Note that $T$ and $J$ have the same eigenvalues.

Special Case: Real Tridiagonal

In the case of a tridiagonal structure with real elements the eigenvalues and eigenvectors can be derived explicity as

${\begin{aligned}\lambda _{k}&=a_{0}+2{\sqrt {a_{1}a_{-1}}}\cos \left({\frac {\pi k}{n+1}}\right)\\v^{k}&=\left(\left({\frac {a_{1}}{a_{-1}}}\right)^{1/2}\sin \left({\frac {1\pi k}{n+1}}\right),\ldots ,\left({\frac {a_{1}}{a_{-1}}}\right)^{n/2}\sin \left({\frac {n\pi k}{n+1}}\right)\right)^{T}\,.\end{aligned}}$

Legendre

Pointwise Evaluations

As shown before the values at the boundary are given by

$P_{n}(1)=1\,,\quad P_{n}(-1)={\begin{cases}1&{\text{for}}\quad n=2m\\-1&{\text{for}}\quad n=2m+1\,.\end{cases}}$

One can show that for $x=0$ the values are given by

$P_{n}(0)={\begin{cases}{\frac {(-1)^{m}}{4^{m}}}{\tbinom {2m}{m}}&{\text{for}}\quad n=2m\\0&{\text{for}}\quad n=2m+1\,.\end{cases}}$

$3$

$1$

Carrier Gen + Recomb

Radiative recombination

During radiative recombination, a form of spontaneous emission, a photon is emitted with the wavelength corresponding to the energy released. This effect is the basis of LEDs. Because the photon carries relatively little momentum, radiative recombination is significant only in direct bandgap materials.

When photons are present in the material, they can either be absorbed, generating a pair of free carriers, or they can stimulate a recombination event, resulting in a generated photon with similar properties to the one responsible for the event. Absorption is the active process in photodiodes, solar cells, and other semiconductor photodetectors, while stimulated emission is responsible for laser action in laser diodes.

In thermal equilibrium the radiative recombination $R_{0}$ and thermal generation rate $G_{0}$ equal each other

$R_{0}=G_{0}=B_{r}n_{0}p_{0}=B_{r}n_{i}^{2}$

where $B_{r}$ is called the radiative capture probability and $n_{i}$ the intrinsic carrier density.

Under steady-state conditions the radiative recombination rate $r$ and resulting net recombination rate $U_{r}$ are

$r=B_{r}np\,,\quad U_{r}=r-G_{0}=B_{r}\left(np-n_{i}^{2}\right)$

where the carrier densities $n,p$ are made up of equilibrium $n_{0},p_{0}$ and excess densities $\Delta n,\Delta p$

$n=n_{0}+\Delta n\,,\quad p=p_{0}+\Delta p\,.$

The radiative lifetime $\tau _{r}$ is given by

$\tau _{r}={\frac {\Delta n}{U_{r}}}={\frac {1}{B_{r}\left(n_{0}+p_{0}+\Delta n\right)}}\,.$

Auger recombination

In Auger recombination the energy is given to a third carrier, which is excited to a higher energy level without moving to another energy band. After the interaction, the third carrier normally loses its excess energy to thermal vibrations. Since this process is a three-particle interaction, it is normally only significant in non-equilibrium conditions when the carrier density is very high. The Auger effect process is not easily produced, because the third particle would have to begin the process in the unstable high-energy state.

The Auger recombination can be calculated from the equation :

$U_{Aug}=\Gamma _{n}\,n(np-n_{i}^{2})+\Gamma _{p}\,p(np-n_{i}^{2})$

In thermal equilibrium the Auger recombination $R_{A}$ and thermal generation rate $G_{0}$ equal each other

$R_{A}=G_{0}=C_{n}n_{0}^{2}p_{0}+C_{p}n_{0}p_{0}^{2}$

where $C_{n},C_{p}$ are the Auger capture probabilities.

The non-equilibrium Auger recombination rate $r_{A}$ and resulting net recombination rate $U_{A}$ under steady-state conditions are

$r_{A}=C_{n}n^{2}p+C_{p}np^{2}\,,\quad U_{A}=r_{A}-G_{0}=C_{n}\left(n^{2}p-n_{0}^{2}p_{0}\right)+C_{p}\left(np^{2}-n_{0}p_{0}^{2}\right)\,.$

The Auger lifetime $\tau _{A}$ is given by

$\tau _{A}={\frac {\Delta n}{U_{A}}}={\frac {1}{n^{2}C_{n}+2n_{i}^{2}(C_{n}+C_{p})+p^{2}C_{p}}}\,.$

Auger recombination in LEDs

The mechanism causing LED efficiency droop was identified in 2007 as Auger recombination, which met with a mixed reaction. In 2013, an experimental study claimed to have identified Auger recombination as the cause of efficiency droop. However, it remains disputed whether the amount of Auger loss found in this study is sufficient to explain the droop. Other frequently quoted evidence against Auger as the main droop causing mechanism is the low-temperature dependence of this mechanism which is opposite to that found for the drop.

${\displaystyle }$

${\displaystyle }$

MINRES

In mathematics, the minimal residual method (MINRES) is an iterative method for the numerical solution of a symmetric but possibly indefinite system of linear equations. The method approximates the solution by the vector in a Krylov subspace with minimal residual. The Lanczos algorithm is used to find this vector.

Introduction

One tries to solve the following square system of linear equations

$Ax=b$

where $x\in \mathbb {R} ^{n}$ is unknown and $A\in \mathbb {R} ^{n\times n}\,,b\in \mathbb {R} ^{n}$ are given.

In the special case of $A$ being symmetric and positive-definite one can use the Conjugate gradient method. For symmetric and possibly indefinite matrices one uses the MINRES method. In the case of unsymmetric and indefinite matrices one needs to fall back to methods such as the GMRES, or Bi-CG. ${\displaystyle }$

The method

Krylov space basis

The matrix $A$ is symmetric and thus one can apply the Lanczos method to find an orthogonal basis for the Krylov subspace ${\displaystyle }$

Denote the Euclidean norm of any vector v by $\|v\|$. Denote the (square) system of linear equations to be solved by

$Ax=b.\,$

The matrix A is assumed to be invertible of size m-by-m. Furthermore, it is assumed that b is normalized, i.e., that $\|b\|=1$.

The n-th Krylov subspace for this problem is

$K_{n}=K_{n}(A,b)=\operatorname {span} \,\{b,Ab,A^{2}b,\ldots ,A^{n-1}b\}.\,$

GMRES approximates the exact solution of $Ax=b$ by the vector $x_{n}\in K_{n}$ that minimizes the Euclidean norm of the residual $r_{n}=Ax_{n}-b$.

The vectors $b,Ab,\ldots A^{n-1}b$ might be close to linearly dependent, so instead of this basis, the Arnoldi iteration is used to find orthonormal vectors $q_{1},q_{2},\ldots ,q_{n}\,$ which form a basis for $K_{n}$. Hence, the vector $x_{n}\in K_{n}$ can be written as $x_{n}=Q_{n}y_{n}$ with $y_{n}\in \mathbb {R} ^{n}$, where $Q_{n}$ is the m-by-n matrix formed by $q_{1},\ldots ,q_{n}$.

The Arnoldi process also produces an ($n+1$)-by-$n$ upper Hessenberg matrix ${\tilde {H}}_{n}$ with

$AQ_{n}=Q_{n+1}{\tilde {H}}_{n}.\,$

Because columns of $Q_{n}$ are orthogonal, we have

$\|Ax_{n}-b\|=\|{\tilde {H}}_{n}y_{n}-Q_{n+1}^{T}b\|=\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|,\,$

where

$e_{1}=(1,0,0,\ldots ,0)^{T}\,$

is the first vector in the standard basis of $\mathbb {R} ^{n+1}$, and

$\beta =\|b-Ax_{0}\|\,,$

$x_{0}$ being the first trial vector (usually zero). Hence, $x_{n}$ can be found by minimizing the Euclidean norm of the residual

$r_{n}={\tilde {H}}_{n}y_{n}-\beta e_{1}.$

This is a linear least squares problem of size n.

This yields the GMRES method. On the $n$-th iteration:

1. calculate $q_{n}$ with the Arnoldi method;
2. find the $y_{n}$ which minimizes $\|r_{n}\|$;
3. compute $x_{n}=Q_{n}y_{n}$;
4. repeat if the residual is not yet small enough.

At every iteration, a matrix-vector product $Aq_{n}$ must be computed. This costs about $2m^{2}$ floating-point operations for general dense matrices of size $m$, but the cost can decrease to $O(m)$ for sparse matrices. In addition to the matrix-vector product, $O(nm)$ floating-point operations must be computed at the n -th iteration.

Convergence

The nth iterate minimizes the residual in the Krylov subspace K_n. Since every subspace is contained in the next subspace, the residual does not increase. After m iterations, where m is the size of the matrix A, the Krylov space K_m is the whole of R^m and hence the GMRES method arrives at the exact solution. However, the idea is that after a small number of iterations (relative to m), the vector x_n is already a good approximation to the exact solution.

This does not happen in general. Indeed, a theorem of Greenbaum, Pták and Strakoš states that for every nonincreasing sequence a₁, …, a_m−1, a_m = 0, one can find a matrix A such that the ||r_n|| = a_n for all n, where r_n is the residual defined above. In particular, it is possible to find a matrix for which the residual stays constant for m − 1 iterations, and only drops to zero at the last iteration.

In practice, though, GMRES often performs well. This can be proven in specific situations. If the symmetric part of A, that is $(A^{T}+A)/2$, is positive definite, then

$\|r_{n}\|\leq \left(1-{\frac {\lambda _{\min }^{2}(1/2(A^{T}+A))}{\lambda _{\max }(A^{T}A)}}\right)^{n/2}\|r_{0}\|,$

where $\lambda _{\mathrm {min} }(M)$ and $\lambda _{\mathrm {max} }(M)$ denote the smallest and largest eigenvalue of the matrix $M$, respectively.

If A is symmetric and positive definite, then we even have

$\|r_{n}\|\leq \left({\frac {\kappa _{2}(A)^{2}-1}{\kappa _{2}(A)^{2}}}\right)^{n/2}\|r_{0}\|.$

where $\kappa _{2}(A)$ denotes the condition number of A in the Euclidean norm.

In the general case, where A is not positive definite, we have

${\frac {\|r_{n}\|}{\|b\|}}\leq \inf _{p\in P_{n}}\|p(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{n}}\max _{\lambda \in \sigma (A)}|p(\lambda )|,\,$

where P_n denotes the set of polynomials of degree at most n with p(0) = 1, V is the matrix appearing in the spectral decomposition of A, and σ(A) is the spectrum of A. Roughly speaking, this says that fast convergence occurs when the eigenvalues of A are clustered away from the origin and A is not too far from normality.

All these inequalities bound only the residuals instead of the actual error, that is, the distance between the current iterate x_n and the exact solution.

Extensions of the method

Like other iterative methods, GMRES is usually combined with a preconditioning method in order to speed up convergence.

The cost of the iterations grow as O(n²), where n is the iteration number. Therefore, the method is sometimes restarted after a number, say k, of iterations, with x_k as initial guess. The resulting method is called GMRES(k) or Restarted GMRES. This methods suffers from stagnation in convergence as the restarted subspace is often close to the earlier subspace.

The shortcomings of GMRES and restarted GMRES are addressed by the recycling of Krylov subspace in the GCRO type methods such as GCROT and GCRODR. Recycling of Krylov subspaces in GMRES can also speed up convergence when sequences of linear systems need to be solved.

Comparison with other solvers

The Arnoldi iteration reduces to the Lanczos iteration for symmetric matrices. The corresponding Krylov subspace method is the minimal residual method (MinRes) of Paige and Saunders. Unlike the unsymmetric case, the MinRes method is given by a three-term recurrence relation. It can be shown that there is no Krylov subspace method for general matrices, which is given by a short recurrence relation and yet minimizes the norms of the residuals, as GMRES does.

Another class of methods builds on the unsymmetric Lanczos iteration, in particular the BiCG method. These use a three-term recurrence relation, but they do not attain the minimum residual, and hence the residual does not decrease monotonically for these methods. Convergence is not even guaranteed.

The third class is formed by methods like CGS and BiCGSTAB. These also work with a three-term recurrence relation (hence, without optimality) and they can even terminate prematurely without achieving convergence. The idea behind these methods is to choose the generating polynomials of the iteration sequence suitably.

None of these three classes is the best for all matrices; there are always examples in which one class outperforms the other. Therefore, multiple solvers are tried in practice to see which one is the best for a given problem.

Solving the least squares problem

One part of the GMRES method is to find the vector $y_{n}$ which minimizes

$\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|.\,$

Note that ${\tilde {H}}_{n}$ is an (n + 1)-by-n matrix, hence it gives an over-constrained linear system of n+1 equations for n unknowns.

The minimum can be computed using a QR decomposition: find an (n + 1)-by-(n + 1) orthogonal matrix Ω_n and an (n + 1)-by-n upper triangular matrix ${\tilde {R}}_{n}$ such that

$\Omega _{n}{\tilde {H}}_{n}={\tilde {R}}_{n}.$

The triangular matrix has one more row than it has columns, so its bottom row consists of zero. Hence, it can be decomposed as

${\tilde {R}}_{n}={\begin{bmatrix}R_{n}\\0\end{bmatrix}},$

where $R_{n}$ is an n-by-n (thus square) triangular matrix.

The QR decomposition can be updated cheaply from one iteration to the next, because the Hessenberg matrices differ only by a row of zeros and a column:

${\tilde {H}}_{n+1}={\begin{bmatrix}{\tilde {H}}_{n}&h_{n+1}\\0&h_{n+2,n+1}\end{bmatrix}},$

where h_n+1 = (h_1,n+1, …, h_n+1,n+1)^T. This implies that premultiplying the Hessenberg matrix with Ω_n, augmented with zeroes and a row with multiplicative identity, yields almost a triangular matrix:

${\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&\rho \\0&\sigma \end{bmatrix}}$

This would be triangular if σ is zero. To remedy this, one needs the Givens rotation

$G_{n}={\begin{bmatrix}I_{n}&0&0\\0&c_{n}&s_{n}\\0&-s_{n}&c_{n}\end{bmatrix}}$

where

$c_{n}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{n}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}.$

With this Givens rotation, we form

$\Omega _{n+1}=G_{n}{\begin{bmatrix}\Omega _{n}&0\\0&1\end{bmatrix}}.$

Indeed,

$\Omega _{n+1}{\tilde {H}}_{n+1}={\begin{bmatrix}R_{n}&r_{n+1}\\0&r_{n+1,n+1}\\0&0\end{bmatrix}}\quad {\text{with}}\quad r_{n+1,n+1}={\sqrt {\rho ^{2}+\sigma ^{2}}}$

is a triangular matrix.

Given the QR decomposition, the minimization problem is easily solved by noting that

$\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|=\|\Omega _{n}({\tilde {H}}_{n}y_{n}-\beta e_{1})\|=\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\|.$

Denoting the vector $\beta \Omega _{n}e_{1}$ by

${\tilde {g}}_{n}={\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}$

with g_n ∈ Rⁿ and γ_n ∈ R, this is

$\|{\tilde {H}}_{n}y_{n}-\beta e_{1}\|=\|{\tilde {R}}_{n}y_{n}-\beta \Omega _{n}e_{1}\|=\left\|{\begin{bmatrix}R_{n}\\0\end{bmatrix}}y_{n}-{\begin{bmatrix}g_{n}\\\gamma _{n}\end{bmatrix}}\right\|.$

The vector y that minimizes this expression is given by

$y_{n}=R_{n}^{-1}g_{n}.$

Again, the vectors $g_{n}$ are easy to update.

Example code

Regular GMRES (MATLAB / GNU Octave)

function [x, e] = gmres( A, b, x, max_iterations, threshold)
  n = length(A);
  m = max_iterations;
  
  %use x as the initial vector
  r=b-A*x;

  b_norm = norm(b);
  error = norm(r)/b_norm;

  %initialize the 1D vectors
  sn = zeros(m,1);
  cs = zeros(m,1);
  e1 = zeros(n,1);
  e1(1) = 1;
  e=[error];
  r_norm=norm(r);
  Q(:,1) = r/r_norm;
  beta = r_norm*e1;
  for k = 1:m                                   
    
    %run arnoldi
    [H(1:k+1,k) Q(:,k+1)] = arnoldi(A, Q, k);
    
    %eliminate the last element in H ith row and update the rotation matrix
    [H(1:k+1,k) cs(k) sn(k)] = apply_givens_rotation(H(1:k+1,k), cs, sn, k);
    
    %update the residual vector
    beta(k+1) = -sn(k)*beta(k);
    beta(k)   = cs(k)*beta(k);
    error  = abs(beta(k+1)) / b_norm;
    
    %save the error
    e=[e; error];
    
    if ( error <= threshold)
      break;
    end
  end

  %calculate the result
  y = H(1:k,1:k) \ beta(1:k);
  x = x + Q(:,1:k)*y; 
end

%----------------------------------------------------%
%                  Arnoldi Function                  %
%----------------------------------------------------%
function [h, q] = arnoldi(A, Q, k)
  q = A*Q(:,k);
  for i = 1:k
    h(i)= q'*Q(:,i);
    q = q - h(i)*Q(:,i);
  end
  h(k+1) = norm(q);
  q = q / h(k+1);
end

%---------------------------------------------------------------------%
%                  Applying Givens Rotation to H col                  %
%---------------------------------------------------------------------%
function [h, cs_k, sn_k] = apply_givens_rotation(h, cs, sn, k)
  %apply for ith column
  for i = 1:k-1                              
    temp     =  cs(i)*h(i) + sn(i)*h(i+1);
    h(i+1) = -sn(i)*h(i) + cs(i)*h(i+1);
    h(i)   = temp;
  end
  
  %update the next sin cos values for rotation
  [cs_k sn_k] = givens_rotation(h(k), h(k+1));
  
  %eliminate H(i+1,i)
  h(k) = cs_k*h(k) + sn_k*h(k+1);
  h(k+1) = 0.0;
end

%%----Calculate the Given rotation matrix----%%
function [cs, sn] = givens_rotation(v1, v2)
  if (v1==0)
    cs = 0;
    sn = 1;
  else
    t=sqrt(v1^2+v2^2);
    cs = abs(v1) / t;
    sn = cs * v2 / v1;
  end
end

See also

• Biconjugate gradient method

Notes

• A. Meister, Numerik linearer Gleichungssysteme, 2nd edition, Vieweg 2005, ISBN 978-3-528-13135-7.
• Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, Society for Industrial and Applied Mathematics, 2003. ISBN 978-0-89871-534-7.
• Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., 7:856–869, 1986. doi:10.1137/0907058.
• S. C. Eisenstat, H.C. Elman and M.H. Schultz, "Variational iterative methods for nonsymmetric systems of linear equations", SIAM Journal on Numerical Analysis, 20(2), 345–357, 1983.
• J. Stoer and R. Bulirsch, Introduction to numerical analysis, 3rd edition, Springer, New York, 2002. ISBN 978-0-387-95452-3.
• Lloyd N. Trefethen and David Bau, III, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, 1997. ISBN 978-0-89871-361-9.
• Dongarra et al. , Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, 2nd Edition, SIAM, Philadelphia, 1994
• Amritkar, Amit; de Sturler, Eric; Świrydowicz, Katarzyna; Tafti, Danesh; Ahuja, Kapil (2015). "Recycling Krylov subspaces for CFD applications and a new hybrid recycling solver". Journal of Computational Physics 303: 222. doi:10.1016/j.jcp.2015.09.040

SOR - Convergence Rate

Convergence

The choice of relaxation factor ω is not necessarily easy, and depends upon the properties of the coefficient matrix. In 1947, Ostrowski proved that if $A$ is symmetric and positive-definite then $\rho (L_{\omega })<1$ for $0<\omega <2$. Thus, convergence of the iteration process follows, but we are generally interested in faster convergence rather than just convergence.

Convergence Rate

The convergence rate for the SOR method can be analytically derived. One needs to assume the following

• the relaxation parameter is appropriate: $\omega \in (0,2)$
• Jacobi's iteration matrix $C_{\text{Jac}}:=I-D^{-1}A$ has only real eigenvalues
• Jacobi's method is convergent: $\mu :=\rho (C_{\text{Jac}})<1$
• a unique solution exists: $\det A\neq 0$.

Then the convergence rate can be expressed as

$\rho (C_{\omega })={\begin{cases}{\frac {1}{4}}\left(\omega \mu +{\sqrt {\omega ^{2}\mu ^{2}-4(\omega -1)}}\right)^{2}\,,&0<\omega \leq \omega _{\text{opt}}\\\omega -1\,,&\omega _{\text{opt}}<\omega <2\end{cases}}$

where the optimal relaxation parameter is given by

$\omega _{\text{opt}}:=1+\left({\frac {\mu }{1+{\sqrt {1-\mu ^{2}}}}}\right)^{2}\,.$

Spectral radius $\rho (C_{\omega })$ of the iteration matrix for the SOR method $C_{\omega }$. The plot shows the dependence on the spectral radius of the Jacobi iteration matrix $\mu :=\rho (C_{\text{Jac}})$.

${\displaystyle }$

${\displaystyle }$

ILU - Stability

Concerning the stability of the ILU the following theorem was proven by Meijerink an van der Vorst.

Let $A$ be an M-matrix, the (complete) LU decomposition given by $A={\hat {L}}{\hat {U}}$, and the ILU by $A=LU-R$. Then

$|L_{ij}|\leq |{\hat {L}}_{ij}|\quad \forall \;i,j$

holds. Thus, the ILU is at least as stable as the (complete) LU decomposition.

ILU - Definition

For a given matrix $A\in \mathbb {R} ^{n\times n}$ one defines the graph $G(A)$ as

$G(A):=\left\lbrace (i,j)\in \mathbb {N} ^{2}:A_{ij}\neq 0\right\rbrace \,,$

which is used to define the conditions a sparsity patterns $S$ needs to fulfill

$S\subset \left\lbrace 1,\dots ,n\right\rbrace ^{2}\,,\quad \left\lbrace (i,i):1\leq i\leq n\right\rbrace \subset S\,,\quad G(A)\subset S\,.$

A decomposition of the form $A=LU-R$ which fulfills

• $L\in \mathbb {R} ^{n\times n}$ is a lower unitriangular matrix
• $U\in \mathbb {R} ^{n\times n}$ is an upper triangular matrix
• $L,U$ are zero outside of the sparsity pattern: $L_{ij}=U_{ij}=0\quad \forall \;(i,j)\notin S$
• $R\in \mathbb {R} ^{n\times n}$ is zero within the sparsity pattern: $R_{ij}=0\quad \forall \;(i,j)\in S$

is called an incomplete LU decomposition (w.r.t. the sparsity pattern $S$).

The sparsity pattern of L and U is often chosen to be the same as the sparsity pattern of the original matrix A. If the underlying matrix structure can be referenced by pointers instead of copied, the only extra memory required is for the entries of L and U. This preconditioner is called ILU(0).

${\displaystyle }$

${\displaystyle }$

CG - Convergence Theorem

Define a subset of polynomials as

$\Pi _{k}^{*}:=\left\lbrace \ p\in \Pi _{k}\ :\ p(0)=1\ \right\rbrace \,,$

where $\Pi _{k}$ is the set of polynomials of maximal degree $k$.

Let $\left(\mathbf {x} _{k}\right)_{k}$ be the iterative approximations of the exact solution $\mathbf {x} _{*}$, and define the errors as $\mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{*}$. Now, the rate of convergence can be approximated as

${\begin{aligned}\left\|\mathbf {e} _{k}\right\|_{\mathbf {A} }&=\min _{p\in \Pi _{k}^{*}}\left\|p(\mathbf {A} )\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq \min _{p\in \Pi _{k}^{*}}\,\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\\&\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\ \left\|\mathbf {e} _{0}\right\|_{\mathbf {A} }\,,\end{aligned}}$

where $\sigma (\mathbf {A} )$ denotes the spectrum, and $\kappa (\mathbf {A} )$ denotes the condition number.

Note, the important limit when $\kappa (\mathbf {A} )$ tends to $\infty$

${\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}}\quad {\text{for}}\quad \kappa (\mathbf {A} )\gg 1\,.$

This limit shows a faster convergence rate compared to the iterative methods of Jacobi or Gauss-Seidel which scale as $\approx 1-{\frac {2}{\kappa (\mathbf {A} )}}$.

${\displaystyle }$

${\displaystyle }$

SOR - Symmetric positive definite case

In case that the system matrix $A$ is of positive definite type one can show convergence.

Let $C=C_{\omega }=I-\left({\frac {1}{\omega }}D+L\right)^{-1}A$ be the iteration matrix. Then, convergence is guarenteed for

$\rho (C_{\omega })<1\quad \Longleftrightarrow \quad \omega \in (0,2)\,.$

Jacobi - Symmetric positive definite case

In case that the system matrix $A$ is of positive definite type one can show convergence.

Let $C=C_{\omega }=I-\omega D^{-1}A$ be the iteration matrix. Then, convergence is guarenteed for

$\rho (C_{\omega })<1\quad \Longleftrightarrow \quad 0<\omega <{\frac {2}{\lambda _{\text{max}}(D^{-1}A)}}\,,$

where $\lambda _{\text{max}}$ is the maximal eigenvalue.

The spectral radius can be minimized for a particular choice of $\omega =\omega _{\text{opt}}$ as follows

$\min _{\omega }\rho (C_{\omega })=\rho (C_{\omega _{\text{opt}}})=1-{\frac {2}{\kappa (D^{-1}A)+1}}\quad {\text{for}}\quad \omega _{\text{opt}}:={\frac {2}{\lambda _{\text{min}}(D^{-1}A)+\lambda _{\text{max}}(D^{-1}A)}}\,,$

where $\kappa$ is the matrix' condition number.

${\displaystyle }$

Hyperbolic system of partial differential equations

The following is a system of $s$ first order partial differential equations for $s$ unknown functions ${\vec {u}}=(u_{1},\ldots ,u_{s})$, ${\vec {u}}={\vec {u}}({\vec {x}},t)$, where ${\vec {x}}\in \mathbb {R} ^{d}$:

$(*)\quad {\frac {\partial {\vec {u}}}{\partial t}}+\sum _{j=1}^{d}{\frac {\partial }{\partial x_{j}}}{\vec {f^{j}}}({\vec {u}})=0,$

where ${\vec {f^{j}}}\in C^{1}(\mathbb {R} ^{s},\mathbb {R} ^{s}),j=1,\ldots ,d$ are once continuously differentiable functions, nonlinear in general.

Next, for each ${\vec {f^{j}}}$ a Jacobian matrix $s\times s$ is defined

$A^{j}:={\begin{pmatrix}{\frac {\partial f_{1}^{j}}{\partial u_{1}}}&\cdots &{\frac {\partial f_{1}^{j}}{\partial u_{s}}}\\\vdots &\ddots &\vdots \\{\frac {\partial f_{s}^{j}}{\partial u_{1}}}&\cdots &{\frac {\partial f_{s}^{j}}{\partial u_{s}}}\end{pmatrix}},{\text{ for }}j=1,\ldots ,d.$

The system $(*)$ is hyperbolic if for all $\alpha _{1},\ldots ,\alpha _{d}\in \mathbb {R}$ the matrix $A:=\alpha _{1}A^{1}+\cdots +\alpha _{d}A^{d}$ has only real eigenvalues and is diagonalizable.

If the matrix $A$ has s distinct real eigenvalues, it follows that it is diagonalizable. In this case the system $(*)$ is called strictly hyperbolic.

If the matrix $A$ is symmetric, it follows that it is diagonalizable and the eigenvalues are real. In this case the system $(*)$ is called symmetric hyperbolic.

Linear system

The case of a linear hyperbolic system of conservation laws (with constant coefficients in one space dimension) is given by

${\begin{aligned}\quad {\frac {\partial {\vec {u}}}{\partial t}}+A{\frac {\partial {\vec {u}}}{\partial x}}&=0\quad {\text{for}}\quad (x,t)\in \mathbb {R} \times (0,\infty )\\{\vec {u}}(x,0)&={\vec {u}}_{0}(x)\quad {\text{for}}\quad x\in \mathbb {R} \,,\end{aligned}}$

where one solves for the unknown function ${\vec {u}}:\mathbb {R} \times [0,\infty )\rightarrow \mathbb {R} ^{s}$ and initial data ${\vec {u}}_{0}:\mathbb {R} \rightarrow \mathbb {R} ^{s}$, and $A\in \mathbb {R} ^{s\times s}$ are given.

A hyperbolic system is real diagonalizable

${\begin{aligned}A&=R\Lambda R^{-1}\\{\text{with}}\quad \Lambda &=\mathrm {diag} (\lambda _{1},\dots ,\lambda _{s})\in \mathbb {R} ^{s\times s}\,,\quad R=({\vec {r}}_{1},\dots ,{\vec {r}}_{s})\in \mathbb {R} ^{s\times s}\,,\quad R^{-1}=({\vec {l}}_{1},\dots ,{\vec {l}}_{s})^{T}\in \mathbb {R} ^{s\times s}\,.\end{aligned}}$

Thus, the conservation law decouples into $s$ independent transport equations

${\begin{aligned}\quad {\frac {\partial {\vec {u}}}{\partial t}}+A{\frac {\partial {\vec {u}}}{\partial x}}&=0\\\Leftrightarrow \quad {\frac {\partial {\vec {v}}}{\partial t}}+\Lambda {\frac {\partial {\vec {v}}}{\partial x}}&=0\,,\quad {\vec {v}}:=R^{-1}{\vec {u}}\\\Leftrightarrow \quad {\frac {\partial v_{i}}{\partial t}}+\lambda _{i}{\frac {\partial v_{i}}{\partial x}}&=0\quad \forall \,i\,.\end{aligned}}$

The general solution is

$v_{i}(x,t)=v_{i}(x-\lambda _{i}t,0)\,,$

and in the original variables for given initial data ${\vec {u}}_{0}$

${\vec {u}}(x,t)=R{\vec {v}}(x,t)=\sum _{i=1}^{s}v_{i}(x-\lambda _{i}t,0){\vec {r}}_{i}=$

${\displaystyle }$

${\displaystyle }$

Example: The Laplace operator

The (continuous) Laplace operator in $n$-dimensions is given by $\Delta u(x)=\sum _{i=1}^{n}\partial _{i}^{2}u(x)$. The discrete Laplace operator $\Delta _{h}u$ depends on the dimension $n$.

In 1D the Laplace operator is approximated as

$\Delta u(x)=u''(x)\approx {\frac {u(x-h)-2u(x)+u(x+h)}{h^{2}}}=:\Delta _{h}u(x)\,.$

This approximation is usually expressed via the following stencil

${\frac {1}{h^{2}}}{\begin{bmatrix}1&-2&1\end{bmatrix}}\,.$

The 2D case shows all the characteristics of the more general nD case. Each second partial derivative needs to be approximated similar to the 1D case

${\begin{aligned}\Delta u(x,y)&=u_{xx}(x,y)+u_{yy}(x,y)\\&\approx {\frac {u(x-h,y)-2u(x)+u(x+h,y)}{h^{2}}}+{\frac {u(x,y-h)-2u(x)+u(x,y+h)}{h^{2}}}\\&={\frac {u(x-h,y)+u(x+h,y)-4u(x)+u(x,y-h)+u(x,y+h)}{h^{2}}}\\&=:\Delta _{h}u(x)\,,\end{aligned}}$

which is usually given by the following stencil

${\frac {1}{h^{2}}}{\begin{bmatrix}&1\\1&-4&1\\&1\end{bmatrix}}\,.$

Consistency

Consistency of the above mentioned approximation can be shown for highly regular functions, such as $u\in C^{4}(\Omega )$. The statement is

$\Delta u-\Delta _{h}u={\mathcal {O}}(h^{2})\,.$

To proof this one needs to substitute Taylor Series expansions up to order 3 into the discrete Laplace operator.

Properties

Subharmonic

Similar to continous subharmonic functions one can define subharmonic functions for finite-difference approximations $u_{h}$

$-\Delta _{h}u_{h}\leq 0\,.$

Mean value

One can define a general stencil of positive type via

${\begin{bmatrix}&\alpha _{N}\\\alpha _{W}&-\alpha _{C}&\alpha _{E}\\&\alpha _{S}\end{bmatrix}}\,,\quad \alpha _{i}>0\,,\quad \alpha _{C}=\sum _{i\in \{N,E,S,W\}}\alpha _{i}\,.$

If $u_{h}$ is (discrete) subharmonic, then the following mean value property holds

$u_{h}(x_{C})\leq {\frac {\sum _{i\in \{N,E,S,W\}}\alpha _{i}u_{h}(x_{i})}{\sum _{i\in \{N,E,S,W\}}\alpha _{i}}}\,,$

where the approximation is evaluated on points of the grid, and the stencil is assumed to be of positive type.

A similar mean value property also holds for the continuous case.

Maximum principle

For a (discrete) subharmonic function $u_{h}$ the following holds

$\max _{\Omega _{h}}u_{h}\leq \max _{\partial \Omega _{h}}u_{h}\,,$

where $\Omega _{h},\partial \Omega _{h}$ are discretizations of the continuous domain $\Omega$, respectively the boundary $\partial \Omega$.

Discontinuous Galerkin Scheme

Scalar hyperbolic conservation law

A scalar hyperbolic conservation law is of the form

${\begin{aligned}\partial _{t}u+\partial _{x}f(u)&=0\quad {\text{for}}\quad t>0,\,x\in \mathbb {R} \\u(0,x)&=u_{0}(x)\,,\end{aligned}}$

where one tries to solve for the unknown scalar function $u\equiv u(t,x)$, and the functions $f,u_{0}$ are typically given.

Space discretization

The $x$-space will be discretized as

$\mathbb {R} =\bigcup _{k}I_{k}\,,\quad I_{k}:=\left(x_{k},x_{k+1}\right)\quad {\text{for}}\quad x_{k}<x_{k+1}\,.$

Furthermore, we need the following definitions

$h_{k}:=|I_{k}|\,,\quad h:=\sup _{k}h_{k}\,,\quad {\hat {x}}_{k}:=x_{k}+{\frac {h_{k}}{2}}\,.$

Basis for function space

We derive the basis representation for the function space of our solution $u$. The function space is defined as

$S_{h}^{p}:=\left\lbrace v\in L^{2}(\mathbb {R} )\;\colon \;{v|}_{I_{k}}\in \Pi _{p}\right\rbrace \quad {\text{for}}\quad p\in \mathbb {N} _{0}\,,$

where ${v|}_{I_{k}}$ denotes the restriction of $v$ onto the interval $I_{k}$, and $\Pi _{p}$ denotes the space of polynomials of maximal degree $p$. The index $h$ should show the relation to an underlying discretization given by $\left(x_{k}\right)_{k}$. Note here that $v$ is not uniquely defined at the intersection points $\left(x_{k}\right)_{k}$.

At first we make use of a specific polynomial basis on the interval $[-1,1]$, the Legendre_polynomials $\left(P_{n}\right)_{n\in \mathbb {N} _{0}}$, i.e.,

$P_{0}(x)=1\,,\quad P_{1}(x)=x\,,\quad P_{2}(x)={\frac {1}{2}}\left(3x^{2}-1\right)\,,\dots$

Note especially the orthogonality relations

$\left\langle P_{i},P_{j}\right\rangle _{L^{2}([-1,1])}={\frac {2}{2i+1}}\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\,.$

Transformation onto the interval $[0,1]$, and normalization is achieved by functions $\left(\phi _{i}\right)_{i}$

$\phi _{i}(x):={\sqrt {2i+1}}P_{i}(2x-1)\quad {\text{for}}\quad x\in [0,1]\,,$

which fulfill the orthonormality relation

$\left\langle \phi _{i},\phi _{j}\right\rangle _{L^{2}([0,1])}=\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\,.$

Transformation onto an interval $I_{k}$ is given by $\left({\bar {\varphi }}_{ki}\right)_{i}$

${\bar {\varphi }}_{ki}:={\frac {1}{\sqrt {h_{k}}}}\phi _{i}\left({\frac {x-x_{k}}{h_{k}}}\right)\quad {\text{for}}\quad x\in I_{k}\,,$

which fulfill

$\left\langle {\bar {\varphi }}_{ki},{\bar {\varphi }}_{kj}\right\rangle _{L^{2}(I_{k})}=\delta _{ij}\quad \forall \,i,j\in \mathbb {N} _{0}\forall \,k\,.$

For $L^{\infty }$-normalization we define $\varphi _{ki}:={\sqrt {h_{k}}}{\bar {\varphi }}_{ki}$, and for $L^{1}$-normalization we define ${\tilde {\varphi }}_{ki}:={\frac {1}{\sqrt {h_{k}}}}{\bar {\varphi }}_{ki}$, s.t.

$\|\varphi _{ki}\|_{L^{\infty }(I_{k})}=\|\phi _{i}\|_{L^{\infty }([0,1])}=:c_{i,\infty }\quad {\text{and}}\quad \|{\tilde {\varphi }}_{ki}\|_{L^{1}(I_{k})}=\|\phi _{i}\|_{L^{1}([0,1])}=:c_{i,1}\,.$

Finally, we can define the basis representation of our solutions $u_{h}$

${\begin{aligned}u_{h}(t,x):=&\sum _{i=0}^{p}u_{ki}(t)\varphi _{ki}(x)\quad {\text{for}}\quad x\in (x_{k},x_{k+1})\\u_{ki}(t)=&\left\langle u_{h}(t,\cdot ),{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}\,.\end{aligned}}$

Note here, that $u_{h}$ is not defined at the interface positions.

DG-Scheme

The conservation law is transformed into its weak form by multiplying with test functions, and integration over test intervals

${\begin{aligned}\partial _{t}u+\partial _{x}f(u)&=0\\\Rightarrow \quad \left\langle \partial _{t}u,v\right\rangle _{L^{2}(I_{k})}+\left\langle \partial _{x}f(u),v\right\rangle _{L^{2}(I_{k})}&=0\quad {\text{for}}\quad \forall \,v\in S_{h}^{p}\\\Leftrightarrow \quad \left\langle \partial _{t}u,{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}+\left\langle \partial _{x}f(u),{\tilde {\varphi }}_{ki}\right\rangle _{L^{2}(I_{k})}&=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}$

By using partial integration one is left with

${\begin{aligned}{\frac {\mathrm {d} }{\mathrm {d} t}}u_{ki}(t)+f(u(t,x_{k+1})){\tilde {\varphi }}_{ki}(x_{k+1})-f(u(t,x_{k})){\tilde {\varphi }}_{ki}(x_{k})-\left\langle f(u(t,\,\cdot \,)),{\tilde {\varphi }}_{ki}'\right\rangle _{L^{2}(I_{k})}=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}$

The fluxes at the interfaces are approximated by numerical fluxes $g$ with

$g_{k}:=g(u_{k}^{-},u_{k}^{+})\,,\quad u_{k}^{\pm }:=u(t,x_{k}^{\pm })\,,$

where $u_{k}^{\pm }$ denotes the left- and right-hand sided limits. Finally, the DG-Scheme can be written as

${\begin{aligned}{\frac {\mathrm {d} }{\mathrm {d} t}}u_{ki}(t)+g_{k+1}{\tilde {\varphi }}_{ki}(x_{k+1})-g_{k}{\tilde {\varphi }}_{ki}(x_{k})-\left\langle f(u(t,\,\cdot \,)),{\tilde {\varphi }}_{ki}'\right\rangle _{L^{2}(I_{k})}=0\quad {\text{for}}\quad \forall \,k\;\forall \,i\leq p\,.\end{aligned}}$