Numerische Modellierung der Diffusion

SW6: Finite-Differenzen Methode FDM

Die Finite Differenzen Methode ist neben dem Finite Volumen Verfahren und dem Finite Elementen Verfahren eine der bekanntesten numerischen Diskretisierungsverfahren für partielle Differentialgleichungen.

Gitter-basierte Verfahren

Diese drei Verfahren gehören zu den Gitter-basierten Verfahren, die die unbekannte Funktion $u(x,t)$ auf endlich vielen Gitterpunkten $P_{i,j}=(x_{i},y_{j}),\ i=1,\ldots N_{x},\ j=1,\ldots N_{y}$ approximieren.

Differenzenquotient

In FDM werden die Ableitungen der gesuchten Funktion durch Differenzenquotienten ersetzt.
Im folgenden wird die Hereitung der Differenzenquotienten für eine Funktion einer Veränderlichen demonstriert.

Taylorpolynom

Approximation von ln(x) durch Taylorpolynome der Grade 1, 2, 3 bzw. 10 um die Entwicklungsstelle 1. Die Polynome konvergieren nur im Intervall (0, 2]. Der Konvergenzradius ist also 1. Author: Georg Johann

Die Taylorentwicklung der Funktion $f:\mathbb {R} \to \mathbb {R}$ im Punkt $x+h$ um den Entwicklungspunkt $x$

$f(x+h)=\sum _{k=0}^{\infty }{\frac {f^{(k)}(x)}{k!}}(h)^{k}=f(x)+f'(x).h+f''(x){\frac {h^{2}}{2!}}+f'''(x).{\frac {h^{3}}{3!}}+\ldots$ .

Die Differenzenquotienten ergeben sich als Approximation der Ableitungen von $f$ nach dem Abbruch der Taylorreihe.

Erster Differenzenquotient

Sei Funktion $f$ hinreichend oft differenzierbar. Aus der Taylorformel mit dem Lagrangeschem Restglied

${\begin{array}{rcl}f(x+h)&=&f(x)+hf'(x)+f''(x){\frac {h^{2}}{2}}+f'''(\theta ){\frac {h^{3}}{3!}}\end{array}}$

mit $\theta \in (x,x+h)$. Damit ergibt sich

• für den ersten (vorwärts-) Differenzenquotient
  

$D_{h}^{+}f(x):={\frac {f(x+h)-f(x)}{h}}=f'(x)+f''(x){\frac {h}{2}}+f'''(\theta _{1}){\frac {h^{2}}{6}},\ \ \theta _{1}\in (x,x+h)$,

• für den ersten (rückwärts-) Differenzenquotient mit Anwendung von $-h$ anstatt von $h$ in der obigen Taylorformel
  

$D_{h}^{-}f(x):={\frac {f(x)-f(x-h)}{h}}=f'(x)-f''(x){\frac {h}{2}}+f'''(\theta _{2}){\frac {h^{2}}{6}},\ \ \theta _{2}\in (x-h,x)$,

• für den ersten zentralen Differenzenquotient durch aus der Summe ${\frac {1}{2}}{\big (}D^{+}f(x)+D^{-}f(x){\big )}$,
  

$D_{h}f(x):={\frac {f(x+h)-f(x-h)}{2h}}=f'(x)+(f'''(\theta _{1})+f'''(\theta _{2})){\frac {h^{2}}{12}},\ \ \theta _{1},\theta _{2}\in (x,x\pm h)$

Güte der Approximation

Alle drei erste Differenzenquotienten approximieren die erste Ableitung:
$D_{h}^{\pm }f(x),D_{h}f(x)\approx f'(x)$ für $h\to 0$ falls $f$ hinreichend oft differenzierbar auf einem Intervall $[a,b]$ ist mit $x,\ x\pm h\in [a,b]$.

Lemma 1:

Ist $f\in C^{2}[a,b]$, dann gilt für alle $x\in [h,b-h]$:
   $|D_{h}^{\pm }f(x)-f'(x)|\leq {\frac {h}{2}}C_{1}$ wobei $C_{1}:=\max _{x\in [a,b]}|f''(x)|$.
Ist $f\in C^{3}[a,b]$, dann gilt  für alle $x\in [h,b-h]$:
   $|D_{h}f(x)-f'(x)|\leq {\frac {h^{2}}{6}}C_{2}$ wobei $C_{2}=\max _{x\in [a,b]}|f'''(x)|$. 

Beweis: Aus den Gleichungen für erste Differenzenquotienten $\ \square$.

Zweiter Differenzenquotient

Der zweite Differenzenquotient approximert die zweite Ableitung $D^{2}f(x)\approx f''(x)$ an der Stelle $x$, falls $f$ hinreichend oft differenzierbar ist.

Die Formel für $D^{2}f(x)$ ergibt sich aus der Differenz der ersten Differenzenquotienten $D_{h}^{+}f(x)-D_{h}^{-}f(x)$ mithilfe der Taylorentwicklung bis zur vierten Ableitung:

$D_{h}^{2}f(x):={\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}=f''(x)+{\big (}f^{(IV)}(\theta _{1})+f^{(IV)}(\theta _{2}){\big )}{\frac {h^{2}}{4!}}.$

Bemerkung

Die Formel für den zweiten Differenzenquotient $D_{h}^{2}f(x)$ kann auch durch sukzessive Anwendung des ersten zentralen Differenzenquotienten mit halber Schrittweite $D_{h/2}f(x)$ hergeleitet werden,
$D_{h/2}{\big (}D_{h/2}f(x){\big )}=D_{h/2}{\Big (}{\frac {f(x+h/2)-f(x-h/2)}{h}}{\Big )}={\frac {1}{h}}{\Big (}{\frac {f(x+h)-f(x)}{h}}-{\frac {f(x)-f(x-h)}{h}}{\Big )}\equiv D_{h}^{2}f(x)$.

Güte der Approximation

Lemma 2:

Ist $f\in C^{4}[a,b]$, dann gilt für alle $x\in [h,b-h]$
  $|D_{h}^{2}f(x)-f''(x)|\leq {\frac {h^{2}}{12}}C_{3}$ wobei $C_{3}:=\max _{x\in [a,b]}|f^{(IV)}(x)|$.

Beweis: Aus der Gleichung für zweiten Differenzenquotient. Der Term ${\Big (}f^{(IV)}(\theta _{1})+f^{(IV)}(\theta _{2}){\Big )}{\frac {h^{2}}{24}}$ im Ausdruck für $D_{h}^{2}f(x)$ trägt zum Approximationsfehler bei $\ \square$.

Zusammenfassung

Die Fehler der Differenzenquotienten kann man mithilfe des Landau-Symbols beschreiben:
$D_{h}f(x)=f'(x)+{\cal {O(h^{2})}}$,
$D_{h}^{\pm }f(x)=f'(x)+{\cal {O(h)}}$,
$D_{h}^{2}f(x)=f''(x)+{\cal {O(h^{2})}}$.

Differenzenquotienten höherer Ordnung

Für die Approximation Ableitungen höherer Ordnung siehe Höhere Differenzenquotienten.

Randwertproblem

Randwertproblem für Diffusionsgleichung beschreibt ein Diffusionsproblem auf einem beschränkten Gebiet $D\subset {\mathbb {R} }^{n}$. Hierbei wird das Verhalten der unbekannten Funktion am Rand des Gebietes $\partial D$ durch eine Bedingung festgelegt. Wir behandeln zuerst die stationäre Poissongleichung.

$$-\Delta u(x)=f(x).$$

Dirichlet Randwertproblem

Ist die gesuchte Funktion $u$ am Rand $\partial D$ durch eine gegebene Funktion festgelegt $g$, erfüllt $u$ die sog.
Dirichlet Randbedingung:

$\displaystyle -\Delta u(x)=f(x)\qquad x\in D$

$\displaystyle \quad u(x)=g(x)\quad \ \ x\in \partial D.$

Neumann Randwertproblem

Richtungsableitung (L.Sittinger, S.Schmitz)

Sind die Richtungsableitungen der Funktion $u$ in der Richtung des äußeren Normalenvektor gegeben, erfüllt $u$ die sog.
Neumann Randbedingung:

$-\Delta u(x)=f(x)\quad x\in D$

${\frac {\partial u(x)}{\partial {\vec {n}}}}\equiv \nabla u\cdot {\vec {n}}=g(x)\quad x\in \partial D.$

Numerische Diskretisierung des Randwertproblems

Dirichletproblem auf dem Intervall (1D)

Sei $x\in [a,b]\subset \mathbb {R}$. Gesucht werden Näherungen für die Funktionswerte von $[a,b]$ auf dem äquidistanten Gitter

$\{x_{i}=a+ih,\ i=1,\ldots ,n\in \mathbb {N} \}\subset [a,b],\quad h={\frac {b-a}{n+1}}$

die Werte $u(a),u(b)$ sind durch die Dirichlet-Randbedingungen (Werte am linken und rechten Rand) bekannt.

Diskretisierung durch Differenzenquotient

In diesem Fall handelt es sich um gewöhnliche Differentialgleichung. Die zweiten Ableitungen werden mit Hilfe des zweiten Differenzenquotienten diskretisiert:

