{ "cells": [ { "cell_type": "markdown", "id": "81b1d2f0", "metadata": {}, "source": [ "---\n", "title: \"M Matrices and Stability in ∞-norm\"\n", "subtitle: \"Lecture 14\"\n", "date: 2026-03-30\n", "abstract-title: Overview\n", "format:\n", " html:\n", " other-links:\n", " - text: This notebook\n", " href: L14.ipynb\n", "---" ] }, { "cell_type": "markdown", "id": "63cf0867", "metadata": {}, "source": [ "::: {.callout-note}\n", "\n", "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.\n", "\n", ":::\n", "\n", "::: {.callout-warning}\n", "\n", "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.\n", "\n", ":::\n", "\n", "::: {.callout-tip}\n", "\n", "This chapter is mostly based on \n", "\n", "* @Burden Chapter 11, \n", "* @FNC Chapter 10\n", "\n", ":::" ] }, { "cell_type": "code", "execution_count": 1, "id": "f8408d86", "metadata": {}, "outputs": [], "source": [ "#| echo: false\n", "\n", "using Plots, LaTeXStrings, LinearAlgebra\n", "\n", "# IF YOU ARE READING THIS, THANK YOU\n", "# FOR DOWNLOADING THE LECTURE NOTES.\n", "# THE FIRST PERSON TO LET ME KNOW,\n", "# WINS SOME CHOCOLATES OF YOUR CHOICE" ] }, { "cell_type": "markdown", "id": "81473756", "metadata": {}, "source": [ "## Recap\n", "\n", "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 \n", "\n", "\\begin{align}\n", " \\mathcal L[u] &= f &&\\text{in } \\Omega \\nonumber\\\\\n", " u &= 0 &&\\text{on } \\partial \\Omega .\n", " \\tag{BVP}\n", "\\end{align}\n", "\n", "This is a problem with *homogeneous Dirichlet boundary conditions* (the values of the solution are prescribed on the boundary (Dirichlet) and are zero (homogeneous)). \n", "\n", "We also consider the specific problem (in 1d):\n", "\n", "::: {#exm-BVPs}\n", "# Linear Advection-Diffusion-Reaction Equations\n", "\n", "Consider\n", "\n", "\\begin{align}\n", " \\mathcal L[u] := -\\Delta u + p \\cdot \\nabla u + q u\n", "\\end{align}\n", "\n", "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: \n", "\n", "* $-\\Delta$ : diffusion; electrons move from areas of high concentration to low, \n", "* $p \\cdot \\nabla u$ : advection (drift); the solution is driven with velocity $p$ (due to e.g. an external electric field), \n", "* $q u$ : reaction; electrons and electron holes are being created or destroyed.\n", "\n", "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, ...... \n", "\n", ":::\n", "\n", "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 \n", "\n", "\\begin{align}\n", " \\begin{split}\n", " &\\mathcal L_h [\\bm U] = \\bm f \\\\\n", " &U_0 = U_n = 0 \n", " \\end{split}\\tag{$\\text{BVP}_h$}\n", "\\end{align}\n", "\n", "where $\\mathcal L_h$ depends on the choice of finite difference scheme.\n", "\n", "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}$:\n", "\n", "\\begin{align}\n", " \\mathcal L_h &= -L + p D^0 + qI \\nonumber\\\\\n", " %\n", " &= -\\frac{1}{h^2} \n", " \\begin{pmatrix}\n", " -2 & 1 & & & \\\\\n", " 1 & -2 & 1 & & \\\\\n", " & 1 & -2 & \\ddots \\\\\n", " & & \\ddots & \\ddots & 1\\\\\n", " & & & 1 & -2 \n", " \\end{pmatrix}\n", " %\n", " + \\frac{p}{2h} \n", " \\begin{pmatrix}\n", " 0 & 1 & & & \\\\\n", " -1 & 0 & 1 & & \\\\\n", " & -1 & 0 & \\ddots \\\\\n", " & & \\ddots & \\ddots & 1\\\\\n", " & & & -1 & 0 \n", " \\end{pmatrix}\n", " % \n", " + q I, \\nonumber\n", "\\end{align}\n", "\n", "a $(n-1)\\times(n-1)$ matrix (acting on the interior points of the grid).\n", "\n", "We saw that this method is 2nd order consistent: \n", "\n", "::: {#def-consistent}\n", "# Consistency\n", "\n", "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)$. \n", "\n", "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$. \n", "\n", ":::\n", "\n", "and, for $q>0$, stable in the following sense:\n", "\n", "::: {#def-stable}\n", "# Stability\n", "\n", "We say that $(\\text{BVP}_h)$ is stable if there exists $C>0$ such that \n", "\n", "\\begin{align}\n", " \\| \\bm U_h - \\bm W_h \\|_{\\infty} \\leq C \\| \\mathcal L_h[\\bm U_h] - \\mathcal L_h[\\bm W_h] \\|_{\\infty}\n", "\\end{align}\n", "\n", "for all $\\bm U_h, \\bm W_h$.\n", "\n", ":::\n", "\n", "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: \n", "\n", "::: {#thm-fundamental}\n", "\n", "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$).\n", "\n", ":::\n", "\n", "::: {.box}\n", "\n", "***Proof.*** \n", "\n", "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 \n", "\n", "\\begin{align}\n", " \\| \\bm U_h - \\bm u \\|_\\infty &\\leq C \\| \\mathcal L_h[\\bm U_h] - \\mathcal L_h[\\bm u] \\|_\\infty \\nonumber\\\\\n", " %\n", " &= C \\| \\bm f - \\mathcal L_h[\\bm u] \\|_\\infty \\nonumber\\\\\n", " %\n", " &= O(h^p)\n", "\\end{align}\n", "\n", "as $h\\to0$. (The final equality is the consistency estimate).\n", "\n", ":::\n", "\n", "To prove the stability estimate, we showed proved the following result: \n", "\n", "::: {#lem-SDD}\n", "\n", "Suppose $A$ is an SDD matrix. Then, \n", "\n", "\\begin{align}\n", " \\|x\\|_\\infty \\leq \n", " %\n", " \\max_i \\left[ |A_{ii}| - \\sum_{j\\not=i} |A_{ij}| \\right]^{-1} \\| Ax\\|_{\\infty}.\n", "\\end{align}\n", "\n", ":::\n", "\n", "::: {#exr-}\n", "\n", "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?\n", "\n", ":::\n", "\n", "It turns out that we can extend these stability estimates to the case where $q=0$:" ] }, { "cell_type": "markdown", "id": "5e412da6", "metadata": {}, "source": [ "## M Matrices\n", "\n", "The aim of this section is to understand and prove the following result: \n", "\n", "::: {#thm-M}\n", "\n", "Suppose that $A$ is an M-matrix with $\\Phi \\succ 0$ and $A\\Phi \\succ 0$.\n", "\n", "Then, \n", "\n", "\\begin{align}\n", " &\\|x\\|_{\\infty} \\leq C_\\Phi \\|A x\\|_{\\infty}, \\nonumber\\\\\n", " %\n", " &\\text{with} \\quad C_\\Phi := \\frac{\\max_i \\Phi_i}{\\min_i (A\\Phi)_i}\n", "\\end{align}\n", "\n", ":::\n", "\n", "Here, we are using the notation \n", "\n", "* $\\Phi \\succ 0$ if and only if $\\Phi_i > 0$ for all $i$, \n", "* Similarly, we use $x \\preccurlyeq y$ or $y \\succcurlyeq x$ to mean $x_i \\leq y_i$ for all $i$,\n", "\n", "and the following definition:\n", "\n", "::: {#def-M}\n", "\n", "A matrix $A$ is an M matrix if \n", "\n", "*(i)* $A_{ii} > 0$ for all $i$, \n", "*(ii)* $A_{ij} \\leq 0$ for all $i \\not= j$, \n", "*(iii)* (Monotonicity) $Ax \\succcurlyeq 0 \\implies x \\succcurlyeq 0$.\n", "\n", ":::\n", "\n", "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:\n", "\n", "::: {#def-IDD}\n", "\n", "A matrix $A$ is *Irreducibly Diagonally Dominant (IDD)* if \n", "\n", "* $|A_{ii}| \\geq \\sum_{j\\not=i} |A_{ij}|$ for all $i$, \n", "* There exists $i^\\star$ such that $|A_{i^\\star i^\\star}| > \\sum_{j\\not=i^\\star} |A_{i^\\star j}|$, \n", "* 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$.\n", "\n", ":::\n", "\n", "It turns out that IDD matrices are invertible. Moreover, we have:\n", "\n", "::: {#lem-IDD-}\n", "\n", "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\n", "\n", ":::\n", "\n", "**Proof.** We suppose that $b := Ax \\succcurlyeq 0$ (we want to show that $x \\succcurlyeq 0$).\n", "\n", "Take $i$ such that $x_i = \\min_j x_j < 0$. We have, \n", "\n", "\\begin{align}\n", " 0 \\leq b_i &= A_{ii} x_i + \\sum_{j\\not= i} A_{ij} x_j \\nonumber\\\\\n", " %\n", " &= |A_{ii}| x_i - \\sum_{j\\not= i} |A_{ij}| x_j \\nonumber\\\\\n", " %\n", " &\\leq \\Big( |A_{ii}| - \\sum_{j\\not= i} |A_{ij}|\\Big) x_i\n", " %\n", " \\leq 0\n", "\\end{align}\n", "\n", "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$. \n", "\n", "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. \n", "\n", "As a result, we have that $x_{i} \\geq 0$ and thus $x \\succcurlyeq 0$, as required. $\\square$\n", "\n", "::: {#exr-}\n", "\n", "For what choice of $p,q$ is $\\mathcal L_h$ an M matrix?\n", "\n", ":::\n", "\n", "::: {#exr-}\n", "\n", "Suppose that $A$ is an M matrix. Show that $[A^{-1}]_{ij} \\geq 0$ for all $i,j$.\n", "\n", ":::\n", "\n", "We now prove Theorem @thm-M:\n", "\n", "**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$.\n", "\n", "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, \n", "\n", "\\begin{align}\n", " \\pm x_i \n", " %\n", " &= \\sum_{j} [A^{-1}]_{ij} (\\pm A x)_j \\nonumber\\\\\n", " %\n", " &\\leq \\sum_{j} [A^{-1}]_{ij} \\|A x\\|_\\infty \\nonumber\\\\\n", " %\n", " &\\leq \\frac{\\Phi_i}{c} \\|A x\\|_\\infty\n", " %\n", " \\leq \\frac{\\max \\Phi_i}{\\min_i (A\\Phi)_i } \\|A x\\|_\\infty,\n", "\\end{align}\n", "\n", "as required. $\\square$\n", "\n", "::: {#exr-}\n", "\n", "*(i)* Notice that $\\Phi_j := \\frac12 jh (1- jh)$ satisfies $\\Phi_j > 0$ for all $j=1,\\dots,n-1$, \n", "*(ii)* Use the fact that $-L$ is second order accurate to show that $-L \\Phi = 1$, \n", "*(iii)* Apply the theorem to with choice of $\\Phi$. What is $C_{\\Phi}$? \n", "*(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.\n", "\n", ":::" ] }, { "cell_type": "markdown", "id": "b94ee96d", "metadata": {}, "source": [ ":::{.callout-tip}\n", "# TO-DO\n", "\n", "* Take a look at the suggestions for [Presentations/Posters/Papers](/Presentations.html) and sign up to present, \n", "* Read through the proofs from this lecture - come with questions on Wednesday, \n", "* (Sign-up to office hours if you would like)\n", "\n", "Previous reading:\n", "\n", "* Chapter 1 reading:\n", " * @Burden sections 5.1 - 5.6, \n", " * @FNC [zero stability](https://fncbook.com/zerostability/) \n", "* Chapter 2 reading:\n", " * Recap matrices: sections 7.1 - 7.2 of @Burden \n", " * Jacobi & Gauss-Seidel: section 7.3 @Burden \n", " * CG: section 7.6 @Burden\n", "* Chapter 3 reading: \n", " * Sections 11.1 - 11.2 of @Burden: shooting methods\n", "\n", ":::" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.11", "language": "julia", "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.11.6" } }, "nbformat": 4, "nbformat_minor": 5 }