14 M Matrices and Stability in ∞-norm
Lecture 14
I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.
These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.
14.1 Recap
We are considering boundary value problems of the following form: Suppose we have a differential operator \mathcal L, a “nice” domain \Omega \subset \mathbb R^d, and data f. Seek u : \Omega \to \mathbb R such that
\begin{align} \mathcal L[u] &= f &&\text{in } \Omega \nonumber\\ u &= 0 &&\text{on } \partial \Omega . \tag{BVP} \end{align}
This is a problem with homogeneous Dirichlet boundary conditions (the values of the solution are prescribed on the boundary (Dirichlet) and are zero (homogeneous)).
We also consider the specific problem (in 1d):
Example 14.1 (Linear Advection-Diffusion-Reaction Equations) Consider
\begin{align} \mathcal L[u] := -\Delta u + p \cdot \nabla u + q u \end{align}
In the context of semiconductor physics, u(x) represents the concentration of charge carriers at a position x \in \Omega. Each of these terms relates to the physics in question:
- -\Delta : diffusion; electrons move from areas of high concentration to low,
- p \cdot \nabla u : advection (drift); the solution is driven with velocity p (due to e.g. an external electric field),
- q u : reaction; electrons and electron holes are being created or destroyed.
The same equation also models e.g pollutant transport in a river, drug delivery in a blood vessel, asset price in finance, heat transfer in fluids, ……
By approximating the continuous problem with finite differences, we obtain a discrete problem: find \bm U \in \mathbb R^{n+1} (with index starting at 0 or equivalently U: \Omega_h \to \mathbb R where \Omega_h := \{x_j\}_{j=0}^n \to \mathbb R with U_j := U(x_j)) such that
\begin{align} \begin{split} &\mathcal L_h [\bm U] = \bm f \\ &U_0 = U_n = 0 \end{split}\tag{$\text{BVP}_h$} \end{align}
where \mathcal L_h depends on the choice of finite difference scheme.
Last time, we approximated the ADR equation with the following finite difference scheme on the mesh \{x_j\}_{j=0}^n with x_j = \frac{j}{n} and mesh size h= \frac{1}{n}:
\begin{align} \mathcal L_h &= -L + p D^0 + qI \nonumber\\ % &= -\frac{1}{h^2} \begin{pmatrix} -2 & 1 & & & \\ 1 & -2 & 1 & & \\ & 1 & -2 & \ddots \\ & & \ddots & \ddots & 1\\ & & & 1 & -2 \end{pmatrix} % + \frac{p}{2h} \begin{pmatrix} 0 & 1 & & & \\ -1 & 0 & 1 & & \\ & -1 & 0 & \ddots \\ & & \ddots & \ddots & 1\\ & & & -1 & 0 \end{pmatrix} % + q I, \nonumber \end{align}
a (n-1)\times(n-1) matrix (acting on the interior points of the grid).
We saw that this method is 2nd order consistent:
Definition 14.1 (Consistency) Suppose that u : \Omega \to \mathbb R is the exact solution to \text{(BVP)} and let \bm u = \{ u(x_j) \}_{j=0}^n be the exact solution evaluated on the mesh. We say the numerical approximation (\text{BVP}_h) is consistent if \| \mathcal L_h[\bm u] - \bm f\|_{\infty} \to 0 as h \to 0 where f_j = f(x_j).
We say the numerical approximation is order-p accurate (or has order of accuracy p) if \| \mathcal L_h[\bm u] - \bm f\|_{\infty} = O(h^p) as h \to 0.
and, for q>0, stable in the following sense:
Definition 14.2 (Stability) We say that (\text{BVP}_h) is stable if there exists C>0 such that
\begin{align} \| \bm U_h - \bm W_h \|_{\infty} \leq C \| \mathcal L_h[\bm U_h] - \mathcal L_h[\bm W_h] \|_{\infty} \end{align}
for all \bm U_h, \bm W_h.
We showed that the stability constant (C in the above defininition) for \mathcal L_h is (at most) C = \frac{1}{q}. Therefore, this finite difference approximation converges:
Theorem 14.1 Suppose that \text{(BVP)} and (\text{BVP}_h) have unique solutions u and \bm U_h. Then, if (\text{BVP}_h) is consistent (with order p) and stable then the method is convergent (with order at least p).
Proof.
Let \bm u be the exact solution to the continuous problem restricted to the mesh and suppose that (\text{BVP}_h) is order p accurate. Then, by stability, we have
\begin{align} \| \bm U_h - \bm u \|_\infty &\leq C \| \mathcal L_h[\bm U_h] - \mathcal L_h[\bm u] \|_\infty \nonumber\\ % &= C \| \bm f - \mathcal L_h[\bm u] \|_\infty \nonumber\\ % &= O(h^p) \end{align}
as h\to0. (The final equality is the consistency estimate).
To prove the stability estimate, we showed proved the following result:
Lemma 14.1 Suppose A is an SDD matrix. Then,
\begin{align} \|x\|_\infty \leq % \max_i \left[ |A_{ii}| - \sum_{j\not=i} |A_{ij}| \right]^{-1} \| Ax\|_{\infty}. \end{align}
Exercise 14.1 Show that, for p \in\mathbb R and q > 0, the approximation scheme (\text{BVP}_h) is stable for all h < \frac{2}{|p|}. What is the corresponding stability constant?
It turns out that we can extend these stability estimates to the case where q=0:
14.2 M Matrices
The aim of this section is to understand and prove the following result:
Theorem 14.2 Suppose that A is an M-matrix with \Phi \succ 0 and A\Phi \succ 0.
Then,
\begin{align} &\|x\|_{\infty} \leq C_\Phi \|A x\|_{\infty}, \nonumber\\ % &\text{with} \quad C_\Phi := \frac{\max_i \Phi_i}{\min_i (A\Phi)_i} \end{align}
Here, we are using the notation
- \Phi \succ 0 if and only if \Phi_i > 0 for all i,
- Similarly, we use x \preccurlyeq y or y \succcurlyeq x to mean x_i \leq y_i for all i,
and the following definition:
Definition 14.3 A matrix A is an M matrix if
(i) A_{ii} > 0 for all i,
(ii) A_{ij} \leq 0 for all i \not= j,
(iii) (Monotonicity) Ax \succcurlyeq 0 \implies x \succcurlyeq 0.
It is quite difficult to check if a matrix is an M matrix directly ((iii) is hard to check). It turns out that all Irreducibly Diagonally Dominant (IDD) matrices that satisfy (i) and (ii) are M matrices:
Definition 14.4 A matrix A is Irreducibly Diagonally Dominant (IDD) if
- |A_{ii}| \geq \sum_{j\not=i} |A_{ij}| for all i,
- There exists i^\star such that |A_{i^\star i^\star}| > \sum_{j\not=i^\star} |A_{i^\star j}|,
- For every i and j there exists k_1,\dots,k_n such that A_{ik_1}, A_{k_1,k_2}, ..., A_{k_n,j} \not= 0.
It turns out that IDD matrices are invertible. Moreover, we have:
Lemma 14.2 Suppose that A is IDD and A_{ii} > 0 for all i and A_{ij} \leq 0 for all i\not= j. Then, A is a M matrix
Proof. We suppose that b := Ax \succcurlyeq 0 (we want to show that x \succcurlyeq 0).
Take i such that x_i = \min_j x_j < 0. We have,
\begin{align} 0 \leq b_i &= A_{ii} x_i + \sum_{j\not= i} A_{ij} x_j \nonumber\\ % &= |A_{ii}| x_i - \sum_{j\not= i} |A_{ij}| x_j \nonumber\\ % &\leq \Big( |A_{ii}| - \sum_{j\not= i} |A_{ij}|\Big) x_i % \leq 0 \end{align}
In the last line, we use the fact that -x_j \leq -x_i (since x_i \leq x_j) and the fact that A is diagonally dominant. That is, all of the above inequalities are actually equalities and thus b_i = 0 and, if A_{ij}\not= 0, then x_i = x_j.
Repeating the exact same argment with j instead of i, we find that b_j = 0 and x_{j} = x_k for all k with A_{jk} \not= 0. Since A is irreducible, we find that this argument “spreads” to all indices and b = 0. However, because there exists a row i^\star for which we have strict diagonal dominance (|A_{i^\star i^\star}| - \sum_{j\not= i^\star} |A_{i^\star j}| >0), we have b_{i^\star} < 0, a contradiction.
As a result, we have that x_{i} \geq 0 and thus x \succcurlyeq 0, as required. \square
Exercise 14.2 For what choice of p,q is \mathcal L_h an M matrix?
Exercise 14.3 Suppose that A is an M matrix. Show that [A^{-1}]_{ij} \geq 0 for all i,j.
We now prove Theorem Theorem thm-M:
Proof. We let A be an M matrix so that A_{ii} > 0 for all i, A_{ij} \leq 0 for all i\not= j and A is monotone: Ax \succcurlyeq 0 implies x \succcurlyeq 0. We also suppose there exists \Phi \succ 0 such that A\Phi \succ 0.
We let c := \min_i (A\Phi)_i so that A\Phi \succcurlyeq c \begin{pmatrix} 1 & 1 & \dots & 1 \end{pmatrix}^\top. Define \Psi := c A^{-1} \begin{pmatrix} 1 & 1 & \dots & 1 \end{pmatrix}^\top and note that A(\Phi - \Psi) \succcurlyeq 0. As a result we have \Phi - \Psi \succcurlyeq 0 and thus \Phi \succcurlyeq \Psi. That is, \sum_j [A^{-1}]_{ij} \leq c^{-1} \Phi_i. Therefore,
\begin{align} \pm x_i % &= \sum_{j} [A^{-1}]_{ij} (\pm A x)_j \nonumber\\ % &\leq \sum_{j} [A^{-1}]_{ij} \|A x\|_\infty \nonumber\\ % &\leq \frac{\Phi_i}{c} \|A x\|_\infty % \leq \frac{\max \Phi_i}{\min_i (A\Phi)_i } \|A x\|_\infty, \end{align}
as required. \square
Exercise 14.4 (i) Notice that \Phi_j := \frac12 jh (1- jh) satisfies \Phi_j > 0 for all j=1,\dots,n-1,
(ii) Use the fact that -L is second order accurate to show that -L \Phi = 1,
(iii) Apply the theorem to with choice of \Phi. What is C_{\Phi}?
(iv) Conclude that approximating -u''(x) = f with u(0) = u(1) = 0 with the discrete problem -L \bm u = \bm f with u_0 = u_n = 0 is 2^\text{nd} order convergent.
- Take a look at the suggestions for Presentations/Posters/Papers and sign up to present,
- Read through the proofs from this lecture - come with questions on Wednesday,
- (Sign-up to office hours if you would like)
Previous reading:
- Chapter 1 reading:
- Burden, Faires, and Burden (2015) sections 5.1 - 5.6,
- Driscoll and Braun (2018) zero stability
- Burden, Faires, and Burden (2015) sections 5.1 - 5.6,
- Chapter 2 reading:
- Chapter 3 reading:
- Sections 11.1 - 11.2 of Burden, Faires, and Burden (2015): shooting methods