$D_{h}^{2}u(x_{i})={\frac {u_{i+1}-2u_{i}+u_{i-1}}{h^{2}}}\approx u''(x_{i})$ für jedes $i=1,\ldots n$.
Hier bezeichnet $u_{i}:=u(x_{i})$.

Lineares Gleichungssystem

Nach der Anwendung der obigen Formel für jedes $i=1,\ldots n$ unter der Beachtung der Randbedingungen $u_{0}=u(a),\ u_{n+1}=u(b)$ in den Formeln für $D_{h}^{2}u(x_{1}),\ D_{h}^{2}u(x_{n})$ für den Vektor den unbekannten Werte ${\vec {u}}_{h}:=(u_{1},\ldots u_{n})^{T}$ ein lineares Gleichungsystem

$$A_{h}{\vec {u}}_{h}={\vec {f}}_{h},$$

mit einer tridiagonaler Systemmatrix $A_{h}$, und dem Vektor der rechten Seite ${\vec {f}}_{h}$,

Systemmatrix und rechte Seite

$A_{h}={\frac {-1}{h^{2}}}.{\begin{pmatrix}-2&1&0&\ldots &0\\1&-2&1&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\0&&\ldots &1&-2\end{pmatrix}},\qquad {\begin{pmatrix}f_{1}+{\frac {u_{0}}{h^{2}}}\\f_{2}\\\vdots \\f_{n-1}\\f_{n}+{\frac {u_{n+1}}{h^{2}}}.\end{pmatrix}}$,

siehe dieses Beispiel.

SW7: Bemerkung: Lösbarkeit des linearen Gleichungssystems

Die Systemmatrix $A_{h},h>0$ ist regulär und damit invertierbar, symmetrisch und positiv definit. Der Lösungsvektor (numerische Lösung) ${\vec {u}}_{h}$ lässt sich durch direkte Methoden für Gleichungssysteme wie Gauß-Eliminierung, oder LU und Choleski-Zerlegung bestimmen.

Außerdem kann man im Falle der streng diagonal-dominanten Matrizen iterative Verfahren wie Jacobi und Gauß-Seidel Verfahren anwenden, siehe Übersicht der Lösungsverfahren für lineare Gleichungssysteme.

SW7: Neumannproblem auf dem Intervall (1D)

Gegeben seien die Randbedingungen für die Ableitungen in den Normalrichtungen $-1,1$:

$$-u'(a)=\alpha ,\ u'(b)=\beta .$$

Um die Ableitungen in den Randpunkten zu approximieren, führen wir neue Unbekannte $u_{0}:=u(a),\ u_{n+1}:=u(b)$ ein.

Approximation erster Ordnung

Verwendet man einseitige Differenzenquotienen für die Approximation der Ableitung an den Rändern,

• den vorwärts -Differenzenquotient $D_{h}^{+}u(a)={\frac {u_{1}-u_{0}}{h}}=u'(a)+{\mathcal {O}}(h)$
• den rückwärts -Differenzenquotient $D_{h}^{-}u(b)={\frac {u_{n+1}-u_{n}}{h}}=u'(b)+{\mathcal {O}}(h)$

erhält man für die neuen Unbekannte $u_{0},\ u_{n+1}$ die 0-te und die ${n+1}$- te Gleichung,
$u_{0}=u_{1}+h\alpha +{\mathcal {O}}(h^{2})$,
$u_{n+1}=u_{n}+h\beta +{\mathcal {O}}(h^{2})$.

Aufgabe Stellen Sie die Systemmatrix ${\tilde {A}}_{h}$ und den Vektor der rechten Seite für diese Approximation der Neumann Randbedingung auf. Ist diese Matrix regulär?

Die Genauigkeit des Approximationsschemas gegeben durch das obere lineare Gleichungssystems ist durch diese Approximation der Ableitungen am Rand des Intervalls reduziert auf die erste Ordnung, der vernachlässigte Fehler ist insgesamt von der Größenordnung ${\mathcal {O}}(h)$.

Approximation zweiter Ordnung

Um die Approximationsgüte der zweiten Differenzquotienten für die Approximation von $u''(x)$ im Inneren des Intervalls ${\mathcal {O}}(h^{2})$, auch für die Randbedingungen zu erhalten, verwendet man den ersten zentralen Differenzenquotient:

$D_{h}u(a)={\frac {u_{1}-u_{-1}}{2h}}=u'(a)+{\mathcal {O}}(h^{2})$

$D_{h}u(b)={\frac {u_{n+2}-u_{n}}{2h}}=u'(b)+{\mathcal {O}}(h^{2})$.

Damit man das lineare Gleichungssystem nicht wieder um neue Gleichungen für $u_{-1},\ u_{n+2}$ ergänzen muss, eliminiert man diese Variablen mithilfe der Randbedingungen. Man erhält aus obigen Gleichungen für $D_{h}$ folgende Gleichungen für $u_{-1},\ u_{n+2}$:

$u_{-1}=u(a-h)=u_{1}-2hu'(a)+{\mathcal {O}}(h^{3})=u_{1}+2h\alpha +{\mathcal {O}}(h^{3})$,

$u_{n+2}=u(b+h)=u_{n}+2hu'(b)+{\mathcal {O}}(h^{3})=u_{n}+2h\beta +{\mathcal {O}}(h^{3})$.

Ersetzt man $u_{-1},\ u_{n+2}$ in den Gleichungen für den zweiten Differenzenquotient $D_{h}^{2}u(x_{0}),\ D_{h}^{2}u(x_{n+1})$ aus obigen Gleichungen, erhält man anstelle der ersten und letzter Gleichung des linearen Systems:
${\frac {2u_{1}-2u_{0}}{h^{2}}}=-2\alpha /h+{\mathcal {O}}(h),\quad {\frac {2u_{n}-2u_{n+1}}{h^{2}}}=-2\beta /h+{\mathcal {O}}(h).$

Basierend auf der Approximation der ersten Ableitungen an den Rändern durch ersten zentralen Differenzenquotient erhält man auf diese Weise insgesamt die Approximation zweiter Ordnung.

Bemerkung: Dass die ersten zwei Gleichungen die Approximation mit einem Fehler insgesamt von der Größenordnung ${\mathcal {O}}(h^{2})$ darstellen, sieht man nach der Multiplikation der obigen zwei Gleichungen mit ${h}/{2}$.

Aufgabe: Stellen Sie die Systemmatrix ${\tilde {A}}_{h}$ und den Vektor der rechten Seite für diese Approximation der Neumann Randbedingung auf.

Bemerkung zur Lösbarkeit des linearen Systems für Neumann-Problem

Die Systemmatrix ${\tilde {A}}_{h}\in \mathbb {R} ^{(n+2)\times (n+2)}$ ergänzt um die neuen Zeilen/Spalten für die Approximation der Neumann Randbedingung ist singulär, denn es gilt ${\tilde {A}}_{h}{\bf {1}}={\bf {0}}$ wobei ${\bf {1}}=(1,1,\ldots ,1),\ {\bf {0}}=(0,0,\ldots ,0).$.

Man kann zeigen (Gauß Eliminierung) $Rang({\tilde {A}}_{h})=n+1,\$.

Damit ist die Lösung dieses linearen System nicht eindeutig.
Dies korrespondiert auch mit der Theorie dieser Differentalgleichung und dem Erhaltungsgesetz:

Nicht-Eindeutigkeit und die notwendige Bedingung für Neumann-Problem

• Wenn die Funktion $u(x)$ eine Lösung des Neumannproblems ist, dann ist diese Lösung nicht eindeutig, auch $u(x)+c$ für jede beliebige Konstante ist eine Lösung.

Die Ableitungen am Rand bestimmen die Steigung der gesuchten Funktion, aber nicht deren Wert. Durch das Festlegen eines Wertes der gesuchten Funktion, z.B. $u_{0}=u(a)$ wird die Lösung eindeutig.

• existiert eine Lösung des Neumann Randwertproblems, dann muss diese nach der Integration erfüllen:
  

$\int _{a}^{b}u''(x)dx=u'(b)-u'(a)=\alpha +\beta =\int _{a}^{b}f(x)dx$
bzw. (mehrdimensional, Anwendung von Satz von Stokes, siehe ):
$\int _{D}\Delta u(x)dx=\int _{\partial D}\nabla u\cdot {\vec {n}}dSx=\int _{\partial D}gdSx=\int _{d}f(x)dx$,
also müssen die Flüsse über den Rand $\alpha$ und $\beta$ , bzw. $g$ mit dem Quellström $f$ im folgenden Sinne korrespondieren:
Im stationären Zustand balanciert der Materialfluss durch die Ränder die Materialquelle im Inneren.

Dirichletproblem in 2 Raumdimensionen

Sei u die gesuchte Funktion definiert auf einem Rechteck D:
$u:\mathbb {R} ^{2}\to \mathbb {R} ,\ {\bf {x}}=(x,y)\in D:=[a,b]\times [c,d]$. Wir betrachten das Dirichlet-Randwertproblem mit homogenen Randbedingungen:

$${\begin{aligned}-\Delta u(x)&=f(x)\ \ {\mbox{in }}\ D\subset \mathbb {R} ^{2},\\u&=0\ \ \ \ \ \ \ \ {\mbox{auf}}\ \ \partial D\end{aligned}}$$

