ETH_NumericalMethodsForPDE

links - NumPDE

website: grausam
script: https://people.math.ethz.ch/~grsam/NUMPDEFL/NUMPDE.pdf
- local copy: NUMPDE.pdf
exercises: grausam/NPDEFL_Problems.pdf
- local copy: NPDEFL_Problems.pdf
exercise repository: https://gitlab.math.ethz.ch/ralfh/NPDERepo
community solutions: exams.vis.ethz.ch

useful links

cpp reference: https://en.cppreference.com/
lehrfm: GitHub / doxygen
eigen documentation: https://eigen.tuxfamily.org/

Info

exam
- NumPDE midterm (@2026-04-17)
  - chapter 1+2 (mainly)
- NumPDE endterm (@2026-05-29)
  - chapter 9 (time-derivative)+12 (??)
practice classes:
- ETH_NumericalMethodsForPDE_Übungen

Info

HPK D 3 (13:45, Thu) https://polybox.ethz.ch/index.php/s/dM5NqE3wqeWYsNp
CAB G 61 (14:00, Thu) https://polybox.ethz.ch/index.php/s/Xen2eXQaoQd6Css

remove pauses in videos:

pip install auto-editor
auto-editor <video> --margin 0.2sec
# if output is black, stream and audio channels are swapped. fix it (before runnign audio-editor):
ffmpeg -i <video> -map 0:v -map 0:a -c copy <output>

Hilfsmittel NumericalMethodsForPartialDifferentialEquations

midterm: none
endterm: probably none

Exercises

-> FS2026_tasks

Week 1 (16.02.-20.02.): 4 h 45 min

Problem 0-1: Setting up implementation environment for NumPDE course (60 min)
Problem 0-2: Using the GDB Debugger (80 min)
Problem 1-1: Quadratic Functionals (100 min) (c) not sure about sums, (b) solved differently)
Problem 1-3: L-infinity Norms are Bounded by H1-seminorms in 1D (45 min) b) solved with Chat

Week 2 (23.02.-27.02.): 6 h 45 min

Problem 1-10: From Boundary Value Problems to Variational Problems (60 min) c) a bit cryptic
Problem 1-6: heat conduction with non-local boundary conditions (75 min) c) important, applies most of the learned norms/inequalities
Problem 1-8: A coupled reaction-diffusion problem (60 min)
Problem 1-5: A second-order elliptic transmission problem in 1D (110 min) do again tomorrow?
Problem 1-2: Linear functionals on Sobolev spaces (100 min)
- solution for h): 20260226_NumPDE_1.2)h)

Week 3 (02.03.-06.03.): 7 h 5 min

Problem 2-3: Pointwise "Exact" Galerkin Solution (65 min)
Problem 2-1: Properties of Galerkin solutions (45 min) much longer than 45min, c), d) enjoyable
Problem 2-2: Transformation of Galerkin Matrices (120 min) skipped d), copied e)
Problem 2-4: Linear finite elements for two-point boundary value problems (195 min) demanding exercise, took ≅3.5h; didn't understood last part (offset function)

Week 4 (09.03.-13.03.): 5 h 45 min

Problem 2-23: Galerkin Matrices (35 min)
Problem 2-20: DofHandler and Assembly (30 min)
Problem 2-9: Handling degrees of freedom (DOFs) in LEHRFEM++ (140 min) last few subtasks done with LM
Problem 2-8: Introduction to local assembly in LEHRFEM++ (140 min) couldn't implement j) myself

Week 5 (16.03.-20.03.): 7 h 5 min

Problem 2-24: Paramatric FEM: Local Computations (20 min) maybe do again
Problem 2-17: Local Computations for Convection Bilinear Form (90 min) solved in <90min, but didn't understand e) completely
Problem 2-21: Blended Parametric Representation of Curvi-Linear triangles (145 min) f) a bit tricky, but mostly understood
Problem 2-13: Local computations for parametric Lagrangian finite elements (170 min) took ~250min, quite difficult (especially writing the tests); I still don't understand the boundary logic

Week 6 (23.03.-27.03.): 5 h 50 min

Problem 2-22: Constrained Neumann Problem (125 min) a)-f) simple, g) uses summation trick that took me quite some time to understand; overall, took ~90min
Problem 2-18: Nitsche's Method for Elliptic BVPs (225 min) d)e)f) (coding) confusing - couldn't find error - should definitely be solved a 2nd time

Week 7 (30.03.-03.04.): 6 h 5 min

Problem 3-6: Projection Onto Constants (70 min) fun task, if you want to practice norms && inequalities
Problem 3-15: Asymptotic Convergence of Finite-Element Discretization Errors and Interpolation Errors (30 min) applying theory, do again
Problem 3-8: Output Functionals for a 2nd-Order Elliptic Boundary Value Problem with Impedance Boundary Conditions (145 min) took much longer, c) useful for practicing LF++, but what is e)-i) about?
Problem 3-2: Debugging Finite Element Codes (120 min) a) is Galerkin assembly (coding)

Week 8 (13.04.-17.04.): 4 h

Problem 9-14: Two-Step Radau Runge-Kutta Single-Step Method for the Heat Equation (50 min) "simple" runge-kutta repetition
Problem 9-20: Method-of-Lines with SDIRK timestepping (60 min)
Problem 9-3: Decaying Method-of-Lines Solution with Implicit-Euler Timestepping (130 min)

Week 9 (20.04.-24.04.): 3 h 30 min

Problem 9-17: Scalar Parabolic Evolution Problems: Continuous and Discrete (75 min) similar to 9.14, 9.20 from last week
Problem 9-11: Gauss-Lobatto IIIC Timestepping (135 min) finally a fun programming execise I did understand

Week 10 (27.04.-01.05.): 4 h 10 min

Problem 9-21: Magnetic Diffusion in a Wire (210 min)
Problem 12-4: FEM for Stokes BvPs (40 min) important basics

Week 11 (04.05.-08.05.): 7 h 15 min

Problem 12-1: Viscous Flow Simulation with Taylor-Hood FEM (435 min) only Hiptmair and god knew what this problem is about. Now, only Hiptmair knows
- this was the most challenging task till now. It takes way longer than 7.5 hours

Week 12 (11.05.-15.05.): 7h 10 min

Problem 10-6: Streamline Upwind Method for Pure Advection Problem (200 min) interesting, nearly succeeded in coding alone -> summary of chapter 2, 10
Problem 12-2: The MINI Element for Stokes (230 min) didn't totally understand b)
- in stokesminielement.h, dimension needs to be changed 9x9 -> 11x11 but is not mentioned in solution (error was extremely confusing)
- in solution of 12-2)i), the norms for 12.2.27 and 12.2.28 are wrong? ask TA?

13. week (18.05.-22.05.): 7 h 10 min

Problem 13-1: Incidence Matrices of a 2D Finite Element Mesh (220 min)
Problem 13-3: Whitney Finite Elements for Magnetostatic BVP in Two Dimensions (210 min)

Videos

0.0 Welcome and introduction (17min)

week (16.02.-20.02.):

1.2.1 Elastic Membranes (17 min)
1.2.2 Electrostatic Fields (6 min)
1.2.3 Quadratic Minimization Problems (21 min)
1.3 Sobolev spaces (22 min)
1.4 Linear Variational Problems (21 min) didn't quite understand

week (23.02.-27.02.):

1.5 Equilibrium Models: Boundary Value Problems (22 min)
1.6 Diffusion Models: Stationary Heat Conduction (5 min)
1.7 Boundary Conditions (07.35 min)
1.8 Second-Order Elliptic Variational Problems (13 min) even solved question :)
1.9 Essential and Natural Boundary Conditions (18 min)

week (02.03.-06.03):

2.2 Principles of Galerkin Discretization (25 min) review questions for 2.2.3 solved -> review questions for this chapter are important !!
2.3 Case Study: Linear FEMfor Two-Point Boundary Value Problems (24 min)
2.4 Case Study: Triangular Linear FEMin Two Dimensions I (28 min)
2.4 Case Study: Triangular Linear FEMin Two Dimensions II (32 min) confusing

week (09.03.-13.03):

2.5 Building Blocks of General Finite Element Methods (22 min) easy :)
2.6 Lagrangian Finite Element Methods (23 min)
2.7.2 Mesh Information and Mesh Data Structures (21 min) first LehrFM introduction
2.7.4 Assembly Algorithms (26 min)
2.7.5 Local Computations (26 min)
2.7.6 Treatment of Essential Boundary Conditions (21 min)

week (16.03.-20.03):

2.8 Parametric Finite Element Methods I (22 min)
2.8 Parametric Finite Element Methods II (27 min) ??
3.1 Abstract Galerkin Error Estimates (20 min)
3.2 Empirical (Asymptotic) Convergence of Lagrangian FEM (26 min) understandable

week (23.03.-27.03.):

3.3 A Priori (Asymptotic) Finite Element Error Estimates I (10 min)
3.3 A Priori (Asymptotic) Finite Element Error Estimates II (18 min)
3.3 A Priori (Asymptotic) Finite Element Error Estimates III (16 min)
3.4 Elliptic Regularity Theory (11 min)
3.5 Variational Crimes (15 min)
3.8 Validation and Debugging of Finite Element Codes (13 min)

week (30.03.-03.04.):

9.2.1 Heat Equation (10 min)
9.2.2 Heat Equation: Spatial Variational Formulation (8 min)
9.2.3 Stability of Parabolic Evaolution Problems (10 min)
9.2.4 Spatial Semi-Discretization: Method of Lines (6 min)
9.2.7 Timestepping for Method-of-Lines ODE (30 min) rewatch
9.2.8 Fully Discrete Method of Lines: Convergence (18 min) first 6min summary of chapter 9)

week (13.04.-17.04.):

10.1.1 Modeling Fluid Flow (7 min)
10.1.2 Heat Convection and Diffusion (5 min)
10.1.3 Incompressible Fluids (10 min)
10.1.4 Time-Dependent (Transient) Heat Flow in a Fluid (4 min)
10.2.1 Singular Perturbation (16 min)
10.2.2 Upwinding (13 min)
10.2.2.1 Upwind Quadrature (17 min)
10.2.2.2 Streamline Diffusion (19 min)
10.3.1 Method of Lines (13 min)

week (20.04.-24.04.):

12.1 Modeling the Flow of a Viscous Fluid (30 min) only 12min ??
12.2.1 The Stokes Equations: Constrained variational formulation (20 min) only 10min ??
12.2.2 The Stokes Equations: Saddle Point Problem Formulation (55 min) (30min), difficult video
12.2.3 Stokes System of Partial Differential Equations (15 min) (7 min)

week (27.04.-01.05.):

12.3.1 Pressure Instability (40 min) (20min)
12.3.2 Stable Galerkin Discretization of Stokes Saddle Point Problem (30 min) (16min)
12.3.3 Convergence of Stable FEM for Stokes (30 min) (17min)
12.3.4 The Taylor-Hood Finite Element Method (15 min) (8m)
12.3.5 The Non-Conforming Crouzeix-Raviart FEM (40 min) (22min)

week (04.05.-08.05.):

13.1.2 Fields and Integral Forms/Currents (28 min) (17min)
13.1.3 Differential Forms (56 min) (32min) ???
13.1.4 Exterior Calculus (70 min) (45min) did not understand a single word
13.1.5 Sobolev Spaces of Differential Forms (18 min) (11min)

week (11.05.-15.05.):

13.2.1 Cochain Calculus (67 min) (40min)
13.2.2 Whitney Forms (70 min) (40min)
13.3.1 Energies and Linear Material Laws (35 min) (21min)
13.3.2.1 Electromagnetic Wave Equations (45 min) (26min)
13.3.2.3 Method of lines for Electromagnetic Wave Equations, Timestepping for Semi-Discrete Electromagnetic Wave Equations (45 min) (30min)

week (18.05.-22.05.):

13.3.3.1 Magnetostatics: Formulations Based on Scalar Potentials (30 min) (18min)
13.3.3.2 Magnetostatics: Vector-Potential-Based Variational Formulations (40 min) (24min)
13.3.3.3-13.3.3.4 Saddle-Point Variational Problems, Discrete LBB-Conditions for a-Based Magnetostatics (60 min) (35min)

week (25.05.-29.05.): no videos

Q&A

#timestamp 2026-03-20

Galerking Orthogonality (3.1.3.2)

$∥ u - u_{h} ∥_{a}^{2} = a (u - u_{h}, u - u_{h}) = a (u - v_{h}, u - u_{h}) + \underset{= 0}{\underset{⏟}{a (v_{h} - u_{h}, u - u_{h})}} \forall v_{h} \in V_{0, h}$

bilinear form is coercive

continuity / boundedness: $| a (u, v) | \leq C ∥ u ∥_{v} ∥ v ∥_{v}$
coercivity: $γ ∥ u ∥_{v}^{2} \leq | a (u, u) |$

from this comes the Poincare inequality and uniform positive definite property

In contrast to Cea's Lemma, coercivity does not depend on the energy norm

∥ \cdot ∥_{a}

Tldr: strong to weak

$- (\partial_{1} \partial_{1} u + \partial_{1} \partial_{2} u + \partial_{2} \partial_{2} u) + ∥ x ∥^{2} u = 0 \in Ω, u = 1 \in \partial Ω$

(simplify strong form, then) convert to weak form

\begin{array}{r} - div ([\begin{array}{c} 1 & \frac{1}{2} \\ \frac{1}{2} & 1 \end{array}] \nabla u) + ∥ x ∥^{2} u = 0 ⟹ - div (A \nabla u) + β u = 0 \\ \int - div (A \nabla u) v + β u v d x = 0 \end{array}

choose suitable test function

e.g. $v \in C_{0}^{\infty} \subset H^{1} (Ω)$

