Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Advanced Topics from Scientific Computing
TU Berlin Winter 2022/23
Notebook 12
Creative Commons Lizenzvertrag Jürgen Fuhrmann

html"""
<b> Advanced Topics from Scientific Computing<br>
TU Berlin Winter 2022/23<br>
Notebook 12<br>
<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons Lizenzvertrag" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/80x15.png" /></a>
Jürgen Fuhrmann
</b> <br>
"""

3.5 ms

Partial Differential Equations

md"""
# Partial Differential Equations
"""
117 ms

Notations

md"""
## Notations
"""

5.9 ms

Given: domain $\Omega\subset\mathbb{R}^d$ ($d=1,2,3 \dots$)

  • Dot product: for $\vec{x},\vec{y}\in \mathbb{R}^d$, $\vec{x}\cdot\vec{y} =\sum_{i=1}^d x_i y_i$

  • Bounded domain $\Omega\subset \mathbb R^d$, with piecewise smooth boundary

  • Scalar function $u: \Omega \to \mathbb R$

  • Vector function $\vec{v}=\begin{pmatrix}v_1\\\vdots\\v_d\end{pmatrix} : \Omega \to \mathbb R^d$

  • Partial derivative $\partial_i u = \frac{\partial u}{\partial x_i}$

  • Second partial derivative $\partial_{ij} u = \frac{\partial^2 u}{\partial x_ix_j}$

md"""
Given: domain ``\Omega\subset\mathbb{R}^d`` (``d=1,2,3 \dots``)
- Dot product: for ``\vec{x},\vec{y}\in \mathbb{R}^d``, ``\vec{x}\cdot\vec{y} =\sum_{i=1}^d x_i y_i``
- Bounded domain ``\Omega\subset \mathbb R^d``, with piecewise smooth boundary
- Scalar function ``u: \Omega \to \mathbb R``
- Vector function ``\vec{v}=\begin{pmatrix}v_1\\\vdots\\v_d\end{pmatrix} : \Omega \to \mathbb R^d``
- Partial derivative ``\partial_i u = \frac{\partial u}{\partial x_i}``
- Second partial derivative ``\partial_{ij} u = \frac{\partial^2 u}{\partial x_ix_j}``
"""

608 ms
  • Gradient of scalar function $u: \Omega \to \mathbb R$:

$$\begin{aligned} \mathrm{grad}=\vec\nabla= \begin{pmatrix} \partial_1 \\ \vdots\\ \partial_d \end{pmatrix}: u \mapsto \vec\nabla u= \begin{pmatrix} \partial_1 u\\ \vdots\\ \partial_d u \end{pmatrix} \end{aligned}$$

  • Divergence of vector function $\vec{v}= \Omega \to \mathbb R^d$:

$$ \begin{aligned} \mathrm{div}=\nabla\cdot:\vec{v}= \begin{pmatrix}v_1 \\\vdots\\v_d\end{pmatrix} \mapsto \nabla\cdot\ \vec{v}= {\partial_1 v_1}+\dots +\partial_d v_d \end{aligned}$$

  • Laplace operator of scalar function $u: \Omega \to \mathbb R$

$$\begin{aligned} \mathrm{div}\cdot\mathrm{grad}&= \nabla\cdot\vec\nabla\\ &=\Delta: u\mapsto \Delta u = \partial_{11}u+\dots +\partial_{dd} u \end{aligned}$$

md"""
- _Gradient_ of scalar function ``u: \Omega \to \mathbb R``:
```math
\begin{aligned}
\mathrm{grad}=\vec\nabla= \begin{pmatrix} \partial_1 \\ \vdots\\ \partial_d \end{pmatrix}:
u \mapsto \vec\nabla u= \begin{pmatrix} \partial_1 u\\ \vdots\\ \partial_d u \end{pmatrix}
\end{aligned}
```
- _Divergence_ of vector function ``\vec{v}= \Omega \to \mathbb R^d``:
```math
\begin{aligned}
\mathrm{div}=\nabla\cdot:\vec{v}= \begin{pmatrix}v_1 \\\vdots\\v_d\end{pmatrix} \mapsto \nabla\cdot\ \vec{v}= {\partial_1 v_1}+\dots +\partial_d v_d
\end{aligned}
```
- _Laplace operator_ of scalar function ``u: \Omega \to \mathbb R``
```math
\begin{aligned}
\mathrm{div}\cdot\mathrm{grad}&= \nabla\cdot\vec\nabla\\
&=\Delta: u\mapsto \Delta u = \partial_{11}u+\dots +\partial_{dd} u
\end{aligned}
```
18.2 ms

Lipschitz domains

md"""
## Lipschitz domains
"""
147 μs

Definition: A connected open subset $\Omega\subset \mathbb{R}^d$ is called domain. If $\Omega$ is a bounded set, the domain is called bounded.

definition"""
__Definition__:
A connected open subset ``\Omega\subset \mathbb{R}^d`` is called _domain_.
If ``\Omega`` is a bounded set, the domain is called _bounded_.