Punktegitter

Der einfachste Fall ist die Verwendung eines äquidistanten Gitters

$\{(x_{i},y_{j}),\ i=1,\ldots n,\,j=1,\ldots m\}$
mit konstanter Gittergröße h:
$h=x_{i}-x_{i-1}=y_{j}-y_{j-1}={\frac {b-a}{n+1}}={\frac {d-c}{m+1}},\ i=1,\ldots n,\,j=1,\ldots m$

Man bezeichne die Approximationen der Lösung in den Gitterpunkten $(x_{i},y_{j}):$

$u_{i,j}:\approx u(x_{i},y_{j}),\ i=1,\ldots n,\ j=1,\ldots m.$

Differenzenquotienten in 2D

Die zweiten Differenzenquotienten $D_{h,x}^{2}u(x_{i},y_{j}),\ D_{h,y}^{2}u(x_{i},y_{j})$ approximieren die zweiten partiellen Ableitungen jeweils in den Richtungen x und y:

${\frac {\partial ^{2}u(x_{i},y_{j})}{\partial x^{2}}}\approx D_{h,x}^{2}u(x_{i},y_{j}):={\frac {u(x_{i+1},y_{j})-2u(x_{i},y_{j})+u(x_{i-1},y_{j})}{h^{2}}}={\frac {u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}}$,
${\frac {\partial ^{2}u(x_{i},y_{j})}{\partial y^{2}}}\approx D_{h,y}^{2}u(x_{i},y_{j}):={\frac {u(x_{i},y_{j+1})-2u(x_{i},y_{j})+u(x_{i},y_{j-1})}{h^{2}}}={\frac {u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}}$
in jedem Gitterpunkt $(x_{i},y_{j})$.

Approximationsschema

Wir erhalten schließlich das Approximationsschema:

$$\Delta u(x_{i},y_{j})\approx {\frac {u_{i+1,j}+u_{i,j+1}-4u_{i,j}+u_{i-1,j}+u_{i,j-1}}{h^{2}}},\qquad i=1,\ldots n,\ j=1,\ldots m.$$

Randwerte

In den Randgitterpunkten $(x_{0},\cdot ),\ (x_{n+1},\cdot ),\ (\cdot ,y_{0}),\ (\cdot ,y_{m+1})$ werden die Werte von $u$ nach der Dirichlet-Randbedingung folgend ersetzt:

$u(x_{0},y_{j})=u(x_{n+1},y_{j})=0,\ j=1,\ldots m$,

$u(x_{i},y_{0})=u(x_{i},y_{m+1})=0,\ i=1,\ldots n$.

Assemblierung

Die unbekannten Approximationen der Lösung in den Gitterpunkten$(x_{i},y_{j}),\ u_{i,j}\approx u(x_{i},y_{j}),\ i=1,\ldots n,\ j=1,\ldots m,$ werden in einen langen Vektor eingeordnet, dabei werden diese Werte in einer bestimmter Reihenfolge in den Vektor aufgestellt.

Beispiel:
Zuerst sammelt man die Unbekannten in der ersten Zeile des Gitters ($y_{1}$ist fest) in den Gitterpunkten von links nach rechts, $x_{1}\to x_{n}$, dann in der zweite Zeile u.s.w.

Die Approximationswerte $u_{i,j}$ sind in dem unbekannten Vektor ${\vec {u}}_{h}$ dann wie folgt geordnet,

$${\vec {u}}_{h}=(u_{1,1},u_{2,1},\ldots u_{n,1}|\ldots |u_{1,m},\ldots u_{n,m})^{T}.$$

In der selben Reihenfolge wird auch der Vektor der rechten Seite aufgestellt,

$${\vec {f}}_{h}=(f_{1,1},f_{2,1},\ldots f_{n,1}|\ldots |f_{1,m},\ldots f_{n,m})^{T},\ f_{i,j}:=f(x_{i},y_{j}).$$

Lineares Gleichungssystem

Nach der oben beschriebener numerischen Diskretisierung und Assemblierung des unbekannten Vektors ergibt sich ein lineares Gleichungssystem $A_{h}{\vec {u}}_{h}={\vec {f}}_{h}$ mit Blocktridiagonaler Matrix $A_{h}\in {\mathbb {R} }^{n.m\times n.m}$

$A_{h}={\frac {-1}{h^{2}}}{\begin{pmatrix}B&I&&\ldots &0\\I&B&I&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\0&&\ldots &I&B\end{pmatrix}}$ mit Diagonalblock $B={\begin{pmatrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{pmatrix}},\quad B\in {\mathbb {R} }^{n\times n}$
und der Einheitsmatrix $I=E_{n}\in {\mathbb {R} }^{n\times n}$ auf der Nebendiagonale.

Inhomogene Randbedingungen

Sind die Randwerte der gesuchten Funktion unterschiedlich von Null, $u|_{\partial D}=g\neq 0$, tragen die bekannte Randwerte der gesuchten Funktion zum Vektor der rechten Seite ${\vec {f}}_{h}$ bei.

Aufgabe: Seien Funktionen $C(y),\ D(y),\ A(x),\ B(x)$ gegeben und die Randwerte der gesuchten Funktion am $\partial D,\ D=[a,b]\times [c,d]$ wie folgt festgelegt:

• Linker/rechter Rand

$u(x_{0},y)=u(a,y)=C(y),$
$u(x_{n+1},y)=u(b,y)=D(y),$

• unterer/oberer Rand

$u(x,y_{0})=u(x,c)=A(x),$
$u(x,y_{m+1})=u(x,d)=B(x).$
Wie verändert sich der Vektor der rechten Seite ${\vec {f}}_{h}$ für diese Randbedingungen ?

Neumann-Problem in 2 Raumdimensionen

Im Neumann-Problem mit äquidistanten Punktegitter wird die Anwendung der zweidimensionaler Differenzenquotienten in dasselbe Stern-Approximationsschema resultieren wie bei Dirichlet-Problem, lediglich wird die Behandlung den Randbedingungen zu veränderter Struktur einiger Blöcke der tridiagonaler Matrix nach der Aufstellung des linearen Gleichungssystem führen.

Neumann-Randbedingungen auf dem Rechteck

Seien die Ableitungen der unbekannten Funktion in den Randpunkten des Rechtecks gegeben:

• linker Rand: ${\frac {\partial u}{\partial {\vec {n}}}}(x_{0},y_{j})=\alpha ,\ \quad {\vec {n}}=(-1,0)^{T},$
• rechter Rand: ${\frac {\partial u}{\partial {\vec {n}}}}(x_{n+1},y_{j})=\beta ,\ \quad {\vec {n}}=(1,0)^{T},\ j=1,\ldots m,$
• unterer Rand: ${\frac {\partial u}{\partial {\vec {n}}}}(x_{i},y_{0})=\gamma ,\ \quad {\vec {n}}=(0,-1)^{T},$
• oberer Rand: ${\frac {\partial u}{\partial {\vec {n}}}}(x_{i},y_{m+1})=\delta ,\ \quad {\vec {n}}=(0,1)^{T},\ i=1,\ldots n,$

wobei ${\frac {\partial u}{\partial {\vec {n}}}}\equiv \nabla u\cdot {\vec {n}}={\frac {\partial u}{\partial x}}n_{1}+{\frac {\partial u}{\partial y}}n_{2}.$

Speziell ${\frac {\partial u}{\partial {\vec {n}}}}=\pm {\frac {\partial u}{\partial x}}$ oder $\pm {\frac {\partial u}{\partial y}}$ am Rand eines Rechtecks.

Approximation erster Ordnung

Bei der Anwendung der einseitigen Differenzenquotienten, analog wie 1-dimensionalem Fall muss der Vektor den Unbekannten ${\vec {u}}_{h}$ um die Randwerte ${u}_{0,j},\ {u}_{n+1,j},\ {u}_{i,0},\ {u}_{i,m+1}$ erweitert werden, diese werden zur Bildung der Differenzenquotienten verwendet:

• linker Rand: den vorwärts -Differenzenquotient $D_{h,x}^{+}u(x_{0},y_{j})={\frac {u_{1,j}-u_{0,j}}{h}}\approx {\frac {\partial u}{\partial x}}(x_{0},y_{j}))\equiv -\alpha$
• rechter Rand: den rückwärts -Differenzenquotient $D_{h,x}^{-}u(x_{n+1},y_{j})={\frac {u_{n+1,j}-u_{n,j}}{h}}\approx {\frac {\partial u}{\partial x}}(x_{n+1},y_{j})\equiv \beta ,\ j=1,\ldots m$
• unterer Rand: den vorwärts -Differenzenquotient $D_{h,y}^{+}u(x_{i},y_{0})={\frac {u_{i,1}-u_{i,0}}{h}}\approx {\frac {\partial u}{\partial y}}(x_{i},y_{0}))\equiv -\gamma$
• oberer Rand: den rückwärts -Differenzenquotient $D_{h,y}^{-}u(x_{i},y_{m+1})={\frac {u_{i,m+1}-u_{i,m}}{h}}\approx {\frac {\partial u}{\partial y}}(x_{i},y_{m+1}))\equiv \delta ,\ i=1,\ldots n.$
  

