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+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$ .

Zweite Differenzenquotient

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 zum vierter 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 Differenzenquotient 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}^{2}f(x)=f''(x)+{\cal {O(h^{2})}}$ ,
$D_{h}^{\pm }f(x)=f'(x)+{\cal {O(h)}}$ .

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änktem 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

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 äquidistantem 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 bekannt.

Diskretisierung durch Differenzenquotient

In diesem Fall handelt es sich um gewöhnliche Differentialgleichung. Die zweiten Ableitungen werden mit Hilfe des zentralen 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 Gleichunssystem

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})

.

Daraus ergeben sich für $u_{-1},\ u_{n+2}$ neue Gleichungen:

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_{1}+2h\beta +{\mathcal {O}}(h^{3})

.

Damit man das lineare Gleichungssystem nicht wieder um neue Gleichungen ergänzen muss, 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.
So erhält man nach der Vernachlässigung der Fehlerterme Anstelle der ersten und letzter Gleichung des linearen Systems:
${\frac {2u_{1}-2u_{0}}{h^{2}}}=-2\alpha /h+{\mathcal {O}}(h^{2}),\quad {\frac {2u_{n}-2u_{n+1}}{h^{2}}}=-2\beta /h+{\mathcal {O}}(h^{2}).$

Der Fehler dieser Approximation ist insgesammt von der Größenordnung ${\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.

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 Neumannproblem 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üße über den Rand $\alpha$ und $\beta$ , bzw. $g$ mit dem Quellström $f$ im folgenden Sinne korrespondieren:
Im stationären Zustand balanziert der Materialfluß 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_{n+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 Gleichunssystem

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_{n+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 im Randpunkten des Rechteckts 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 Folgendem N die Gesamtanzahl der Gitterpunkte bezeichnet (in fereinfachter 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

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

$-c\Delta {\vec {u}}={\vec {f}}$ veträglich (Konsistenzbegriff).

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 $u$ des Originalproblems ein $C\in {\mathbb {R} }_{+}$ existiert

\displaystyle \|A_{h}{\vec {u}}-{\vec {f}}\|=\|A_{h}{\vec {u}}+c\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,\partial _{x_{2}}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}}=-c\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

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$ .
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 }=c\max _{i}\left|-D_{h}^{2}u(x_{i})+u''(x_{i})\right|\leq c{\frac {h^{2}}{12}}\max _{x\in D}|{u}^{(IV)}(x)|=c{\frac {h^{2}}{12}}\max _{x\in D}|{f}''(x)|.

Damit ist die Konstante $C$ aus der Definition der Konsistenz im eindimensionalem Fall $C={\frac {c}{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.

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)_{ij}^{-1}\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 ^[1] $\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 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 Vorausstetzung 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-diagonaldominan:
${\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}$ untersuchen wir in der Maximumnorm $\||\cdot \||_{\infty }$ (Zeilensummennorm), die von der Maximum-Vektornom $\|\cdot \|_{\infty }$ induziert ist, siehe . ^[2]

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:
Wäre $w\leq 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}$

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 Maximumnorm folg, 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, erhlaten 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,

\||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 die Systemmatrix  $L_{h}$   des Finite-Differenzen-Verfahren  $L_{h}{\vec {u}}_{h}={\vec {f}}$  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 Interval [a,b]. Um die Konvergenz des FDM- Verfahrens auf einem Rechteck-Gebiet in 2D nachzuweisen, muss mann lediglich die Monotonie-Eigenschaften der Blocktridiagonaler Systemmatrix untersuchen.

SW10: Numerische Diskretisierung der Reaktionsdiffusiongleichung 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äumliche 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 Aproximationsschema 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, Gesammtinfizierte

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 kummulierte 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 Befö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 Infektinsverbreitung bei einer Kontaktrate: Wahrscheinlichkeit des Kontakts eines Infizierten mit einem Anfälligen.

Diskretisierung des Reaktionsterms, aktuell Infizierte

Im Falle der Modellierung der Dichte der aktuellen Infizierten $u_{I}$ in Wechselwirkung mit Susceptibles $u_{S}$ nach dem Prinzip des SIR-Kompartimentmodells werden zwei einander gekopelte Reaktionsdiffusionsgleichungen für die Dichtefunktionen $u_{I},\ u_{S}$ gelöst. Die entsprechende Reaktionsterme $f_{S}(u_{S},u_{I}),f_{I}(u_{S},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),

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

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 ergiben 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 neuem Zeitschritt ${\vec {u}}_{h}(t_{k+1})$ erhält man iterativ durch einsetzen der Lösung aus vorherrigen Zeitschritt ${\vec {u}}_{h}(t_{k})$ auf 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 nnumerische 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 in neuem Zeitschritt ${\vec {u}}_{h}(t_{k+1})$ erhält man nach dem Lösen des obigen System von nichlinearem Gleichungen für ${\vec {u}}_{h}(t_{k+1})$ und benötigt weitere numerische Approximationsverfahren für nichlineare 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 Diffusionskoffizientes

Um die räumlichen Epidemieausbreitung regional zu differenzieren, nehmen wir an, dass der Diffusionskoeffizienzt 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 Faktor $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 linkem und rechtem 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 unterem und oberem 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 Differenzenqutiont 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 Differenzenqutiont 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 )}

