Numerische Modellierung der Diffusion

Aus Wikiversity
Zur Navigation springen Zur Suche springen

SW6: Finite-Differenzen Methode FDM[Bearbeiten]

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[Bearbeiten]

FDM mesh consisting of grid points, adaptation of https://commons.wikimedia.org/wiki/File:Digital_Image_ACC_Coordinate_Assignment.PNG
Diese drei Verfahren gehören zu den Gitter-basierten Verfahren, die die unbekannte Funktion auf endlich vielen Gitterpunkten approximieren.

Differenzenquotient[Bearbeiten]

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[Bearbeiten]

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 im Punkt um den Entwicklungspunkt

.

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

Erster Differenzenquotient[Bearbeiten]

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

mit . Damit ergibt sich

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

,

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

,

  • für den ersten zentralen Differenzenquotient durch aus der Summe ,

Güte der Approximation[Bearbeiten]

Alle drei erste Differenzenquotienten approximieren die erste Ableitung:
für falls hinreichend oft differenzierbar auf einem Intervall ist mit .

Lemma 1:

Ist , dann gilt für alle :
    wobei .
Ist , dann gilt  für alle :
    wobei . 

Beweis: Aus den Gleichungen für erste Differenzenquotienten .


Zweite Differenzenquotient[Bearbeiten]

Zweite Differenzenquotient approximert die zweite Ableitung an der Stelle falls hinreichend oft differenzierbar ist.

Die Formel für ergibt sich aus der Differenz der ersten Differenzenquotienten mithilfe der Taylorentwicklung bis zum vierter Ableitung:

Bemerkung[Bearbeiten]

Die Formel für den zweiten Differenzenquotient kann auch durch sukzessive Anwendung des ersten zentralen Differenzenquotient mit halber Schrittweite hergeleitet werden,
.

Güte der Approximation[Bearbeiten]

Lemma 2:

Ist , dann gilt für alle 
   wobei .

Beweis: Aus der Gleichung für zweiten Differenzenquotient. Der Term im Ausdruck für trägt zum Approximationsfehler bei .

Zusammenfassung[Bearbeiten]

Die Fehler der Differenzenquotienten kann man mithilfe des Landau-Symbols beschreiben:
,
,
.

Differenzenquotienten höherer Ordnung[Bearbeiten]

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



Randwertproblem[Bearbeiten]

Randwertproblem für Diffusionsgleichung beschreibt ein Diffusionsproblem auf einem beschränktem Gebiet . Hierbei wird das Verhalten der unbekannten Funktion am Rand des Gebietes durch eine Bedingung festgelegt. Wir behandeln zuerst die stationäre Poissongleichung.


Dirichlet Randwertproblem[Bearbeiten]

Ist die gesuchte Funktion am Rand durch eine gegebene Funktion festgelegt , erfüllt die sog.
Dirichlet Randbedingung:

Neumann Randwertproblem[Bearbeiten]

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

Numerische Diskretisierung des Randwertproblems[Bearbeiten]

Dirichletproblem auf dem Intervall (1D)[Bearbeiten]

Sei . Gesucht werden Näherungen für die Funktionswerte von auf dem äquidistantem Gitter

die Werte sind bekannt.

Diskretisierung durch Differenzenquotient[Bearbeiten]

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

für jedes .
Hier bezeichnet .

Lineares Gleichunssystem[Bearbeiten]

Nach der Anwendung der obigen Formel für jedes unter der Beachtung der Randbedingungen in den Formeln für für den Vektor den unbekannten Werte ein lineares Gleichungsystem

mit einer tridiagonaler Systemmatrix , und dem Vektor der rechten Seite ,

Systemmatrix und rechte Seite[Bearbeiten]

,

siehe dieses Beispiel.

SW7: Bemerkung: Lösbarkeit des linearen Gleichungssystems[Bearbeiten]

Die Systemmatrix ist regulär und damit invertierbar, symmetrisch und positiv definit. Der Lösungsvektor (numerische Lösung) 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) [Bearbeiten]

Gegeben seien die Randbedingungen für die Ableitungen in den Normalrichtungen :