(Hier wurde direkt die obige Neumann-Randbedingung für ein Rechteckgebiet $\ D=[a,b]\times [c,d]$ eingesetzt.)

Lineares Gleichungssystem

Nach der Einordnung der unbekannter Werte in den Lösungsvektor wie obenbeschrieben (siehe Assemblierung ) ergibt sich wie zuvor eine erweiterte Blocktridiagonale Systemmatrix ${\tilde {A}}_{h}\in {\mathbb {R} }^{(n+2)(m+2)\times (n+2)(m+2)}$:

• erweitert um die 0-te und m+1-te Blockzeile für die Approximation der Ableitung nach y in den Gitterpunkten am oberem und unterem Rand,
• die Blöcke der Matrix: ${B,\ I},\ 0$ sind um die 0-te und n+1-te Zeile für die Approximation der Ableitung nach x in den Gitterpunkten am linkem und rechtem Rand erweitert, der Block $B$ ist mit der Approximationsmatrix ${\tilde {A}}_{h}$ des Neumann-Problem in einer Raumdimension identisch.

Aufgabe: Stellen Sie die Systemmatrix ${\tilde {A}}_{h}$ und die rechte Seite des linearen System des Approximationschemas für den oben definerten Neumann-Randwertproblem auf dem Rechteck D auf.

SW 8: Konvergenz und Stabilität

Wir untersuchen, unter welcher Vorausstetzungen und mit welcher Approximationsgüte
die numerische Lösung für das Randwertproblem für Poissongleichung, $u_{i}$ (in 1D), ggf. $u_{i,j}$ (auf dem Recheteck in 2D), repräsentiert durch den Vektor ${\vec {u}}_{h}$

zu der exakten Lösung in den Gitterpunkten $u(x_{i})$, ggf. $u(x_{i},y_{j})$ konvergiert.

Im folgenden bezeichnen wir mit ${\vec {u}}$ den Vektor der Restriktion der Funktion der exakten Lösung an die Gitterpunkte $x_{i}\ i=1,\ldots N$, ggf $(x_{i},y_{j})\ i=1,\ldots N,\ j=1,\ldots M$ (nach der Assemblierung).

Konvergenz der numerischen Lösung

Wir untersuchen den Unterschied der exakten und der numerischen Lösung in einer geeigneter Vektornorm auf ${\mathbb {R} }^{N}$, wobei im Folgenden N die Gesamtanzahl der Gitterpunkte bezeichnet (in vereinfachter Indexierung der Gitterpunkte).

$\displaystyle \|{\vec {u}}-{\vec {u}}_{h}\|,\quad {\vec {u}}=(u(x_{1}),\ldots u(x_{N}))^{T},\ {\vec {u}}_{h}=(u_{1},\ldots u_{N})^{T},\ u_{i}\approx u(x_{i}).$

Im folgenden werden wir zeigen, dass die Konvergenz $\|{\vec {u}}-{\vec {u}}_{h}\|\to 0$ für Gittergröße $h\to 0$ von

• den Eigenschaften der Systemmatrix $A_{h}$
• der Approximationsgüte des Approximationschema (der Differenzenquotienten)

abhängig ist.

Voraussetzungen der Konvergenz

Sei Matrix $A_{h}$ invertierbar und ${\vec {u}}_{h}$ die Lösung des linearen Systems $A_{h}{\vec {u}}_{h}={\vec {f}}_{h}=:{\vec {f}}=(f(x_{1}),\ldots ,(x_{N}))^{T}$.

Unter Anwendung der Verträglichkeit der Matrixnorm folgt

$\displaystyle \|{\vec {u}}-{\vec {u}}_{h}\|=\|A_{h}^{-1}A_{h}{\vec {u}}-A_{h}^{-1}{\vec {f}}\|=\|A_{h}^{-1}(A_{h}{\vec {u}}-{\vec {f}})\|\leq \||A_{h}^{-1}\||.\|A_{h}{\vec {u}}-{\vec {f}}\|,$

wobei $\||\cdot \||$ die Vektornorm-induzierte (natürliche) Matrixnorm bezeichnet.

Folglich ist die Konvergenz $\|{\vec {u}}-{\vec {u}}_{h}\|\to 0$ gegeben falls

1. die Norm von $A_{h}^{-1}$ gleichmäßig (bez. $h$) beschränkt ist: $\||A_{h}^{-1}\||\leq r,\ r\in {\mathbb {R} }$ - (Stabilität)
2. $\|A_{h}{\vec {u}}-{\vec {f}}\|\to 0$, für $h\to 0$, d.h., das Approximationsschema $A_{h}{\vec {u}}_{h}={\vec {f}}$ ist mit dem Originalproblem, $-\Delta {\vec {u}}={\vec {f}}$ veträglich (Konsistenz).

Der zweite Punkt beschreibt, dass Restriktion der exakten Lösung auf die Gitterpunkte das Approximationschema annäherungsweise erfüllt.

Definition der Konsistenz

Wir setzen voraus, dass die Funktion der rechten Seite $f$ stetig ist.

Definition: (Konsistenzordnung)
Das Finite Differenzenverfahren für die Poisson Randwertaufgabe hat bezüglich einer Vektornorm in ${\mathbb {R} }^{N}$ die Konsistenzordnung $p$ wenn für jede hinreichend oft stetig differenzierbare Lösung ${\vec {u}}$ des Originalproblems ein $C\in {\mathbb {R} }_{+}$ existiert

$\displaystyle \|A_{h}{\vec {u}}-{\vec {f}}\|=\|A_{h}{\vec {u}}+\Delta {\vec {u}}\|\leq C.h^{p}$

Konsistenz elliptischer Randwertprobleme

Den Begriff der Konsistenz, Stabilität und schließlich der Beweis der Konvergenz für das FDM- Verfahren kann man auf elliptische Randwertprobleme mit elliptischen Differenzialoperator $L:u\to L(u)$ für die Funktion $u$ ausbreiten.

Elliptisches Randwertproblem

Vorausgesetzt $f,c:D\to {\mathbb {R} },\ b:D\to {\mathbb {R} }^{n},$ stetige Funktionen auf ${\mathbb {R} }^{n}$ sind gegeben und $c\geq 0$.

$\displaystyle L(u):=-\Delta u(x)+b(x)\cdot \nabla u(x)+c(x)u(x)=f(x),\ \ x\in D\subset {\mathbb {R} }^{n},$

$u(x)=0,\ \ x\in \partial D$

Approximationschema des elliptischen Operators

Das entsprechende Approximationsschema $L_{h}(u)$ zum Differenzialoperator $L(u)$ enthält im Fall vom Rechteckgebiet $D\in {\mathbb {R} }^{2}$ das Stern-Approximationsschema für das Laplace Operator anhand der zweiten Differenzenquotienten $D_{h}^{2}$ (Matrix $A_{h}$) und ein weiteres Beitrag durch die Approximation des Gradienten $\nabla u=(\partial _{x_{1}}u,\ldots ,\partial _{x_{n}}u)^{T}$ mit einseitigen, oder zentralen Differenzenquotienten $D_{h}^{\pm },\ D_{h}$, siehe erste Differenzenquotient.
Die Konsistenzdefinition für den elliptischen Operator lautet:

$\displaystyle \|L_{h}{\vec {u}}-{\vec {f}}\|=\|L_{h}{\vec {u}}-L({\vec {u}})\|\leq C.h^{p}.$

Satz: Konsistenz des FDM Verfahren

Sei $u\in C^{4}(D)$ die exakte Lösung des obigen elliptischen Randwertproblems $L(u)=f$ mit homogenen Randbedingungen. 

Dann hat das Finite-Differenzen-Verfahren mit dem Approximationsschema $L_{h}$
i) die Konsistenzordnung $p=1$ unter der Anwendung von einseitigen Differenzenquotienten $D_{h}^{\pm }$
ii) die Konsistenzordnung $p=2$ unter der Anwendung von zentralen Differenzenquotienten $D_{h}$
für den Term $\nabla u$. $\quad \diamond$

Beweis: basiert auf der Approximationsgüte der ersten und der zweiten Differenzenquotienten (Lemma 1 und 2). Für das Approximationsschema $A_{h}$ ergibt sich daraus $A_{h}{\vec {u}}=-\Delta {\vec {u}}+{\cal {O}}(h^{2})$, für das Approximationsschema $L_{h}$ dagegen entweder $L_{h}{\vec {u}}=L({\vec {u}})+{\cal {O}}(h^{2})+{\cal {O}}(h)=L({\vec {u}})+{\cal {O}}(h)$ oder $L_{h}{\vec {u}}=L({\vec {u}})+{\cal {O}}(h^{2})$ $\square$.

Folgerungen

1. Im Fall $b=c=0,\ D\in {\mathbb {R} }^{2}-$ Rechteckgebiet, erhalten wir die Randwertaufgabe für Poissongleichung mit homogenen Dirichletrandbedingungen. In diesem Fall ist das Approximationsschema $L_{h}=A_{h}$ von der Konsistenzordnung $p=2$.
   
2. Im Fall dass $D\in {\mathbb {R} }$ ein geschlossenes Intervall ist, erhalten wir aus Lemma 2 in der Maximum-Vektornorm $\|\cdot \|_{\infty }$ für ${\vec {u}},\ {\vec {f}}\in {\mathbb {R} }^{N}$ konkret

