{grid}
```{grid-item-card}
:class-header: bg-light
Voraussetzungen
^^^
- mathematische Grundlagen von PDEs und der Linienmethode
- Kenntnis des expliziten Euler-Verfahrens und des Crank-Nicolson-Verfahrens
- Partielle Differentialgleichungen mit Randbedingungen
```
```{grid-item-card}
:class-header: bg-light
Lerninhalte
^^^
- Umsetzung der Verfahrensvorschriften zur Umsetzung der Linienmethode
- Einschätzung der Auswirkungen numerischer Ungenauigkeiten
```
(conv_diff_pde)=
{image}
:alt: Aedes Albopictus
:width: 512px
:align: center
Wir wollen simulieren, wie sich ein Insektengift in einer Blutbahn ausbreitet. Zum einen wird sich durch ein Diffusionsprozess das Insektengift mit dem Blut vermischen. Zum anderen wird es mit der Fließgeschwindigkeit des Blutes transportiert (Konvektion). Es ergibt sich ein sogenannter Konvektions-Diffusionsprozess. $\rho(t,x)$ sei die Konzentration des Insektengiftes entlang der Blutbahn.
Die Kontinuitätsgleichung in einer Dimension lautet
$$ 0 = \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial x} J $$mit einem Massestrom $J$. Dieser besteht aus einem Diffusionsterm $-a \frac{\partial \rho}{\partial x}$ der den Konzentrationsgradienten des Insektengiftes ausgleicht, sowie einem Geschwindigkeitsterm $c\rho$,
$$ J = -a \frac{\partial \rho}{\partial x} + c(t) \rho, $$wobei $c(t)$ die Fließgeschwindigkeit des Blutes zum Zeitpunkt $t$ bezeichnet. Als Vereinfachung gehen wir von einer zeitlich konstanten Fließgeschwindigkeit $c(t)=c=\text{const}$ aus. Es sollte nicht allzu schwer sein, das Modell zu einem späteren Zeitpunkt durch eine pulsierende Fließgeschwindigkeit zu erweitern.
Setzen wir den Massestrom in die Kontinuitätsgleichung ein, erhalten wir die Konvektions-Diffusionsgleichung
{math}
:label: convection-diffusion
\frac{\partial \rho}{\partial t} + c \frac{\partial \rho}{\partial x} = a \frac{\partial^2 \rho}{\partial x^2}
Führen Sie eine Diskretisierung des Problems durch (örtlich und zeitlich), um eine Vorschrift zur numerischen Lösung der Gleichung zu ermitteln. Gehen Sie von zeitlich konstanten Dirichlet-Randbedingungen aus. Verwenden Sie für die örtlichen Ableitungen zentrale Differenzenquotienten und für die zeitliche
Gehen Sie wie folgt vor:
Diskretisieren Sie zunächst im Ort, so dass Sie $N$ neue Funktionen $y_i(t) \approx \rho(x_i, t), i=1,...,N$ erhalten. Hierbei bezeichnet $N$ die Anzahl der Diskretisierungspunkte und $x_i$ bezeichnet den $i$-ten Diskretisierungspunkt Ihres Intervals.
Schreiben Sie die approximierte Lösung $y_i(t_j)$ der inneren Punkte $i=2,...,N-1$ in einen Vektor
$$ \mathbf{y}^{(j)} = \begin{bmatrix} y_2(t_j) \\ \vdots \\ y_{N-1}(t_j) \end{bmatrix} = \begin{bmatrix} y_2^{(j)} \\ \vdots \\ y_{N-1}^{(j)} \end{bmatrix} $$Bringen Sie die Ortsableitungen der DGL auf die rechte Seite, ersetzen Sie die partiellen Ortsableitungen durch Differenzenquotienten, so dass Sie für die inneren Punkte Folgendes schreiben können:
$$ \frac{\partial \mathbf{y}(t_j)}{\partial t} = A \cdot \mathbf{y}^{j} + \mathbf{b}. $$Dabei ergeben sich $A \in \mathbb{R}^{N-1 \times N-1}$ und $\mathbf{b} \in \mathbb{R}^{N-1}$ aus den verwendeten Differenzenquotienten sowie den Randwerten $\rho(x_0, t) = \rho_a, \rho(x_N, t) = \rho_b$.
Diskretisieren Sie nun in der Zeit, indem Sie eine Vorschrift
$$ \mathbf{y}^{(j)} = \mathbf{y}^{(j-1)} + \Delta t \phi(\mathbf{y}^{(j)}, \mathbf{y}^{(j-1)}) $$finden. Je nachdem ob Sie das explizite Euler-Verfahren oder das Crank-Nicolson-Verfahren verwenden, erhalten Sie eine anderen Gestalt für $\phi$.
Schreiben Sie ein Matlab-Programm, das die Gleichung {eq}convection-diffusion
für ein beliebiges Zeitinterval $[t_0, t_1]$ und Ortsintervall $[x_0, x_1]$ mit zeitlich konstanten Dirichlet-Randbedingungen $\rho(x_0, t) = \rho_a, \rho(x_N, t) = \rho_b$ löst. Implementieren Sie sowohl das explizite Euler-Verfahren als auch das Crank-Nicolson-Verfahren.
Visualisieren Sie die Lösung mit Hilfe einer Animation. Unten sehen Sie eine Beispiellösung.
{admonition}
:class: warning
Beachten Sie bei der Wahl der Diskretisierung die [Courant-Friedrichs-Lewy-Zahl](https://de.wikipedia.org/wiki/CFL-Zahl).
Prüfen Sie, für welche CFL-Zahl das explizite Euler-Verfahren nicht mehr funktioniert. Stellen Sie durch Ausprobieren auf 2 Nachkommastellen fest, wo der Grenzwert ist. Der Wert liegt vermutlich zwischen 0.1 und 1. Stellen Sie Ihren Code so ein, dass die CFL-Zahl direkt von Ihnen - auf einen klein genugen Wert - festgelegt werden kann.
Erweitern Sie Ihr Programm, so dass die Nutzenden eine zeitabhängige (z.B. pulsierende) Fließgeschwindigkeit $c(t)$ vorgeben können. Passen Sie gegebenenfalls die Schrittweite an, um die CFL-Zahl einzuhalten.
Ihre Lösungen können beispielsweise für das Parameterset
so aussehen (nur jeder 50. Zeitschritt wurde ausgegeben):
{image}
:alt: "blood_animated_cos4pit"
:width: 800px
:align: center
{admonition}
Die Parameter sind so gewählt, dass die physikalischen Effekte in der graphischen Ausgabe wahrnehmbar sind. Realistischer wären Parameter wie $c = 22.88 + 1.58\cdot\sin(2\cdot\pi \,t/\text{s}) \,\frac{\text{cm}}{\text{s}}$ (vgl. {cite}`Mohammadkarim2017`), $a$ in der Größenordnung $10^{-9}$ (vgl. [Wikipedia](https://de.wikipedia.org/wiki/Diffusionskoeffizient)), Strecken in der Größenordnung $\text{cm}$ und ein Puls von 60-80, statt wie hier 120.
Für $a \ll 1$ (hier $a=0.01$) hängen sich unabhängig von der Lösungsmethode hinter dem "Hauptstrom" Oszillationen an:
{image}
:alt: "blood_animated_a001"
:width: 800px
:align: center
Der Ursprung dieser Oszillationen liegt in der Numerik: Die Steilheit der linken Flanke führt dazu, dass die Dichte unter- und später überschätzt wird. Dies lässt sich auch mit einem Differenzenquotienten höherer Ordnung für den Diffusionsanteil nicht lösen. Eine dagegen einfache Abhilfe ist die Verwendung von mehr Ortsschritten (hier $x_n = 500$ statt $100$):
{image}
:alt: "blood_animated_a001_xn500"
:width: 800px
:align: center
{admonition}
Weiterführend können Sie sich über ENO (Essentially Non-Oscillatory) Methoden informieren. Diese sind speziell dafür entwickelt, numerische Oszillationen zu vermeiden. Ein Beispiel für eine ENO Methode finden Sie in dem Kapitel {ref}`traffic_pde`.
{bibliography}
:filter: docname in docnames