Um die Ableitungen in den Randpunkten zu approximieren, führen wir neue Unbekannte ein.

Approximation erster Ordnung[Bearbeiten]

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

  • den vorwärts -Differenzenquotient
  • den rückwärts -Differenzenquotient

erhält man für die neuen Unbekannte die 0-te und die - te Gleichung,
,
.

[Bearbeiten]

Aufgabe Stellen Sie die Systemmatrix 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 .

Approximation zweiter Ordnung[Bearbeiten]


Um die Approximationsgüte der zweiten Differenzquotienten für die Approximation von im Inneren des Intervalls() auch für die Randbedingungen zu erhalten, verwendet man den ersten zentralen Differenzenquotient:

.

Daraus ergeben sich für neue Gleichungen:

,
.
[Bearbeiten]

Damit man das lineare Gleichungssystem nicht wieder um neue Gleichungen ergänzen muss, ersetzt man in den Gleichungen für den zweiten Differenzenquotient aus obigen Gleichungen.
So erhält man nach der Vernachlässigung der Fehlerterme Anstelle der ersten und letzter Gleichung des linearen Systems:

Der Fehler dieser Approximation ist insgesammt von der Größenordnung .

Aufgabe: Stellen Sie die Systemmatrix 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[Bearbeiten]


Die Systemmatrix ergänzt um die neuen Zeilen/Spalten für die Approximation der Neumann Randbedingung ist singulär, denn es gilt wobei .

Man kann zeigen (Gauß Eliminierung) .

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[Bearbeiten]
  • Wenn die Funktion eine Lösung des Neumannproblem ist, dann ist diese Lösung nicht eindeutig, auch 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. wird die Lösung eindeutig.


bzw. (mehrdimensional, Anwendung von Satz von Stokes, siehe ):
,
also müssen die Flüße über den Rand und , bzw. mit dem Quellström im folgenden Sinne korrespondieren:
Im stationären Zustand balanziert der Materialfluß durch die Ränder die Materialquelle im Inneren.

Dirichletproblem in 2 Raumdimensionen[Bearbeiten]



Sei u die gesuchte Funktion definiert auf einem Rechteck D:
. Wir betrachten das Dirichlet-Randwertproblem mit homogenen Randbedingungen:

Punktegitter[Bearbeiten]

Der einfachste Fall ist die Verwendung eines äquidistanten Gitters


mit konstanter Gittergröße h:


Man bezeichne die Approximationen der Lösung in den Gitterpunkten

Differenzenquotienten in 2D[Bearbeiten]


Die zweiten Differenzenquotienten approximieren die zweiten partiellen Ableitungen jeweils in den Richtungen x und y:

,

in jedem Gitterpunkt .

Approximationsschema[Bearbeiten]

Wir erhalten schließlich das Approximationsschema:


Randwerte[Bearbeiten]

In den Randgitterpunkten werden die Werte von nach der Dirichlet-Randbedingung folgend ersetzt:

,

.

Assemblierung[Bearbeiten]


Die unbekannten Approximationen der Lösung in den Gitterpunkten 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 (ist fest) in den Gitterpunkten von links nach rechts, , dann in der zweite Zeile u.s.w.

Die Approximationswerte sind in dem unbekannten Vektor dann wie folgt geordnet,


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

Lineares Gleichunssystem[Bearbeiten]

Nach der oben beschriebener numerischen Diskretisierung und Assemblierung des unbekannten Vektors ergibt sich ein lineares Gleichungssystem mit Blocktridiagonaler Matrix

mit Diagonalblock
und der Einheitsmatrix auf der Nebendiagonale.

Inhomogene Randbedingungen[Bearbeiten]




Sind die Randwerte der gesuchten Funktion unterschiedlich von Null, , tragen die bekannte Randwerte der gesuchten Funktion zum Vektor der rechten Seite bei.

Aufgabe: Seien Funktionen gegeben und die Randwerte der gesuchten Funktion am wie folgt festgelegt:

  • Linker/rechter Rand



  • unterer/oberer Rand



Wie verändert sich der Vektor der rechten Seite für diese Randbedingungen ?

Neumann-Problem in 2 Raumdimensionen[Bearbeiten]



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[Bearbeiten]