$\|A_{h}{\vec {u}}-{\vec {f}}\|_{\infty }=\|A_{h}{\vec {u}}+{\vec {u}}''\|_{\infty }=\max _{i}\left|-D_{h}^{2}u(x_{i})+u''(x_{i})\right|\leq {\frac {h^{2}}{12}}\max _{x\in D}|{u}^{(IV)}(x)|={\frac {h^{2}}{12}}\max _{x\in D}|{f}''(x)|.$

Damit ist die Konstante $C$ aus der Definition der Konsistenz im eindimensionalen Fall $C={\frac {1}{12}}\max _{x\in D}|{f}''(x)|.$

Definition der Stabilität

Definition: (Stabilität)
Das Finite-Differenzen-Verfahren für das elliptische Randwertproblem $\displaystyle L(u(x))=f(x),\ \ x\in D\subset {\mathbb {R} }^{n},$ $u(x)=0,\ \ x\in \partial D$ heißt stabil (bezüglich einer Vektornorm $\|\cdot \|$),
wenn $C,h_{0}$ existieren, $C,h_{0}\in {\mathbb {R} }_{+}$, sodass die Matrix (Approximationsschema) $L_{h}$ invertierbar ist und es gilt $\||L_{h}^{-1}\||\leq C$ für $0<h<h_{0}$.

Die Untersuchung der Stabilität anhand der Konstruktion der inversen Matrix $L_{h}^{-1}$ ist aufwändig. Deswegen wird die Stabilität anhand anderer Kriterien untersucht, der Monotonie Eigenschaften der Matrix $L_{h}$.

Bemerkung:

Die Stabilitätsbedingung garantiert auch die Beschränkheit der numerischen Lösung als Lösung des linearen Systems $L_{h}{\vec {u}}_{h}={\vec {f}}$, denn $\|{\vec {u}}_{h}\|=\|L_{h}^{-1}{\vec {f}}\|\leq \||L_{h}^{-1}\||.\|{\vec {f}}\|$ in einer Vektornorm. Hier $\||L_{h}^{-1}\||$ ist eine natürliche Matrixnorm.

Monotonie Eigenschaften

Definiton (M-Matrix)

Eine reguläre Matrix $A\in {\mathbb {R} }^{N\times N}$ mit $a_{ij}\leq 0$ für $i\neq j$

und $(A^{-1})_{ij}\geq 0$, kurz $A^{-1}\geq 0$, heißt M-Matrix, oder Monotonie-Matrix.

Folgerung: Monotonie Eigenschaft

Für eine M-Matrix $A\in {\mathbb {R} }^{N\times N}$ und Vektoren ${\vec {x}},{\vec {y}}\in {\mathbb {R} }^{N}$gilt:

${\vec {x}}\leq {\vec {y}}\implies A^{-1}{\vec {x}}\leq A^{-1}{\vec {y}}$

Beweis: Ist ${x}_{j}\leq {y}_{j}\ \ \forall j=1,\ldots N\ \implies \ \sum _{j=1}^{N}(A^{-1})_{ij}x_{j}\leq \sum _{j=1}^{N}(A^{-1})_{ij}y_{j},\$ da die Elemente $(A^{-1})_{ij}\geq 0\ \ \square .$

SW9 Kriterien für M-Matrix

In diesem Abschnitt werden die Kriterien für Tridiagonalmatrizen wie die Systemmatrix $A_{h}$ des Laplace-Approximationsschemas in 1D und die Matrix $L_{h}$ des elliptischen Randwertproblems in einer Raumdimension $x\in {\mathbb {R} }$ untersucht.

Satz: Monotonie tridiagonaler Matrizen

Jede irreduzibel diagonaldominante  Tridiagonalmatrix mit positiven Diagonalelementen und negativen Nebendiagonalelementen ist eine M-Matrix. $\quad \diamond$

Beweis: siehe $\square$

Bemerkung:
Tridiagonale Matrizen sind irreduzibel falls alle Sub- und Superdiagonalelemente ungleich Null sind.

Folgerung 1

Die Matrix   $A_{h}$ des Laplace-Approximationsschemas in 1D für das Dirichlet Randwertproblem is eine M-Matrix.

Beweis:
$A_{h}$ ist tridiagonal, hat alle Sub- und Superdiagonalelementen ungleich Null und für erste und letzte Zeile gilt $|a_{i,i}|>|a_{i,i\pm 1}|$. Damit ist $A_{h}$ eine M-Matrix $\square$.

Folgerung 2

Die Matrix   $L_{h}$ des elliptischen Randwertproblems in 1D unter der Verwendung des  ersten zentralen Differenzenquotienten $D_{h}$   für $u'(x_{i}),\ i=1,\ldots N$ ist  für 
$0<h<h_{0}\equiv {\frac {2}{\max _{x\in D}|b(x)|}}$   eine M-Matrix.

Beweis:

$L_{h}={\frac {1}{h^{2}}}{\begin{pmatrix}d_{1}&s_{1}&0&\ldots &0\\r_{2}&d_{2}&s_{2}&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\0&&\ldots &r_{N}&d_{N}\end{pmatrix}}$ ist tridiagonal, mit

${\begin{matrix}d_{i}=2+h^{2}c(x_{i}),\\r_{i}=-1-{\frac {h}{2}}b(x_{i}),\\s_{i}=-1+{\frac {h}{2}}b(x_{i}).\end{matrix}}$
Nach der Voraussetzung ist $\left|{\frac {h}{2}}b(x_{i})\right|<1,\ i=1,\ldots N$ und somit sind $r_{i},\ s_{i}\neq 0$ und negativ.
$\implies L_{h}$ ist irreduzibel.
$L_{h}$ ist auch schwach-diagonaldominant:
${\frac {1}{h^{2}}}(|r_{i}|+|s_{i}|)={\frac {2}{h^{2}}}\leq {\frac {d_{i}}{h^{2}}}$ für $c(x)\geq 0,\ i=2,\ldots N-1,$ wobei für $i=1,N$ die strikte Ungleichung gilt $\ \square$.

Monotonie der Systemmatrix und die Stabilität

Hier befassen wir ist mit der Frage wie aus der Monotonie Eigenschaft die Stabilität des Approximationsschemas folgt.

Die Beschränkheit der Matrix $L_{h}^{-1}$ untersuchen wir in der Maximumnorm $\||\cdot \||_{\infty }$ (Zeilensummennorm), die von der Maximum-Vektornom $\|\cdot \|_{\infty }$ induziert ist, siehe .

Hilfsproblem

Wir untersuchen folgendes Randwertproblem: gegeben sei $b\in C^{2}[0,1]$:

$-w''(x)+b(x)w'(x)=1,$

$\ w(0)=w(1)=0.$

Sei die Lösungsfunktion $w\in C^{4}[0,1]$.

Lemma 3:

Die Lösung $w$ des obigen Randwertproblem ist $w\geq 0\ \forall x\in [0,1]$.

Beweis:
Angenommen $w\neq 0$ (Fall $w=0$ können wir ausschließen).

Wäre $w<0$ mindesten in einem Punkt, hätte $w$ ein lokales Minimum in einem $x_{0}\in (0,1)$. d.h.,
$w'(x_{0})=0$ und $w''(x_{0})\geq 0$.

Damit erhalten wir aus der obigen Differenzialgleichung $-w''(x_{0})+b(x_{0})w'(x_{0})=-w''(x_{0})=1\implies w''(x_{0})\leq 0$, was im Widerspruch zu der Bedingung des lokalen Minima steht. $\square$

Beschränkheit der Matrix $L_{h}^{-1}$

Sei zusätzlich eine stetige nichtnegative Funktion $c\in C[0,1],\ c\geq 0$ gegeben.

Für den elliptischen Operator mit Funktion $b$ aus obigen Hilfsproblem gilt nach Lemma 3:

$L(w)=-w''+bw'+cw\geq 1,$

Aus dem Satz über die Konsistenz des Approximationsschema für elliptische Operatoren erhalten wir für das Approximationsschema $L_{h}$ zum Operator $L$ (unter Verwendung der zentralen Differenzenquotioenten) und für den Vektor ${\vec {w}}$ der Restriktion der Lösung ${w}$ auf die Gitterpunkte,

$\|L_{h}{\vec {w}}-L({\vec {w}})\|_{\infty }\leq C.h^{2}.$

Aus dieser Abschätzung in Maximumnorm folgt, dass auch gilt $|(L_{h}{\vec {w}}-L({\vec {w}}))_{i}|\leq Ch^{2}$ , wobei $i=1\ldots ,N$ und schließlich

$L({\vec {w}})+Ch^{2}{\bf {1}}\geq L_{h}{\vec {w}}\geq L({\vec {w}})-Ch^{2}{\bf {1}},\quad {\bf {1}}=(1,\ldots ,1)^{T}.$

Da $L(w)\geq {\bf {1}},$ s. oben, erhalten wir aus der rechen Ungleichung $L_{h}{\vec {w}}\geq {\bf {1}}-Ch^{2}{\bf {1}}.$

Folglich gilt für ausreichend kleine $h<h_{0}$, $h_{0}$ aus Folgerung 2

${\frac {1}{2}}{\bf {1}}\leq L_{h}{\vec {w}}.$