.

Insgesammt 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 versetzen:

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 brauch 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 obige die 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 allgemeiner 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

↑ M. Hanke Burgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens, Kap. XV, Springer Verlag
↑ Natürliche Matrixnormen Wikiversity Artikel zu Matrixnom, abgerufen am 17.06.2020

This article is issued from Wikiversity. The text is licensed under Creative Commons - Attribution - Sharealike. Additional terms may apply for the media files.

[hb-1] M. Hanke Burgeois: Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens, Kap. XV, Springer Verlag

[wikiMatrixnorm-2] Natürliche Matrixnormen Wikiversity Artikel zu Matrixnom, abgerufen am 17.06.2020

[1]

[2]

SW6: Finite-Differenzen Methode FDM

Gitter-basierte Verfahren

Differenzenquotient

Taylorpolynom

Erster Differenzenquotient

Güte der Approximation

Zweite Differenzenquotient

Bemerkung

Güte der Approximation

Zusammenfassung

Differenzenquotienten höherer Ordnung

Randwertproblem

Dirichlet Randwertproblem

Neumann Randwertproblem

Numerische Diskretisierung des Randwertproblems

Dirichletproblem auf dem Intervall (1D)

Diskretisierung durch Differenzenquotient

Lineares Gleichunssystem

Systemmatrix und rechte Seite

SW7: Bemerkung: Lösbarkeit des linearen Gleichungssystems

SW7: Neumannproblem auf dem Intervall (1D)

Approximation erster Ordnung

Approximation zweiter Ordnung

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

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

Dirichletproblem in 2 Raumdimensionen

Punktegitter

Differenzenquotienten in 2D

Approximationsschema

Randwerte

Assemblierung

Lineares Gleichunssystem

Inhomogene Randbedingungen

Neumann-Problem in 2 Raumdimensionen

Neumann Randbedingungen auf dem Rechteck

Approximation erster Ordnung

Lineares Gleichungssystem

SW 8: Konvergenz und Stabilität

Konvergenz der numerischen Lösung

Voraussetzungen der Konvergenz

Definition der Konsistenz

Konsistenz elliptischer Randwertprobleme

Elliptisches Randwertproblem

Approximationschema des elliptischen Operators

Satz: Konsistenz des FDM Verfahren

Folgerungen

Definition der Stabilität

Monotonie Eigenschaften

Definiton (M-Matrix)

Folgerung: Monotonie Eigenschaft

SW9 Kriterien für M-Matrix

Satz: Monotonie tridiagonaler Matrizen

Folgerung 1

Folgerung 2

Monotonie der Systemmatrix und die Stabilität

Hilfsproblem

Beschränkheit der Matrix L h {\displaystyle L_{h}}

Satz (Stabilität des Approximationschema L h {\displaystyle L_{h}} )

Hauptsatz: Zweite Konvergenzordnung des FDM Verfahrens

SW10: Numerische Diskretisierung der Reaktionsdiffusiongleichung mit FDM

Räumliche Semidiskretisierung für die Diffusiongleichung

System von gewöhnlichen Differentialgleichungen

Räumliche Semidiskretisierung für die Reaktionsdiffusiongleichung

Diskretisierung des Reaktionsterms, Gesammtinfizierte

Diskretisierung des Reaktionsterms, aktuell Infizierte

System von gewöhnlichen Differentialgleichungen

Zeitdiskretisierung der Reaktion-Diffusionsgleichung

Zeitdiskretisierung mit explizitem Eulerverfahren

Zeitdiskretisierung mit implizitem Eulerverfahren

Weitere Zeitdiskretisierungsverfahren für gewöhnliche Differentialgleichungen

Animation

Raumdiskretisierung des regionaldifferenzierten Diffusionskoffizientes

Diskretisierung der Randbedingungen

Diskretisierung des Diffusionsterms

Approximationsschema

Bemerkung

Seiteninformation

Literatur

Beschränkheit der Matrix $L_{h}$

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