Seien die Ableitungen der unbekannten Funktion im Randpunkten des Rechteckts gegeben:

  • linker Rand:
  • rechter Rand:
  • unterer Rand:
  • oberer Rand:

wobei

Speziell oder am Rand eines Rechtecks.

Approximation erster Ordnung[Bearbeiten]



Bei der Anwendung der einseitigen Differenzenquotienten, analog wie 1-dimensionalem Fall muss der Vektor den Unbekannten um die Randwerte erweitert werden, diese werden zur Bildung der Differenzenquotienten verwendet:

  • linker Rand: den vorwärts -Differenzenquotient
  • rechter Rand: den rückwärts -Differenzenquotient
  • unterer Rand: den vorwärts -Differenzenquotient
  • oberer Rand: den rückwärts -Differenzenquotient

(Hier wurde direkt die obige Neumann-Randbedingung für ein Rechteckgebiet eingesetzt.)

Lineares Gleichungssystem[Bearbeiten]

Nach der Einordnung der unbekannter Werte in den Lösungsvektor wie obenbeschrieben (siehe Assemblierung ) ergibt sich wie zuvor eine erweiterte Blocktridiagonale Systemmatrix :

  • 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: 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 ist mit der Approximationsmatrix des Neumann-Problem in einer Raumdimension identisch.


Aufgabe: Stellen Sie die Systemmatrix 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[Bearbeiten]

Wir untersuchen, unter welcher Vorausstetzungen und mit welcher Approximationsgüte
die numerische Lösung für das Randwertproblem für Poissongleichung, (in 1D), ggf. (auf dem Recheteck in 2D), repräsentiert durch den Vektor

zu der exakten Lösung in den Gitterpunkten , ggf. konvergiert.

Im folgenden bezeichnen wir mit den Vektor der Restriktion der Funktion der exakten Lösung an die Gitterpunkte , ggf (nach der Assemblierung).


Konvergenz der numerischen Lösung[Bearbeiten]



Wir untersuchen den Unterschied der exakten und der numerischen Lösung in einer geeigneter Vektornorm auf , wobei im Folgendem N die Gesamtanzahl der Gitterpunkte bezeichnet (in fereinfachter Indexierung der Gitterpunkte).


Im folgenden werden wir zeigen, dass die Konvergenz für Gittergröße von

  • den Eigenschaften der Systemmatrix
  • der Approximationsgüte des Approximationschema (der Differenzenquotienten)

abhängig ist.



Voraussetzungen der Konvergenz[Bearbeiten]




Sei Matrix invertierbar und die Lösung des linearen Systems .

Unter der Anwednung der Verträglichkeit der Matrixnorm erfolgt

wobei die Vektornorm-induzierte (natürliche) Matrixnorm bezeichnet.

Folglich ist die Konvergenz gegeben falls

  1. die Norm von gleichmäßig (bez. ) beschränkt ist: - (Stabilitätsbegriff)
  2. , für das Approximationsschema ist mit dem Originalproblem,

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[Bearbeiten]



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

Definition: (Konsistenzordnung)
Das Finite Differenzenverfahren für die Poisson Randwertaufgabe hat bezüglich einer Vektornorm in die Konsistenzordnung wenn für jede hinreichend oft stetig differenzierbare Lösung des Originalproblems ein existiert


Konsistenz elliptischer Randwertprobleme[Bearbeiten]



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 für die Funktion ausbreiten.

Elliptisches Randwertproblem[Bearbeiten]

Vorausgesetzt stetige Funktionen auf sind gegeben und .


Approximationschema des elliptischen Operators[Bearbeiten]

Das entsprechende Approximationsschema zum Differenzialoperator enthält im Fall vom Rechteckgebiet das Stern-Approximationsschema für das Laplace Operator anhand der zweiten Differenzenquotienten (Matrix ) und ein weiteres Beitrag durch die Approximation des Gradienten mit einseitigen, oder zentralen Differenzenquotienten , siehe erste Differenzenquotient.
Die Konsistenzdefinition für den elliptischen Operator lautet:

Satz: Konsistenz des FDM Verfahren[Bearbeiten]