Da $L_{h}$ gleichzeitig eine M-Matrix ist, erhalten wir aus der Eigenschaft der Monotonie der Matrix $L_{h}$

${\frac {1}{2}}L_{h}^{-1}{\bf {1}}\leq {\vec {w}}.$

Weil $L_{h}^{-1}$ nur nichtnegative Elemente enthält ist $\||L_{h}^{-1}\||_{\infty }=\|L_{h}^{-1}{\bf {1}}\|_{\infty }$ und somit nach obiger Ungleichung $\||L_{h}^{-1}\||_{\infty }\leq 2\|{\vec {w}}\|_{\infty }\leq 2\max _{x\in [0,1]}|w|.$

Nach der Voraussetzung (für die zweite Konsistenzordnung) ist $w$ viermal stetig differenzierbar, also auch stetig und damit ist $\max _{x\in [0,1]}|w|$ beschränkt, folglich

$\||L_{h}^{-1}\||_{\infty }\leq K\in {\mathbb {R} }_{+}.$

Wir haben bewiesen:

Satz (Stabilität des Approximationschema $L_{h}$)

Sei $b\in C^{2}[0,1],\ c\in C[0,1],\ c\geq 0$.
Unter der Verwendung des zentralen Differenzenquotienten $D_{h}u(x)\approx u'(x)$ ist das Finite-Differenzen-Verfahren $L_{h}{\vec {u}}_{h}={\vec {f}}$  unter der Anwendung der Systemmatrix $L_{h}$ für elliptische Randwertprobleme stabil. $\quad \diamond$

Hauptsatz: Zweite Konvergenzordnung des FDM-Verfahrens

Sei $u\in C^{4}[0,1]$ die exakte Lösung des  elliptischen Randwertproblems,  ${\vec {u}}$ die Restriktion auf die Gitterpunkte und sei ${\vec {u}}_{h}$ die numerische Lösung, d.h. die Lösung von $L_{h}{\vec {u}}_{h}={\vec {f}}:={\vec {f}}_{h}$, wobei $L_{h}$   das Approximationschema mit zentralen Differenzenquotienten $D_{h}u(x)\approx u'(x)$ ist.
 

Dann gilt für eine hinreichend kleine Schrittweite $h$ die zweite Konvergenzordnung: $\|{\vec {u}}-{\vec {u}}_{h}\|_{\infty }\leq Ch^{2}$ mit einer Konstante $C>0$. $\quad \diamond$

Dieses Ergebniss lässt sich verallgemeinern auf Intervall [a,b]. Um die Konvergenz des FDM-Verfahrens auf einem Rechteck-Gebiet in 2D nachzuweisen, muss man lediglich die Monotonie-Eigenschaften der Blocktridiagonaler Systemmatrix untersuchen.

SW10: Numerische Diskretisierung der Reaktionsdiffusionsgleichung mit FDM

Im folgenden wir die numerische Diskretisierung des Reaktionsdiffusionsprozesses in der Populationsdynamik beschrieben.

$$\partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=f(u(x,t)).$$

Hier beschreibt die unbekannte Funktion $u$ die Populationsdichte der (aktuell oder kummulierten) Infizierten Individuen, wobei der Reaktionsterm $f=f(u)$ die zeitliche Dynamik der aktuell oder der kummulierten Infizierten steuert. Der Reaktionsterm stellt den Zusammenhang der räumlichen Infektionsverbreitung als Diffusionsprozess und der dynamischen Kompartimentmodellen dar.

Alternativ kann $f=f(u)$ auch den unbeschränkten Epideme-Ausbruch modellieren, siehe unbeschränktes Wachstum.

Räumliche Semidiskretisierung für die Diffusiongleichung

Zuerst wird Neumann Randwertproblem für die zeitabhängige homogene Diffusionsgleichung auf einem Rechteckgebiet D betrachtet.

$\quad \partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=0,\ x\in \ D\subset \mathbb {R} ^{2},\ t\in (t_{0},T),$

$\quad {\frac {\partial u(x)}{\partial {\vec {n}}}}\equiv \nabla u\cdot {\vec {n}}=g(x)\quad x\in \partial D.$

Im folgenden nehmen wir an, dass der Diffusionskoeffizient ist konstant, $c(x)=c$. Somit erhalten wir die homogene Diffusionsgleichung in der Form $\partial _{t}u(x,t)-c\Delta u(x,t)=0$.

System von gewöhnlichen Differentialgleichungen

Mithilfe der räumlichen Semidiskretisierung des Laplace Operators mit Finite-Differenzenverfahren unter Beachtung der homogenen Neumann Randbedingungen $\alpha =\beta =\gamma =\delta =0$ und Approximation der Richtungsableitungen erster Ordnung erhält man für ein festes $t\in (t_{0},T)$ nach dem Stern-Approximationsschema

$-c\Delta {\vec {u}}(t)\approx A_{h}{\vec {u}}_{h}$,

hierbei ist ${\vec {u}}(t)=(u(x_{0},y_{0},t),\ldots ,u(x_{n+1},y_{m+1},t))^{T}$ der Vektor der Restriktionen der exakten Lösung an die Gitterpunkte $(x_{i},y_{j}),\ i=0,\ldots ,n+1,\ j=0,\ldots ,m+1$ in Zeitpunkt $t$ und ${\vec {u}}_{h}(t)=(u_{0,0}(t),\ldots ,u_{n+1,m+1}(t))^{T}$ der Vektor der numerischen Lösung in den Gitterpunkten in Zeitpunkt $t$.

Die Matrix des Approximationsschema für Neumann Randbedingungen hat die Form:

$A_{h}={\frac {-c}{h^{2}}}{\begin{pmatrix}-hI&hI&&\ldots &0\\{\tilde {I}}&B&{\tilde {I}}&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\0&&\ldots &hI&-hI\end{pmatrix}},$

wobei $B={\begin{pmatrix}-h&h&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &\\0&\ldots &1&-4&1&\\0&\ldots &0&h&-h\end{pmatrix}},\quad {\tilde {I}}={\begin{pmatrix}0&&\ldots &&0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\ldots &0&1&0\\0&&\ldots &&0\end{pmatrix}}\quad {\tilde {I}},B\in {\mathbb {R} }^{(n+2)\times (n+2)},$
und der Einheitsmatrix $I=E_{n}\in {\mathbb {R} }^{(n+2)\times (n+2)}$.

Schließlich ergibt sich nach der räumlichen Semidiskretisierung folgendes System der gewöhnlichen Differentialgleichungen:

${\frac {d}{dt}}{\vec {u}}_{h}(t)=-A_{h}{\vec {u}}_{h}(t)$,

Der Vektor der numerischen Lösung ${\vec {u}}_{h}$ wird hier als Funktion von Zeit betrachtet.

Räumliche Semidiskretisierung für die Reaktionsdiffusiongleichung

Diskretisierung des Reaktionsterms, Gesamtinfizierte

Der Reaktionsterm $f(u(x,y,t))$ an den Gitterpunkten wird Mithilfe der Restriktionen der exakten Lösung ${\vec {u}}(t)$formuliert, für die kumulierte Infizierte (Logistisches Wachstum):

$f({\vec {u}}(t))\approx F({\vec {u}}_{h}(t)):={\frac {c}{\vec {B}}}*{\vec {u}}_{h}(t)*({\vec {B}}-{\vec {u}}_{h}(t))\equiv c{\begin{pmatrix}{\frac {u_{0,0}}{B_{0,0}}}(B_{0,0}-u_{0,0})\\\vdots \\{\frac {u_{n+1,m+1}}{B_{n+1,m+1}}}(B_{n+1,m+1}-u_{n+1,m+1})\end{pmatrix}}(t),$

wobei Vektor ${\vec {B}}=(B_{0,0},\ldots ,B_{n+1,m+1})^{T}$ der nach der Assemblierung aus der Matrix der Bevölkerungsdichte $(B)_{i,j}=B(x_{i},y_{j})$ entstanden ist und $*$ ist komponentenweise Multiplikation.

Bemerkung: die logistische Infektionsrate $k:={\frac {c}{\vec {B}}}$, vergleiche Modifiziertes SIR Modell berechnet sich als Rate der Infektionsverbreitung bei einer Kontaktrate: Wahrscheinlichkeit des Kontakts eines Infizierten mit einem Anfälligen.

Diskretisierung des Reaktionsterms, aktuell Infizierte

Im Falle der räumlichen Modellierung der Dichte der aktuellen Infizierten $u_{I}$ in Wechselwirkung mit Susceptibles $u_{S}$ nach dem Prinzip des SIR-Kompartimentmodells werden drei einander gekopelte Reaktionsdiffusionsgleichungen für die Dichtefunktionen $u_{I},\ u_{S},\ u_{R}$ gelöst.

$$\partial _{t}u_{S}(x,t)-div\ (c(x)\nabla u_{S}(x,t))=f_{S}(u_{S}(x,t),u_{I}(x,t)).$$

$$\partial _{t}u_{I}(x,t)-div\ (c(x)\nabla u_{I}(x,t))=f_{I}(u_{S}(x,t),u_{I}(x,t)).$$

$$\partial _{t}u_{R}(x,t)-div\ (c(x)\nabla u_{R}(x,t))=f_{R}(u_{I}(x,t)).$$