IBP (use Green's formula)

\int_{Ω} A \nabla u \nabla v + β u v d v - \underset{0 by b.c.}{\underset{⏟}{\int_{\partial Ω} A \nabla u \cdot n v d l}} = 0

Incorporate b.c. (if needed)

for Dirichlet (essential) -> nothing to do

choose suitable fct. space

minimal suitable fct. space (often $H^{1}$ )
here, we need at list one derivative for $\nabla u$

\int A \nabla u \nabla v + β u v d x = 0 \forall v \in H_{0}^{1} (Ω)

#timestamp 2026-05-15

endterm

chapter 9 will be about discretising Method of Lines (butcher tableau)

for forms: see https://www.math.ucla.edu/~tao/preprints/forms.pdf

Vorlesung

1. Second-Order Scalar Elliptic Boundary Value Problems

Video 1.2.1: Elastic Membranes

Goal: compute shape of string (1D)/membrane (2D) under vertical loading

Notion 1.2.1.1. Configuration space

The configuration space of a mathematical model is a set, each of whose elements completely describe all relevant aspects of a state of the modeled physical system.

string must be continuous $⟺$ string not ripped apart
b.c.: u(a),u(b) given

Pasted image 20260216223006.png
Pasted image 20260216223153.png
Pasted image 20260216223325.png
-> but need formula for potential energy
Pasted image 20260216223452.png

Video 1.2.2 Eloectrostatic Fields (6min)

Pasted image 20260217102058.png
Pasted image 20260217102244.png
Note: $\hat{V} \subset C^{0} (\overset{―}{Ω})$

Equilibrium condition for electrostatics

The physical field will attain minimal electromagnetic field energy (1.2.2.16)

u_{*} = \underset{u \in {\hat{V}}_{E}}{argmin} J_{E} (u) .

Video 1.2.3 Quadratic Minimization Problems (21 min)

Pasted image 20260217102711.png

Definition 1.2.3.13. Quadratic minimization problem

A minimization problem

w_{*} = \underset{w \in V_{0}}{argmin} J (w)

is called a quadratic minimization problem, if $J$ is a =quadratic functional on a real vector space $V_{0}$ .

extension to affine space:

\hat{V} = \underset{"offset function"}{\underset{⏟}{μ_{0}}} + V_{0} \subset V

Pasted image 20260217103305.png
Pasted image 20260217103546.png
Recall from LA:

A = A^{⊤} s.p.d ⟹ (\vec{μ}, \vec{η}) \to {\vec{μ}}^{⊤} A \vec{η} \hat{=} inner product

Pasted image 20260217104011.png
Pasted image 20260217104355.png

Tldr

we are trying to find the most stable/lowest energy state of the system. The idea is to find a general formula for the structurally similar string/membrane/electricity problem, which can be described by a Quadratic functional:

J (v) = \underset{internal energy}{\underset{⏟}{\frac{1}{2} a (v, v)}} - \underset{external load}{\underset{⏟}{l (v)}} + c = \frac{1}{2} x^{⊤} A x - b^{⊤} x + c (a bilinear + symmetric, l linear)

However, we can only solve the problem if the energy formula has a minimum (imagine the energy formula as a paraboloid "bowl"):
Pasted image 20260217110348.png|300
Thus, for $\exists$ of solution, the problem must fulfill:

$a (v, v) > 0 \forall v \neq 0$ -> bilinear form must be at least semi-positive definite
$| l (v) | \leq C | | v | |_{a}$ , where $| | v | |_{a} = \sqrt{a (v, v)}$ energy norm. -> linear form $l (v)$ must be bounded/continuous with respect to energy norm. This ensures the "load" isn't so strong that it rips the system apart.

#todo add some integrals to the tldr

Video 1.3 Sobolev spaces (22 min)

have a bad reputation as having a useless mathematical sophistication

Pasted image 20260217113247.png
choose

V_{0} := {functions v on Ω : a (v, v) < \infty}

=> energy space generated by $a (\cdot, \cdot)$
Pasted image 20260217113422.png
Pasted image 20260217113818.png
in contrast to $L^{2}$ -Space, here the boundary conditions are not dropped; the reason for this are the norms:
Pasted image 20260217114017.png
Sobolev space without boundary conditions:
Pasted image 20260217114124.png
the $L^{2}$ can be bounded by the size of the domain and the energy norm:

Theorem 1.3.4.17: First Poincaré-Friedrichs inequality

If $Ω \subset R^{d}$ , $d \in N$ , is bounded, then

∥ u ∥_{0} \leq diam (Ω) ∥ grad u ∥_{0} \forall u \in H_{0}^{1} (Ω) .

$⟹ | ℓ (v) | \leq \underset{C}{\underset{⏟}{∥ f ∥_{0} diam (Ω)}} ∥ v ∥_{a} \forall v \in H_{0}^{1} (Ω)$

Pasted image 20260217114451.png

How to "work with" Sobolev spaces

Most concrete results about Sobolev spaces boil down to relationships between their norms. The spaces themselves remain intangible, but the norms are very concrete and can be computed and manipulated as demonstrated above.

Tip

Do not be afraid of Sobolev spaces!
It is only the norms that matter for us, the 'spaces' are irrelevant!

Corollary 1.3.4.26.

H^{1}

-norm of piecewise smooth functions

Under the assumptions of Thm. 1.3.4.23 we have for a continuous, piecewise smooth function $u \in C^{0} (Ω)$

{| u |}_{H^{1} (Ω)}^{2} = | u |_{H^{1} (Ω_{1})}^{2} + | u |_{H^{1} (Ω_{2})}^{2} = \int_{Ω_{1}} | grad u (x) |^{2} d x + \int_{Ω_{2}} | grad u (x) |^{2} d x .

If $u \in C^{0} (Ω) \Rightarrow$ compute $H^{1}$ -norm piecewise

Tldr

$L^{2}$ space - space of all functions where $\int f^{2} < \infty$ . Here, b.c. are meaningless, since changing the function at a single point costs zero energy
$H^{1}$ space - "energy space" for membranes; includes functions where derivatives are square-integrable. Unlike $L^{2}$ , the b.c. cannot be ignored, since changing a point blows up the derivative
$H_{0}^{1}$ -space - $H^{1}$ with zero b.c.
- one can think of $H^{1}$ as the space of peicewise smooth, globally continuous functions

Tldr

Functions of an $L^{2}$ space are square-integrable on $Ω$ :

\int_{Ω} | u (x) |^{2} d x < \infty

Functions of the $H^{1}$ space are in $L^{2}$ , too. Additionally, the first (hence ^1) derivative is in $L^{2}$ , too:

H^{1} (Ω) = {u \in L^{2} (Ω) | \nabla u \in L^{2} (Ω)}

Note that $H^{1}$ functions do not need to be in $C^{1}$ or even $C^{0}$ ! A counterexample is

u (x, y) = {\begin{cases} 0, & x \leq \frac{1}{2} \\ 1, & x > \frac{1}{2} \end{cases} in Ω = [0, 1]^{2}

Video 1.4 Linear Variational Problems (21 min)

Pasted image 20260219081701.png
Pasted image 20260219082140.png
Pasted image 20260219082308.png
Poincare-friedrich?
Pasted image 20260219082838.png

Video 1.5 Equilibrium Models: Boundary Value Problems (22 min)

Pasted image 20260224174310.png
Pasted image 20260224175303.png
=> 2-pt BVP (a Dirichlet problem) in strong form (vs. "weak form" = LVP)

Pasted image 20260224175546.png
Pasted image 20260224175607.png

Lemma 1.5.3.4: Fundamental lemma of calculus of variations in higher dimensions

If $f \in L^{2} (Ω)$ satisfies

\int_{Ω} f (x) v (x) d x = 0 \forall v \in C_{0}^{\infty} (Ω),

then $f \equiv 0$ can be concluded.

Pasted image 20260224180021.png
Pasted image 20260224180509.png

Video 1.6 Diffusion Models: Stationary Heat Conduction (5 min)

Pasted image 20260224182905.png
Pasted image 20260224183014.png
by combining those two using Gauss theorem, we get:
Pasted image 20260224183120.png
=> we reach at the same PDE as for the membrane model

Video 1.7 Boundary Conditions (07.35 min)

Pasted image 20260224183718.png
example of a mixed BVP:
Pasted image 20260224184006.png

Boundary conditions for second-order elliptic BVPs

For second order elliptic boundary value problems exactly one boundary condition is needed on every part of $\partial Ω$ .

Video 1.8 Second-Order Elliptic Variational Problems (13 min)

Pasted image 20260224190120.png
for Dirichlet BVP:
Pasted image 20260224213109.png
for Radiation BVP:
Pasted image 20260224213302.png
for Neumann BVP:
Pasted image 20260224214942.png
Pasted image 20260224215319.png
Pasted image 20260224215332.png

Theorem 1.8.0.20: Second Poincaré-Friedrichs inequality

If $Ω \subset R^{d}$ , $d \in N$ , is bounded, then

\exists C = C (Ω) > 0 : ∥ u ∥_{0} \leq C diam (Ω) ∥ grad u ∥_{0} \forall u \in H_{*}^{1} (Ω) .

Video 1.9 Essential and Natural Boundary Conditions (18 min)

summary of 1.4 - 1.9:
Pasted image 20260225081041.png
Pasted image 20260225081059.png
Pasted image 20260225081424.png
Pasted image 20260225082319.png
Pasted image 20260225082332.png

Tldr

For a summary of this week, see the notes from the exercise session:
-> ETH_NumericalMethodsForPDE_Übungen#^75b794

2. Finite Element Methods (FEM)

Video 2.2 Principles of Galerkin Discretization (25 min)

graph LR
	A[Variational boundary
value problems] -- discretization --> B["System of a finite number of equations for (real) unknowns"]

Pasted image 20260301110814.png

garlekin discretization in two steps:

$V_{0} \to V_{0, h} \subset V, \dim V_{0, h} < \infty$

Galerkin discretization (step 1)

Replace the infinite-dimensional function space $V_{0}$ in the linear variational problem (2.2.0.2) with a finite-dimensional subspace $V_{0, h} \subset V_{0}$ .

-> "h" $\hat{=}$ tag for s.th. finite

Pasted image 20260301111127.png
Pasted image 20260301111209.png

Choose basis, then plug basis expansions int DVP, and use (bi-linearity)

Pasted image 20260301111355.png
Pasted image 20260301111722.png|697

Theorem 2.2.2.7: Independence of Galerkin solution of choice of basis

The choice of the basis $B_{h}$ has no impact on the (set of) Galerkin solutions $u_{h}$ of (2.2.1.1).

Pasted image 20260301111839.png
Pasted image 20260301112119.png

Properties all Garlekin matrices for a DVP (all matrices from a congruence clas) share:

regularity (invertibility)
symmetry
pos.dev.

Garleking pos. dev / symm. ⟺ a (\cdot, \cdot) s.p.d on V_{0, h}

Pasted image 20260301112433.png

Video 2.3 Case Study: Linear FEMfor Two-Point Boundary Value Problems (24 min)

u^{″} = f in [a, b] u (a) = u (b) = 0

NOTE: $f \in L^{2} (Ω)$ is given inprocedural form, $f \in C^{0} ([a, b])$
Pasted image 20260301112821.png

Idea underlying the finite element method: approximate $u$ by means of continuous, piecewise polynomial functions on a mesh:

Pasted image 20260301113228.png
Pasted image 20260301113412.png
Pasted image 20260301114959.png

Video 2.4 Case Study: Triangular Linear FEMin Two Dimensions I (28 min)

In 1D, we partition the domain $Ω$ into an 1D Mesh
In 2D, we use triangulation

Pasted image 20260301182551.png

1. means: no hanging nodes (nodes that are on the edge of another triangle)

Pasted image 20260301183038.png

$ϵ (M) =$ set of edges
$M = {K_{1}, \dots, K_{M}}, M \in N$
$v (M) = {{\vec{x}}^{1}, \dots, {\vec{x}}^{N}}, N \in N$

test space in 2 (space of M-p.w. linear continuous function on $Ω$ ):
Pasted image 20260301183355.png

tent function alternative in 2D:
Pasted image 20260301184029.png

In 2D, the Galerking matrix does not have a nice structure like in 1D; however, it is still sparse, since nodes unconnected by edges result in a 0 in the Matrix:
Pasted image 20260301184700.png

Video 2.4 Case Study: Triangular Linear FEMin Two Dimensions II (32 min)

Galerkin discretization based on a nodal basis of tent functions and a triangular mesh:
Pasted image 20260301203218.png
Pasted image 20260301203646.png
there are two approaches to creating the Galerkin-matrix:

assemble matrix from all vertices

Problems:

Edges not available in data structure (see first 2.4 video)
Triangles adjacent to nodes not available

In my opinion, it is not fair to compare it with the cell-oriented approach, when the data class is also cell-oriented

cell-oriented approach, by "distribute & add"

for finite-element algorithms, always do it cell-oriented

Pasted image 20260301205155.png
dof: degrees of freedom
Pasted image 20260301205231.png

now, we assemble everything:
Pasted image 20260301205651.png
for the RHS, we can compute the integral using composite trapezoidal rule (we are working with triangles anyways):
Pasted image 20260301205803.png

Tldr

We are trying to get variational problems from infinite-dimensional function spaces to something discrete by

choosing a finite-dimensional subspace $V_{0, h} \subset V_{0}$ , where the bilinear and linear form $a (\cdot, \cdot), l (\cdot)$ remain well-defined
choosing an ordered basis ${b_{1}, \dots, b_{N}}$ for the subspace $V_{0, h}, \dim (V_{0, h}) = N$ , and using the basis to convert problem into a linear system

In those two steps, the problem can be transformed into a LSE

A μ = ϕ

Tldr - Step 1

In the first step of the discretization, our goal is to replace space $V_{0}, \dim (V_{0}) = \infty$ with finite $V_{0, h} \subset V_{0}$ . The space should preserve the bilinear form $a (\cdot, \cdot)$ and linear form $l (\cdot)$ . The new problem is:

find u_{h} \in V_{0, h} s.t. a (u_{h}, v_{h}) = l (v_{h}) \forall v_{h} \in V_{0, h}

define Mesh $M$ . Normally, we use triangulation - pay attention that the mesh does not have "hanging nodes" (2.4.1)
choose a "building block" function to use for each cell. Normally, pecewise linear functions are used, evtl. quadratic or cubic polynomials

ensure $C^{0}$ -> do not use piecewise constants, since then $V_{0, h} ⊄ V_{0}$ , e.g. if $V_{0} = H^{1}$

incoroporate boundary conditions, in case of Dirichlet by discarding any basis functions with nodes on the boundary

Tldr - Step 2

4. choose a basis $B_{h} = {b_{1}^{h}, \dots, b_{N}^{h}}$ , then write every $v_{h} \in V_{0, h}$ as linear combination of basis functions:

u_{h} = \sum_{j = 1}^{N} μ_{j} b_{j}^{h}, v_{h} = \sum_{i = 1}^{N} ν_{i} b_{i}^{h}

convert to a linear system by exploiting bilinearity of $a (\cdot, \cdot)$ and linearity of $l (\cdot)$ :

\sum_{j = 1}^{N} μ_{j} a (b_{j}^{h}, b_{i}^{h}) = l (b_{i}^{h}) for i = 1, \dots, N

⟺ A μ = ϕ

$A$ : Galerkin/stiffness matrix, $A_{i j} = a (b_{j}^{h}, b_{i}^{h})$ , how basis functions interact through the bilinear form
$ϕ$ : load vector, $ϕ_{i} = l (b_{i}^{h})$
$μ$ : solution vector; the final Galerkin solution $u_{h} = \sum b_{i}^{h} \cdot μ_{i}$ (basis, weighted by solution vector)

Basis choice

no impact on Galerking solution; however:
basis functions with small local support (e.g. tent function) results in a sparse matrix -> desired outcome
basis choice determines condition number of matrix -> sensitivity to rounding errors

Video 2.5 Building Blocks of General Finite Element Methods (22 min)

meshes
2D: {cells, edges, nodes}=(mesh/geometric) entities

polynomials
Pasted image 20260306093818.png
1st option: e.g.

p = 2 : P_{2} (R^{2}) = span {1, x, y, x^{2}, x y, y^{2}}

Definition 2.5.2.7: Tensor product polynomials

Space of tensor product polynomials of degree $p \in N$ in each coordinate direction

Q_{p} (R^{d}) := {x \mapsto \sum_{ℓ_{1} = 0}^{p} \dots \sum_{ℓ_{d} = 0}^{p} c_{ℓ_{1}, \dots, ℓ_{d}} x_{1}^{ℓ_{1}} \dots x_{d}^{ℓ_{d}}, c_{ℓ_{1}, \dots, ℓ_{d}} \in R}

= Span {x \mapsto p_{1} (x_{1}) \dots p_{d} (x_{d}), p_{i} \in P_{p} (R), i = 1, \dots, d} .

Lemma 2.5.2.8: Dimension of spaces of tensor product polynomials

$\dim Q_{p} (R^{d}) = (p + 1)^{d} for all p \in N_{0}, d \in N$

2nd option: e.g.

Q_{2} (R^{2}) = span {\begin{matrix} 1, x, x^{2} \\ y, y x, y x^{2} \\ y^{2}, y^{2} x, y^{2} x^{2} \end{matrix}}

basis functions (also called global shape functions (GSF), degrees of freedom)

Locally supported basis functions

Third main ingredient of FEM: locally supported basis functions (see Section 2.2 for role of bases in Galerkin discretization)

Basis functions $b_{h}^{1}, \dots, b_{h}^{N}$ for a finite element trial/test space $V_{0, h}$ built on a mesh $M$ must satisfy:

( $B_{1}$ ) $B_{h} := {b_{h}^{1}, \dots, b_{h}^{N}}$ is basis of $V_{0, h} ⟹ N = \dim V_{0, h}$ ,
( $B_{2}$ ) each $b_{h}^{i}$ is associated with a single geometric entity (cell/edge/face/vertex) of $M$ ,
( $B_{3}$ ) $supp (b_{h}^{i}) = ⋃ {\overset{―}{K} : K \in M, p \subset \overset{―}{K}}$ , if $b_{h}^{i}$ associated with cell/edge/face/vertex $p$ .

we want local support:

if $supp (b_{h}^{i})$ , $supp (b_{h}^{j})$ don't overlap, entries in Galerkin matrix $(A)_{i j} = 0$ => sparse matrix

generalization of barycentric coordinate functions:
Pasted image 20260306094819.png

Video 2.6 Lagrangian Finite Element Methods (23 min)

Pasted image 20260310115750.png
one triangle, with Kardinal property:
Pasted image 20260310115836.png
gluing triangles (interpolation nodes):
Pasted image 20260310115922.png
-> globally continuous piecewise function

for more dimensions, we "subdivide"; then, we might also have interpolation nodes inside the shape
Pasted image 20260310120326.png
for hybrid mesh:
Pasted image 20260310120410.png
works, becuase functions $\in Q_{1} / P_{1}$ are all linear when restricted to straigt edge (which we have)

works for any degree!

Video 2.7.2 Mesh Information and Mesh Data Structures (21 min)

problems with FEM:
Pasted image 20260310121851.png

// factory object creates mesg entities auto of 2D mesh data structure from information contained in mesh_file
auto factory = std::make_unique<lf::mesg::hybrid2d::MeshFactory>(2);
lf::io::GmshReader readermove(factory), mesh_file.string();

we need

container (tranversal + ordering/indexing of entities)
topology (incidence/adjecency)
geometry (node locations, entity shapes)

Pasted image 20260310122228.png

mesh

index access
loops:

for(const lf::mesh::Entity &entity : mesh.Entitiy(codimension));

topology layer
Pasted image 20260310122516.png
Pasted image 20260310122650.png
access:
Pasted image 20260310122726.png

auto sub_entities_range = entity.SubEntities(sub_codimension);
for(const: lf::mesh::Entity &subentity : sub_entities_range)
    mesh.Index(subentity);

geometry layer

geometry object attached for every entity interface

const lf::geometry::Geometry *geo_ptr = entity.Geometry();
Eigen::MatrixXd corners = lf::geometry::Corners(*geo_ptr);

Video 2.7.4 Assembly Algorithms (26 min)

continuation of 2.4

assembly algorithms: building of FE Galerkin amtrix/r.h.s. vector from local contributions

we can split the integrals from the weak form into sums of local contributions (see exercise lession, week3, thu, hongg)
Pasted image 20260310132555.png
Pasted image 20260310133311.png
Pasted image 20260310133322.png
LehrFM numbering rules:
Pasted image 20260310133410.png
Pasted image 20260310133459.png
Pasted image 20260310133904.png

dimension of Space
Number of local Dofs (number of local shape functions)
array of locglobmap for all local indices of entity
number of local/global shape functions associatet with entity
global indices of associated local shape functions of indices (not discussed until now)
global index of entity associated with dof

Initalize dynamic Dof Handler:
Pasted image 20260310134015.png
for lagr. FE: uniform no of GSFs associated with entities of same type -> initalization can be simplified:

using dof_map_t = std::map<lf::base::RefEl, base::size_type>;
UniformFEDofHandler(std::shared_ptr<const lf::mesh::Mesh> mesh_p,
                    dof_map_t dofmap);

example:
Pasted image 20260310134056.png
distribute pattern used for assembly

loops over entities + DofHandler queries + local operation

not copied: code 2.7.4.23: Assembly function of LehrFEM++

Video 2.7.5 Local Computations (26 min)

either by analytic computation (only if we have polynomials -> locally constant functions on each cell of the mesh) or

analytic computation:

Lemma 2.7.5.5. Integration of powers of barycentric coordinate functions

For any non-degenerate $d$ -simplex $K$ with barycentric coordinate functions $λ_{1}, \dots, λ_{d + 1}$ and exponents $α_{j} \in N, j = 1, \dots, d + 1$ ,

\begin{matrix} (2.7.5.6) & \int_{K} λ_{1}^{α_{1}} \cdot \dots \cdot λ_{d + 1}^{α_{d + 1}} d x = d! | K | \frac{α_{1}! α_{2}! \cdot \dots \cdot α_{d + 1}!}{(α_{1} + α_{2} + \dots + α_{d + 1} + d)!} \forall α \in N_{0}^{d + 1} . \end{matrix}

$↪$ depends only on area of $K$

e.g.
Pasted image 20260310135246.png
local quadrature

in 1D: affine mapping $[a, b] \to [- 1, 1]$ + transformation formula for integrals
local quadrature rules used in finite element methods are botained by transformation from (a few) local quadrature rules defined on reference elements.
in 2D:

(affine mapping maps lines to lines)

general formula:
Pasted image 20260310135653.png
example on $K \in M$ :
Pasted image 20260310135938.png
however, if $ϕ_{K}$ affine => $det (D ϕ_{K}) \equiv c o n s t .$ , and we get:
Pasted image 20260310140013.png
in LehrFEM, we can use:
Pasted image 20260310145530.png

Definition 2.7.5.29: Order of a local quadrature rule

A local quadrature rule according to Def. 2.7.5.9 is said to be order $q \in N$ , if

for a simplex $K$ (triangle, tetrahedron) it is exact for all polynomials $f \in P_{q - 1} (R^{d})$ ,
for a tensor product element $K$ (rectangle, brick) it is exact for all tensor product polynomials $f \in Q_{q - 1} (R^{d})$ .

"order $↑ ⟹$ no. of nodes $↑$ "

transformation on affine maps preserves polynomials

-> integral transformation doesn't change order of polynomials

Pasted image 20260310145821.png
quadrature rules in LehrFEM: reference object shape + desired order
Pasted image 20260310145901.png

not copied: code 2.7.5.41: Entity-based composite numerical quadrature

Video 2.7.6 Treatment of Essential Boundary Conditions (21 min)

#todo understand offset function:

Pasted image 20260310150708.png
Pasted image 20260310151152.png
Pasted image 20260310152333.png

maybe missed something relevant...

Pasted image 20260310152411.png

using LehrFM:
Pasted image 20260310152445.png

not copied code 2.7.6.17: Transformation function, which generates the block matrix from above

Tldr: LehrFEM++ example for dirichlet b.c. problem

Problem: $- Δ u = f$
Weak form: $\int_{Ω} \nabla u \nabla v d x = \int_{Ω} f v d x, \forall v \in H_{0}^{1} (Ω)$

Boundary conditions: Dirichlet (essential)

// 1. Load Mesh and Setup DofHandler for Linear FE (1 DOF per node)
auto factory = std::make_unique<lf::mesh::hybrid2d::MeshFactory>(2);
lf::io::GmshReader readermove(factory), "mesh.msh";
auto mesh_p = reader.mesh();
lf::assemble::UniformFEDofHandler dofh(mesh_p, {{lf::base::RefEl::kPoint(), 1}});

// 2. Prepare Global System Containers (N x N)
lf::assemble::COOMatrix<double> A_triplets(dofh.NumDofs(), dofh.NumDofs());
Eigen::VectorXd rhs = Eigen::VectorXd::Zero(dofh.NumDofs());

// 3. Assemble Stiffness Matrix (using built-in Laplace Provider)
lf::uscalfe::LinearFELaplaceElementMatrix elmat_builder;
lf::assemble::AssembleMatrixLocally(0, dofh, dofh, elmat_builder, A_triplets);

// 4. Handle Dirichlet Boundary Conditions (Offset Function)
auto bd_flags = lf::mesh::utils::flagEntitiesOnBoundary(mesh_p, 2); // Flag nodes
auto selector = [&](unsigned int idx) -> std::pair<bool, double> {
 auto node = dofh.Entity(idx);
 return {bd_flags(*node), 0.0}; // Fix boundary nodes to 0.0
};
lf::assemble::FixFlaggedSolutionCompAlt(selector, A_triplets, rhs);

// 5. Solve
Eigen::SparseMatrix<double> A = A_triplets.makeSparse();
Eigen::SparseLU<Eigen::SparseMatrix<double>> solver(A);
if(solver.info() != Eigen::Success) abort();
Eigen::VectorXd solution = solver.solve(rhs);

codim: 0->cells, 1->edges, 2->nodes
COOMatrex avoids memory shifts during assembly
AssembleMatrixLocally is templated -> works independent of test/trial space

Video 2.8 Parametric Finite Element Methods I (22 min)

remaining challenge: what to do with general mesh (curved edges, triangles + quadrilaterals)

Affine equivalence
Pasted image 20260314201011.png
affine pullback preserves polynomials:

Lemma 2.8.1.3: Preservation of polynomials under affine pullback

If $Φ : R^{d} \mapsto R^{d}$ is an affine (linear) transformation ( $\to$ Def. 2.7.5.13), then

Φ^{*} (P_{p} (R^{d})) = P_{p} (R^{d}) and Φ^{*} (Q_{p} (R^{d})) = Q_{p} (R^{d}) .

general quadrilateral lagrangian FE

so far: GSF -> LSF #todo: look this up again
now: start from LSF ${\hat{b}}^{j}$ on reference element $\hat{K}$ -> inverse pullback:

b_{K}^{j} := (ϕ_{K}^{- 1})^{⋆} {\hat{b}}^{j}

-> gluing: GSF

bilinear transformation: componentwise (tensor product) polynomial $\deg = 2$
Pasted image 20260314202207.png

gluing:
-> local shape functions are no longer bilinear, since $ϕ_{K}^{- 1}$ is no longer polynomial (rational function)

Non-polynomial nature of bilinear shape functions

The local shape functions $b_{K}^{i}$ defined by (2.8.2.4), where $Φ_{K}$ is a bilinear transformation and ${\hat{b}}^{i}$ are the bilinear local shape functions on the unit square, are not polynomial in general.

Video2.8 Parametric Finite Element Methods II (27 min)

Pasted image 20260314211639.png
Note that possibility of gluing needs to be checked!

#todo: what are we doing here?

[...]
Pasted image 20260314213531.png
Pasted image 20260314213641.png

Information defining a

H^{1} (Ω)

-conforming parametric finite element

In order to define a $H^{1} (Ω)$ -conforming parametric finite element in a finite element code we have to implement for all reference elements and $i = 1, \dots, Q$

a function realizing $\hat{x} \mapsto {\hat{b}}^{i} (\hat{x}), x \in \hat{K}$ ,
and the evaluation $\hat{x} \mapsto grad {\hat{b}}^{i} (\hat{x}), x \in \hat{K}$ .

// interface in LF++:
size_type ScalarReferenceFiniteElement::NumRefShapeFunctions() const;
size_type ScalarReferenceFiniteElement::NumRefShapeFunctions(dim_t codim) const;
size_type ScalarReferenceFiniteElement::NumRefShapeFunctions(dim_t codim, sub_idx_t subidx) const;

Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
ScalarReferenceFiniteElement::EvalReferenceShapeFunctions(
    const Eigen::MatrixXd& refcoords) const;

Eigen::Matrix<SCALAR, Eigen::Dynamic, Eigen::Dynamic>
ScalarReferenceFiniteElement::GradientsReferenceShapeFunctions(
    const Eigen::MatrixXd& refcoords) const;

didn't copy code 2.8.3.29: Eval() method for parametric element matrix provider

Boundary Approximation by Parametric FEM:
Pasted image 20260314214831.png
"6-node triangle" in LehrFEM++:
Pasted image 20260314215207.png

Degeneracy of the Transformation

$Φ_{K}$ may be "degenerate" $⟺ D Φ_{K} (\hat{x})$ singular for some $\hat{x} \in \hat{K}$ .

Question: why is this better than simply splitting the curve into multiple triangles? The approximation wouldn't be much worse than by this approach.

Tldr: Pullback

Let $K$ be a cell, $\hat{K}$ a reference cell. Then we define the bijection $Φ_{K}$

\begin{aligned} Φ_{K} : \hat{K} & \to K \\ Φ_{K} (\hat{x}) & = x for \hat{x} \in \hat{K}, x \in K \\ Φ_{K}^{- 1} (x) & = \hat{x} (the inverse mapping) \\ u (x) & = (Φ_{K}^{- 1})^{*} \hat{u} (x) = \hat{u} (Φ_{K}^{- 1} (x)) \\ \hat{u} (\hat{x}) & = (Φ^{⋆} u) (\hat{x}) = u (Φ (\hat{x})) = u (x) \end{aligned}

For Gradients,

\begin{aligned} (2.8.3.11) & {grad}_{\hat{x}} \hat{u} (\hat{x}) & = (D Φ (\hat{x}))^{⊤} ({grad}_{x} u) (Φ (\hat{x})) \\ {grad}_{x} u (x) & = (D Φ (\hat{x}))^{- ⊤} {grad}_{\hat{x}} \hat{u} (\hat{x}) \end{aligned}

#todo what is the difference between $Φ^{- 1}$ and $Φ^{⋆}$ ? Is this correct:

Notation	Name	Acts On	Direction
$Φ$	Mapping	Points (Coordinates)	$\hat{K} \to K$
$Φ^{- 1}$	Inverse Mapping	Points (Coordinates)	$K \to \hat{K}$
*$Φ^{}$**	Pullback	Functions / Forms	$K \to \hat{K}$
*$(Φ^{- 1})^{}$**	Inverse Pullback	Functions / Forms	$\hat{K} \to K$

3. FEM: Convergence and Accuracy

Video 3.1 Abstract Galerkin Error Estimates (20 min)

convergence (how fast does $u_{h}$ get closer to $μ$ ?) and accuracy (how close is $u_{h}$ to $μ$ ?)

Discretization error $u - u_{h} \in V_{0}$ (a function in BVP context!)

Our goal:

Bound relevant norm of discretization error $u - u_{h}$

A relevant norm:

LVP $⟺ u = \underset{v \in V_{0}}{argmin} J (v), J (v) := \frac{1}{2} a (v, v) - ℓ (v)$ (energy function)

Lemma 3.1.2.3. Energy norm and energy deviation

Let $u \in V_{0}$ solve the linear variational prolem (2.2.0.2),

\begin{matrix} (2.2.0.2) & u \in V_{0} : a (u, v) = ℓ (v) \forall v \in V_{0} \end{matrix}

and let Ass. 3.1.1.2–Ass. 3.1.1.4 be satisfied. Then, with $J (v) := \frac{1}{2} a (v, v) - ℓ (v)$ ,

\begin{matrix} (3.1.2.4) & J (w) - J (u) = ∥ w - u ∥_{a}^{2} \forall w \in V_{0} \end{matrix}

galerkin orthogonality:
Pasted image 20260315120409.png
from this follows Pythagoras theorem:
Pasted image 20260315121030.png
Pasted image 20260315121044.png
or shorter:

Theorem 3.1.3.7. Cea's lemma

Under Ass. 3.1.1.2–Ass. 3.1.1.4 the energy norm of the Galerkin discretization error for (2.2.0.2) satisfies

∥ u - u_{h} ∥_{a} = inf_{v_{h} \in V_{0, h}} ∥ u - v_{h} ∥_{a} .

-> Galerking solution is optimal (in energy)

minimizes $J (v)$ over $V_{0, h}$

Optimality of Galerkin solutions:

$\begin{matrix} (3.1.3.9) & \underset{((energy) norm of) discretization error}{\underset{⏟}{∥ u - u_{h} ∥_{a}}} = \underset{best approximation error in V_{0, h}}{\underset{⏟}{inf_{v_{h} \in V_{0, h}} ∥ u - v_{h} ∥_{a}}}, \end{matrix}$

-> As regards the energy norm, the Galerkin solution is the best possible solution we can obtain in a given trial space.

Galerking error estimates in energy form approximation error estimates
we can improve the result by enlarging the space (refinement)

Pasted image 20260315121738.png

2d triangular: regular refinement by splitting triangle into 4 congruent triangles
use lf::refinement in LF++

Video 3.2 Empirical (Asymptotic) Convergence of Lagrangian FEM (26 min)

Asymptotic Convergence

Convergence: asymptotic perspective

Crucial: our notion of convergence is asymptotic!
sequence of discrete models $⟹$ sequence of approximate solutions $(u_{h}^{(i)})_{i \in N}$
$↕$ $⟹$ study sequence $(∥ u_{h}^{(i)} - u ∥)_{i \in N}$ as $i \to \infty$
created by variation of a discretization parameter.

Definition 3.2.1.4. Mesh width

Given a mesh $M = {K}$ , its mesh width $h_{M}$ is defined as

h_{M} := max {diam K : K \in M}, diam K := max {| p - q | : p, q \in K} .

limit $h_{m}$ -> 0
for p-refinement, discretication DP $\hat{=}$ degree p

cost might depend on dimension $V_{0, h}$ (or $n n z (A)$ )

algebraic and exponential convergence

Definition 3.2.2.1. Types of convergence

\to

$\begin{aligned} ∥ u - u_{N} ∥ & = O (N^{- α}), α > 0 & ⟺ algebraic convergence with rate α \\ ∥ u - u_{N} ∥ & = O (\exp (- γ N^{δ})), with γ, δ > 0 & ⟺ exponential convergence \end{aligned}$

(asymptotically for $N \to \infty$ )

log-log plot with straigt line: algebraic convergence
lin-log plor with line: exponential convergence

-> (integral) error norms can usually be computed only by numerical quadrature (/sampling) of sufficiently high order

Example:

alg. cvg., but slow convergence for non-smooth solutions (h-refinement)
- $L^{2}$ cvg. seems to be faster than $∥ \cdot ∥_{a}$ -cvg (?)
exp. cvg.(p-refinement)

-> there might be pre-asymptotic behaviour (no convergence until mesh fine / polynomial degree high)

Tldr: h- vs p- refinement

Feature	h-Refinement	p-Refinement
Strategy	Subdivide mesh cells ( $h \to 0$ )	Increase local polynomial degree ( $p \to \infty$ )
Typical Convergence	Algebraic: $O (h^{p})$	Exponential: $O (\exp (- γ p^{δ}))$ for analytic solutions
Limiting Factor	Polynomial degree $p$ and smoothness $k$	Smoothness $k$ of the exact solution
Space Property	Often results in nested spaces (e.g., via regular refinement)	Inherently leads to nested finite element spaces

Video 3.3 A Priori (Asymptotic) Finite Element Error Estimates I (10 min)

Foundation: Cea's lemma

∥ u - u_{h} ∥_{a} = inf_{v_{h} \in V_{0, h}} ∥ u - v_{h} ∥ \leq ∥ u - w_{h} ∥_{a} \forall w_{h} \in V_{0, h}

a suitable "candidate function" $w_{h}$ obtained by interpolation

Note: For 2nd-order elliptic Dirichlet BVP: $∥ v ∥_{a} \sim | v |_{H^{1} (Ω)}$

Global estimates for norms

$\begin{matrix} (3.3.1.14) & | u - I_{1} u |_{H^{1} ([a, b])} \leq h_{M} ∥ u^{″} ∥_{L^{2} ([a, b])} \end{matrix}$

-> algebraic cvg., rate 1

\begin{matrix} (3.3.1.11) & ∥ u - I_{1} u ∥_{L^{2} ([a, b])} \leq h_{M}^{2} ∥ u^{″} ∥_{L^{2} ([a, b])} \end{matrix}

-> algebraic cvg. rate 2

\begin{matrix} (3.3.1.4) & ∥ u - I_{1} u ∥_{L^{\infty} ([a, b])} \leq \frac{1}{4} h_{M}^{2} \underset{smoothness requirement}{\underset{⏟}{∥ u^{″} ∥_{L^{\infty} ([a, b])}}} \end{matrix}

-> algebraic cvg. rate 2

Video 3.3 A Priori (Asymptotic) Finite Element Error Estimates II (18 min)

interpolation errors in order to bound the energy norm of the discretization error for lagrangian FEM
still looking for piecewise linear FE, but now in 2D

Definition 3.3.2.1. Linear interpolation in 2D

The linear interpolation operator $I_{1} : C^{0} (\bar{Ω}) \mapsto S_{1}^{0} (M)$ is defined by

I_{1} u \in S_{1}^{0} (M), I_{1} u (p) = u (p) \forall p \in V (M) .

uniquely fixes element of piecweise FE space

Pasted image 20260321201604.png

Definition 3.3.2.20. Shape regularity measure

For a simplex $K \in R^{d}$ we define its shape regularity measure as the ratio

ρ_{K} := h_{K}^{d} : | K |

and the shape regularity measure of a simplicial mesh $M = {K}$

ρ_{M} := max_{K \in M} ρ_{K} .

$ρ_{K}$ is similarity-invariant (rotation, scaling, translation)
depends only on it's smallest angle
expresses distortion compared to equilateral triangle $△$ :

Pasted image 20260321202155.png

Theorem 3.3.2.21. Error estimate for piecewise linear interpolation

For any $u \in C^{2} (\bar{Ω})$ and 2D piecewise linear interpolation $I_{1} : C^{0} (\bar{Ω}) \to S_{1}^{0} (M)$ , $M$ a triangular mesh, holds

∥ u - I_{1} u ∥_{L^{2} (Ω)} \leq \sqrt{\frac{3}{8}} h_{M}^{2} ∥ | D^{2} u |_{F} ∥_{L^{2} (Ω)},

-> alg. cvg. rate 2

∥ grad (u - I_{1} u) ∥_{L^{2} (Ω)} \leq \sqrt{\frac{3}{24}} ρ_{M} h_{M} ∥ | D^{2} u |_{F} ∥_{L^{2} (Ω)} .

-> alg. cvg. rate 1
where $h_{M}$ denotes the mesh width ( $\to$ Def. 3.2.1.4) and $ρ_{M}$ the shape regularity measure ( $\to$ Def. 3.3.2.20) of $M$ .

#todo learn multi-index notation

Definition 3.3.3.1. Higher order Sobolev spaces/norms

The $m$ -th order Sobolev norm, $m \in N_{0}$ , for $u : Ω \subset R^{d} \mapsto R$ (sufficiently smooth) is defined by

∥ u ∥_{H^{m} (Ω)}^{2} := \sum_{k = 0}^{m} \sum_{α \in N^{d}, | α | = k} \int_{Ω} | D^{α} u |^{2} d x, where D^{α} u := \frac{\partial^{| α |} u}{\partial x_{1}^{α_{1}} \dots \partial x_{d}^{α_{d}}} .

Sobolev space $H^{m} (Ω) := {v : Ω \mapsto R :== ∥ v ∥_{H^{m} (Ω)} < \infty ==}$ .

$L^{2} (Ω) = H^{0} (Ω) \supset H^{1} (Ω) \supset H^{2} (Ω) \supset H^{3} (Ω) \supset \dots$
$\overset{}{\to}$
increasing regularity of functions

$∥ \cdot ∥_{H^{m}}$ measures "ruggedness" of functions

Definition 3.3.3.3. Higher order Sobolev semi-norms

The $m$ -th order Sobolev semi-norm, $m \in N$ , for sufficiently smooth $u : Ω \mapsto R$ is defined by

| u |_{H^{m} (Ω)}^{2} := \sum_{α \in N^{d}, | α | = m} \int_{Ω} | D^{α} u |^{2} d x .

rewriting error bounds using Sobolev norms, we get:

Corollary 3.3.3.4. Error estimate for piecewise linear interpolation in 2D

Under the assumptions/with notations of Thm. 3.3.2.21 ( $\forall u \in H^{2} (Ω)$ )

\begin{aligned} ∥ u - I_{1} u ∥_{L^{2} (Ω)} & \leq \sqrt{\frac{3}{8}} h_{M}^{2} | u |_{H^{2} (Ω)} \\ | u - I_{1} u |_{H^{1} (Ω)} & \leq \sqrt{\frac{3}{24}} ρ_{M} h_{M} | u |_{H^{2} (Ω)} \end{aligned}

Theorem 3.3.3.8. Sobolev embedding theorem

$m > \frac{d}{2} \Rightarrow H^{m} (Ω) \subset C^{0} (\bar{Ω}) \land \exists C = C (Ω) > 0 : ∥ u ∥_{\infty} \leq C ∥ u ∥_{H^{m} (Ω)} \forall u \in H^{m} (Ω)$

Theorem answers:

"How many derivatives must a function have (in a weak, square-integrable sense) before we can guarantee it is actually a continuous, point-wise defined function?"

e.g. $d = 2$ : m>1 => $H^{2} (Ω)$ functions are guaranteed to be continuous

Video 3.3 A Priori (Asymptotic) Finite Element Error Estimates III (16 min)

Theorem 3.3.5.6. Best approximation error estimates for Lagrangian finite elements

Let $Ω \subset R^{d}, d = 1, 2, 3$ , be a bounded polygonal/polyhedral domain equipped with a mesh $M$ consisting of simplices or parallelepipeds. Then, for each $k \in N$ , there is a constant $C > 0$ depending only on $k$ and the shape regularity measure $ρ_{M}$ such that

inf_{v_{h} \in S_{p}^{0} (M)} ∥ u - v_{h} ∥_{H^{1} (Ω)} \leq C \underset{discretization params.}{\underset{⏟}{{(\frac{h_{M}}{p})}^{min {p + 1, k} - 1}}} \underset{smoothness reqirement}{\underset{⏟}{∥ u ∥_{H^{k} (Ω)}}} \forall u \in H^{k} (Ω) (3.3 .5 .7)

best approximation error $\to$ bounds for FE discretization error for 2nd-order elliptic BVPs (there: $u \hat{=}$ exact solution of BVP) in energy norm
"generic constant":
- actual value unknown
- dependence known

If we assume uniform shape-regularity and local quasi-uniformity of the mesh, then, from theorem 3.3.5.6 we get:

\begin{matrix} (3.3.5.19) & inf_{v_{h} \in S_{p}^{0} (M)} | u - v_{h} |_{H^{1} (Ω)} \leq C N^{- \frac{min p, k - 1}{d}} | u |_{H^{k} (Ω)} \end{matrix}

=> Algebraic convergence of discretization error with respect to N (DOFs), rate

\frac{1}{d} min (p, k - 1)

rate limited by
- polynomial degree $p$
- smoothness
rate deteriorates for larger $d$

By what extra cost can we buy a certain improvement of the FEM solution? (3.3.5.23)

Assuming the a priori error estimate is sharp (the actual error follows the theoretical trend), the relationship between error reduction and computational cost is:

The Goal (Gain): Reduce the discretization error (energy norm) by a factor of $ρ > 1$ .
The Price (Cost): Increase the number of unknowns ( $N$ ) by a factor of:

Cost Factor \approx ρ^{\frac{d}{min p, k - 1}}

takeaways

The Curse of Dimension ( $d$ ): As the spatial dimension $d$ increases, the cost of reducing error grows exponentially. Accuracy is much "cheaper" in 1D than in 3D.
Optimal Efficiency ( $p = k - 1$ ): For best gain/cost, polynomial degree $p$ should match solution's smoothness $k$ ( $H^{k}$ )
Smoothness Matters: For very smooth solutions ( $k$ is large), high-order elements ( $p > 1$ ) reduce the required increase in $N$ significantly compared to linear elements ( $p = 1$ ).
The "Slow" Limit: If the solution is not smooth (small $k$ ), increasing the polynomial degree $p$ does not help efficiency, because the cost factor will be stuck at $ρ^{d / (k - 1)}$ regardless of how high you choose $p$ .

type of refinement	alg. cvg. rate
h-refinement ( $p$ fixed & small: $p \leq k - 1$ )	$- p / d$
h-refinement (low regularity of u: $k \leq p + 1$ )	$- (k - 1) / d$
p-refinement ( $h_{m}$ fixed, limited regularity $u$ )	$- (k - 1) / d$
p-refinement (perfectly smooth (analytic) solution) (because $C$ depends on $k$ in unspecified way)	theorem doesn't work

Here is a compressed summary of Unit 3.3 for your markdown callout, based on the course videos and documents:

Tldr 3.3: A Priori Finite Element Error Estimates

Cea's Lemma: energy norm of discretization $\hat{=}$ best approximation error of FE space
accuracy is governed by the mesh width ( $h_{M}$ ), the polynomial degree ( $p$ ), and the Sobolev smoothness ( $u \in H^{k}$ ) of the exact solution.

General $H^{1}$ -Error Estimate:

inf_{v_{h} \in S_{p}^{0} (M)} | u - v_{h} |_{H^{1} (Ω)} \leq C {(\frac{h_{M}}{p})}^{min {p, k - 1}} | u |_{H^{k} (Ω)}

where the generic constant $C$ depends on the shape regularity ( $ρ_{M}$ ) of the mesh but is independent of $h_{M}$ .

Algebraic Convergence Trends ( $h$ -refinement): On quasi-uniform meshes, the error decays relative to the number of unknowns ( $N$ ) as:

Error Norm (N) \approx C N^{- α}, with rate α = \frac{min {p, k - 1}}{d}

Curse of Dimension: Convergence speed deteriorates as the spatial dimension $d$ increases.
Rate Limits: The rate is limited by either the polynomial degree $p$ (for smooth solutions) or the solution's lack of regularity $k$ (for rough solutions).

Asymptotic Efficiency (Cost-Benefit): Reducing the error by a factor $ρ > 1$ requires increasing the problem size $N$ by a factor of $ρ^{d / min p, k - 1}$ . Efficiency is theoretically optimal when $p = k - 1$ .

#todo Learn more about Sobolev regularity (solution $u \subset H^{k} (Ω)$ , which $k$ ?)

Video 3.4 Elliptic Regularity Theory (11 min)

WHat is the Sobolev regularity of solutions of elliptic BVP?

Theorem 3.4.0.2. Smooth elliptic lifting theorem

If $\partial Ω$ is $C^{\infty}$ -smooth, ie. possesses a local parameterization by $C^{\infty}$ -functions, and $σ \in C^{\infty} (\overset{―}{Ω})$ , then, for any $k \in N$ ,

\begin{aligned} Dirichlet problem & u \in H_{0}^{1} (Ω) and - div (σ grad u) \in H^{k} (Ω) \\ Neumann prob. & u \in H^{1} (Ω), - div (σ grad u) \in H^{k} (Ω) and grad u \cdot n = 0 on $ \partial Ω \end{aligned}

$\Rightarrow u \in H^{k + 2} (Ω)$ .

In addition, for such $u$ there is $C = C (k, Ω, σ)$ such that

∥ u ∥_{H^{k + 2} (Ω)} \leq C ∥ div (σ grad u) ∥_{H^{k} (Ω)}

probably irrelevant:

Theorem 3.4.0.6. Corner singular function decomposition

Let $Ω \subset R^{2}$ be a polygon with $J$ corners $c^{i}$ . Denote the polar coordinates in the corner $c^{i}$ by $(r_{i}, φ_{i})$ and the inner angle at the corner $c^{i}$ by $ω_{i}$ . Additionally, let $f \in H^{l} (Ω)$ with $l \in N_{0}$ and $l \neq λ_{i k} - 1$ , where the $λ_{i k}$ are given by the singular exponents

\begin{matrix} (3.4.0.7) & λ_{i k} = \frac{k π}{ω_{i}} for k \in N . \end{matrix}

Then $u \in H_{0}^{1} (Ω)$ with = $- Δ u = f$ = in $Ω$ (Dirichlet prob.) can be decomposed

\begin{matrix} (3.4.0.8) & u = u^{0} + \sum_{i = 1}^{J} ψ (r_{i}) \sum_{λ_{i k} < l + 1} κ_{i k} s_{i k} (r_{i}, φ_{i}), κ_{i k} \in R, \end{matrix}

with regular part $u^{0} \in H^{l + 2} (Ω)$ , cut-off functions $ψ \in C^{\infty} (R^{+})$ ( $ψ \equiv 1$ in a neighborhood of $0$ ), and corner singular functions

\begin{matrix} (3.4.0.9) & \begin{aligned} λ_{i k} \notin N : & s_{i k} (r, φ) = r^{λ_{i k}} \sin (λ_{i k} φ), \\ λ_{i k} \in N : & s_{i k} (r, φ) = r^{λ_{i k}} (\ln r) \sin (λ_{i k} φ) . \end{aligned} \end{matrix}

=> reentrant corners lead to very poor regularity

Intermediate

$Ω \subset R^{2}$ has re-entrant corners $⟹$ if $u$ solves $Δ u = f$ in $Ω$ , $u = 0$ on $\partial Ω$ , then $u \notin H^{2} (Ω)$ in general.

however, for convex problems, the dirichlet solution is in $H^{2}$ (BVP is 2-regular)

Theorem 3.4.0.10. Elliptic lifting theorem on convex domains

If $Ω \subset R^{d}$ convex, $u \in H_{0}^{1} (Ω), Δ u \in L^{2} (Ω) \Rightarrow u \in H^{2} (Ω)$ .

Video 3.5 Variational Crimes (15 min)

"theory is incredibly difficult"

Variational crime: instead of solving exact discrete variational problem, we solve a perturbed one:

a_{h} ({\tilde{u}}_{h}, v_{h}) = l_{h} (v_{h})

-> inevitable, since some functions are only given in procedural form

Guideline for acceptable variational crimes

Variational crimes must not affect (type and rate) of asymptotic convergence

Impact of Numerical Quadrature

Note: To show that QR is of order 2 -> show that it integrates polybomials of order 1 exactly (n-1) -> integrates $λ_{1}, λ_{2}, λ_{3}$ exactly

Intermediate

Finite element theory tells us that the above guideline can be met, if the local numerical quadrature rule has sufficiently high order. The quantitative results can be condensed into the following rules of thumb:

Target Accuracy	Required Quadrature Order
$\| u - u_{h} \|_{1} = O (h_{M}^{p})$	Order $2 p - 1$ sufficient for $f_{h}$
$\| u - u_{h} \|_{1} = O (h_{M}^{p})$	Order $2 p - 1$ sufficient for $a_{h}$

Example of crime: Using midpoint rule (order 2) for quadratic elements ( $p = 2$ ) will cause convergence rate in terms of $N$ to drop.

Impact of boundary approximation

Rule of thumb deduced from sophisticated finite element theory:

If $V_{0, h} = S_{p}^{0} (M)$ , use boundary fitting with polynomials of degree $p$ .

for quadratic FEM, use quadratic boundary approximation (see chapter 2)

Example of crime: When $Ω$ has curved boundary, approximation with polygon (linear) is a crime (perturbed domain $Ω_{h}$ )

When using LehrFEM to refine meshes, the boundary stays the same -> round boundary does not get "rounded" when refining

Video 3.8 Validation and Debugging of Finite Element Codes (13 min)

matching theoretical rates are a strong circumstatial evidence of correctness

-> if code converges mroe slowly than theory predicts, there is almost certainly a bug

3 main methods:

observing convergence between mesh levels

test code by comparing solutions on a sequence of nested meshes created by regular refinement

refine mesh -> $| u_{k} - u_{k - 1} |$ should shrink
norm of the difference, $| | u_{k} - u_{k - 1} | |$ should display the same algebraic convergence rate ( $α$ ) as the actual discretization error
shortcut: compare the difference of their norms ( $| | | u_{k} | | - | | u_{k - 1} | | |$ ), should also converge at rate $α$

method of manufactured solutions

9. Second-Order Linear Evolution Problems

Video 9.2.1 Heat Equation (10 min)

linear evolution: time-dependent / transient
posed on space-time cylinder

9.2 Parabolic Initial-Boundary Value Problems

balance law:
Pasted image 20260328150854.png
using localization and fourier law, we get:

\frac{\partial}{\partial t} (ρ u) - div (κ (x) grad u) = f in \tilde{Ω} := Ω \times] 0, T [

spatial b.c.:

For second order parabolic evolutions we can/must use the same spatial boundary conditions as for stationary second order elliptic boundary value problems

(although the boundary conditions are now dependent from x and t)

initial conditions (temporal b.c.):

u (x, 0) = u_{0} (x) BVP

combined:

\begin{aligned} (7.1.1.5) & \frac{\partial}{\partial t} (ρ (x) u) - div (κ (x) grad u) & = f in \tilde{Ω} := Ω \times] 0, T [, \\ (7.1.1.7) & boundary condition: u (x, t) & = g (x, t) for (x, t) \in \partial Ω \times] 0, T [, \\ (7.1.1.8) & initial condition: u (x, 0) & = u_{0} (x) for all x \in Ω . \end{aligned}

Video 9.2.2 Heat Equation: Spatial Variational Formulation (8 min)

Spatial variational form of (7.1.1.5)–(7.1.1.7): seek $t \in] 0, T [\mapsto u (t) \in H_{0}^{1} (Ω)$

\begin{aligned} (6.2.2.4) & \underset{=: m (\dot{u}, v)}{\underset{⏟}{\int_{Ω} ρ (x) \dot{u} (t) v d x}} + \underset{=: a (u (t), v)}{\underset{⏟}{\int_{Ω} κ (x) \nabla u (t) \cdot \nabla v d x}} & = \underset{=: l (t) (v)}{\underset{⏟}{\int_{Ω} f (x, t) v (x) d x}} \forall v \in H_{0}^{1} (Ω), \\ (6.2.2.5) & initial condition ⟶ u (0) & = u_{0} \in H_{0}^{1} (Ω) \end{aligned}

abstract linear evolution problem:

\begin{array}{r} (6.2.2.8) & t \in] 0, T [\mapsto u (t) \in V_{0} : {\begin{cases} m (\dot{u} (t), v) + a (u (t), v) = ℓ (t) (v) \forall v \in V_{0}, \\ u (0) = u_{0} \in V_{0} . \end{cases} \end{array}

a, m are independent of time and s.p.d
solution will depend linearly on f, $u_{0}$

Video 9.2.3 Stability of Parabolic Evaolution Problems (10 min)

#todo sould I add solution of 1D heat equation?
#todo TA prioritized a different part of the viedo/chapter

It is always simpler to look at the square of the norms

reminder: Poincare-Friedrichs inequality:

γ \underset{= L^{2} - norm}{\underset{⏟}{∥ v ∥_{m}^{2}}} \leq \underset{\sim H^{1} - seminorm}{\underset{⏟}{∥ v ∥_{a}^{2}}} \forall v \in V_{0}

Lemma 6.2.3.3. Decay of solutions of parabolic evolutions

For $f \equiv 0$ the solution $u (t)$ of (7.1.2.4) satisfies

∥ u (t) ∥_{m} \leq e^{- γ t} ∥ u_{0} ∥_{m}, ∥ u (t) ∥_{a} \leq e^{- γ t} ∥ u_{0} ∥_{a} \forall t \in] 0, T [

where $γ > 0$ is the constant from (7.1.3.1), and $∥ \cdot ∥_{a}$ , $∥ \cdot ∥_{m}$ stand for the energy norms induced by $a (\cdot, \cdot)$ and $m (\cdot, \cdot)$ , respectively.

Dissipation of energy in parabolic evolutions

Exponential decay of energy during parabolic evolution without excitation ("Parabolic evolutions dissipate energy")

Video 9.2.4 Spatial Semi-Discretization: Method of Lines (6 min)

-> spatial Galerkin discretization of 6.2.2.8 yields:

\sum_{i = 1}^{N} {\dot{μ}}_{i} (t) m (b_{h}^{i}, b_{h}^{j}) + \sum μ_{i} (t) a (b_{h}^{i}, b_{h}^{j}) = ℓ (t) (b_{h}^{j}) \forall v_{h} \in V_{0, h}, j = 1, \dots, N

Method-of-lines ordinary differential equation

Combining (7.1.4.1) and (7.1.4.2) we obtain

\begin{matrix} (6.2.4.4) & (7.1 .4 .1) \Rightarrow {\begin{cases} M {\frac{d}{d t} \vec{μ} (t)} + A \vec{μ} (t) = \vec{φ} (t) \\ \vec{μ} (0) = {\vec{μ}}_{0} \end{cases} for 0 < t < T, \end{matrix}

with

s.p.d. stiffness matrix $A \in R^{N, N}$ , $(A)_{i j} := a (b_{h}^{j}, b_{h}^{i})$ (independent of time),
s.p.d. mass matrix $M \in R^{N, N}$ , $(M)_{i j} := m (b_{h}^{j}, b_{h}^{i})$ (independent of time),
source (load) vector $\vec{φ} (t) \in R^{N}$ , $(\vec{φ} (t))_{i} := ℓ (t) (b_{h}^{i})$ (time-dependent),
${\vec{μ}}_{0} ≜$ coefficient vector of a projection of $u_{0}$ onto $V_{0, h}$ .

Note

(7.1.4.4) is an ordinary differential equation (ODE) for $t \mapsto \vec{μ} (t) \in R^{N}$

but only semi-discretization, since $\vec{μ} (t)$ is still function of continuous time

=> MOL-ODE is semi-discrete

Video 9.2.7 Timestepping for Method-of-Lines ODE (30 min)

\begin{aligned} (7.1.7.2) & Expl. Eul. (7.1.7.2): & {\vec{μ}}^{(j)} = {\vec{μ}}^{(j - 1)} + τ_{j} M^{- 1} (\vec{φ} (t_{j - 1}) - A {\vec{μ}}^{(j - 1)}) \\ (7.1.7.3) & Impl. Eul. (7.1.7.3): & {\vec{μ}}^{(j)} = (τ_{j} A + M)^{- 1} (M {\vec{μ}}^{(j - 1)} + τ_{j} \vec{φ} (t_{j})) \end{aligned}

since s.p.d., can be invertedy -> use Cholesky-decomposition ( #todo what is that?)

MOL-ODE is stiff!

Explicit Euler timestepping is not stable (exponential blow-up) in case of large timestep combined with fine mesh.

Lemma 6.2.7.32. Behavior of generalized eigenvalues

Let $M$ be a simplicial mesh and $A$ , $M$ denote the Galerkin matrices for the bilinear forms $a (u, v) = \int_{Ω} grad u \cdot grad v d x$ and $m (u, v) = \int_{Ω} u (x) v (x) d x$ , respectively, and $V_{0, h} := S_{p, 0}^{0} (M)$ . Then the smallest and largest generalized eigenvalues of $A \vec{μ} = λ M \vec{μ}$ , denoted by $λ_{min}$ and $λ_{max}$ , satisfy

\frac{1}{diam (Ω)^{2}} \leq λ_{min} \leq C, λ_{max} \geq C h_{M}^{- 2},

where the "generic constants" ( $\to$ Rem. 3.3.5.10) depend only on the polynomial degree $p$ , the domain $Ω$ , and the shape regularity measure $ρ_{M}$ .

-> timestep constraint: $τ \leq C h_{m}^{2}$ (similar to ode45)

diagonalization of RK-SSM ( $\vec{φ} \equiv 0$ )
Pasted image 20260328180532.png

Unconditional stability of single step methods

Necessary condition for $universal unconditional stability$ of a single step method for semi-discrete parabolic evolution problem (7.1.4.4) ("method of lines"):

The discrete evolution $Ψ_{λ}^{τ} : R \mapsto R$ of the single step method applied to the scalar ODE $\dot{u} = - λ u$ satisfies

\begin{matrix} (6.2.7.45) & λ > 0 ⟹ lim_{j \to \infty} (Ψ_{λ}^{τ})^{j} u_{0} = 0 \forall u_{0} \in R, \forall τ > 0. \end{matrix}

RK-SSM for $\dot{u} = - λ u \Rightarrow u^{(j)} = \underset{stability functions}{\underset{⏟}{S (- λ τ)}} u^{(j - 1)}$

\begin{array}{r} \begin{array}{cc} c & A \\ b^{⊤} \end{array} \\ (6.2.7.47a) & Ψ_{λ}^{t, t + τ} u = S (- λ τ) u, \end{array}

with rational stability function

\begin{matrix} (6.2.7.47b) & S (z) = 1 + z b^{⊤} (I - z A)^{- 1} 1 = \frac{det (I - z A + z 1 b^{⊤})}{det (I - z A)}, 1 = [1, \dots, 1]^{⊤} \in R^{s} . \end{matrix}

$\to$ a rational function in $z$

RK-SSM is unconditionally stable for MOL-ODE if $| S (z) | \leq 1 \forall z < 0$

evne desirable $| S (z) | ≪ 1$ for $z < 0$

Definition 6.2.7.48. L(

π

)-stability

A Runge-Kutta single-step method satisfying (7.1.7.45) is called $L (π) -stable$ , if its stability function $S (z)$ according to (7.1.7.47b) satisfies

$| S (z) | < 1$ for all $z < 0$ , and
" $S (- \infty)$ " $:= lim_{z \in R \to - \infty} S (z) = 0$ .

Pasted image 20260328181937.png
necessarily implicit

Video 9.2.8 Fully Discrete Method of Lines: Convergence (18 min)

The Parabolic Discretization Path

Heat Equation (Strong Form): transient heat equation on a space-time cylinder $Ω \times] 0, T [$ , subject to initial and boundary conditions: $\frac{\partial u}{\partial t} - Δ u = f (x, t) in Ω, u (x, t) = 0 on \partial Ω$
Spatial Variational Formulation: By treating time $t$ as a passive parameter and applying integration by parts in the spatial domain, the PDE is converted to a weak form:

\int_{Ω} \frac{\partial u}{\partial t} v d x + \int_{Ω} grad u \cdot grad v d x = \int_{Ω} f (t) v d x \forall v \in H_{0}^{1} (Ω)

Abstract Formulation: The weak form is generalized into an abstract parabolic evolution problem using the mass bilinear form $m (\cdot, \cdot)$ and the stiffness bilinear form $a (\cdot, \cdot)$ :

m (\dot{u}, v) + a (u, v) = ℓ (t) (v) \forall v \in V_{0}

Method of Lines ODE: Applying Galerkin discretization in space (e.g., using finite elements) transforms the abstract evolution into a system of ordinary differential equations for the coefficient vector $\vec{μ} (t)$ :

M \dot{\vec{μ}} (t) + A \vec{μ} (t) = \vec{φ} (t)

where $M$ is the mass matrix and $A$ is the stiffness matrix.
5. Runge-Kutta Method (Full Discretization): Finally, the semi-discrete ODE is integrated in time using a single-step method (SSM). For parabolic problems, L-stable Runge-Kutta methods (such as Radau methods) are preferred to ensure stability:

{\vec{μ}}^{(j)} = Ψ_{τ, t_{j - 1}} {\vec{μ}}^{(j - 1)}

resulting in a fully discrete scheme where the total error is the combined effect of spatial and temporal discretization.

Pasted image 20260328195207.png

Refinement for fully discrete parabolic evolution problems

Guideline: spatial and temporal resolution have to be adjusted in tandem

for a conditionally stable explicit method, we often have timestep constraint $τ \leq C h^{2}$

-> if we reduce mesh for accuracy, we have to reduce the timestep more severly than accuarcy requires -> inefficiency "overkill constraint"

-> use unconditionally stable implicit methods instead

10. Modeling Fluid Flow

Video 10.1.1 Modeling Fluid Flow (7 min)

flow field $v : Ω \to R^{d} \hat{=}$ velocity field

needs to be at least Lipschitz-continuous
if $v \cdot \vec{n} = 0$ -> no in/outflow -> streamlines exists for all times

streamline ODE

$\dot{y} = v (y)$

flow map

Φ : {\begin{cases} R \times Ω ⟶ Ω \\ (t, y_{0}) ⟶ Φ^{t} y_{0} := y (t) \end{cases}

where $t \mapsto y (t)$ solves IVP for streamline with $y (0) = y_{0}$
properties:

\begin{aligned} Φ^{0} y_{0} & = y_{0} \forall y_{0} \in Ω \\ Φ^{s + t} & = Φ^{s} \circ Φ^{t} \end{aligned}

Video 10.1.2 Heat Convection and Diffusion (5 min)

\underset{balance law}{\underset{⏟}{div j = f}} + \underset{generalized fourier’s law}{\underset{⏟}{j (x) = - κ grad u (x) + v (x) ρ u (x)}}

=> linear scalar convection-diffusion equation for energy (CDE) (1.6.0.3) + (10.1.2.4):

\begin{matrix} (10.1.2.8) & - div (κ grad u) + div (ρ v (x) u) = f in Ω \end{matrix}

\begin{aligned} - div (κ grad u) & + div (ρ v (x) u) = f \\ ↓ & ↓ \\ diffusive term & convective term \\ (2nd-order) & (1st-order) \end{aligned}

Intermediate

Suitable boundary conditions fitting a PDE are determined by the highest-order term.

Here: diffusion term

Video 10.1.3 Incompressible Fluids (10 min)

Definition 10.1.3.1. Incompressible flow field

A fluid flow is called incompressible, if the associated flow map $Φ^{t}$ is $v o l u m e p r e s e r v i n g$ ,

| Φ^{t} (V) | = | Φ^{0} (V) | = | V |

for all sufficiently small $t > 0$ , and for all control volumes $V$ .

mathematically:

| Φ^{t} (V) | = \int_{Φ^{t} (V)} 1 d x = \int_{V} | det D_{x} Φ^{t} (\hat{x}) | d \hat{x}

\frac{d}{d t} | Φ^{t} (V) | = 0 ⟺ \int_{V} \frac{\partial}{\partial t} | det D_{x} Φ^{t} (x) | d x = \frac{\partial}{\partial t} det D_{x} Φ^{t} (x) = 0 \forall V

Theorem 10.1.3.7. Differentiation formula for determinants

Let $S : I \subset R \mapsto R^{n, n}$ be a smooth matrix-valued function. If $S (t_{0})$ is regular for some $t_{0} \in I$ , then

\frac{d}{d t} (det \circ S) (t_{0}) = det (S (t_{0})) tr (\frac{d S}{d t} (t_{0}) S^{- 1} (t_{0})),

where $det : R^{n, n} \to R$ is the matrix determinant and $tr$ stands for the trace of a matrix.

Theorem 10.1.3.8. Divergence-free velocity fields for incompressible flows

A stationary fluid flow in $Ω \subset R^{d}$ is incompressible ( $\to$ Def. 10.1.3.1), if and only if its associated velocity field $v = [v_{1}, \dots, v_{d}]^{⊤} : Ω \to R^{d}$ satisfies $div v = \sum_{j = 1}^{d} \frac{\partial v_{j}}{\partial x_{j}} = 0$ everywhere in $Ω$ .

scalar convection-diffusion equation for incompressible flow:

\begin{matrix} (10.1.3.11) & - κ Δ u + ρ v \cdot grad u = f in Ω \end{matrix}

where $κ = c o n s t ., ρ = c o n s t .$

Theorem 10.1.3.13. Maximum principle for scalar 2nd-order convection diffusion equations

Let $v : Ω \mapsto R^{d}$ be a continuously differentiable vector field and $u \in C^{0} (\overset{―}{Ω}) \cap C^{2} (Ω)$ . Then there holds the maximum principle

\begin{array}{r} - Δ u + v \cdot grad u \underset{heat sources only}{\underset{⏟}{\geq}} 0 ⟹ min_{x \in \partial Ω} u (x) = min_{x \in Ω} u (x), \\ - Δ u + v \cdot grad u \leq 0 ⟹ max_{x \in \partial Ω} u (x) = max_{x \in Ω} u (x) . \end{array}

Video 10.1.4 Time-Dependent (Transient) Heat Flow in a Fluid (4 min)

new: time-dependent velocity field $v = v (x, t)$

now, conservation of energy also contains change in time:

\underset{energy stored in V}{\underset{⏟}{\frac{d}{d t} \int_{V} ρ u d x}} + \underset{power flux through \partial V}{\underset{⏟}{\int_{\partial V} j \cdot n d S}} = \underset{heat generation in V}{\underset{⏟}{\int_{V} f d x}} for all "control volumes" V . (9.2 .1 .3)

yields transient CDE:

\frac{\partial}{\partial t} (ρ u) - div (κ grad u) + div (ρ v (x, t) u) = f (x, t) in \tilde{Ω} := Ω \times] 0, T [. (10.1 .4 .1)

$+$ "elliptic" boundary conditions
$+$ initial condition : $u (x, 0) = u_{0} (x)$

for incompressible fluid ${div}_{x} v (x, t) = 0$ :

\begin{matrix} (10.1.4.2) & \frac{\partial}{\partial t} (ρ u) - κ Δ u + ρ v (x, t) \cdot grad u = f (x, t) in \tilde{Ω} := Ω \times] 0, T [ \end{matrix}

Video 10.2.1 Singular Perturbation (16 min)

linear variational problem ( $u, w \in H_{0}^{1} (Ω)$ ):

\begin{matrix} (10.2.0 .3) & \underset{bilinear form a (u, w)}{\underset{⏟}{ϵ \int_{Ω} grad u \cdot grad w d x + \int_{Ω} (v \cdot grad u) w d x}} = \underset{linear form ℓ (w)}{\underset{⏟}{\int_{Ω} f (x) w (x) d x}} \end{matrix}

$a (\cdot, \cdot)$ not symmetric $⟹$ no energy norm
$a (\cdot, \cdot)$ continuous on $H^{1} (Ω)$ , since $v$ bounded
$w \in H_{0}^{1} (Ω) : a (w, w) = ϵ \int_{Ω} ∥ grad w ∥^{2} d x > 0 w \neq 0$ (2nd part of bilinear form vanishes)
- $⟹ a (\cdot, \cdot)$ is positive definite

#todo maybe some notes on the derivation of notion?

use method of characteristics (MoC) to compute $u (y (t))$ :

\begin{matrix} (10.2.1.8) & u (y (t)) = u (x_{0}) + \int_{0}^{t} f (y (s)) d s \end{matrix}

Notion 10.2.1.9. Singularly perturbed problem

A boundary value problem depending on parameter $ϵ \approx ϵ_{0}$ is called singularly perturbed, if the limit problem for $ϵ \to ϵ_{0}$ is not compatible with the boundary conditions.

For the pure transport problem $v \cdot grad u = f$ Dirichlet boundary conditions can be imposed only on the inflow boundary part

\begin{matrix} (10.2.1.11) & Γ_{in} := {x \in \partial Ω : v (x) \cdot n (x) < 0} . \end{matrix}

but not on its complement in $\partial Ω$ , the outflow boundary part

\begin{matrix} (10.2.1.12) & Γ_{out} := {x \in \partial Ω : v (x) \cdot n (x) > 0} . \end{matrix}

Pasted image 20260421090327.png|400

$ϵ > 0$ : Dirichlet b.c. can be imposed everywhere on $\partial Ω$
$ϵ = 0$ : Dirichlet b.c. can be imposed only on $Γ_{in}$

=> singular perturbation

Pasted image 20260421090250.png|400

On a closed streamline we cannot find a point $x_{0} \in Γ_{in}$ , which renders the solution formula (10.2.1.8) meaningless.

Video 10.2.2 Upwinding (13 min)

$Ω \subset R^{d} \hat{=}$ bounded computational domain

\begin{matrix} (10.2.0.1) & \underset{\begin{matrix} diffusive term \\ (2nd-order term) \end{matrix}}{\underset{⏟}{- ϵ Δ u}} + \underset{\begin{matrix} convective term \\ (1st-order term) \end{matrix}}{\underset{⏟}{v (x) \cdot grad u}} = f in Ω, u = 0 on \partial Ω . \end{matrix}

singularly perturbed for $ϵ \to 0$
$ϵ = 0$ : Dirichlet b.c. can only be imposed on inflow boundary

#todo look at derivation of 10.2.2.2

for $ϵ ≪ 1$ spurious oscillation due to even-odd decoupling

Alrebraic sign conditions

(Linearly interpolated) discrete solution satisfies maximum principle (3.7.2.1). $⟺$ system matrix complies with sign-conditions

(3.7.2.8): positive diagonal entries, $(A)_{i i} > 0$ ,
(3.7.2.9): non-positive off-diagonal entries, $(A)_{i j} \leq 0, if i \neq j$ ,
“(3.7.2.10)”: diagonal dominance, $\sum_{j} (A)_{i j} \geq 0$ .

Tldr: Upwinding: idea

Goal: $ϵ$ -robust method (satisfaction of maximum principle)
Mechanism of Upwinding: use one-sided difference quotients in direction of velocity ("upwind" direction)
- In 1D: use backwards diff. quotient to satisfy algebraic sign conditions regardless of $h, ϵ$ (forwards diff. quotient only works for $h < 2 ϵ$ )
- 2D/3D: standard Galerking methods produces "terrible oscillation" -> upwinding needed

"tradeoff between cvg. speed and accuracy ()"

Video 10.2.2.1 Upwind Quadrature (17 min)

different derivation of upwinding for 1D: apply global composite trapezoidal rule to 2nd term of convertcion-diffusion 2-pt BVP:

\begin{matrix} (10.2.2.1) & - ϵ \frac{d^{2} u}{d x^{2}} + \frac{d u}{d x} = f (x) in Ω, u (0) = 0, u (1) = 0 \end{matrix}

$u \in H_{0}^{1} (] 0, 1 [) :$

\underset{=: a (u, v)}{\underset{⏟}{ϵ \int_{0}^{1} \frac{d u}{d x} (x) \frac{d v}{d x} (x) d x + \int_{0}^{1} \frac{d u}{d x} (x) v (x) d x}} = \underset{=: ℓ (v)}{\underset{⏟}{\int_{0}^{1} f (x) v (x) d x}} \forall v \in H_{0}^{1} (] 0, 1 [)

Pasted image 20260421110419.png
upwind quadrature rule approximates the convective bilinear form as a weighted sum over all mesh nodes $p \in V (M)$ :

\int_{Ω} (v \cdot grad u_{h}) w_{h} d x \approx \sum_{p \in V (M)} m (p) \cdot (v (p) \cdot (grad u_{h})_{u p} (p)) w_{h} (p)

$m (p)$ (Vertex Mass): This is the weight associated with node $p$ , representing the area of the node patch surrounding it. It is calculated as one-third of the sum of the areas of all triangles $K$ sharing vertex $p$ :

m (p) = \frac{1}{3} \sum_{K \in U_{p}} | K |

$w_{h} (p)$ : The value of the test function at the node (which is $1$ if $w_{h}$ is the basis function $b_{p}$ and $0$ otherwise).

p.w. linear functions: $grad u_{h}$ is constant on each triangle but discontinuous across edges and at nodes

=> idea: use values from upstream triangle $K_{u}$
Pasted image 20260421110516.png|400

K_{u} (p) := {K \in M : p \in K and \exists s > 0 such that p - s v (p) \in K}

The gradient value used in the quadrature is then defined by the limit:

v (p) \cdot grad u_{h} (p) := lim_{δ \to 0} v (p) \cdot grad u_{h} (p - δ v (p))

using "tent function" basis ${b_{i}}_{i = 1}^{N}$ , entries of the convective part of the Galerkin matrix $B$ are given by:

B_{j, i} = m (a_{j}) (v (a_{j}) \cdot (grad b_{i}) |_{K_{u} (a_{j})})

=> fulfills algebraic sign conditions => $ϵ$ -robust

#todo not happy with this summary

Always be sure in which direction the velocity is going!

Video 10.2.2.2 Streamline Diffusion (19 min)

Artificial diffusion/viscosity CFD

upwinding = $h$ -dependent enhancement of diffusive term

problem:
Pasted image 20260421111344.png
enhancing diffusion $⟹$ smearing of internal layers and ridges

we are no longer solving the right problem!

Heuristics of streamline upwinding (SV)

Since the solution is smooth along streamlines, then adding diffusion in the direction of streamlines cannot do much harm.

solution: anisotropic artificial diffusion in streamline direction

on cell $K$ replace:

ϵ \leftarrow \underset{new diffusion tensor}{\underset{⏟}{ϵ I + = δ_{K} v_{K} v_{K}^{T} =}} \in R^{2, 2}

$v_{K} \hat{=}$ local velocity (e.g., obtained by averaging)
$δ_{K} > 0 \hat{=}$ method parameter controlling the strength of anisotropic diffusion ( $ε_{K}, h_{K}$ - dependent

$⟹$ SU-enhanced variational problem

\begin{matrix} (10.2.2.30) & \int_{Ω} (ϵ I + = δ_{K} v_{K} v_{K}^{T} =) grad u \cdot grad w + v (x) \cdot grad u w d x = \int_{Ω} f w d x \forall w \in H_{0}^{1} (Ω) . \end{matrix}

This tampering affects the solution $u$ : (solution of (10.2.2.30) $\neq$ solution of (10.2.0.1))

we want the original solution $u$ still solves the stabilized VP -> consistent modification

Definition 10.2.2.31. Consistent modifications of variational problems

A variational problem is called a consistent modification of another, if both possess the same (unique) solution(s): If $u \in V$ solves

a (u, v) = ℓ (v) \forall v \in V_{0}

then, it should also solve modified LVP:

\tilde{a} (u, v) = \tilde{ℓ} (v) \forall v \in V_{0}

to satisfy this, we add a residual (term that vanishes for exact solution)

\begin{aligned} (10.2.2.35) & \int_{Ω} ϵ grad u \cdot grad w + (v (x) \cdot grad u) w d x \\ + & \underset{stabilization term = 0, if u is exact solution}{\underset{⏟}{\sum_{K \in M} δ_{K} \int_{K} \underset{―}{(- ϵ Δ u + v \cdot grad u - f)} \cdot (v \cdot grad w) d x}} = \int_{Ω} f w d x \forall w \in H_{0}^{1} (Ω) \end{aligned}

Choice of $δ_{K}$ : (for p.w. linear FE)

δ_{K} := {\begin{cases} ϵ^{- 1} h_{K}^{2} & , if \frac{∥ v ∥_{K, \infty} h_{K}}{2 ϵ} \leq 1, \leftarrow diffusion dominated \\ h_{K} & , if \frac{∥ v ∥_{K, \infty} h_{K}}{2 ϵ} > 1. \leftarrow convection dominated \end{cases}

pros: sharp jump -> no smoothing (compared to upwind quadrature)
pros: faster convergence than upwind quadrature (upwind quadrature sacrifices cvg.)
cons: ("small") overshoots
cons: violation of max principle

"Meta theorem" :

No $ϵ$ -robust linear method that respects the maximum principle strictly can be more than 1st-order convergent in $L^{2}$ !

Video 10.3.1 Method of Lines (13 min)

Method of Lines (MOL) for transient convection-diffusion problems: "half-way" discretization strategy that converts a time-dependent partial differential equation (PDE) into a system of ordinary differential equations (ODEs) by first discretizing only the spatial variables:

transient convection-diffusion equation with Dirichlet boundary conditions:

\begin{matrix} (10.3.1.1) & {\begin{cases} \frac{\partial u}{\partial t} - ϵ Δ u + v (x, t) \cdot grad u = f & in \tilde{Ω} := Ω \times] 0, T [, \\ u (x, t) = g (x, t) \forall x \in \partial Ω, 0 < t < T, & u (x, 0) = u_{0} (x) \forall x \in Ω . \end{cases} \end{matrix}

Spatial Discretization: $t$ is treated as parameter, spatial domain is discretized (e.g., using finite elements on a fixed mesh) => system of ODEs for the time-dependent basis expansion coefficients

\begin{matrix} (10.3.1.2) & M \frac{d \vec{μ}}{d t} (t) + ϵ A \vec{μ} (t) + B (t) \vec{μ} (t) = \vec{φ} (t) \end{matrix}

* $\vec{μ} (t)$ : The time-dependent vector of basis expansion coefficients representing $u_{h} (x, t)$ .
* $M$ : mass matrix
* $A$ : FE Galerkin matrix for $- Δ$
* $B (t)$ : transport matrix (discretization of the convective term $v \cdot grad u$ )
* $\vec{φ} (t)$ : load vector (depending on $f$ and $g$ -> time-independet)
2. Temporal Discretization: solve ODE system using numerical integrator (timestepping method) to obtain a fully discrete solution.

robustness:

spatial discretization: use $ϵ$ -robust discretization for convective term (e.g. upwind schema)
timestepping stability: use e.g. implicit euler/SDIRK-2

12. Finite Elements for the Stokes Equations

Video 12.1 Modeling the Flow of a Viscous Fluid (30 min)

(substantial modelling component in this chapter)

thermodynamic perspective:
- $div \vec{v} = 0 in Ω$ (incompressible)
assumptions:
- small Reynolds number -> slow flows / small $Ω$ / sticky fluid
- energy dissipation mainly/only occurs in eddies -> curl
  - dissipation: kinetic energy -> heat
boundary conditions:
- no-slip: $v = 0$ on $\partial Ω$

with these assumptions, we get the space:

\begin{matrix} (12.1.0.6) & V := {v : \overset{―}{Ω} \mapsto R^{d} continuous, div v = 0, v_{| \partial Ω} = 0} . \end{matrix}

balance between 2nd law of T.D. (entropy maximized) and 1st law of T.D. (conservation of energy):

\begin{aligned} (12.1.0.12) & \underset{dissipated energy}{\underset{⏟}{\int_{Ω} μ ∥ curl v (x) ∥^{2} d x}} = \underset{energy injected through forces}{\underset{⏟}{\int_{Ω} f \cdot v d x}} \\ (12.1.0.13) & v^{*} = argmax & {\int_{Ω} μ ∥ curl v (x) ∥^{2} d x : v \in V, v satisfies (12.1.0.12)} \end{aligned}

converted to quadratic minimization problem:

\begin{matrix} (12.1.0.18) & v^{*} = \underset{v \in V}{argmin} \frac{1}{2} \int_{Ω} μ ∥ curl v (x) ∥^{2} d x - \int_{Ω} f \cdot v d x . \end{matrix}

#todo review questions

Video 12.2.1 The Stokes Equations: Constrained variational formulation (20 min)

Lemma 12.2.1.2.

- Δ = curl curl - grad div

For $v \in C^{2} (\overset{―}{Ω}), v_{| \partial Ω} = 0$ , holds ( $v = [v_{1}, \dots, v_{d}]^{⊤}$ )

\int_{Ω} ∥ curl v ∥^{2} d x + \int_{Ω} | div v |^{2} d x = \int_{Ω} ∥ D v ∥_{F}^{2} d x = \int_{Ω} \overset{Jacobian}{\overset{⏞}{D v}} : D v d x = \sum_{i = 1}^{d} \int_{Ω} ∥ grad v_{i} ∥^{2} d x

where $D v : D v$ is the Euclidean inner product of matrices

using lemma 12.2.1.2 + divergence constraint ( $div v = 0$ ), we get:

\begin{aligned} (12.2.1.4) & v^{*} & = \underset{v \in V}{argmin} (\frac{1}{2} a (v, v) - ℓ (v)) \\ (12.2.1.12) & a (v, w) & := \int_{Ω} μ D v : D w d x = μ \sum_{i = 1}^{d} \int_{Ω} grad v_{i} \cdot grad w_{i} d x \\ ℓ (w) & := \int_{Ω} f \cdot w d x (f external stirring force field) \end{aligned}

on the space

\begin{matrix} (12.2.1.13) & V = H_{0}^{1} (div 0, Ω) := {v \in (H_{0}^{1} (Ω))^{d} : div v = 0} \end{matrix}

s.p.d. $a$ (since $a (v, v) = 0 ⟹$ all gradients of components vanish)
not decoupled (since $div \vec{v} = 0$ couples components of velocity field, unfortunately)

No simple FE subspace of

H_{0}^{1} (div 0, Ω)

-> variational form useless

#todo review questions

Video 12.2.2 The Stokes Equations: Saddle Point Problem Formulation (55 min)

[AI] We want to solve a constrained minimization problem. We introduce a state space $U$ and a multiplier space $Q$ . For a convex differentiable functional $J (v)$ subject to a linear constraint $B v = 0$ , we define the Lagrangian functional $L (v, p) := J (v) + (p, B v)_{Q}$ , where $p \in Q$ is the Lagrange multiplier. The solution is found by solving the min-max problem $v^{*} = {argmin}_{v \in U} sup_{p \in Q} L (v, p)$ .

Lemma 12.2.2.4. Necessary conditions for existence of solution of saddle point problems

Any solution $v^{*}$ of (12.2.2.3) will be the first component of a zero $(v^{*}, p^{*})$ of the derivative ("gradient") of the Lagrangian functional $L$ : $= D L (v^{*}, p^{*}) = 0$

$\hat{=}$ solutions of min-max problems correspond to saddle points:
Pasted image 20260422220206.png
[AI] To find the saddle point, we evaluate the partial derivatives of $L$ :

\begin{aligned} \frac{\partial L}{\partial v} (v, p) (w) & = D J (v) (w) + (p, B w)_{Q} = 0 \forall w \in U \\ \frac{\partial L}{\partial p} (v, p) (q) & = (q, B v)_{Q} = 0 \forall q \in Q \end{aligned}

Setting the derivatives to zero yields a system of variational equations. If $J (v)$ is quadratic, we replace $D J (v) (w)$ with a bilinear form $a (v, w)$ and a linear form $ℓ (w)$ .

...we get to the variational saddle point problem (SPP)

\begin{matrix} (12.2.2.9) & \begin{aligned} \underset{a (v^{⋆}, w)}{\underset{⏟}{⟨ D J (v^{*}), w ⟩}} + (p^{*}, B w)_{Q} & = \underset{ℓ (w)}{\underset{⏟}{0}} \forall w \in U, \\ (q, B v^{*})_{Q} & = 0 \forall q \in Q . \end{aligned} \end{matrix}

$p$ is the Lagrange multiplier / "pressusre"
$a$ is also called "elliptic"

Seek the velocity field $v \in (H_{0}^{1} (Ω))^{d}$ and a Lagrange multiplier $p \in L^{2} (Ω)$ such that

\begin{matrix} (12.2.2.16) & \begin{aligned} \underset{a (v^{*}, w)}{\underset{⏟}{\int_{Ω} μ D v : D w d x}} + \underset{b (w, p)}{\underset{⏟}{\int_{Ω} div w p d x}} & = \underset{ℓ (w)}{\underset{⏟}{\int_{Ω} f \cdot w d x}} & \forall w \in (H_{0}^{1} (Ω))^{d} \\ \underset{b (v, q)}{\underset{⏟}{\int_{Ω} div v q d x}} & = 0 & \forall q \in L^{2} (Ω) \end{aligned} \end{matrix}

[AI] Ensuring uniqueness of pressure: By Gauss' theorem, $\int_{Ω} div v d x = \int_{\partial Ω} v \cdot n d S = 0$ since $v = 0$ on $\partial Ω$ . We infer that the pressure solution $p$ is only unique up to a constant (adding any constant $c$ to $p$ still retains a valid solution). To restore uniqueness, we impose a vanishing average condition on the pressure. We replace $L^{2} (Ω)$ with the constrained space $L_{*}^{2} (Ω) := {q \in L^{2} (Ω) : \int_{Ω} q d x = 0}$ .

Theorem 12.2.2.23. Ladyzhenskaya–Babuška–Brezzi conditions (LBB-conditions)

If the following two conditions are satisfied,
(i) the bilinear form $a$ is $N (b)$ -elliptic (ellipticity on the kernel)

\begin{matrix} (LBB1) & \exists α > 0 : | a (v, v) | \geq α ∥ v ∥_{U}^{2} \forall v \in N (b), \end{matrix}

(ii) the bilinear form $b$ satisfies the inf-sup condition

\begin{matrix} (LBB2) & \exists β > 0 : sup_{v \in U} \frac{| b (v, q) |}{∥ v ∥_{U}} \geq β ∥ q ∥_{Q} \forall q \in Q, \end{matrix}

then the variational saddle point problem (12.2.2.50) possesses a unique solution $(v, p) \in U \times Q$ , which satisfies

\begin{matrix} (12.2.2.24) & ∥ v ∥_{U} + ∥ p ∥_{Q} \leq C (sup_{w \in U} \frac{| ℓ (w) |}{∥ w ∥_{U}} + sup_{q \in Q} \frac{| g (q) |}{∥ q ∥_{Q}}), \end{matrix}

with $C > 0$ independent of $ℓ$ and $g$ .

Pasted image 20260422222216.png

LBB2 $⟹ N (B^{⊤}) = 0$ (trivial kernel) $⟹ B$ has full rank $⟹ n \leq m ⟹ R (B) = R^{n} ⟹ A |_{N (B)}$ is regular
LBB1 $⟹ A$ is symmetric positive definite (elliptic) on the null space of $B$ $⟹$ the saddle point system matrix is invertible, giving a unique solution for $\vec{μ}$ and $\vec{π}$ $⟹ \exists$ unique solution

Applied to stokes SPP:

\begin{aligned} \int_{Ω} μ D v : D w d x + \int_{Ω} div w p d x & = \int_{Ω} f \cdot w d x \forall w \in (H_{0}^{1} (Ω))^{d}, \\ (12.2.2.19) & \int_{Ω} div v q d x & = 0 \forall q \in L_{*}^{2} (Ω), \end{aligned}

\begin{array}{r} (12.2.2.21) & \begin{array}{c} v \in U \\ p \in Q \end{array} : {\begin{cases} a (v, w) + b (w, p) = ℓ (w) & \forall w \in U \\ b (v, q) = g (q) & \forall q \in Q \end{cases} \end{array}

where

$U$ and $Q$ are Hilbert spaces, with norms $∥ \cdot ∥_{U}$ and $∥ \cdot ∥_{Q}$ , respectively,
$a : U \times U \to R$ and $b : U \times Q \to R$ are continuous/bounded bilinear forms,
$ℓ : U \to R$ and $g : Q \to R$ are continuous/bounded linear forms.

note: complex derivation using Theorem (12.2.2.38): existence of stable velocity potentials can be looked up in the script; however, we get:

Theorem 12.2.2.40. Existence and uniqueness of weak solutions of Stokes problem ("

stability estimate

The linear variational saddle point problem (12.2.2.19) ("Stokes problem") has a unique solution $(v, p) \in H_{0}^{1} (div 0, Ω) \times L_{*}^{2} (Ω)$ , which satisfies

\begin{matrix} (12.2.2.41) & \exists C = C (Ω) > 0 : ∥ v ∥_{H^{1} (Ω)} + ∥ p ∥_{L^{2} (Ω)} \leq C ∥ f ∥_{L^{2} (Ω)} . \end{matrix}

Video 12.2.3 Stokes System of Partial Differential Equations (15 min)

Pasted image 20260428103723.png

12.2.3.1. Extra somoothness of velocity and pressure

For the solution $(v, p)$ of (12.2.2.19) we assume

v = (H^{2} (Ω))^{d} and u \in H^{1} (Ω)

Stokes equation:

\begin{matrix} (12.2.3.5) & \begin{aligned} - μ \underset{―}{Δ} u + \nabla p & = f & momentum balance \\ div u & = 0 & in Ω & incompressibility \\ \int_{Ω} p d x & = 0 & zero mean for pressure \\ u & = 0 & on \partial Ω & no-slip b.c. \end{aligned} \end{matrix}

pressure poisson: If $v, p$ sufficiently smooth, we have $div (\underset{―}{Δ} u) = \underset{―}{Δ} (div u)$ , and we get:

Δ p = div f in Ω

-> cannot be solved, missing b.c. for $p$

#todo review questions

12.3 Galerking Discretization of the Stokes Saddle Point Problem

Video 12.3.1 Pressure Instability (40 min)

applying galerking discretization to the problem from chapter 12.1/.2:

Step 1 (space discretization)

\begin{array}{r} M ⟶ M_{n} \subset M \\ Q ⟶ Q_{n} \subset Q \end{array} in (2.19)

N := \dim M_{n}, M := \dim Q_{n} < \infty

DVP:

\begin{matrix} (12.3.0.4) & \begin{array}{l} v_{h} \in U_{h} \\ p_{h} \in Q_{h} \end{array} : \begin{aligned} a (v_{h}, w_{h}) + b (w_{h}, p_{h}) & = ℓ (w_{h}) \forall w_{h} \in U_{h}, \\ b (v_{h}, q_{h}) 0 & = 0 \forall q_{h} \in Q_{h} . \end{aligned} \end{matrix}

Step 2 (choosing basis functions) ${b_{h}^{i}}_{i = 1}^{N}$ of $M_{h}$ , ${β_{h}^{j}}_{j = 1}^{M}$ of $Q_{h}$
- + expansion:

{\underset{―}{v}}_{h} = \sum_{i = 1}^{N} v_{i} b_{h}^{i} p_{h} = \sum_{j = 1}^{M} π_{j} β_{h}^{j}

leads to saddle point LSE:

\begin{matrix} (12.3.0.7) & [\begin{matrix} A & B^{T} \\ B & 0 \end{matrix}] [\begin{matrix} \vec{v} \\ \vec{π} \end{matrix}] = [\begin{matrix} \vec{ϕ} \\ 0 \end{matrix}], \end{matrix}

with matrix entries:

\begin{aligned} (12.3.0.8a) & A & := {[a (b_{h}^{j}, b_{h}^{i})]}_{i, j = 1}^{N} = {[\int_{Ω} μ D b_{h}^{j} (x) : D b_{h}^{i} (x) d x]}_{i, j = 1}^{N} \in R^{N, N}, \\ (12.3.0.8b) & B & := {[b (b_{h}^{j}, β_{h}^{i})]}_{\begin{array}{c} 1 \leq i \leq M \\ 1 \leq j \leq N \end{array}} = {[\int_{Ω} div b_{h}^{j} (x) β_{h}^{i} (x) d x]}_{\begin{array}{c} 1 \leq i \leq M \\ 1 \leq j \leq N \end{array}} \in R^{M, N}, \\ (12.3.0.8c) & \vec{ϕ} & := {[ℓ (b_{h}^{j})]}_{j = 1}^{N} = {[\int_{Ω} f (x) \cdot b_{h}^{j} (x) d x]}_{j = 1}^{N} \in R^{N} . \end{aligned}

Note that saddle point matrix is symmetric, but not s.p.d.

pressure instability
$B$ has to be chosen carefully to prevent spurious oscillations

$B^{T} \in R^{N \times M}$ -> only if $M > N ⟹ N (B^{⊤}) = {0}$ -> $\vec{π}$ unique
- only if $B$ has trivial kernel, the solution is unique
however, this is only a necessary (not sufficient!) condition -> we might get oscillations anyway (BAD NEWS)

Intermediate

$\dim U_{h} \geq \dim Q_{h}$ is a necessary condition for the uniqueness of the pressure solution $p_{h}$ of the discrete saddle point variational problem (12.3.0.4).

Video 12.3.2 Stable Galerkin Discretization of Stokes Saddle Point Problem (30 min)

$div v_{h}$ must be "large enough to fix the pressure" $p_{h} \in Q_{h}$

#todo why, exactly?

example: by rising velocity field space $P^{1} \to P^{2}$ , checkerboard modes instability (see previous video) is cured

Definition 12.3.2.7. Stable finite element pair

A pair of finite element spaces $U_{h} \subset H_{0}^{1} (Ω)$ , $Q_{h} \subset L_{*}^{2} (Ω)$ is a stable finite element pair, if the solution $(v_{h}, p_{h})$ of the discrete saddle point problem (12.3.0.4) satisfies

\exists C > 0 : ∥ v_{h} ∥_{a} + ∥ p_{h} ∥_{L^{2} (Ω)} \leq C sup_{w \in (H_{0}^{1} (Ω))^{d}} \frac{| \int_{Ω} f (x) w (x) d x |}{∥ w ∥_{a}} \forall f \in L^{2} (Ω),

where the constant $C > 0$ may depend only on $Ω$ , the coefficient $μ$ , and the shape regularity measure ( $\to$ Def. 3.3.2.20) of $M$ .

$C$ must be independent of $h_{M}$ !

velocity component: unconditional stability
pressure component:

Theorem 12.3.2.15. inf-sup condition

The finite element spaces $U_{h} \subset H_{0}^{1} (Ω)$ , $Q_{h} \subset L_{*}^{2} (Ω)$ provide a stable finite element pair ( $\to$ Def. 12.3.2.7) for the Stokes problem (12.2.2.19)/(12.3.0.2) if there is a constant $β > 0$ depending only on $Ω$ and the shape regularity measure ( $\to$ Def. 3.3.2.20) of $M$ such that

inf-sup cond.: $sup_{w_{h} \in U_{h}} \frac{| b (w_{h}, q_{h}) |}{∥ w_{h} ∥_{a}} \geq β ∥ q_{h} ∥_{L^{2} (Ω)} \forall q_{h} \in Q_{h} .$

LBB1 is directly satisfied by the fact that $a (\cdot, \cdot)$ is p.d.
LBB2 comes from 12.3.2.15 (above)

-> #todo conclusion?

Video 12.3.3 Convergence of Stable FEM for Stokes (30 min)

abstract lax equivalence principle

Assumption 12.3.3.5. Stability of the discrete variational problem

$\begin{matrix} (12.3.3.6) & \exists C_{s} > 0 : ∥ u_{h} ∥ \leq C_{s} sup_{w_{h} \in H_{h}} \frac{| ℓ (w_{h}) |}{∥ w_{h} ∥} \forall ℓ, \end{matrix}$

using residual equation + assumption + $△$ -inequality:

\underset{discretizat. error}{\underset{⏟}{∥ u - u_{h} ∥}} \leq \underset{independent of u}{\underset{⏟}{(1 + C_{s} C_{c})}} \underset{best approx. error}{\underset{⏟}{inf_{v_{h} \in H_{h}} ∥ u - v_{h} ∥}}

-> $u$ is quasi-optimal (optimal (see Cea's lemma) + some constant)

Quasi-optimality: $C > 0$ independent of data, discretization parameters

Discrete stability

⟹

quasi-optimality

If a discrete linear variational problem is stable according to Ass. 12.3.3.5, then Galerkin solutions are quasi-optimal.

which leads to:

Theorem 12.3.3.13. Convergence of stable FE for Stokes problem

If $U_{h}, Q_{h}$ is =a stable finite element pair= ( $\to$ Def. 12.3.2.7) for the Stokes variational saddle point problem (12.2.2.19), then the corresponding finite element Galerkin solution $(v_{h}, p_{h})$ satisfies

∥ v - v_{h} ∥_{H^{1} (Ω)} + ∥ p - p_{h} ∥_{L^{2} (Ω)} \leq C \underset{sum of both best approximation errors}{\underset{⏟}{(inf_{w_{h} \in U_{h}} ∥ v - w_{h} ∥_{H^{1} (Ω)} + inf_{q_{h} \in Q_{h}} ∥ p - q_{h} ∥_{L^{2} (Ω)})}},

with a constant $C > 0$ that depends only on $Ω, μ$ , and the shape regularity of the finite element mesh.

Note: norm of velocity error also depends best-approx. error of pressure space
slower convergence error matters more!

but we don't want a poor pressure space to affect the convergence of velocity solution:

Definition 12.3.3.16. Pressure-robust Galerkin discretization of the Stokes problem

A pair of finite element spaces $U_{h} \subset (H_{0}^{1} (Ω))^{d}$ , $Q_{h} \subset L^{2} (Ω)$ for the Galerkin discretization of the Stokes variational saddle point problem (12.2.2.19) is called pressure-robust, if the velocity component $v_{h}$ of the Galerkin solution satisfies

∥ v - v_{h} ∥_{H^{1} (Ω)} \leq C inf_{w_{h} \in U_{h}} ∥ v - w_{h} ∥_{H^{1} (Ω)}

with $C > 0$ independent of the driving force field $f$ .

Lemma 12.3.3.17. Criterion for pressure robustness

A pair of finite element spaces $U_{h} \subset (H_{0}^{1} (Ω))^{d}, Q_{h} \subset L^{2} (Ω)$ for the Galerkin discretization of the Stokes variational saddle point problem (12.2.2.19) is pressure-robust in the sense of Def. 12.3.3.16, if

div U_{h} \subset Q_{h} .

If we assume this, from

\begin{matrix} (12.3.3.18) & \begin{array}{l} v_{h} \in U_{h}, \\ p_{h} \in Q_{h} \end{array} : \begin{aligned} \int_{Ω} μ D v_{h} : D w_{h} d x + \int_{Ω} div w_{h} p_{h} d x & = \int_{Ω} f \cdot w_{h} d x \forall w_{h} \in U_{h}, \\ \int_{Ω} div v_{h} q_{h} d x & = 0 \forall q_{h} \in Q_{h} . \end{aligned} \end{matrix}

we get:

\int_{Ω} μ D (v - v_{h}) : D w_{h} d x = 0 \forall w_{h} \in U_{h} (div 0)

and since this is a version of Galerkin orthogonality, (similar to proof for Cea's lemma)

∥ \underset{―}{v} - {\underset{―}{v}}_{h} ∥_{a} \leq inf_{{\underset{―}{w}}_{h} \in M_{h} (div 0)} ∥ \underset{―}{v} - {\underset{―}{w}}_{h} ∥_{a}

-> optimality, if approximating function in discrite space of functions with vanishing div
-> however, no practical spaces have been found that satisfy this criterion -> currently useless

#todo review questions

Video 12.3.4 The Taylor-Hood Finite Element Method (15 min)

stable FE are balanced, if convergence (#^12-3-3-13) has same order for smooth $v, p$
requirements for Stokes FEM:

Intermediate

The pair $(U_{h}, Q_{h})$ of finite element spaces must be stable ( $\to$ Def. 12.3.2.7)
The velocity finite element space $U_{h}$ should provide the same rate of algebraic convergence of the $H^{1} (Ω)$ -best approximation error w.r.t. $h_{M} \to 0$ as the pressure space $Q_{h}$ in $L^{2} (Ω)$ .
The velocity finite element space $U_{h}$ should guarantee 1 and 2 with as few degrees of freedom as possible.

-> 2. and 3. concern efficiency. A stable way to satisfy the requirements is the Taylor-Hood (PS-P1) finite element method for Stokes problem:

a triangular/tetrahedral or rectangular/hexahedral mesh $M$ of $Ω$ , which may even be hybrid, see Section 2.5.1,
the velocity space: $U_{h} := (S_{2, 0}^{0} (M))^{2} \subset (H_{0}^{1} (Ω))^{d}$ , and
the pressure space: $Q_{h} := S_{1}^{0} (M)$ , which means a continuous pressure approximation.

$U_{h}$ and $Q_{h}$ are a stable pair!
algebraic convergence, order 2 for
- velocity, $H^{1}$ -norm
- pressure, $L^{2}$ -norm

Video 12.3.5 The Non-Conforming Crouzeix-Raviart FEM (40 min)

non-conforming: $M_{h} ⊄ (H_{0}^{1} (Ω))^{d}$

Definition 12.3.5.2. 2D scalar Crouzeix-Raviart finite element space

The Crouzeix-Raviart finite element space on a triangular mesh $M$ is

CR (M) := {v \in L^{2} (Ω) : v_{| K} \in P_{1} (K) \forall K \in M,

$v continuous at midpoints of interior edges of M} .$

global shape functions:
local shape functions (LSF) $b_{K}^{1}, b_{K}^{2}, b_{K}^{3}$ associated with edges of triangle $K \in M$ :

CR-P0 FEM for Stokes:

\begin{aligned} (12.3.5.8) & {CR}_{0} (M) & := {v_{h} \in CR (M) : v_{h} (m_{e}) = 0 \forall e \subset \partial Ω} \\ U_{h} & := ({CR}_{0} (M))^{2} ⊄ (H_{0}^{1} (Ω))^{2} \\ Q_{h} & := S_{0, *}^{- 1} (M) \end{aligned}

Pasted image 20260428145349.png
DVP: seek $v_{h} \in (CR (M))^{2}, p_{h} \in S_{0, *}^{- 1} (M)$ such that

\begin{matrix} (12.3.5.16) & \begin{aligned} \sum_{K \in M} \int_{K} μ D {v_{h} |}_{K} : D {w_{h} |}_{K} d x + \sum_{K \in M} \int_{K} div {w_{h} |}_{K} p_{h} d x & = \int_{Ω} f \cdot w_{h} d x, \\ \sum_{K \in M} \int_{K} div {v_{h} |}_{K} q_{h} d x & = 0, \end{aligned} \end{matrix}

for all $w_{h} \in ({CR}_{0} (M))^{2}, q_{h} \in S_{0, *}^{- 1} (M)$ .
derivatives restricted to cell $K$ , since not known outside

spaces are balanced (even optimal)
stable, if:

Lemma 12.3.5.19. Inf-sup condition for

({CR}_{0} (M))^{2}

S_{0}^{- 1} (M)

-pair

With a constant $C > 0$ depending only on $Ω$ and the shape-regularity measure of $M$

\begin{matrix} (12.3.5.20) & sup_{w_{h} \in ({CR}_{0} (M))^{2}} \frac{1}{∥ w_{h} ∥_{H^{1} (M)}} | \sum_{K \in M} \int_{K} div {w_{h} |}_{K} q_{h} d x | \geq \frac{1}{C} ∥ q_{h} ∥_{L^{2} (Ω)} \forall q_{h} \in S_{0, *}^{- 1} (M) \end{matrix}

Note: here, we refer to the broken $H^{1}$ norm and $b$ bilinear form:

∥ \underset{―}{w} ∥_{H^{1} (M)}^{2} := \sum_{K \in M} ∥ \underset{―}{w} |_{K} ∥_{H^{1} (K)}^{2} ("broken" H^{1} -norm)

for proof, see script #todo see script
- for proof, we use the edge average interpolation operator:

Lemma 12.3.5.22. Continuity of

I_{CR}

With a constant $C_{I} > 0$ depending only on the shape-regularity measure of $M$

\begin{matrix} (12.3.5.23) & ∥ I_{CR} v ∥_{H^{1} (M)} \leq C_{I} ∥ v ∥_{H^{1} (Ω)} \forall v \in H^{1} (Ω) \end{matrix}

intermediate result:

Lemma 12.3.5.25. Fortin-projector property of

I_{CR}

$\begin{matrix} (12.3.5.26) & \sum_{K \in M} \int_{K} div ({I_{CR} v |}_{K}) q_{h} d x = \int_{Ω} div v q_{h} d x \forall v \in (H^{1} (Ω))^{2}, q_{h} \in S_{0}^{- 1} (M) \end{matrix}$

and remember:

Theorem 12.2.2.38. Existence of stable velocity potentials

$\exists C = C (Ω) > 0 : \forall q \in L_{*}^{2} (Ω) : \exists v \in (H_{0}^{1} (Ω))^{d} : q = div v \land ∥ v ∥_{H^{1} (Ω)} \leq C ∥ q ∥_{L^{2} (Ω)} .$

#todo review questions
#todo tldr not useful

Tldr: Stokes saddle point variational problem

$v \in (H_{0}^{1} (Ω))^{d}, p \in L_{*}^{2} (Ω) s.t.$

\begin{matrix} (12.2.2.19) & \begin{aligned} \int_{Ω} μ D v : D w d x + \int_{Ω} div w p d x & = \int_{Ω} f \cdot w d x \forall w \in (H_{0}^{1} (Ω))^{d}, \\ \int_{Ω} div v q d x & = 0 \forall q \in L_{*}^{2} (Ω) . \end{aligned} \end{matrix}

with a, b form:

\begin{matrix} (12.3.0.2) & \begin{matrix} v \in U := (H_{0}^{1} (Ω))^{d} \\ p \in Q := L_{*}^{2} (Ω) \end{matrix} : \begin{aligned} a (v, w) + b (w, p) & = ℓ (w) \forall w \in U, \\ b (v, q) & = 0 \forall q \in Q . \end{aligned} \end{matrix}

13. Finite-Element Exterior Calculus (FEEC)

13.1 (Differential) Forms

Video 13.1.2 Fields and Integral Forms/Currents (28 min)

Electric Field ( $\underset{―}{e}$ ): Measured via work along directed paths $\to$ Fundamental nature is a 1-form.
Magnetic Induction ( $\underset{―}{b}$ ): Measured via flux through oriented surfaces $\to$ Fundamental nature is a 2-form.

-> $\underset{―}{e}$ and $\underset{―}{b}$ have totally different physical natures and should be discretized differently.

Orientation:
Pasted image 20260510135853.png

S_{ℓ} (Ω)

$S_{ℓ} (Ω) ≜$ set of all $ℓ$ -dimensional oriented "surfaces" (sub-manifolds) in $Ω \subset R^{d}$

Notion 13.1.2.7. Integral

ℓ

-forms (IF)/

ℓ

-currents

An (integral) $ℓ$ -form/ $ℓ$ -current $\underset{―}{ω}$ , $0 \leq ℓ \leq n$ , on $\underset{ambient domain}{\underset{⏟}{Σ}}$ is a continuous $(*)$ and additive $(* *)$ mapping $\underset{―}{ω} : ⋃_{k \leq ℓ} S_{k} (Σ) \to R$ satisfying

\begin{matrix} (13.1.2.8) & \begin{aligned} \underset{―}{ω} (- Γ) = - \underset{―}{ω} (Γ) & \forall Γ \in S_{ℓ} (Σ), \\ \underset{―}{ω} (γ) = 0 & \forall γ \in S_{k} (Σ), k < ℓ, \end{aligned} \end{matrix}

where $- Γ$ denotes the $ℓ$ -dimensional surface/sub-manifold with flipped orientation.

on domain $Ω \subset R^{3}$ :
Pasted image 20260510135621.png
Notation: $I^{ℓ} (Ω) ≜$ vector space of $ℓ$ -forms

$ω \in I^{ℓ} (Ω) : ⟨ ω, T ⟩ := ω (T)$

Field Type	Form Degree	Integrated Over...	Examples
Scalar Field	0-form	Points	Temperature, Electric Potential
Circulation	1-form	Curves	Electric Field ( $\underset{―}{e}$ ), Magnetic Field ( $\underset{―}{h}$ )
Flux	2-form	Surfaces	Magnetic Induction ( $\underset{―}{b}$ ), Heat Flux, Current Density ( $\underset{―}{j}$ )
Density	3-form	Volumes	Mass Density, Charge Density

Integral

ℓ

-form /

ℓ

-current,

I^{ℓ} (Ω)

An integral $ℓ$ -form $\underset{―}{ω}$ is a mapping $\underset{―}{ω} : S_{ℓ} (Ω) \to R$ that is:

Continuous: Small surface deformations lead to small changes in value.
Additive: $\underset{―}{ω} (Γ_{1} \cup Γ_{2}) = \underset{―}{ω} (Γ_{1}) + \underset{―}{ω} (Γ_{2})$ for disjoint $Γ_{1}, Γ_{2}$ .
Orientation-reversing: $\underset{―}{ω} (- Γ) = - \underset{―}{ω} (Γ)$ .
Dimension-consistent: $\underset{―}{ω} (γ) = 0$ for surfaces $γ$ with $\dim (γ) < ℓ$ .

$I^{ℓ} (Ω)$ denotes the vector space of these forms.

Notation: as a duality pairing: $⟨ \underset{―}{ω}, Γ ⟩ := \underset{―}{ω} (Γ)$ .

Video 13.1.3 Differential Forms (56 min)

we now go from global measurements of fields (integral forms) to local characterization:

$e (\vec{x}, t)$ can be regarded as a linear mapping from displacements into $R$
$b (\vec{x}, t)$ should be read as an anti-symmetric bilinear form $(\vec{β φ}, \vec{m}) \to R$

Definition 13.1.3.4. Alternating

ℓ

-linear form

An alternating $ℓ$ -linear form $μ$ , $ℓ \in N$ , on a real vector space $V$ is a mapping

μ : \underset{ℓ times}{\underset{⏟}{V \times V \times \dots \times V}} \to R,

which

is linear in each of its $ℓ$ arguments, and which
changes sign, when any two arguments are swapped ("alternating").
$⟺ = 0$ when any two arg. vectors are the same
$⟺ = 0$ when args. are linearly depependent

it returns the signed volume of the parallelotope spanned by its arguments

Notation: $\overset{ℓ}{⋀} (V) ≜$ vector space

$\overset{0}{⋀} (V) = R$
$\overset{1}{⋀} (V) ≜$ linear forms on $V$
dimension formula: ( $n := dim (V)$ , is $# ℓ$ -subsets of a basis of $V$ )

\overset{ℓ}{⋀} (V) = (\binom{n}{ℓ}), 0 if ℓ > n

Vector proxies: differential forms in $R^{3}$ are identified with classical fields via proxies

0-form $ω$ : Proxy is a scalar function $u (x)$ .
1-form $ω$ : Proxy is a vector field $u (x)$ such that $ω (x) (v) = u (x) \cdot v$ .
2-form $ω$ : Proxy is a vector field $u (x)$ such that $ω (x) (v_{1}, v_{2}) = u (x) \cdot (v_{1} \times v_{2})$ .
3-form $ω$ : Proxy is a scalar function $u (x)$ such that $ω (x) (v_{1}, v_{2}, v_{3}) = u (x) det (v_{1}, v_{2}, v_{3})$ .

Integration over manifolds: integral of a diff. form $ω$ over an oriented $ℓ$ -dimensional manifold $Γ$ :

\begin{matrix} (13.1.3.34) & \int_{Γ} ω = {\begin{cases} ω (x) & for ℓ = 0, Γ = {x}, x \in Ω, \\ \int_{Γ} v.p. (ω) (x) \cdot d \vec{s} (x) & for ℓ = 1, \\ \int_{Γ} v.p. (ω) (x) \cdot n (x) d S (x) & for ℓ = 2, flux integral \\ \int_{Γ} v.p. (ω) (x) d x & for ℓ = 3. \end{cases} \end{matrix}

Intermediate

Continuous differential $ℓ$ -forms $ω_{1} \in C^{0} Λ^{ℓ} ({\overset{―}{Ω}}_{1})$ , $ω_{2} \in C^{0} Λ^{ℓ} ({\overset{―}{Ω}}_{2})$ , give rise to a valid integral $ℓ$ -form on $Ω$ by piecewise integration, if and only if

\begin{matrix} (13.1.3.37) & ω_{1} (x) (t_{1}, \dots, t_{ℓ}) = ω_{2} (x) (t_{1}, \dots, t_{ℓ}) \forall t_{1}, \dots, t_{ℓ} \in T_{x} (Σ) \forall x \in Σ . \end{matrix}

"For piecewise smooth forms to define a valid global integral form, they must agree on the tangent space of the interface $Σ$ "

0-forms: Must be globally continuous.
1-forms: Must have tangential continuity (jump in $u \times n = 0$ ).
2-forms: Must have normal continuity (jump in $u \cdot n = 0$ ).
3-forms: No continuity required.

Video 13.1.4 Exterior Calculus (70 min)

pullbacks:
Pasted image 20260512143308.png

Definition 13.1.4.4. Pullback/transformation for differential forms

Given a $C^{1}$ -mapping $Φ : \hat{Ω} \to Ω$ , $\hat{Ω}, Ω \subset R^{d}$ domains, the pullback $Φ^{*} : Λ^{ℓ} (Ω) \to Λ^{ℓ} (\hat{Ω})$ , $0 \leq ℓ \leq d$ , is defined as

(Φ^{*} ω) (\hat{x}) (v_{1}, \dots, v_{ℓ}) = ω (Φ (\hat{x})) (D Φ (\hat{x}) v_{1}, \dots, D Φ (\hat{x}) v_{ℓ}) \forall v_{i} \in R^{d} \forall \hat{x} \in \hat{Ω} .

Theorem 13.1.4.5. Invariance of integrals under pullback

Using the notations of Def. 13.1.4.4, the pullback $Φ^{*} : Λ^{ℓ} (Ω) \to Λ^{ℓ} (\hat{Ω})$ satisfies

\begin{matrix} (13.1.4.6) & \int_{\hat{Γ}} (Φ^{*} ω) = \int_{Φ (\hat{Γ})} ω \forall \hat{Γ} \in S_{ℓ} (\hat{Ω}), \end{matrix}

which ensures compatibility with the pullback of integral forms.

Vector Proxy Formulas (in 3D):

0-forms: $\hat{u} = u (Φ (\hat{x}))$ (Standard function pullback).
1-forms: $\hat{u} = (D Φ)^{⊤} u$ (Covariant transformation).
2-forms: $\hat{u} = det (D Φ) (D Φ)^{- 1} u$ (Piola transformation).

orientation of $\partial \sum (t)$ = induced orientation
Pasted image 20260512144020.png

Definition 13.1.4.27. Exterior derivative of integral forms

Let $Σ$ be an $n$ -dimensional manifold/domain, $n \in N$ . For $0 \leq ℓ < n$ the exterior derivative operators (for integral forms)

{\underset{―}{d}}_{ℓ} : I^{ℓ} (Σ) \to I^{ℓ + 1} (Σ)

are defined as

\begin{matrix} (13.1.4.28) & ⟨ {\underset{―}{d}}_{ℓ} \underset{―}{ω}, Γ ⟩ = ⟨ \underset{―}{ω}, \partial Γ ⟩ \forall \underset{―}{ω} \in I^{ℓ} (Σ) \forall Γ \in S_{ℓ + 1} (Σ), \end{matrix}

where the boundary $\partial Γ$ is endowed with the induced orientation.

^{"} ⟨ {\underset{―}{d}}_{ℓ} \cdot, \cdot ⟩ = ⟨ \cdot, \partial \cdot ⟩^{"}

fundamental identity

\begin{aligned} \partial \partial Γ & = \emptyset & ⟹ & {\underset{―}{d}}_{ℓ + 1} \circ {\underset{―}{d}}_{ℓ} = 0 \\ \partial Φ (Γ) & = Φ (\partial Γ) & ⟹ & Φ^{*} \circ {\underset{―}{d}}_{ℓ} = {\underset{―}{d}}_{ℓ} \circ Φ^{*} \end{aligned}

Definition 13.1.4.38. Exterior derivative for

C^{1}

differential forms

The exterior derivative of a continuously differentiable differential form $ω \in C^{1} Λ^{ℓ} (Ω)$ of order $ℓ \in {0, \dots, d - 1}$ on a domain $Ω \subset R^{d}$ is the continuous differential $ℓ + 1$ -form $d_{ℓ} ω$ defined as $^{a}$

\begin{matrix} (13.1.4.39) & d_{ℓ} ω (x) (v_{1}, \dots, v_{ℓ + 1}) := \sum_{k = 1}^{ℓ + 1} (- 1)^{k - 1} D ω (x) v_{k} (v_{1}, \dots, {\overset{ˇ}{v}}_{k}, \dots, v_{ℓ + 1}) \end{matrix}

for all $x \in Ω$ and "tangent vectors" $v_{i} \in R^{d}$ .

$^{a}$ In this formula $D ω : Ω \to L (R^{d}, Λ^{ℓ} (R^{d}))$ is the (Fréchet) derivative of $ω$ , which maps into the space of linear mappings $R^{d} \to Λ^{ℓ} (R^{d})$ .

Pasted image 20260512145150.png

Theorem 13.1.4.57. Generalized Stokes theorem

For every domain $Ω \subset R^{d}$ and $ℓ \in {0, \dots, d - 1}$ holds true

\begin{matrix} (13.1.4.58) & \int_{Γ} d_{ℓ} ω = \int_{\partial Γ} ω \forall ω \in C^{1} Λ^{ℓ} (Σ), \forall Γ \in S_{ℓ + 1} (Ω) . \end{matrix}

3D Vector Proxy Identifications:

$d_{0} \leftrightarrow grad$ (0-form $\to$ 1-form).
$d_{1} \leftrightarrow curl$ (1-form $\to$ 2-form).
$d_{2} \leftrightarrow div$ (2-form $\to$ 3-form).

Theorem 13.1.4.64.

d \circ d = 0

For every domain $Ω \subset R^{d}$ and $ℓ \in {0, \dots, d - 2}$ holds true

\begin{matrix} (13.1.4.65) & d_{ℓ + 1} \circ d_{ℓ} = 0 on C^{2} Λ^{ℓ} (Ω) \Leftrightarrow R (d_{ℓ}) \subset N (d_{ℓ + 1}), \end{matrix}

that is, the image of $d_{ℓ}$ lies in the kernel of $d_{ℓ + 1}$ .

curl \circ grad = 0 div \circ curl = 0

Theorem 13.1.4.81. Existence of potentials

Let $Ω \subset R^{d}$ be a domain and $1 \leq ℓ \leq d$ . Topological assumption: If every oriented, compact and closed $^{⋆}$ $ℓ$ -dimensional sub-manifold $\in S_{ℓ} (Ω)$ is the boundary of an $ℓ + 1$ -dimensional sub-manifold $\in S_{ℓ + 1} (Ω)$ , then

\begin{matrix} (13.1.4.82) & ω \in C^{1} Λ^{ℓ} (Ω), d_{ℓ} ω = 0 \Rightarrow \exists η \in C^{1} Λ^{ℓ - 1} (Ω) such that \underset{a potential}{\underset{⏟}{d_{ℓ - 1} η = ω}} . \end{matrix}

$^{⋆}$ We call $Γ \in S_{ℓ} (Ω)$ closed, if $\partial Γ = \emptyset$

Definition 13.1.4.102. Exterior product/wedge product of alternating multilinear forms

For a vector space $V$ the exterior product (also called wedge product)

\land : Λ^{ℓ} (V) \times Λ^{k} (V) \to Λ^{ℓ + k} (V)

( $Λ^{n} (V)$ is the vector space of alternating $n$ -linear forms on $V$ , see Def. 13.1.3.4) is defined as

(α \land β) (v_{1}, \dots, v_{ℓ + k}) = \frac{1}{ℓ! k!} \sum_{σ \in Π_{ℓ + k}} \underset{Permutations of (1, \dots, ℓ + k)}{\underset{⏟}{sgn (σ) α (v_{σ (1)}, \dots, v_{σ (ℓ)}) β (v_{σ (ℓ + 1)}, \dots, v_{σ (ℓ + k)})}}

for all $v_{1}, \dots, v_{ℓ + k} \in V$ .

Exterior Product $\land$ is associative, bilinear & (anti-)commutative

\begin{matrix} (13.1.4.113) & v.p. (ω \land η) = {\begin{cases} v.p. (ω) v.p. (η) & for ω \in C^{0} Λ^{0} (Ω), η \in C^{0} Λ^{k} (Ω), k = 0, 1, 2, 3, \\ v.p. (ω) \times v.p. (η) & for ω \in C^{0} Λ^{1} (Ω), η \in C^{0} Λ^{1} (Ω), \\ v.p. (ω) \cdot v.p. (η) & for ω \in C^{0} Λ^{1} (Ω), η \in C^{0} Λ^{2} (Ω) . \end{cases} \end{matrix}

Theorem 13.1.4.117. Pullback and exterior product/wedge product

Let $\hat{Ω}, Ω \subset R^{d}$ be two domains connected by a $C^{1}$ -mapping $Φ : \hat{Ω} \to Ω$ . Then the pullback $Φ^{*}$ of differential forms commutes with the exterior product of differential forms:

\begin{matrix} (13.1.4.118) & Φ^{*} (ω \land η) = (Φ^{*} ω) \land (Φ^{*} η) \forall ω \in C^{0} Λ^{ℓ} (Ω), η \in C^{0} Λ^{k} (Ω), ℓ, k \in N_{0} . \end{matrix}

Theorem 13.1.4.119. Product rule for exterior derivatives

For any domain $Ω \subset R^{d}$ and $ℓ, k \in {0, \dots, d}$

\begin{matrix} (13.1.4.120) & d_{ℓ + k} (ω \land η) = d_{ℓ} ω \land η + (- 1)^{ℓ} ω \land d_{k} η \forall ω \in C^{1} Λ^{ℓ} (Ω), η \in C^{1} Λ^{k} (Ω) . \end{matrix}

example from video: 20260513_NumPDE_video13.1.5_maxwell

Video 13.1.5 Sobolev Spaces of Differential Forms (18 min)

Definition 13.1.5.1. Space of square-integrable differential forms

Let $u = [u_{1} \dots u_{m}]^{⊤} =: Ω \to R^{(\binom{d}{ℓ})}$ by any "component representation" of a differential $ℓ$ -form $ω \in Λ^{ℓ} (Ω), ℓ \in {0, \dots, d}$ .
Then we define the Hilbert space of square-integrable differential $ℓ$ -forms

L^{2} Λ^{ℓ} (Ω) := {ω \in Λ^{ℓ} (Ω) : u_{i} \in L^{2} (Ω)} .

Definition 13.1.5.2. Sobolev space of differentiable differential forms

For a domain $Ω \subset R^{d}$ and $ℓ \in {0, \dots, d - 1}$ we define

\underset{"finite-energy fields"}{\underset{⏟}{H Λ^{ℓ} (d, Ω)}} := {ω \in L^{2} Λ^{ℓ} (Ω) : d_{ℓ} ω \in L^{2} Λ^{ℓ + 1} (Ω)},

which becomes a Hilbert space with the basis-dependent norm

∥ ω ∥_{H Λ^{ℓ} (d, Ω)}^{2} := ∥ ω ∥_{L^{2} Λ^{ℓ} (Ω)}^{2} + ∥ d_{ℓ} ω ∥_{L^{2} Λ^{ℓ + 1} (Ω)}^{2} .

3. Vector Proxy Identifications (in 3D)

Form Degree	Sobolev Space	Vector Analytic Proxy	Property
0-form	$H Λ^{0} (d, Ω)$	$H^{1} (Ω)$	$ω, grad (ω) \in L^{2}$
1-form	$H Λ^{1} (d, Ω)$	$H (curl, Ω)$	$u, curl (u) \in L^{2}$
2-form	$H Λ^{2} (d, Ω)$	$H (div, Ω)$	$u, div (u) \in L^{2}$

Pasted image 20260512152537.png|300

Theorem 13.1.5.10. Transmission condition for

H Λ^{ℓ} (d, Ω)

Let $ω_{1} \in C^{1} Λ^{ℓ} ({\overset{―}{Ω}}_{1})$ and $ω_{2} \in C^{1} Λ^{ℓ} ({\overset{―}{Ω}}_{2}), ℓ \in {0, \dots, d - 1}$ , be two continuously differentiable $ℓ$ -forms that are continuous up to boundary of $Ω_{i}$ and set

ω := {\begin{cases} ω_{1} & in Ω_{1}, \\ ω_{2} & in Ω_{2} . \end{cases}

Then

ω \in H Λ^{ℓ} (d, Ω) \Leftrightarrow \underset{space of tangent vectors}{\underset{⏟}{ω_{1} (x) (t_{1}, \dots, t_{ℓ}) = ω_{2} (x) (t_{1}, \dots, t_{ℓ})}} \forall t_{i} \in T_{x} (Σ), x \in Σ :

$ω \in H Λ^{ℓ} (d, Ω)$ , if and only if the restrictions of $ω_{1}$ and $ω_{2}$ to $Σ$ coincide as differential $ℓ$ -forms on $Σ$ .

Practical 3D Continuity:

0-forms: Global continuity.
1-forms: Tangential continuity ( $u_{1} \times n = u_{2} \times n$ ).
2-forms: Normal continuity ( $u_{1} \cdot n = u_{2} \cdot n$ ).

For

d > 1

, integration over sub-manifolds is not a bounded linear functional on these Sobolev spaces (similar to how point evaluation is not bounded in

H^{1}

)

Tldr: Sobolev Spaces of Differential Forms

$L^{2}$ -Space

coefficient functions $u_{i} \in L^{2} (Ω)$
intrinsic space (basis-independent), norm basis-dependent

Sobolev Space of "Finite Energy" Fields

H Λ^{ℓ} (d, Ω) := {ω \in L^{2} Λ^{ℓ} (Ω) : d_{ℓ} ω \in L^{2} Λ^{ℓ + 1} (Ω)}

Norm: $∥ ω ∥_{H Λ^{ℓ}}^{2} := ∥ ω ∥_{L^{2}}^{2} + ∥ d_{ℓ} ω ∥_{L^{2}}^{2}$ .
The exterior derivative $d_{ℓ}$ is a bounded linear operator mapping $H Λ^{ℓ} \to H Λ^{ℓ + 1}$ .

Zero Boundary Conditions ( $H_{0} Λ^{ℓ}$ )

forms with vanishing tangential trace on $\partial Ω$
example: $e \times n = 0$ (PEC boundary conditions for the electric field, 13.1.3.43)

13.2 Discrete Differential Forms

13.2.1 Cochain Calculus (67 min)

Introduce discrete integral $ℓ$ -forms as mappings from a finite (and suitable: meaningful $\partial$ ) subset of $S_{ℓ} (Ω)$ into $R$ .

reminder: Integral form: assigns value -> oriented surface

Cochain ("discrete integral form") Assigns values to specific $ℓ$ -facets of a mesh

$ℓ$ -Cochain ( $\vec{ω}$ ): A mapping from the set of $ℓ$ -facets $F_{ℓ} (M)$ of a mesh to $R$ .

in 3D:

Cochain	assignes values to...
0-cochain	nodes
1-cochain	edges
2-cochain	faces
3-cochain	cells

Because there are a finite number of facets ( $N_{ℓ}$ ), an $ℓ$ -cochain $\vec{ω}$ can be identified simply as a vector in $R^{N_{ℓ}}$ .

Definition 13.2.1.1. Generalized oriented meshes

A generalized oriented mesh $M$ of a domain $Ω \subset R^{d}$ is a collection of sets of $ℓ$ -facets,

\begin{matrix} (13.2.1.1) & M := {(F_{ℓ} (M))}_{ℓ = 0}^{d} \end{matrix}

such that

$F_{ℓ} (M) \subset S_{ℓ} (Ω)$ : every $ℓ$ -facet $f \in F_{ℓ} (M)$ is an oriented (*), relatively open, bounded, $ℓ$ -dimensional, non-degenerate $^{a}$ $C_{pw}^{1}$ -surface $\subset R^{d}$ ,
the facets form a partition of $\overset{―}{Ω}$ :

\overset{―}{Ω} = ⋃_{ℓ = 0}^{d} ⋃_{F \in F_{ℓ} (M)} F, f \cap f^{'} = \emptyset \forall f \in F_{j} (M), f^{'} \in F_{k} (M), k, j \in {0, \dots, d}

for every $F \in F_{ℓ} (M), 0 < ℓ \leq d$ , its boundary $\partial F$ is the union of the closures of finitely many $ℓ - 1$ -facets:

\forall F \in F_{ℓ} (M) : \exists {f_{1}, \dots, f_{m}} \subset F_{ℓ - 1} (M) such that \partial F = {\overset{―}{f}}_{1} \cup \dots \cup {\overset{―}{f}}_{m}

the intersection of the closures of any two $ℓ$ -facets is an $ℓ - 1$ -facet,

\forall ℓ \in {1, \dots, d}, \forall f, f^{'} \in F_{ℓ} (M) : \overset{―}{f} \cap {\overset{―}{f}}^{'} \in F_{ℓ - 1} (M)

for every $f \in F_{ℓ} (M), 0 \leq ℓ < d$ , there is a $F \in F_{ℓ + 1} (M)$ such that $f \subset \partial F$ .

$^{a}$ An $ℓ$ -dimensional sub-manifold $\subset R^{d}$ is called non-degenerate, if it has non-zero $ℓ$ -dimensional volume.

For cochain to define values consistently, the mesh must be oriented. Each facet has an intrinsic orientation (e.g., the direction of an edge). When we look at the boundary of a face, the face's orientation induces a direction on its edges.

Definition 13.2.1.13. Relative orientation of mesh facets

Given $F \in F_{ℓ} (M), 0 < ℓ \leq d$ , and $f \in F_{ℓ - 1} (M)$ , we define their relative orientation as

o_{r} (f, F) := {\begin{cases} 1 & , if f \subset F and the orientations of f and \partial F agree, \\ - 1 & , if f \subset F and the orientations of f and \partial F are opposite, \\ 0 & , if f ⊄ \partial F . \end{cases}

example:
Pasted image 20260513172938.png|300

Definition 13.2.1.16. Exterior derivative for cochains

Given a generalized oriented mesh $M$ of $Ω \subset R^{d}$ and $ℓ \in {0, \dots, d - 1}$ , the exterior derivative $\underset{linear operator}{\underset{⏟}{{\vec{d}}_{ℓ} : C^{ℓ} (M) \to C^{ℓ + 1} (M)}}$ is defined by

\begin{matrix} (13.2.1.17) & ⟨ {\vec{d}}_{ℓ} \vec{ω}, F ⟩ = \sum_{f \in F_{ℓ} (M)} o_{r} (f, F) ⟨ \vec{ω}, f ⟩ \forall F \in F_{ℓ + 1} (M) . \end{matrix}

incidence matrix representation $D_{ℓ}$ :

matrix representation of the discrete exterior derivative, $(D_{ℓ})_{i, j} = o r (f_{j}, F_{i})$

Matrix	Calculus Proxy	Maps...
$D_{0}$	Discrete Gradient	Nodes $\to$ Edges
$D_{1}$	Discrete Curl	Edges $\to$ Faces
$D_{2}$	Discrete Divergence	Faces $\to$ Cells

Theorem 13.2.1.21.

\vec{d} \circ \vec{d} = 0

For any generalized mesh $M$ of a domain $Ω \subset R^{d}$ and $ℓ \in {0, \dots, d - 2}$ holds true

{\vec{d}}_{ℓ + 1} \circ {\vec{d}}_{ℓ} = 0 ⟺ D_{ℓ + 1} D_{ℓ} = O ⟺ R (D_{ℓ}) \subset N (D_{ℓ + 1}) .

the boundary of a boundary is empty

Theorem 13.1.4.81. Existence of potentials

Let $Ω \subset R^{d}$ be a domain and $1 \leq ℓ \leq d$ . If

(topological assumption:) every oriented, compact and closed $ℓ$ -dimensional sub-manifold of $\in S_{ℓ} (Ω)$ is the boundary of an $ℓ + 1$ -dimensional sub-manifold $\in S_{ℓ + 1} (Ω)$

then

\begin{matrix} (13.1.4.82) & ω \in C^{1} Λ^{ℓ} (Ω), d_{ℓ} ω = 0 ⟹ \exists η \in C^{1} Λ^{ℓ - 1} (Ω) such that d_{ℓ - 1} η = ω . \end{matrix}

Theorem 13.2.1.22. Existence of cochain potentials

Let the domain $Ω \subset R^{d}$ satisfy the assumptions of Thm. 13.1.4.81. Then for any generalized oriented mesh $M$ of $Ω$ and for any $1 \leq ℓ \leq d$ we have

\vec{ω} \in C^{ℓ} (M) : {\vec{d}}_{ℓ} \vec{ω} = 0 ⟹ \exists \vec{η} \in C^{ℓ - 1} (M) such that {\vec{d}}_{ℓ - 1} \vec{η} = \vec{ω}

\begin{matrix} (13.2.1.23) & \Leftrightarrow R ({\vec{d}}_{ℓ - 1}) = N ({\vec{d}}_{ℓ}) \Leftrightarrow R (D_{ℓ - 1}) = N (D_{ℓ}) . \end{matrix}

From now on,

C^{0}

is not a continuous function but an 0-cochain!

Summary for 3D till now

Category	Degree $ℓ = 0$	Degree $ℓ = 1$	Degree $ℓ = 2$	Degree $ℓ = 3$
Physical Nature	Scalar Value	Circulation	Flux	Density
Mathematical Nature	Real Number	Linear Form (Functional)	Alternating Bilinear Form	Alternating Trilinear Form
Local Action $ω (x) (\cdot)$	Scalar value $c \in R$	Maps 1 vector: $v \mapsto R$	Maps 2 vectors: $(v_{1}, v_{2}) \mapsto R$	Maps 3 vectors: $(v_{1}, v_{2}, v_{3}) \mapsto R$
Integral Form ( $I^{ℓ}$ )	Value at a point	Work along a curve	Flow through a surface	Amount in a volume
Differential Form ( $Λ^{ℓ}$ )	Scalar function $u (x)$	1-linear form (Linear functional)	Alternating 2-linear form	Alternating 3-linear form
Alternating Property	N/A	(Linear by default)	$ω (v_{1}, v_{2}) = - ω (v_{2}, v_{1})$	Swapping any two vectors flips the sign
3D Vector Proxy ( $v.p.$ )	Scalar function $u (x)$	Vector field $u (x)$	Vector field $u (x)$	Scalar function $u (x)$
Proxy Evaluation	$u (x)$	Dot product: $u (x) \cdot v$	Cross product: $u (x) \cdot (v_{1} \times v_{2})$	Determinant: $u (x) det (v_{1}, v_{2}, v_{3})$
Integration (Eq. 13.1.3.34)	Value at point $x$	Path integral: $\int_{Γ} u \cdot d s$	Flux integral: $\int_{Γ} u \cdot n d S$	Volume integral: $\int_{Γ} u d x$
Exterior Deriv. ( $d_{ℓ}$ )	Gradient	Curl	Divergence	Zero ( $d \circ d = 0$ )
Sobolev Space	$H^{1} (Ω)$	$H (curl, Ω)$	$H (div, Ω)$	$L^{2} (Ω)$
Mesh Entity (Facet)	Node (Vertex)	Edge	Face	Cell
Cochain ( $C^{ℓ}$ )	Values on nodes	Values on edges	Values on faces	Values on cells

13.2.2 Whitney Forms (70 min)

Goal:

\begin{matrix} ℓ -cochain \vec{ω} \\ (mapping F_{ℓ} (M) \to R) \end{matrix} ⟶ \begin{matrix} matching integral ℓ -form \underset{―}{ω} \\ (mapping S_{ℓ} (Ω) \to R) \end{matrix}

interpolation condition: $⟨ \underset{―}{ω}, f ⟩ = ⟨ \vec{ω}, f ⟩ \forall f \in F_{ℓ} (M)$

Plaintext: goal is to turn cochain values into integral form that can be localized to a differential form

Local Whitney Forms (Basis Functions) On a simplex $K$ with barycentric coordinates $λ_{i}$ , the local Whitney $ℓ$ -forms are the basis for reconstruction:

$ℓ = 0$ : $β_{i}^{0} = λ_{i}$ (Standard "tent" or P1 functions).
$ℓ = 1$ : $β_{i, j}^{1} = λ_{i} d λ_{j} - λ_{j} d λ_{i}$ (Edge-associated basis).
Key Property: They satisfy the cardinal basis property: $\int_{f^{'}} β_{f}^{ℓ} = δ_{f f^{'}}$ .

general formula for d-simplex $K$ :

\begin{matrix} (13.2.2.17) & (W_{ℓ, K} \vec{ω}) (x) = ℓ! \sum_{f \in F_{ℓ} (K)} \underset{β_{f, K}^{ℓ}}{\underset{⏟}{\sum_{j = 0}^{ℓ} (- 1)^{j} (λ_{i_{j}} d_{0} λ_{i_{0}} \land \dots \land d_{0} λ_{i_{j}} \land \dots \land d_{0} λ_{i_{ℓ}})}} \cdot ⟨ \vec{ω}, f ⟩, \end{matrix}

Lemma 13.2.2.2.22. Local interpolation commutes with exterior derivative

The local interpolant of the exterior derivative of an $ℓ$ -cochain is the same as the exterior derivative of the interpolant of that cochain,

d_{ℓ} (W_{ℓ, K} \vec{ω}) = W_{ℓ + 1, K} (\vec{d} \vec{ω}) \begin{matrix} \forall \vec{ω} \in C^{ℓ} (M), \\ \forall K \in F_{d} (M) \end{matrix} \Leftrightarrow \begin{matrix} C^{ℓ} (M) & \overset{W_{ℓ, K}}{\to} & C^{\infty} Λ^{ℓ} (K) \\ {\vec{d}}_{ℓ} ↓ & ↓ d_{ℓ} \\ C^{ℓ + 1} (M) & \to_{W_{ℓ + 1, K}}^{} & C^{\infty} Λ^{ℓ + 1} (K) \end{matrix} .

d \circ W_{ℓ} = W_{ℓ + 1} \circ d

-> Whitney forms are not just any interpolant; they are structure-preserving

taking the derivative of the reconstructed field is the same as reconstructing the field from the discrete derivative of the cochain

\begin{matrix} (13.2.2.37) & span {v. p. (β_{f, K}^{2})}_{f \in F_{2} (K)} = {x \in K \mapsto a x + b, a \in R, b \in R^{3}} \subset (P_{1} (K))^{3} . \end{matrix}

vector proxy formulas for the spaces spanned by local Whitney $ℓ$ -forms on $K$

Degree ℓ	span{v. p.(βf,Kℓ)}f∈Fℓ(K)	Cochain coeff.
$ℓ = 0$	${x \in K \mapsto a + b \cdot x, a \in R, b \in R^{3}}$	$u \mapsto u (a_{i})$
$ℓ = 1$	${x \in K \mapsto a + b \times x, a, b \in R^{3}}$	$u \mapsto \int_{Σ [a_{i}, a_{j}]} u \cdot d \vec{s}$
$ℓ = 2$	${x \in K \mapsto a x + b, a \in R, b \in R^{3}}$	$u \mapsto \int_{Σ [a_{i}, a_{j}, a_{k}]} u \cdot n dS$
$ℓ = 3$	${x \mapsto c, c \in R}$	$u \mapsto \int_{K} u dx$

Definition 13.2.2.43. Generalized interpolation operator for cochains

Given an $ℓ$ -cochain $\vec{ω} \in C^{ℓ} (M)$ , we define the generalized interpolation operators,

{\underset{―}{W}}_{ℓ} : C^{ℓ} (M) \to I^{ℓ} (Ω), 0 \leq ℓ \leq d,

\begin{matrix} (13.2.2.44) & ({\underset{―}{W}}_{ℓ} \vec{ω}) (x) := (W_{ℓ, K} \vec{ω}) (x) for all x \in K, K \in F_{d} (M), \end{matrix}

with $W_{ℓ, K}$ . In particular, ${\underset{―}{W}}_{ℓ} \vec{ω}$ is a $M$ -piecewise smooth differential $ℓ$ -form.

Continuity: Inherently satisfy the tangential continuity required for membership in $H \overset{ℓ}{⋀} (d, Ω)$
Association: Each global basis function is associated with exactly one mesh facet (node, edge, or face) and is locally supported on adjacent cells.
Boundary Conditions: W0ℓ(M) is spanned by basis forms associated only with interior facets.

We call the integral $ℓ$ -form

\begin{matrix} (13.2.2.52) & β_{f}^{ℓ} := {\underset{―}{W}}_{ℓ} {\vec{ω}}_{f} \in I^{ℓ} (Ω) \end{matrix}

the $\underset{A global shape function !}{\underset{⏟}{Whitney ℓ -form}}$ associated with the $ℓ$ -facet $f$ .

Lemma

$\begin{matrix} (13.2.2.53) & supp (β_{f}^{ℓ}) = ⋃ {K \in F_{d} (M) : f \subset \overset{―}{K}} . \end{matrix}$

Intermediate

given a simplicial mesh $M$ of $Ω$ we have identified the Whitney finite-element space of degree $ℓ \in {0, \dots, d}$ ,

W^{ℓ} (M) := {\underset{―}{W}}_{ℓ} (C^{ℓ} (M)) = span ({β_{f}^{ℓ}}_{f \in F_{ℓ} (M)}) \subset H Λ^{ℓ} (d, Ω),

as a finite-element subspace of the Sobolev space $H Λ^{ℓ} (d, Ω)$ of (differential) $ℓ$ -forms on $Ω$ .

Convergence

reminder:if a finite-element space $V_{h}$ locally contains all polynomials of degree $p$ , then, for "sufficiently smooth" $u$ , and on uniformly shape regular meshes with meshwidth $h_{M}$ ,

inf_{v_{h} \in V_{h}} ∥ u - v_{h} ∥_{L^{2} (Ω)} = O (h_{M}^{p + 1}) for h_{M} \to 0,

here $p = 0$ , and thus algebraic cvg. of order 1

Theorem 13.2.2.59.

L^{2}

-norm best approximation estimates for

W^{ℓ} (M)

If $ω \in H^{1} Λ^{ℓ} (Ω)$ then

inf_{v_{h} \in W^{ℓ} (M)} ∥ ω - v_{h} ∥_{L^{2} Λ^{ℓ} (Ω)} \leq C_{ℓ} h_{M} ∥ ω ∥_{H^{1} Λ^{ℓ} (Ω)},

with a constant $C_{ℓ} > 0$ that depends only on $ℓ$ and the shape regularity measure (Def. 3.3.2.20) of $M$ .

13.3 Whitney FEM for Computational Electromagnetism

13.3.1 Energies and Linear Material Laws (35 min)

Maxwells equations are underdetermined -> need "closure equations"

We have two 1-forms ( $e, h$ ) and two 2-forms ( $d, b$ ), plus the current $j$ .
Material Laws provide the required link: they map 1-forms to 2-forms.

Field energy functionals

Electric Energy: $E_{e l} (e) = \frac{1}{2} m_{ϵ} (e, e)$ .
Magnetic Energy: $E_{m a g} (h) = \frac{1}{2} m_{μ} (h, h)$ .

Assumption 13.3.1.4. Linear materials

Both electric and magnetic field energies are homogeneous quadratic functionals, that is, there are bilinear forms $m_{ϵ} : L^{2} Λ^{1} (Ω) \times L^{2} Λ^{1} (Ω) \to R$ and $m_{μ} : L^{2} Λ^{1} (Ω) \times L^{2} Λ^{1} (Ω) \to R$ such that

\begin{matrix} (13.3.1.5) & E_{el} (e (t)) = \frac{1}{2} m_{ϵ} (e (t), e (t)) and E_{mag} (h (t)) = \frac{1}{2} m_{μ} (h (t), h (t)) . \end{matrix}

Moreover we expect $m_{ϵ}$ and $m_{μ}$ to be bounded from above and below

\begin{matrix} (13.3.1.6) & \exists ϵ^{-}, ϵ^{+} > 0 : ϵ^{-} ∥ ω ∥_{L^{2} Λ^{1} (Ω)}^{2} \leq m_{ϵ} (ω, ω) \leq ϵ^{+} ∥ ω ∥_{L^{2} Λ^{1} (Ω)}^{2} \forall ω \in L^{2} Λ^{1} (Ω), \end{matrix}

\begin{matrix} (13.3.1.7) & \exists μ^{-}, μ^{+} > 0 : μ^{-} ∥ ω ∥_{L^{2} Λ^{1} (Ω)}^{2} \leq m_{μ} (ω, ω) \leq μ^{+} ∥ ω ∥_{L^{2} Λ^{1} (Ω)}^{2} \forall ω \in L^{2} Λ^{1} (Ω) . \end{matrix}

-> we want $m_{ϵ}$ , $m_{μ}$ to be symmetric bilinear, bounded, uniformly positive definite

3D: For local linear materials, the energy is the integral of a local quadratic form:

Permittivity Tensor ( $ϵ$ ): $d (x) = ϵ (x) e (x)$ .
Permeability Tensor ( $μ$ ): $b (x) = μ (x) h (x)$ .
Nature: $ϵ (x)$ and $μ (x)$ are symmetric positive definite $3 \times 3$ matrices.

using the wedge product, we get weak formulation

\begin{aligned} \int_{Ω} d (t) \land e^{'} & = m_{ϵ} (e (t), e^{'}) \forall e^{'} \in L^{2} Λ^{1} (Ω) \\ \int_{Ω} b (t) \land h^{'} & = m_{μ} (h (t), h^{'}) \forall h^{'} \in L^{2} Λ^{1} (Ω) . \end{aligned}

this is a bijection

inserting this, we get:
Combining the expressions from image_1e0f33.png and image_1e0f17.png for the electromagnetic energy balance, the relationship is given by:

\begin{aligned} \underset{EM energy balance}{\underset{⏟}{\frac{d E^{⋆}}{d t} (t)}} & = - \int_{Ω} e (t) \land j (t) + \underset{Poynting vector}{\underset{⏟}{\int_{\partial Ω} h (t) \land e (t)}} \\ = power delivered by j + power flux through \partial Ω \\ \hat{=} Poynting’s theorem \end{aligned}

13.3.2.1 Electromagnetic Wave Equations (45 min)

IBVP for EM fields:

we have Maxwell's equations, Material Laws (weak form) and b.c.

\begin{aligned} ME: & d_{1} e (t) = - \frac{d}{d t} b (t) \in L^{2} Λ^{2} (Ω), d_{1} h (t) = \frac{d}{d t} d (t) + j (t) \in L^{2} Λ^{2} (Ω) \\ ML: & \int_{Ω} d (t) \land e^{'} = m_{ϵ} (e (t), e^{'}) \forall e^{'} \in L^{2} Λ^{1} (Ω) \\ \int_{Ω} b (t) \land h^{'} = m_{μ} (h (t), h^{'}) \forall h^{'} \in L^{2} Λ^{1} (Ω) \\ b.c.: & (artificial / idealization) \end{aligned}

In variational form:

derivation: time differentiation, substitution of time derivatives with maxwell's equations, integration by parts

Intermediate

Seek $e : [0, T] \to \underset{PEC b.c. are essential b.c.}{\underset{⏟}{H_{0} Λ^{1} (d, Ω)}}$ and $h : [0, T] \to L^{2} Λ^{1} (Ω)$ such that

\begin{matrix} (13.3.2.6) & \begin{aligned} m_{ϵ} (\frac{d e}{d t} (t), e^{'}) - \int_{Ω} h (t) \land d_{1} e^{'} & = - \int_{Ω} j (t) \land e^{'} \forall e^{'} \in H_{0} Λ^{1} (d, Ω), \\ m_{μ} (\frac{d h}{d t} (t), h^{'}) + \int_{Ω} d_{1} e (t) \land h^{'} & = 0 \forall h^{'} \in L^{2} Λ^{1} (Ω) . \end{aligned} \end{matrix}

when we combine the terms, we get:

\begin{aligned} (13.3.2.14) & m_{ϵ} (\frac{d e}{d t} (t), e (t)) + m_{μ} (\frac{d h}{d t} (t), h (t)) & = - \int_{Ω} j (t) \land e (t) . \\ ⟹ \frac{d}{d t} E_{t o t} (t) & = \underset{power injected by src. current}{\underset{⏟}{- \int_{Ω} j (t) \land e (t)}} \end{aligned}

compared to Poynting's theorem, power flux through $\partial Ω$ is missing because of the b.c. -> energy change is only caused by work of source current $j$

use material law, differentiate in time, test with exterior derivative, 2nd term vanish (both for 1st and 2nd formula), (integration by parts), extract PDE:

#todo missed some derivations that seem interesting

We get to

Seek $e :] 0, T [\to H_{0} (curl, Ω)$ such that

\int_{Ω} (ϵ (x) \frac{d^{2} \overset{ˇ}{e}}{d t^{2}} (x, t)) \cdot {\overset{ˇ}{e}}^{'} (x) d x + \int_{Ω} (μ^{- 1} (x) curl \overset{ˇ}{e} (x, t)) \cdot curl {\overset{ˇ}{e}}^{'} (x) d x = - \int_{Ω} \frac{\partial \overset{ˇ}{j}}{\partial t} (x, t) \cdot {\overset{ˇ}{e}}^{'} (x) d x

for all ${\overset{ˇ}{e}}^{'} \in H_{0} (curl, Ω)$ (13.3.2.32)

Note

Which is a form we've already seen in chapter 9.2:

\begin{matrix} (9.3.1.8) & u (t) \in V_{0} : m (\frac{d^{2} u}{d t^{2}}, v) + a (u, v) = ℓ (v) \forall v \in V_{0} . \end{matrix}

with function space $V_{0} := H_{0} (curl, Ω),$

-> electromagnetic wave equations have a lot in common with scalar wave equations

finite speed of propagation
expands with at most max speed (for $supp e_{0}, h_{0} \subset Ω_{0}$ ):

supp e (t), supp h (t) \subset {x \in R^{3} : dist (Ω_{0}, x) \leq t \cdot s_{max}}

domain of influence and dependence

13.3.2.3 Method of lines for Electromagnetic Wave Equations, Timestepping for Semi-Discrete Electromagnetic Wave Equations (45 min)

we start from #^13-3-2-6:

Seek $e : [0, T] \to H_{0} Λ^{1} (d, Ω)$ and $h : [0, T] \to L^{2} Λ^{1} (Ω)$ such that
$\begin{matrix} (13.3.2.6) & \begin{aligned} m_{ϵ} (\frac{d e}{d t} (t), e^{'}) - \int_{Ω} h (t) \land d_{1} e^{'} & = - \int_{Ω} j (t) \land e^{'} \forall e^{'} \in H_{0} Λ^{1} (d, Ω), \\ m_{μ} (\frac{d h}{d t} (t), h^{'}) + \int_{Ω} d_{1} e (t) \land h^{'} & = 0 \forall h^{'} \in L^{2} Λ^{1} (Ω) . \end{aligned} \end{matrix}$

method lines (MOL):

stage 1: spatial galerkin discretization

Trial (& test spaces) -> fin.dim. subspaces independent of time

For $H_{0} Λ^{1} (d, Ω)$ use $W_{0}^{1} (M)$
For $L^{2} Λ^{1} (Ω)$ use $\underset{M -p.w. constant 1-forms}{\underset{⏟}{P_{0} Λ^{1} (M)}}$

choose ordered bases

$W_{0}^{1} (M)$ : use Whitney 1-forms
$P_{0} Λ^{1} (M)$ : 3 basis 1-forms per cell, v.p. $[\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}], [\begin{matrix} 0 \\ 1 \\ 0 \end{matrix}], [\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}]$

MOL-ODE: Seek time-dependent basis expansion coefficient vectors $\vec{e} : [0, T] \to R^{n}$ and $\vec{h} : [0, T] \to R^{m}$ such that
$\begin{matrix} (13.3.2.36) & \begin{aligned} M_{ϵ} \frac{d \vec{e}}{d t} (t) - B^{⊤} \vec{h} (t) & = \vec{ψ} (t), \\ M_{μ} \frac{d \vec{h}}{d t} (t) + B \vec{e} (t) & = 0, \end{aligned} \end{matrix}$

preserves constraits
- $div (b) = 0$
- continuity equation
structure defines discrete Hamiltonian system -> reason why semi-descrete system (and Carank-Nicolson discretization, see later) conserve discrete $E_{t o t}$

some other notes:

Recall:

\begin{aligned} (13.3.2.21) & \int_{Ω} b (t) \land d_{0} v^{'} & = \underset{\frac{d}{d t} d_{2} b (t) = 0}{\underset{⏟}{\int_{Ω} b_{0} \land d_{0} v^{'}}} \forall v^{'} \in H Λ^{0} (d, Ω) and all times t, \\ (13.3.2.27) & \int_{Ω} \frac{d d}{d t} (t) \land d_{0} v^{'} & = \underset{continuity equ.}{\underset{⏟}{- \int_{Ω} j (t) \land d_{0} v^{'}}} \forall v^{'} \in H_{0} Λ^{0} (d, Ω) and for all times t . \end{aligned}

Observe:
- $d_{0} W_{0} Λ^{0} (M) \subset W_{0} Λ^{1} (M)$
- $d_{0} W Λ^{0} (M) \subset P_{0} Λ^{1} (M)$

stage 2: timestepping -> see chapter 9.3.4 (which was not part of the course)

temporal grid $0 = t_{0} < t < \dots < t_{M} = T$ , $T > 0$ final time
Crank-Nicolson timestepping (RK-SSM/implicit midpoint SSM)
- replace time derivatives with diff. quotients, averages

\begin{matrix} (13.3.2.51) & [\begin{matrix} M_{ϵ} & - \frac{1}{2} τ_{j} B^{⊤} \\ \frac{1}{2} τ_{j} B & M_{μ} \end{matrix}] [\begin{matrix} {\vec{e}}^{(j)} \\ {\vec{h}}^{(j)} \end{matrix}] = [\begin{matrix} M_{ϵ} & τ_{j} \frac{1}{2} B^{⊤} \\ - \frac{1}{2} τ_{j} B & M_{μ} \end{matrix}] [\begin{matrix} {\vec{e}}^{(j - 1)} \\ {\vec{h}}^{(j - 1)} \end{matrix}] + [\begin{matrix} \vec{ψ} (t_{j - \frac{1}{2}}) \\ 0 \end{matrix}] . \end{matrix}

fully implicit method (have to calculate LSE for every timestep)
discrete energy balance

\begin{aligned} (13.3.2.47) & E_{tot}^{(j)} & = E_{tot}^{(j - 1)} + τ_{j} \frac{1}{2} ({\vec{e}}^{(j)} + {\vec{e}}^{(j - 1)})^{⊤} \vec{ψ} (t_{j - \frac{1}{2}}), \\ E_{tot}^{(k)} & := \frac{1}{2} ({\vec{e}}^{(k)})^{⊤} M_{ϵ} {\vec{e}}^{(k)} + \frac{1}{2} ({\vec{h}}^{(k)})^{⊤} M_{μ} {\vec{h}}^{(k)} \\ = 1 / 2 m_{ϵ} (e_{h}^{(k)}, e_{h}^{(k)}) + 1 / 2 m_{μ} (h_{h}^{(k)}, h_{h}^{(k)}) \end{aligned}

timestepping alternative: Leapfrog
- use staggered temporal grid
only partially implicit (but also only conditionally stable, not discussed in this video)
- you only need to solve linear systems for $M_{ϵ}, M_{μ}$ separately

\begin{aligned} M_{ϵ} {\vec{e}}^{(j)} & = M_{ϵ} {\vec{e}}^{(j - 1)} + τ B^{⊤} {\vec{h}}^{(j - \frac{1}{2})} + τ \vec{ψ} (t_{j - \frac{1}{2}}), & (I) \\ M_{μ} {\vec{h}}^{(j + \frac{1}{2})} & = M_{μ} {\vec{h}}^{(j - \frac{1}{2})} - τ B {\vec{e}}^{(j)} . & (II) \end{aligned}

we cannot use Mass lumping, but not available for general for Whitney FE

boundary conditions:

space $W_{0}^{1} (M)$ imposes Perfect Electrical Conductor (PEC) b.c.

e \times n = 0 on \partial Ω

13.3.3 Magnetostatic

13.3.3.1 Magnetostatics: Formulations Based on Scalar Potentials (30 min)

static settings: fields do not depend on time, $\frac{d}{d t} = 0$

built-in constraints are no longer automatically enforced by reduced Maxwell's equations, need additional demands:

\begin{aligned} (13.3.3.3) & electrostatics: (FL): & d_{1} e = 0, & \underset{constraints in transient setting}{\underset{⏟}{d_{2} d = q}}, \\ (13.3.3.4) & magnetostatics : (AL): & \underset{necessary: d_{2} j = 0}{\underset{⏟}{d_{1} h = j}}, & d_{2} b = 0. \end{aligned}

$+$ demand: $supp (j) \subset Ω ⟹ \underset{[v.p.: j \cdot n = 0 on \partial Ω]}{\underset{⏟}{j |_{\partial Ω} = 0}}$

Magnetostatic BVP:

\begin{matrix} (13.3.3.9) & d_{1} h = j, d_{2} b = 0 in Ω \subset R^{3}, \underset{[P M C]}{\underset{⏟}{h |_{\partial Ω} = 0}} on \partial Ω, \end{matrix}

weak MLs:

\int_{Ω} b \land h^{'} = m_{μ} (h, h^{'}) \forall h^{'} \in L^{2} Λ^{1} (Ω), \int_{Ω} h \land b^{'} = m_{μ}^{*} (b, b^{'}) \forall b^{'} \in L^{2} Λ^{2} (Ω) .

two ways to convert to variational form:

Formulations based on Scalar Potential

Assumption 13.3.3.7. Simple topology of

Ω

We assume that the assumptions of Thm. 13.1.4.81 on $Ω$ for $ℓ = 1, 2$ are satisfied:

Every oriented closed 2-surface $Ω$ is the boundary of a sub-domain $\subset Ω$ .
Every closed directed curve in $Ω$ is the boundary of an oriented surface $\subset Ω$ .

and we get (after some transformations, i.b.p.):

Intermediate

Seek $v \in H_{0} Λ^{0} (d, Ω)$ such that

\begin{matrix} (13.3.3.12) & m_{μ} (d_{0} v, d_{0} v^{'}) = m_{μ} (h_{s}, d_{0} v^{'}) \forall v^{'} \in H_{0} Λ^{0} (d, Ω) \end{matrix}

For local, linear ML, v.p., $v \in H_{0}^{1} (Ω) :$

\begin{matrix} (13.3.3.13) & \int_{Ω} (μ (x) grad v (x)) \cdot grad v^{'} (x) d x = \int_{Ω} (μ (x) h_{s} (x)) \cdot grad v^{'} (x) d x \forall v^{'} \in H_{0}^{1} (Ω) . \end{matrix}

$\to$ 2nd-order elliptic BVP in weak form

now, we can replace

$H_{0}^{1} (Ω) \leftarrow W_{0}^{0} (M) = S_{1, 0}^{0} (M)$
seek $h_{s} \in H_{0} \overset{1}{⋀} (d_{1}, Ω)$ in $h_{s} \in W_{0}^{1} (M)$

option 2 for conversion to variational form not explained

13.3.3.2 Magnetostatics: Vector-Potential-Based Variational Formulations (40 min)

\begin{matrix} (13.3.3.21) & a \in H Λ^{1} (d, Ω) : m_{μ}^{*} (d_{1} a, d_{1} a^{'}) = \int_{Ω} j \land a^{'} \forall a^{'} \in H Λ^{1} (d, Ω) . \end{matrix}

LLML & vector proxy formulation:

\begin{aligned} \overset{ˇ}{a} \in H (curl, Ω) : & \int_{Ω} μ (x)^{- 1} curl \overset{ˇ}{a} (x) \cdot curl {\overset{ˇ}{a}}^{'} (x) d x \\ (13.3.3.22) & = & \int_{Ω} j (x) \cdot {\overset{ˇ}{a}}^{'} (x) d x \forall {\overset{ˇ}{a}}^{'} \in H (curl, Ω) . \end{aligned}

with b.c. hidden in formulation

#todo didn't understand this part

minute 13: great explanation of Lagrange multiplier

Augmented LSE :

\begin{matrix} (13.3.3.26) & \underset{regular square matrix}{\underset{⏟}{[\begin{matrix} A & Z \\ Z^{⊤} & O \end{matrix}]}} [\begin{matrix} \vec{ζ} \\ \vec{π} \end{matrix}] = [\begin{matrix} \vec{φ} \\ \vec{0} \end{matrix}] . \end{matrix}

applied to our variational form:
Seek $a \in H Λ^{1} (d, Ω)$ , $p \in H_{*} Λ^{0} (d, Ω) := {q \in H Λ^{0} (d, Ω) : \int_{Ω} v . p . (q) (x) d x = 0}$

\begin{matrix} (13.3.3.28) & \begin{aligned} m_{μ}^{*} (d_{1} a, d_{1} a^{'}) + \overset{p = 0!}{\overset{⏞}{m (a^{'}, d_{0} p)}} & = \int_{Ω} j \land a^{'} \forall a^{'} \in H Λ^{1} (d, Ω), \\ \underset{enforces m (\cdot, \cdot) -orthogonality to N (d_{1})}{\underset{⏟}{m (a, d_{0} p^{'})}} & = 0 \forall p^{'} \in H_{*} Λ^{0} (d, Ω) . \end{aligned} \end{matrix}

notice that $N (d_{0}) = P_{0} \overset{0}{⋀} (Ω)$

or LLML, vector proxy notation:
Seek $\overset{˘}{a} \in H (curl, Ω), \overset{˘}{p} \in H_{*}^{1} (Ω) := {v \in H^{1} (Ω) : \int_{Ω} v (x) d x = 0}$ such that

\begin{matrix} (13.3.3.29) & \begin{aligned} \int_{Ω} μ (x)^{- 1} curl \overset{˘}{a} (x) \cdot curl {\overset{˘}{a}}^{'} (x) d x + \int_{Ω} {\overset{˘}{a}}^{'} (x) \cdot grad \overset{˘}{p} (x) d x & = \int_{Ω} \overset{˘}{j} (x) \cdot {\overset{˘}{a}}^{'} (x) d x, \\ \int_{Ω} \overset{˘}{a} (x) \cdot grad {\overset{˘}{p}}^{'} (x) d x & = 0. \end{aligned} \end{matrix}

for all ${\overset{˘}{a}}^{'} \in H (curl, Ω), {\overset{˘}{p}}^{'} \in H_{*}^{1} (Ω)$ .

Galerking discretization

choose f.d. subspaces -> Whitney FE spaces
- also add vanishing mean constraint

Seek $a_{h} \in W^{1} (M)$ , $p_{h} \in W^{0} (M)$ , and $λ \in R$ such that
$\begin{matrix} (13.3.3.33) & \begin{aligned} m_{μ}^{*} (d_{1} a_{h}, d_{1} a_{h}^{'}) + m (a_{h}^{'}, d_{0} p_{h}) & = \int_{Ω} j \land a_{h}^{'}, & \forall a^{'} \in W^{1} (M) \\ \overset{orthogonality to (discrete) kernel}{\overset{⏞}{m (a_{h}, d_{0} p_{h}^{'}) + \underset{λ = 0}{\underset{⏟}{λ \int_{Ω} p_{h} (x)^{'} d x}}}} & = 0, & \forall p_{h}^{'} \in W^{0} (M) \cap H_{⋆} \land^{0} (d, Ω) \\ \underset{"orthogonality to const."}{\underset{⏟}{\int_{Ω} p_{h} (x) d x}} & = 0. \end{aligned} \end{matrix}$
for all ${\overset{˘}{a}}_{h}^{'} \in W^{1} (M)$ , $p_{h}^{'} \in W^{0} (M)$ .

ordered bases (whitney 1/O-forms)

$\begin{matrix} (13.3.3.34) & [\begin{matrix} A & B^{⊤} & 0 \\ B & O & c \\ 0 & c^{⊤} & 0 \end{matrix}] [\begin{matrix} {\vec{a}}_{h} \\ {\vec{p}}_{h} \\ λ \end{matrix}] = [\begin{matrix} \vec{ψ} \\ 0 \\ 0 \end{matrix}], \end{matrix}$
with ${\vec{p}}_{h}, λ$ dummy unknowns

13.3.3.3-13.3.3.4 Saddle-Point Variational Problems, Discrete LBB-Conditions for a-Based Magnetostatics (60 min)

start: #^13-3-3-28 Is a saddle point variational problem

we define the kernel of $b$ :

N (b) = {u \in M : b (u, q) = 0 \forall q \in Q} \subset M

In Euclidean setting, we get a saddle point LSE:

\begin{matrix} (12.2.2.32) & \begin{matrix} \vec{v} \in R^{m} \\ \vec{π} \in R^{n} \end{matrix} : \begin{aligned} A \vec{v} + B^{⊤} \vec{π} & = \vec{φ}, \\ B \vec{v} & = \vec{γ}, \end{aligned} \end{matrix}

Pasted image 20260514135747.png|500

checking #^LBB conditions:
- LBB2: $N (B^{⊤}) = {0} ⟹ Rank (B^{⊤}) = n ⟹ k = m - n$
- LBB2: $m \geq n$
  - $R (B) = N (B^{⊤})^{⊥} = {0}^{⊥} = R^{n}$
- LBB1: $\exists \vec{π}$

Galerkin discretization:

step 1

\begin{aligned} U_{h} \subset U & , \dim U_{h} = m \\ Q_{h} \subset Q & , \dim Q_{h} = n \end{aligned}

LBB2 in discrete setting $\Rightarrow m \geq n$
Discrete LBB conditions $\Rightarrow quasi-optimality$

Theorem 13.3.3.40. Galerkin discretization error estimate for discrete variational saddle point problems

Let the assumptions of Thm. 12.2.2.23 be satisfied. Also assume discrete counterparts of (LBB1) and (LBB2), namely the discrete LBB-conditions

\begin{aligned} (LBB 1_{h}) & \exists α_{h} > 0 : & | a (v_{h}, v_{h}) | \geq α_{h} ∥ v_{h} ∥_{U}^{2} \forall v_{h} \in \underset{}{\underset{⏟}{N_{h} (b)}}, \\ (LBB 2_{h}) & \exists β_{h} > 0 : & sup_{v_{h} \in U_{h}} \frac{| b (v_{h}, q_{h}) |}{∥ v_{h} ∥_{U}} \geq β_{h} ∥ q_{h} ∥_{Q} \forall q_{h} \in Q_{h} . \end{aligned}

then the discrete variational saddle point problem (12.2.2.49) possesses a unique solution $(v_{h}, p_{h}) \in U_{h} \times Q_{h}$ , which satisfies the Galerkin discretization error estimate

\begin{matrix} (13.3.3.41) & ∥ v - v_{h} ∥_{U} + ∥ p - p_{h} ∥_{Q} \leq C (inf_{u_{h} \in U_{h}} ∥ v - u_{h} ∥_{U} + inf_{q_{h} \in Q_{h}} ∥ p - q_{h} ∥_{Q}), \end{matrix}

where $(v, p) \in U \times Q$ solves (12.2.2.49) and with a constant $C > 0$ depending only on $C_{a}, C_{b}$ and $α_{h}, β_{h}$ }.

Tldr: discrete #^LBB

ellipticity on kernel of $b$
in-sup condiction

-> existence, uniqueness, stability of solution + quasi-optimality

quasi-optimality: $∥ discretization error ∥ \leq ∥ best approx. error | ∥$

step 2 of Galerkin discretization:

ordered bases $\begin{matrix} B_{U} := {b_{h}^{1}, \dots, b_{h}^{N}}, \\ B_{Q} := {β_{h}^{1}, \dots, β_{h}^{M}} \end{matrix}$ of $\begin{matrix} U_{h}, & N := \dim U_{h}, \\ Q_{h}, & M := \dim Q_{h}, \end{matrix}$

saddle point LSE:

$\begin{matrix} (12.3.0.7) & [\begin{matrix} A & B^{T} \\ B & 0 \end{matrix}] [\begin{matrix} \vec{v} \\ \vec{π} \end{matrix}] = [\begin{matrix} \vec{ψ} \\ \vec{φ} \end{matrix}], \end{matrix}$

note, that for our problem, while LBB2 always holds, LBB1 only holds if we assume simple topology (#^simple-topology), which we can prove using

Theorem 13.3.3.58. Poincaré-Friedrichs-type inequality for

d_{1}

Under Ass. 13.3.3.7 and with a constant $C_{1} > 0$ depending only on $Ω$ (and on the definition of the $L^{2}$ -norm) holds true

\begin{matrix} (13.3.3.59) & ∥ v ∥_{L^{2} Λ^{1} (Ω)} \leq C_{1} ∥ d_{1} v ∥_{L^{2} Λ^{2} (Ω)} \forall v \in N (b) . \end{matrix}

Theorem 13.3.3.64. Discrete Poincaré-Friedrichs-type inequality for

d_{1}

Taking for granted Ass. 13.3.3.7 (#^simple-topology) we have

\begin{matrix} (13.3.3.65) & m_{μ}^{*} (d_{1} a_{h}^{'}, d_{1} a_{h}^{'}) \geq α ∥ a_{h}^{'} ∥_{L^{2} Λ^{1} (Ω)}^{2} \forall a_{h}^{'} \in N_{h} (b), \end{matrix}

with a constant $α > 0$ depending only on $Ω, μ_{*}^{\pm}$ , and the shape regularity measure of $M$ (and on the definition of the $L^{2}$ -norm).

since LBB1 & LBB2 holds, we get quasi-optimality (for uniformly shape regular families of meshes)

\begin{matrix} (13.3.3.70) & \begin{aligned} ∥ a - a_{h} ∥ & _{H Λ^{1} (d, Ω)} + ∥ p - p_{h} ∥_{H Λ^{0} (d, Ω)} \\ \leq C (\underset{= O (h_{m})}{\underset{⏟}{inf_{v_{h} \in W^{1} (M)} ∥ a - v_{h} ∥_{H Λ^{1} (d, Ω)}}} + \underset{= O (h_{m})}{\underset{⏟}{inf_{q_{h} \in W^{0} (M)} ∥ p - q_{h} ∥_{H Λ^{0} (d, Ω)}}}), \end{aligned} \end{matrix}

Corollary 13.3.3.71. Asymptotic Whitney finite element error estimate

Provided that $a \in H^{1} Λ^{1} (Ω)$ , $d_{1} a \in H^{1} Λ^{2} (Ω)$ , $p \in H^{2} Λ^{0} (Ω)$ , the asymptotic discretization error estimate

\begin{matrix} (13.3.3.72) & ∥ a - a_{h} ∥_{H Λ^{1} (d, Ω)} + ∥ p - p_{h} ∥_{H Λ^{0} (d, Ω)} = O (h_{M}) for h_{M} \to 0 \end{matrix}

holds true for uniformly shape regular families of simplicial meshes $M$ with meshwidths $h_{M}$ .

algebraic cvg with rate 1!

comparing the two different discretizations:
#todo

Formatting

We're converting image into markdown + latex formulas. For normal text, normal markdown. For boxes, use callouts. The format of callouts is the following:

> [!name] <title>
> content

Note that there is no empty line between the callout header and the content!

orange boxes: convert to !intermediate
lemmas, theorems (yellow): convert to !lemma, !theorem
definitions (red): convert to !definition

For inline latex, the format is the following: $< c o n t e n t >$

For one-line latex, the format is the following:

< c o n t e n t >

note the newlines in between $$ and the content
for latex blocks in callouts, please do not write newlines: if possible, the whole latex block should be on one line

For multiline latex, the format is the following:

\begin{array}{r} < c o n t e n t > \end{array}

Highlighted text can be highlighted in markdown, too: highlight. This is not standard commonmark, but if you see highlighted text in the image, highlight it when converting, too. Same for bold: if a text is emphasised, make it bold.

Note that highlighting can only be used for markdown, not for latex
If something cannot be converted to pure markdown, warn me instead of attempting to use non-standard syntax.
For this task, you may omit any custom formats specified earlier.
Your answer should consist solely of the markdown (and the warning if something cannot be converted). No follow-up suggestions, questions, clarifications.
for colors in latex, use \textcolor{color}{content} instead of \color{color}{content}
whenever possible (i.e. for latex blocks (not inline)), use \tag{...} to write the number