"""
8.7 ms

Definition:

  • Let $D \subset \mathbb{R}^n$. A function $f: D\to \mathbb{R}^m$ is called Lipschitz continuous if there exists $c>0$ such that $||f(x)-f(y)|| \leq c ||x-y||$ for any $x,y\in D$

  • A hypersurface in $\mathbb{R}^n$ is a graph if for some $k$ it can be represented on some domain $D\subset \mathbb{R}^{n-1}$ as

$$ x_k=f(x_1,\dots,x_{k-1},x_{k+1},\dots, x_n)$$

  • A domain $\Omega\subset \mathbb{R}^n$ is a Lipschitz domain if for all $x\in \partial\Omega$, there exists a neigborhood of $x$ on $\partial \Omega$ which can be represented as the graph of a Lipschitz continuous function.

definition"""
__Definition__:

- Let ``D \subset \mathbb{R}^n``. A function ``f: D\to \mathbb{R}^m`` is called _Lipschitz continuous_ if
there exists ``c>0`` such that ``||f(x)-f(y)|| \leq c ||x-y||`` for any ``x,y\in D``
- A hypersurface in ``\mathbb{R}^n`` is a _graph_ if for some ``k`` it can be represented on some domain ``D\subset \mathbb{R}^{n-1}`` as
```math
x_k=f(x_1,\dots,x_{k-1},x_{k+1},\dots, x_n)
```
- A domain ``\Omega\subset \mathbb{R}^n`` is a _Lipschitz domain_ if for all ``x\in \partial\Omega``,
there exists a neigborhood of ``x`` on ``\partial \Omega`` which can be represented as the graph of
a Lipschitz continuous function.

"""

2.4 ms

Standard PDE calculus happens in Lipschitz domains

important"""
Standard PDE calculus happens in Lipschitz domains
"""
2.0 ms
let
PyPlot.clf()
fig=PyPlot.figure(10)
PyPlot.subplots(nrows=1, ncols=3)
PyPlot.subplot(1,3,1)
# ax=PyPlot.axes(aspect=1.0)
T=collect(0:0.01:2π)
plot(sin.(T),0.1.*sin.(4*T).+cos.(T),linewidth=3)
PyPlot.axis("off")
gca().get_xaxis().set_visible(false)
gca().axes.get_yaxis().set_visible(false)

PyPlot.subplot(1,3,2)
# ax=PyPlot.axes(aspect=1.0)
T=collect(0:0.75:2π)
push!(T,T[1])
plot(1.5.*sin.(T),cos.(T),linewidth=3)
PyPlot.axis("off")
gca().get_xaxis().set_visible(false)
gca().axes.get_yaxis().set_visible(false)

PyPlot.subplot(1,3,3)
X=collect(-1:0.01:1)
plot(X,sqrt.(abs.(X)),color=:blue,linewidth=3)
xtop=[-1,-1,1,1]
ytop=[1,2,2,1]
1.2 s
  • Boundaries of Lipschitz domains are continuous

  • Polygonal domains are Lipschitz

  • Boundaries of Lipschitz domains have no cusps (e.g. the graph of $y=\sqrt{|x|}$ has a cusp at $x=0$)

statement"""
- Boundaries of Lipschitz domains are continuous
- Polygonal domains are Lipschitz
- Boundaries of Lipschitz domains have no cusps
(e.g. the graph of ``y=\sqrt{|x|}`` has a cusp at ``x=0``)