Diese partielle Differentialgleichungen werden durch die reche Seiten $f_{S},\ f_{I},\ f_{S}$ gekoppelt. Die entsprechende Reaktionsterme $f_{S}(u_{S},u_{I}),f_{I}(u_{S},u_{I})$ und $f_{R}(u_{I})$ werden an den Gitterpunkten ausgewertet:

$f_{S}({\vec {u}}(t))\approx F_{S}({\vec {u}}_{h}(t)):=-{\frac {c}{\vec {B}}}*{\vec {u_{I}}}_{h}(t)*{\vec {u_{S}}}_{h}(t)$

$f_{I}({\vec {u}}(t))\approx F_{I}({\vec {u}}_{h}(t)):={\frac {c}{\vec {B}}}*{\vec {u_{I}}}_{h}(t)*{\vec {u_{S}}}_{h}(t)-w{\vec {u_{I}}}_{h}(t),$

$f_{R}({\vec {u}}(t))\approx F_{R}({\vec {u}}_{h}(t)):=w{\vec {u_{I}}}_{h}(t),$

wobei $w$ die Wechselrate zwischen den Infizierten und Genesenen ist.

Da die erste und zweite partielle Differentialgleichung für $u_{I},\ u_{S}$ nicht von $u_{R}$ abhängig ist, kann man diese zwei Gleichungen separat lösen und $u_{R}$ dann im Nachgang berechnen.)

System von gewöhnlichen Differentialgleichungen

Nach der Zunahme des diskretisierten Reaktionstern $F({\vec {u}}_{h}(t))$ zum System von gewöhnlichen Differentialgleichungen des räumlich diskretisierten Neumann Randwertproblems für die Diffusion ergibt sich folgendes System der gewöhnlichen Differentialgleichungen (GDGL) für die Reaktionsdiffusionsgleichung:

${\frac {d}{dt}}{\vec {u}}_{h}(t)=-A_{h}{\vec {u}}_{h}(t)+F({\vec {u}}_{h}(t))$

mit einem linearen Anteil $-A_{h}{\vec {u}}_{h}(t)$ und einem nichtlinearen Anteil $F({\vec {u}}_{h}(t))$.

Dieses GDGL System, ergänzt um die Anfangsbedingung gegeben durch eine bekannte Funktion $u^{0}(x,y)$

${\vec {u}}(t_{0})={\vec {u}}_{h}(t_{0})={\vec {u}}^{0}=(u^{0}(x_{0},y_{0}),\ldots ,u^{0}(x_{n+1},y_{m+1}))^{T}$

definiert ein Anfangswertproblem für Reaktionsdiffusionsgleichung.

Zeitdiskretisierung der Reaktion-Diffusionsgleichung

Unter der Zeitdiskretisierung versteht man ein approximatives Vorgehen zur Bestimmung der Lösung eines Anfangswertproblems in diskreten Zeitpunkten. Der Zeitintervall $(t_{0},T)$ wird auf äquidistante Teilintervalle $(t_{k},t_{k+1})$ der Länge $\Delta t=t_{k+1}-t_{k},\quad k=0,1,\ldots K,\ T=t_{0}+K\Delta t$ geteilt.

Das einfachste numerische Verfahren, das Euler- oder Polygonzugverfahren entsteht durch das Ersetzten der Zeitableitung ${\frac {d{\vec {u}}_{h}}{dt}}$ mit Vorwärts- oder Rückwardsdifferenzenquotienten

$D_{\Delta t}^{\pm }{\vec {u}}_{h}(t)={\frac {{\vec {u}}_{h}(t\pm \Delta t)-{\vec {u}}_{h}(t)}{\Delta t}}\approx {\frac {d{\vec {u}}_{h}(t)}{dt}}$.

Dieses Verfahren hat die Konvergenzordnung 1, siehe Approximationsgüte der Differenzenquotienten.

Nach dem Anwenden des ersten Differenzenquotientes im obigen System von GDL ergeben sich folgende Approximationsschemas erster Ordnung:

Zeitdiskretisierung mit explizitem Eulerverfahren

Unter Verwendung von $D_{\Delta t}^{+}$ erhält man die Berechnungsformel

${\vec {u}}_{h}(t_{k+1})={\vec {u}}_{h}(t_{k})+\Delta t{\Big (}-A_{h}{\vec {u}}_{h}(t_{k})+F({\vec {u}}_{h}(t_{k})){\Big )},\quad k=0,1,\ldots K,\ T=t_{0}+K\Delta t.$

mit dem Startvektor $u^{0}(x,y)$

${\vec {u}}_{h}(t_{0})={\vec {u}}^{0}=(u^{0}(x_{0},y_{0}),\ldots ,u^{0}(x_{n+1},y_{m+1}))^{T}$

Die numerische Lösung in einem neuen Zeitschritt ${\vec {u}}_{h}(t_{k+1})$ erhält man iterativ durch Einsetzen der Lösung aus dem vorherigen Zeitschritt ${\vec {u}}_{h}(t_{k})$ in die rechte Seite der obigen Formel. Dieses Verfahren hat beschränkten Bereich der zulässigen Schrittweiten $\Delta t$ im Hinblick auf die Stabilität der numerischen Lösung, die Schrittweite $\Delta t$ muss aus diesem Grund ausreichend klein gewählt werden, ansonsten kann die numerische Lösung nach wenigen Schritten instabil werden und 'explodieren'.

Zeitdiskretisierung mit implizitem Eulerverfahren

Unter Verwendung von $D_{\Delta t}^{-}$ erhält man die implizite Berechnungsformel

${\vec {u}}_{h}(t_{k+1})={\vec {u}}_{h}(t_{k})+\Delta t{\Big (}-A_{h}{\vec {u}}_{h}(t_{k+1})+F({\vec {u}}_{h}(t_{k+1})){\Big )},\quad k=0,1,\ldots K,\ T=t_{0}+K\Delta t.$

mit dem Startvektor $u^{0}(x,y)$

${\vec {u}}_{h}(t_{0})={\vec {u}}^{0}=(u^{0}(x_{0},y_{0}),\ldots ,u^{0}(x_{n+1},y_{m+1}))^{T}$

Die numerische Lösung im neuen Zeitschritt ${\vec {u}}_{h}(t_{k+1})$ erhält man nach dem Lösen des obigen Systems von nichlinearen Gleichungen für ${\vec {u}}_{h}(t_{k+1})$ und benötigt weitere numerische Approximationsverfahren für nichtlineare Gleichungssysteme. Dieses Verfahren hat gute Stabilitätseigenschaften in Hinblick auf die Schrittweite $\Delta t$ (es ist A-stabil), die Schrittweite $\Delta t$ ist stabilitätsbedingt nicht beschränkt.

Weitere Zeitdiskretisierungsverfahren für gewöhnliche Differentialgleichungen

Sei eine gewöhnliche Differentialgleichung oder ein System von DGL gegeben allgemein als

${\bf {\dot {y}}}(t)=f(t,{\bf {y}}(t)),\quad {\bf {y}}:{\mathbb {R} }\to {\mathbb {R} }^{d},\ d=1,2,,\ldots .$

Unter Verwendung des zentralen Differenzenquotienten $D_{\Delta t/2}$ erhält man
die explizite Mittelpunktregel (verbessertes Eulerverfahren) zweiter Ordnug:

${\bf {y}}_{k+1}={\bf {y}}_{k}+\Delta t{\Big (}f(t_{k+1/{2}},{\bf {y}}_{k+1/{2}}){\Big )},\quad {\bf {y}}_{k+1/{2}}={\bf {y}}_{k}+{\frac {\Delta t}{2}}f(t_{k},{\bf {y}}_{k})$.

Ein weiteres Verfahren zweiter Ordnung ist die Trapezregel (Newmark scheme):

${\bf {y}}_{k+1}={\bf {y}}_{k}+{\frac {\Delta t}{2}}{\Big (}f(t_{k},{\bf {y}}_{k})+f(t_{k+1},{\bf {y}}_{k+1}){\Big )}$.

Für weitere Diskretisierungsverfahren für die Anfangswertprobleme ggf. höherer Konvergenzordnungen siehe Runge-Kutta und Mehrschrittverfahren.

Animation

Inhomogene Populationsdichte in unserem Beispiel, angepasst von Monitor Einwohnerzahl 2011

In folgender Animation wurde die Reaktionsdiffusiongleichung für kummulierte Infizierte mit explizitem Eulerverfahren auf einem Gebet von 15 x 18 km berechnet.

Animation der Epidemieausbreitung mit geographisch differenzierter Populationsdichte, Anzahl der Gesammtinfizierten pro Rasterzelle (220mx220m), konstanter Diffusionskoeffizient c=0.1, Infektionsrate 0.3.

Animation mit konstantem Diffusionskoeffizient $c=0.1$ und inhomogener Populationsdichte $B(x,y)$ mit $\max _{x\in D}B(x)=800/[km^{2}]$, siehe Bevölkerungsdichte Deutschland 2011.

Raumdiskretisierung des regionaldifferenzierten Diffusionskoffizienten

Um die räumlichen Epidemieausbreitung regional zu differenzieren, nehmen wir an, dass der Diffusionskoeffizient von der Populationsdichte abhängig ist, siehe nicht-konstanter Diffusionkoeffizient