Sei  die exakte Lösung des obigen elliptischen Randwertproblems  mit homogenen Randbedingungen. 

Dann hat das Finite-Differenzen-Verfahren mit dem Approximationsschema 
i) die Konsistenzordnung  unter der Anwendung von einseitigen Differenzenquotienten 
ii) die Konsistenzordnung  unter der Anwendung von zentralen Differenzenquotienten 
für den Term . 

Beweis: basiert auf der Approximationsgüte der ersten und der zweiten Differenzenquotienten (Lemma 1 und 2). Für das Approximationsschema ergibt sich daraus , für das Approximationsschema dagegen entweder oder .

Folgerungen[Bearbeiten]



  1. Im Fall Rechteckgebiet, erhalten wir die Randwertaufgabe für Poissongleichung mit homogenen Dirichletrandbedingungen. In diesem Fall ist das Approximationsschema von der Konsistenzordnung .
  2. Im Fall dass ein geschlossenes Intervall ist erhalten wir aus Lemma 2 in der Maximum-Vektornorm für konkret

Damit ist die Konstante aus der Definition der Konsistenz im eindimensionalem Fall

Definition der Stabilität[Bearbeiten]



Definition: (Stabilität)
Das Finite-Differenzen-Verfahren für das elliptische Randwertproblem heißt stabil (bezüglich einer Vektornorm ),
wenn existieren, , sodass die Matrix (Approximationsschema) invertierbar ist und es gilt für .

Die Untersuchung der Stabilität anhand der Konstruktion der inversen Matrix ist aufwändig. Deswegen wird die Stabilität anhand anderer Kriterien untersucht, der Monotonie Eigenschaften der Matrix .

Bemerkung:

Die Stabilitätsbedingung garantiert auch die Beschränkheit der numerischen Lösung als Lösung des linearen Systems , denn in einer Vektornorm.

Monotonie Eigenschaften[Bearbeiten]



Definiton (M-Matrix)[Bearbeiten]

Eine reguläre Matrix mit für

und , kurz , heißt M-Matrix, oder Monotonie-Matrix.

Folgerung: Monotonie Eigenschaft[Bearbeiten]

Für eine M-Matrix und Vektoren gilt:


Beweis: Ist da die Elemente

SW9 Kriterien für M-Matrix[Bearbeiten]



In diesem Abschnitt werden die Kriterien für Tridiagonalmatrizen wie die Systemmatrix des Laplace-Approximationsschemas in 1D und die Matrix des elliptischen Randwertproblems in einer Raumdimension untersucht.

Satz: Monotonie tridiagonaler Matrizen[Bearbeiten]
Jede irreduzibel diagonaldominante  Tridiagonalmatrix mit positiven Diagonalelementen und negativen Nebendiagonalelementen ist eine M-Matrix. 


Beweis: siehe [1]

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

Folgerung 1[Bearbeiten]
Die Matrix    des Laplace-Approximationsschemas in 1D eine M-Matrix.

Beweis:
ist tridiagonal, hat alle Sub- und Superdiagonalelementen ungleich Null und für erste und letzte Zeile gilt . Damit ist eine M-Matrix .

Folgerung 2[Bearbeiten]



Die Matrix    des elliptischen Randwertproblems in 1D unter der Verwendung des  ersten zentralen Differenzenquotienten    für  ist  für 
   eine M-Matrix.

Beweis:

ist tridiagonal, mit


Nach der Vorausstetzung ist und somit sind und negativ.
ist irreduzibel.
ist auch schwach-diagonaldominan:
für wobei für die strikte Ungleichung gilt .

Monotonie der Systemmatrix und die Stabilität[Bearbeiten]



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

Die Beschränkheit der Matrix untersuchen wir in der Maximumnorm (Zeilensummennorm), die von der Maximum-Vektornom induziert ist, siehe . [2]


Hilfsproblem[Bearbeiten]

Wir untersuchen folgendes Randwertproblem: gegeben sei :

Sei die Lösungsfunktion .

Lemma 3:

Die Lösung  des obigen Randwertproblem ist .

Beweis:
Wäre mindesten in einem Punkt, hätte ein lokales Minimum in einem . d.h.,
und .