"""

3.0 ms

Divergence theorem (Gauss' theorem)

md"""
## Divergence theorem (Gauss' theorem)

"""

139 μs

Theorem: Let $\Omega\subset \RR^d$ be a bounded Lipschitz domain and $\vv: \Omega\to \RR^d$ be a continuously differentiable vector function. Let $\vn$ be the outward normal to $\Omega$. Then,

$$ \int_\Omega \nabla\cdot \vv \,d\vx = \int_{\partial\Omega} \vv\cdot\vn\,ds.$$

statement"""
__Theorem__:
Let $\Omega\subset \RR^d$ be a bounded Lipschitz domain and $\vv: \Omega\to \RR^d$ be a
continuously differentiable vector function. Let $\vn$ be the outward normal to $\Omega$.
Then,
```math
\int_\Omega \nabla\cdot \vv \,d\vx = \int_{\partial\Omega} \vv\cdot\vn\,ds.
```
"""
87.9 μs

This is a generalization of the Newton-Leibniz rule of calculus:

Let $d=1$, $\Omega=(a,b)$. Then:

  • $n_a=(-1)$

  • $n_b=(1)$

  • $\nabla\cdot v=v'$

$$\int_\Omega \nabla\cdot \vec{v} \,d\vec{x} = \int_a^b v'(x) \,dx=v(b)-v(a) = v(a)n_a+v(b)n_b$$

md"""

This is a generalization of the Newton-Leibniz rule of calculus:
Let ``d=1``, ``\Omega=(a,b)``. Then:
- ``n_a=(-1)``
- ``n_b=(1)``
- ``\nabla\cdot v=v'``


```math
\int_\Omega \nabla\cdot \vec{v} \,d\vec{x} = \int_a^b v'(x) \,dx=v(b)-v(a) = v(a)n_a+v(b)n_b
```

"""

3.0 ms

Species evolution in a domain Ω

md"""
## Species evolution in a domain Ω
"""

143 μs

Let

  • $\Omega$: domain, $(0,T)$: time evolution interval

  • $u(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}$: time dependent local amount of species (aka species concentration)

  • $f(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}$: species sources/sinks

  • $\vec{j}(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}^d$: vector field of the species flux

md"""
Let
- ``\Omega``: domain, ``(0,T)``: time evolution interval
- ``u(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}``: time dependent _local amount of species_ (aka species concentration)
- ``f(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}``: species _sources/sinks_
- ``\vec{j}(\vec{x},t): \Omega \times [0,T] \to \mathbb{R}^d``: vector field of the species _flux_

"""

467 μs

Representative Elementary Volume (REV)

Let $\omega\subset \Omega$: be a representative elementary volume (REV) Define averages:

  • $J(t)=\int_{\partial\omega}\vec{j}(\vec{x},t) \cdot\vec{n}\;ds$: flux of species trough $\partial \omega$ at moment $t$

  • $U(t)=\int_\omega u(\vec{x},t) \;d\vec{x}$: amount of species in $\omega$ at moment $t$

  • $F(t)=\int_\omega f(\vec{x},t) \;d\vec{x}$: rate of creation/destruction at moment $t$

md"""
### Representative Elementary Volume (REV)
Let ``\omega\subset \Omega``: be a _representative elementary volume (REV)_
Define averages:
- ``J(t)=\int_{\partial\omega}\vec{j}(\vec{x},t) \cdot\vec{n}\;ds``: flux of species trough ``\partial \omega`` at moment ``t``
- ``U(t)=\int_\omega u(\vec{x},t) \;d\vec{x}``: amount of species in ``\omega`` at moment ``t``
- ``F(t)=\int_\omega f(\vec{x},t) \;d\vec{x}``: rate of creation/destruction at moment ``t``
"""
6.5 ms

Species conservation

Let $(t_0,t_1)\subset (0,T)$. The Change of the amount of species in $\omega$ during $(t_0,t_1)$ is proportional to the sum of the amount transported through boundary and the amount created/destroyed

$$U(t_1)-U(t_0)+ \int_{t_0}^{t_1} J(t)\;dt = \int_{t_0}^{t_1} F(t)\;dt $$

Using the definitions of U,F,J, we get

$$\int_{\omega} (u(\vx,t_1)- u(\vx,t_0))\,d\vx + \int_{t_0}^{t_1} \int_{\partial\omega} \vj(\vx,t) \cdot\vn\,ds\,dt=\int_{t_0}^{t_1}\int_{\omega} f(\vx,t)\,ds$$

Gauss' theorem gives

$$\int_{t_0}^{t_1} \int_{\omega} \partial_t u(\vx,t)\,d\vx\,dt + \int_{t_0}^{t_1} \int_{\omega} \nabla\cdot \vj(\vx,t) \,d\vx\,dt =\int_{t_0}^{t_1}\int_{\omega} f(\vx,t)\,ds$$

md"""
### Species conservation