$c(x,y)={\tilde {c}}(B(x,y))$.

Für diesen Fall muss Diffusionsmatrix (das Approximationsschema) $A_{h}$ entsprechend angepasst werden, hierbei werden anstatt eines konstanten Faktors $c/h^{2}$ vor der Matrix entsprechende Gitterwerte der Funktion $c(x,y)$ in der Matrix $A_{h}$ figurieren.

Diskretisierung der Randbedingungen

Die homogene Neumann-Randbedingung in diesem Fall lautet:

${\frac {\partial (cu)}{\partial {\vec {n}}}}\equiv c(x,y)\nabla u(x,y)\cdot {\vec {n}}$.

Mit Anwendung der einseitigen Differenzenquotienten und der Bezeichnung $c(x_{i},y_{j})=:c_{i,j}$ erhält man

• am linken und rechten Rand des Rechtecks:

$c_{0,j}{\frac {u_{1,j}-u_{0,j}}{h}}.(-1)=0,\quad c_{n+1,j}{\frac {u_{n+1,j}-u_{n,j}}{h}}.(1)=0,\quad j=0,1,\ldots m+1$,

• am unteren und oberen Rand des Rechtecks:

$c_{i,0}{\frac {u_{i,1}-u_{i,0}}{h}}.(-1)=0,\quad c_{i,m+1}{\frac {u_{i,m+1}-u_{i,m}}{h}}.(1)=0,\quad i=0,1,\ldots n+1$.

Diskretisierung des Diffusionsterms

Das Approximationsschema für den Diffusionsterm

$div{\Big (}c(x,y)\nabla u(x,y){\Big )}=\partial _{x}{\Big (}c(x,y)\partial _{x}u(x,y){\Big )}+\partial _{y}{\Big (}c(x,y)\partial _{y}u(x,y){\Big )}$

kann mit zweifacher Anwendung des zentralen Differenzenquotienten im obigen Ausdruck in den inneren Gitterpunkten $(x_{i},y_{j}),i=1,\ldots n,\ j=1,\ldots m$ hergeleitet werden:

$div{\Big (}c(x,y)\nabla u(x,y){\Big )}\approx D_{h/2,x}{\Big (}c(x_{i},y_{j})D_{h/2,x}u(x_{i},y_{j}){\Big )}+D_{h/2,y}{\Big (}c(x_{i},y_{j})D_{h/2,y}u(x_{i},y_{j}){\Big )}$ .

Herleitung des zweiten Differenzenquotienten bez. x:

$D_{h/2,x}{\Big (}c(x_{i},y_{j})D_{h/2,x}u(x_{i},y_{j}){\Big )}=D_{h/2,x}{\Big (}c_{i,j}{\frac {u_{i+1/2,j}-u_{i-1/2,j}}{h}}{\Big )}$,

$\quad ={\frac {1}{h}}{\Big (}{\frac {c_{i+1/2,j}u_{i+1,j}-c_{i-1/2,j}u_{i,j}}{h}}-{\frac {c_{i+1/2,j}u_{i,j}-c_{i-1/2,j}u_{i-1,j}}{h}}{\Big )}={\frac {1}{h^{2}}}{\Big (}u_{i+1,j}c_{i+1/2,j}-u_{i,j}(c_{i-1/2,j}+c_{i+1/2,j})+u_{i-1,j}c_{i-1/2,j}{\Big )}$,

hier bezeichnet $u_{i+1/2,j}:=u{\big (}x_{i}+{h}/{2}{\big )}$ die Werte der Funktion zwischen den Gitterpunkten $(x_{i},y_{j}),(x_{i+1},y_{j})$.

Analog kann man für den zweiten Differenzenquotienten bez. y zeigen:

$D_{h/2,y}{\Big (}c(x_{i},y_{j})D_{h/2,y}u(x_{i},y_{j}){\Big )}={\frac {1}{h^{2}}}{\Big (}u_{i,j+1}c_{i,j+1/2}-u_{i,j}(c_{i,j-1/2}+c_{i,j+1/2})+u_{i,j-1}c_{i,j-1/2}{\Big )}$.

Insgesamt gilt für den Diffusionsterm:

• $div{\Big (}c(x,y)\nabla u(x,y){\Big )}\approx {\frac {1}{h^{2}}}{\Big (}u_{i+1,j}c_{i+1/2,j}+u_{i,j+1}c_{i,j+1/2}-u_{i,j}(c_{i-1/2,j}+c_{i+1/2,j}+c_{i,j-1/2}+c_{i,j+1/2})+u_{i-1,j}c_{i-1/2,j}+u_{i,j-1}c_{i,j-1/2}{\Big )}.$

Approximationsschema

Unter der Anwendung der Approximation der Randbedingungen und des Diffusionsterms erhalten wir folgende Systemmatrix

$A_{h}={\frac {-1}{h^{2}}}{\begin{pmatrix}-hI_{0}&hI_{0}&&\ldots &0\\{\tilde {I}}_{\ell ,1}&B_{1}&{\tilde {I}}_{r,1}&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\0&&\ldots &hI_{m+1}&-hI_{m+1}\end{pmatrix}},$

wobei

$I_{0}=diag({\vec {c}}_{0}),\ I_{m+1}=diag({\vec {c}}_{m+1}),\qquad {\vec {c}}_{0}=(c_{0,0},\ldots ,c_{n+1,0}),\quad {\vec {c}}_{m+1}=(c_{0,m+1},\ldots ,c_{n+1,m+1}),$

$B_{j}={\begin{pmatrix}-hc_{0,j}&hc_{0,j}&0&\ldots &0\\c_{1/2,j}&-(c_{{\frac {1}{2}},j}+c_{{\frac {3}{2}},j}+c_{1,j+{\frac {1}{2}}}+c_{1,j-{\frac {1}{2}}})&c_{3/2,j}&\ldots &0\\\vdots &\ddots &\ddots &\ddots &\\0&\ldots &c_{n-1/2,j}&-(c_{n-{\frac {1}{2}},j}+c_{n+{\frac {1}{2}},j}+c_{n,j+{\frac {1}{2}}}+c_{n,j-{\frac {1}{2}}})&c_{n+1/2,j}&\\0&\ldots &0&hc_{n+1,j}&-hc_{n+1,j}\end{pmatrix}},$

${\tilde {I}}_{\ell ,j}={\begin{pmatrix}0&&\ldots &&0\\0&c_{1,j-{\frac {1}{2}}}&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\ldots &0&c_{n,j-{\frac {1}{2}}}&0\\0&&\ldots &&0\end{pmatrix}},\quad {\tilde {I}}_{r,j}={\begin{pmatrix}0&&\ldots &&0\\0&c_{1,j+{\frac {1}{2}}}&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &\vdots \\0&\ldots &0&c_{n,j+{\frac {1}{2}}}&0\\0&&\ldots &&0\end{pmatrix}},\quad j=1,\ldots m,\ \quad {\tilde {I}}_{\ell /r,j},B_{j}\in {\mathbb {R} }^{(n+2)\times (n+2)}$ und $I_{0},\ I_{m+1}=\in {\mathbb {R} }^{(n+2)\times (n+2)}$.

Bemerkung

In den obigen Matrizen kann man die Zwischenwerte der Funktion c durch die Mittelwerte ersetzen:

$c_{i\pm {\frac {1}{2}},j}\approx {\frac {1}{2}}(c_{i\pm 1,j}+c_{i,j}),$

$c_{i,j\pm {\frac {1}{2}}}\approx {\frac {1}{2}}(c_{i,j\pm 1}+c_{i,j}).$

In diesem Fall braucht man für den raumabhängigen Diffusionsokeffizient nur die Matrix der Werte an den Gitterpunkten $c_{i,j}$ speichern.

Für stückweise lineare Funktion $c$ sind die obigen Mittelwerte identisch mit den Zwischenwerten, d.h in den obigen Formeln für $c_{i\pm {\frac {1}{2}},j},\ c_{i,j\pm {\frac {1}{2}}}$ gilt Gleichheit '='.

Im Fall einer allgemeinen Funktion ist die Approximation durch die Mittelwerte von zweiter Ordnung, was nach der Addition der Taylorpolynome für $c{\Big (}x_{i}\pm {\frac {h}{2}},y_{j}{\Big )},$ bzw. $c{\Big (}x_{i},y_{j}\pm {\frac {h}{2}}{\Big )}$ folgt.

Seiteninformation

Dieser Wiki2Reveal Foliensatz wurde für den Lerneinheit Kurs:Räumliche Modellbildung' erstellt der Link für die Wiki2Reveal-Folien wurde mit dem Wiki2Reveal-Linkgenerator erstellt.

• Die Seite wurde als Dokumententyp PanDocElectron-SLIDE erstellt.

• siehe auch weitere Informationen zu Wiki2Reveal und unter Wiki2Reveal-Linkgenerator.

Literatur

• Wikiversity Beitrag zur Modellierung von COVID-19
• Wikipedia Artikel Taylorreihe
• Wikipedia Artikel zu Finiten-Differenzen-Methode
• Numerische Verfahren für Randwertaufgaben Technische Universität Hamburg-Harburg, Lehrmaterial