Damit erhalten wir aus der obigen Differenzialgleichung , was im Widerspruch zu der Bedingung des lokalen Minima steht.

Beschränkheit der Matrix [Bearbeiten]



Sei zusätzlich eine stetige nichtnegative Funktion gegeben.

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

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


Aus dieser Abschätzung Maximumnorm folg, dass auch gilt , wobei und schließlich



Da s. oben, erhalten wir aus der rechen Ungleichung

Folglich gilt für ausreichend kleine , aus Folgerung 2

Da gleichzeitig eine M-Matrix ist, erhlaten wir aus der Eigenschaft der Monotonie der Matrix

Weil nur nichtnegative Elemente enthält ist und somit nach obiger Ungleichung

Nach der Voraussetzung (für die zweite Konsistenzordnung) ist viermal stetig differenzierbar, also auch stetig und damit ist beschränkt,

Wir haben bewiesen:

Satz (Stabilität des Approximationschema )[Bearbeiten]

Sei .
Unter der Verwendung des zentralen Differenzenquotienten ist die Systemmatrix des Finite-Differenzen-Verfahren für elliptische Randwertprobleme stabil.

Hauptsatz: Zweite Konvergenzordnung des FDM Verfahrens[Bearbeiten]



Sei  die exakte Lösung des  elliptischen Randwertproblems,   die Restriktion auf die Gitterpunkte und sei  die numerische Lösung, d.h. die Lösung von , wobei    das Approximationschema mit zentralen Differenzenquotienten  ist. 
Dann gilt für eine hinreichend kleine Schrittweite die zweite Konvergenzordnung: mit einer Konstante .


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[Bearbeiten]

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

Hier beschreibt die unbekannte Funktion die Populationsdichte der (aktuell oder kummulierten) Infizierten Individuen, wobei der Reaktionsterm 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 auch den unbeschränkten Epideme-Ausbruch modellieren, siehe unbeschränktes Wachstum.

Räumliche Semidiskretisierung für die Diffusiongleichung[Bearbeiten]

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

Im folgenden nehmen wir an, dass der Diffusionskoeffizient ist konstant, . Somit erhalten wir die homogene Diffusionsgleichung in der Form .

System von gewöhnlichen Differentialgleichungen[Bearbeiten]

Mithilfe der räumliche Semidiskretisierung des Laplace Operators mit Finite-Differenzenverfahren unter Beachtung der homogenen Neumann Randbedingungen und Approximation der Richtungsableitungen erster Ordnung erhält man für ein festes nach dem Stern-Approximationsschema

,

hierbei ist der Vektor der Restriktionen der exakten Lösung an die Gitterpunkte in Zeitpunkt und der Vektor der numerischen Lösung in den Gitterpunkten in Zeitpunkt .

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




wobei
und der Einheitsmatrix .


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

,

Der Vektor der numerischen Lösung wird hier als Funktion von Zeit betrachtet.

Räumliche Semidiskretisierung für die Reaktionsdiffusiongleichung[Bearbeiten]

Diskretisierung des Reaktionsterms, Gesammtinfizierte[Bearbeiten]

Der Reaktionsterm an den Gitterpunkten wird Mithilfe der Restriktionen der exakten Lösung formuliert, für die kummulierte Infizierte (Logistisches Wachstum):

wobei Vektor der nach der Assemblierung aus der Matrix der Befölkerungsdichte entstanden ist und ist komponentenweise Multiplikation.

Bemerkung: die logistische Infektionsrate , 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[Bearbeiten]

Im Falle der Modellierung der Dichte der aktuellen Infizierten in Wechselwirkung mit Susceptibles nach dem Prinzip des SIR-Kompartimentmodells werden zwei einander gekopelte Reaktionsdiffusionsgleichungen für die Dichtefunktionen gelöst. Die entsprechende Reaktionsterme werden an den Gitterpunkten ausgewertet:

wobei die Wechselrate zwischen den Infizierten und Genesenen ist.

System von gewöhnlichen Differentialgleichungen[Bearbeiten]

Nach der Zunahme des diskretisierten Reaktionstern 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:

mit einem linearen Anteil und einem nichtlinearen Anteil .

Dieses GDGL System, ergänzt um die Anfangsbedingung gegeben durch eine bekannte Funktion

definiert ein Anfangswertproblem für Reaktionsdiffusionsgleichung.

Zeitdiskretisierung der Reaktion-Diffusionsgleichung[Bearbeiten]

Unter der Zeitdiskretisierung versteht man ein approximatives Vorgehen zur Bestimmung der Lösung eines Anfangswertproblems in diskreten Zeitpunkten. Der Zeitintervall wird auf äquidistante Teilintervalle der Länge geteilt.

Das einfachste numerische Verfahren, das Euler- oder Polygonzugverfahren entsteht durch das Ersetzten der Zeitableitung mit Vorwärts- oder Rückwardsdifferenzenquotienten

.

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 [Bearbeiten]

Unter Verwendung von erhält man die Berechnungsformel

mit dem Startvektor

Die numerische Lösung in neuem Zeitschritt erhält man iterativ durch einsetzen der Lösung aus vorherrigen Zeitschritt auf die rechte Seite der obigen Formel. Dieses Verfahren hat beschränkten Bereich der zulässigen Schrittweiten im Hinblick auf die Stabilität der numerischen Lösung, die Schrittweite 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 [Bearbeiten]

Unter Verwendung von erhält man die implizite Berechnungsformel

mit dem Startvektor

Die numerische Lösung in neuem Zeitschritt erhält man nach dem Lösen des obigen System von nichlinearem Gleichungen für und benötigt weitere numerische Approximationsverfahren für nichlineare Gleichungssysteme. Dieses Verfahren hat gute Stabilitätseigenschaften in Hinblick auf die Schrittweite (es ist A-stabil), die Schrittweite ist Stabilitätsbedingt nicht beschränkt.

Weitere Zeitdiskretisierungsverfahren für gewöhnliche Differentialgleichungen[Bearbeiten]

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

Unter Verwendung des zentralen Differenzenquotienten erhält man
die explizite Mittelpunktregel (verbessertes Eulerverfahren) zweiter Ordnug:

.

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

.


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

Animation[Bearbeiten]

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 und inhomogener Populationsdichte mit , siehe Bevölkerungsdichte Deutschland 2011.



Raumdiskretisierung des regionaldifferenzierten Diffusionskoffizientes[Bearbeiten]

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

.

Für diesen Fall muss Diffusionsmatrix (das Approximationsschema) entsprechend angepasst werden, hierbei werden anstatt eines konstanten Faktor vor der Matrix entsprechende Gitterwerte der Funktion in der Matrix figurieren.

Diskretisierung der Randbedingungen[Bearbeiten]

Die homogene Neumann Randbedingung in diesem Fall lautet:

.

Mit Anwendung der einseitigen Differenzenquotienten und der Bezeichnung erhält man

  • am linkem und rechtem Rand des Rechtecks:
,
  • am unterem und oberem Rand des Rechtecks:
.

Diskretisierung des Diffusionsterms[Bearbeiten]

Das Approximationsschema für den Diffusionsterm

kann mit zweifacher Anwendung des zentralen Differenzenquotienten im obigen Ausdruck in den inneren Gitterpunkten hergeleitet werden:

.

Herleitung des zweiten Differenzenqutiont bez. x:

,
,

hier bezeichnet die Werte der Funktion zwischen den Gitterpunkten .

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

.

Insgesammt gilt für den Diffusionsterm:

Approximationsschema[Bearbeiten]

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



wobei






und .

Bemerkung[Bearbeiten]

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

In diesem Fall brauch man für den raumabhängigen Diffusionsokeffizient nur die Matrix der Werte an den Gitterpunkten speichern.

Für stückweise lineare Funktion sind obige die Mittelwerte identisch mit den Zwischenwerten, d.h in den obigen Formeln für gilt Gleichheit '='.

Im Fall einer allgemeiner Funktion ist die Approximation durch die Mittelwerte von zweiter Ordnung, was nach der Addition der Taylorpolynome für bzw. folgt.

Seiteninformation[Bearbeiten]

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.



Literatur[Bearbeiten]

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