Let ``(t_0,t_1)\subset (0,T)``.
The Change of the amount of species in $\omega$ during $(t_0,t_1)$ is proportional to the sum of the amount transported through boundary and the amount created/destroyed
```math
U(t_1)-U(t_0)+ \int_{t_0}^{t_1} J(t)\;dt = \int_{t_0}^{t_1} F(t)\;dt
```
Using the definitions of U,F,J, we get

```math
\int_{\omega} (u(\vx,t_1)- u(\vx,t_0))\,d\vx + \int_{t_0}^{t_1} \int_{\partial\omega} \vj(\vx,t) \cdot\vn\,ds\,dt=\int_{t_0}^{t_1}\int_{\omega} f(\vx,t)\,ds
```

Gauss' theorem gives
```math
\int_{t_0}^{t_1} \int_{\omega} \partial_t u(\vx,t)\,d\vx\,dt + \int_{t_0}^{t_1} \int_{\omega} \nabla\cdot \vj(\vx,t) \,d\vx\,dt
=\int_{t_0}^{t_1}\int_{\omega} f(\vx,t)\,ds
```
"""

3.8 ms

Continuity equation

The above is true for all $\omega\subset \Omega$, $(t_0,t_1)\subset (0,T)$ $\Rightarrow$

$$\partial_t u(\vx,t) + \nabla\cdot \vj(\vx,t)=f(\vx,t) \quad \text{in}\; \Omega \times [0,T]$$

  • While this sounds obvious, mathematical reasoning about this is more complex

  • Whenever one encounters the divergence operator, chances are that it describes a conservation law for certain species. This physical meaning is very concrete and, if possible should be preserved during the process of discretizing PDEs.

md"""
### Continuity equation

The above is true for all $\omega\subset \Omega$, $(t_0,t_1)\subset (0,T)$ $\Rightarrow$
```math
\partial_t u(\vx,t) + \nabla\cdot \vj(\vx,t)=f(\vx,t) \quad \text{in}\; \Omega \times [0,T]
```
- While this sounds obvious, mathematical reasoning about this is more complex
- Whenever one encounters the divergence operator, chances are that it describes a conservation law for certain species. This physical meaning is very concrete and, if possible should be preserved during the process of discretizing PDEs.
"""
351 μs

Flux expressions

As a rule, species flux is proportional to the negative gradient of the species concentration: $\vec{j}(\vec{x},t) \sim -\vec\nabla u(\vec{x},t)$. This corresponds to the direction of steepest descend.

Therefore we set $\vec{j}=-\delta \vec\nabla u$, where $\delta>0$ can be constant, space/time dependent or even depend on $u$. For simplicity, we assume $\delta$ to be constant, unless stated otherwise.

md"""
## Flux expressions


As a rule, species flux is proportional to the negative gradient of the species concentration: ``\vec{j}(\vec{x},t) \sim -\vec\nabla u(\vec{x},t)``. This corresponds to the direction of steepest descend.

Therefore we set ``\vec{j}=-\delta \vec\nabla u``, where ``\delta>0`` can be constant,
space/time dependent or even depend on ``u``. For simplicity, we assume ``\delta`` to be constant, unless stated otherwise.

"""
271 μs

Heat conduction

  • $u=T$: temperature

  • $\delta=\lambda$: heat conduction coefficient

  • $f$: heat source

  • $\vec{j}=-\lambda\vec\nabla T$: Fourier law

md"""
### Heat conduction

- ``u=T``: temperature
- ``\delta=\lambda``: heat conduction coefficient
- ``f``: heat source
- ``\vec{j}=-\lambda\vec\nabla T``: _Fourier law_

"""
411 μs
Loading...