Kurs:Räumliche Modellbildung/Gruppe 16

Gruppenseite - (EGL)

Diese Seite enthält die gesammelten Materialien der Gruppe 16 für die Portfolioprüfung.

• COVID-19/Mathematische Modellbildung

Teilnehmer-innen

• Ellerwald, Jonas
• von Gerichten, Marco
• Lapport, Felix

Diskrete Modellierung

Infektionsepidemiologie

Unter der Infektionsepidemiologie versteht man eine medizinische Fachdisziplin und ein Teilgebiet der Epidemiologie, welche sich mit der Ausbreitung und der Eindämmung von Infektions- und parasitären Krankheiten, unter Einbezug ihres räumlichen und zeitlichen Auftreten, beschäftigt. Dabei können die meisten Infektionskrankheiten mathematisch modelliert werden, um eben das bereits angesprochene zeitliche und räumliche Verhalten zu untersuchen und zu prognostizieren.

Zyklus 1: SIR-Modell

Erste Gedankenspiele

Wir wollten eine Graphische Darstellung des SIR Modells anhand der Zahlen der Stadt Kaiserslautern. Die Einwohnerzahl von Kaiserslautern beträgt ca. 100.000. Zusätlich haben wir die Infektionsrate auf 30% täglich festgelegt.

Um uns die Funktionen der Infizierten, Empfänglichen und Genesenen deutlich zu machen, begannen wir mit einer Exeltabelle

Anhand dieser Tabelle und der Bildlichen Programierung der Zellen entwickelten wir das Pendant in Octave:

S = 100885 ;  %Möglichen Infizierten
I = 0  ; %Infizierten
R = 0  ;  %Genesenen

k = 0.3 ;   %Infektionsrisiko
w = 1/14;   %Genesungszeit
Zeit = 30;  %Laufzeit des Modells in Tagen

Erstellung von Vectoren von S, I, R in der länge der zu berechnenden Zeit

Sneu = [0:Zeit];
Ineu = [0:Zeit];
Rneu = [0:Zeit];

Vectoren Nullsetzen

Sneu(2:Zeit) = 0;
Ineu(2:Zeit) = 0;
Rneu(1:Zeit) = 0;

Eintragung des Startwerts an der ersten Stelle des Vectors

Sneu(1) = S;
Ineu(1) = I;
Rneu(1) = R;

Da uns die Genesungsrate von 1/14 nicht zusagte, da dies nicht wirklich der konkreten Darstellung des Genesungsprozesses abbildet, entschieden wir uns für die Genesung nach 14 Tagen. Wegen der Genesungszeit von 14 Tagen, mussten wir für den Zähler eine if Abfrage einbauen, um den Fehler des Programms bezüglich der Vectoren abfrage zu vermeiden.

for t = 1:Zeit
 
 if t >= 15
   Rneu(t+1) = Ineu(t-14);
 else
   Rneu(t+1) = 0; 
 endif

 Ineu(t+1) = k*(Sneu(t))+Ineu(t);
 Sneu(t+1) = Sneu(1) - Ineu(t+1); 
 
endfor

Die ersten Gedankenspiele unserer weiteren epidemiologischen Untersuchung haben wir nun zunächst unabhängig von Wechselwirkungen der Orte für die Ortschaften Kaiserlautern, Neustadt an der Weinstraße, Landau, Pirmasens, Impflingen und Rumbach ausgeweitet Dementsprechend haben wir die Excel-Programmierung bzw. die Libre Office-Programmierung auch durch diese Orte erweitert.

Fazit & weitere Optimierungsmöglichkeiten

Der Zyklus 1 ist eine zunächst noch sehr einfaches Modell, das das Infektionsgeschehen in unseren gewählten Knoten ohne Bewegung der Bevölkerung modelliert. Die Knoten selbst liegen zumindest in Form der Städte Kaiserlautern, Neustadt, Landau und Pirmasens entlang der Schnell- und Kraftfahrtstraßen und bilden somit die Knotenpunkte eines Wegenetzes. Die Modellierung in Excel und Octave des SIR-Modells können für unsere weiteren Zyklen erneut genutzt werden und bilden die Basis für das weitere Vorgehen. Letztendlich fehlt es dem Zyklus jedoch an den Verbindungen zwischen den Knoten und somit auch an Personenverkehr zwischen unseren gewählten Orten. Somit ergeben sich bereits weitere Möglichkeiten einer Optimierung des Modells. So werden im nächsten Schritt und dementsprechend im zweiten Zyklus die Verbindungen zwischen den Knoten ergänzt und der Personenverkehr könnte anhand einer Transportmatrix simuliert werden.

Zyklus 2: SIR-Modell mit Transportmatrix

In Zyklus zwei wurde dann der von uns ausgewählte Raum durch eine endliche Menge von Knoten und Verbindungen mittels eines ungerichteten Graphen diskretisiert. Der Graph kodiert dabei das Wegenetz von Impflingen über Landau, Neustadt, Kaiserslautern, Pirmasens nach Rumbach. Es ergibt sich ein ungerichteter Graph mit 6 Knoten wie der Abbildung zu entnehmen ist. Mit der Menge der Knoten $V=\{Kaiserslautern,Neustadt,Landau,Pirmasens,Impflingen,Rumbach\}$ und der Menge $E$ der ungerichteten Knotenverbindungen: $E=\{\{Impflingen,Landau\},\{Landau,Kaiserslautern\},\{Landau,Neustadt\},\{Neustadt,Kaiserslautern\},\{Kaiserslautern,Pirmasens\},\{Pirmasens,Rumbach\}\}$

Anzahl_Orte = 6
vec_u0 = [100000 53148 46677 40403 805 472]  ;
Matrix_Tag_0 = [100000 1 0; 53148 0 0; 46677 0 0; 40403 0 0; 805 0 0; 472 0 0];
r = 0.03 ;
Tage = 5;
w = 1/14;

Die Verbindungen und Übergänge zwischen den jeweiligen Ortsbewohnern haben wir mit Hilfe sogenannter Übergangsmatrizen in Octave modelliert. Unter Beachtung der vorhandenen Verbindungen aus Menge $E$ haben wir eine Transportmatrix erstellt.

#Implementieren der Transportmatrix für 6 Orte
Transportmatrix = [0.6 0.3 0.3 0.3 0 0; 0.2 0.5 0.2 0 0 0; 0.1 0.2 0.3 0.2 0.2 0; 0.1 0 0.18 0.48 0 0.15; 0 0 0.02 0 0.8 0; 0 0 0 0.02 0 0.85 ]

Transportmatrix

Ort/Stadt             Kaiserslautern| Neustadt a. d. W. | Landau i. d. P. | Pirmasens | Impflingen | Rumbach 
---------------------------------------------------------------------------------------------------------------
Kaiserslautern    |       0.6       |        0.3        |       0.3       |     0.3   |      0     |    0 
Neustadt a. d. W. |       0.2       |        0.5        |       0.2       |      0    |      0     |    0
Landau i. d. P.   |       0.1       |        0.2        |       0.3       |     0.2   |     0.2    |    0
Pirmasens         |       0.1       |         0         |       0.18      |     0.48  |      0     |   0.15
Impflingen        |        0        |         0         |       0.02      |      0    |     0.8    |    0
Rumbach           |        0        |         0         |         0       |     0.02  |      0     |   0.85

Grundsätzlich gibt diese uns Auskunft über den Personenverkehr zwischen den jeweiligen Verbindungen. Die Einträge mit Null kennzeichnen die nicht vorhandenen Verbindungen zwischen den Knoten. Somit kann aus epidemiologischer Sicht kein Austausch von Personen stattfinden, da diese keine epidemiologischen Nachbarn sind. Mit einer for-Schleife haben wir die täglichen Bewegungen zwischen den Knoten simuliert und für eine endliche Anzahl von Tagen laufen lassen.

 for i = 1:Tage
   %Impffaktor wird festgelegt
   Impffaktor = 0;
   %Ab 4tem Tag wird geimpft (mit dem Faktor 0.0005)
   if i >= 4
     Impffaktor = 0.005;
   endif
   %Iteration läuft über alle Orte
   for j = 1:Anzahl_Orte
     %Epidemiologischer Prozess
     Ergebnismatrix(j,:) = SIR_Impf_Vector(Ergebnismatrix(j,1),Ergebnismatrix(j,2),Ergebnismatrix(j,3),r,w,Impffaktor);
     %Räumlicher Prozess
     %if j == 1 || j == 3 
       %Die Infizierten bleiben an ihrem Ort (j=1 -> S ; j=3 -> R) für die die Transportmatrix angewendet wird
       %Ergebnismatrix(:,j) = Transportmatrix*Ergebnismatrix(:,j);
     %endif
 
   endfor
 Ergebnismatrix
 endfor

Fazit & Ausblick

Durch das Bilden eines ungerichteten Graphen konnten wir unsere Knoten verbinden und das Wegenetz und somit die Transportprozesse zwischen den jeweiligen Knoten darstellen. Über den Transport und das Anwenden des SIR-Modells kann das Pandemiegeschehen in den jeweiligen Orten modelliert werden. Dies zeigt wie stark der Einfluss der Mobilität auf das jeweilige Geschehen ist. Für den weiteren Zyklus wäre es nun eine Optimierung das SIR-Modell durch das modifizierte SIR-Modell zu ersetzen.

Zyklus 3: Modifiziertes SIR-Modell mit Transportmatrix

Zunächst haben wir wieder unseren ungerichteten Graphen mit unseren 6 Knoten bzw. unseren 6 gewählten Orten in Octave implementiert. Der Vektor vec_u0 gibt die jeweilige Bevölkerungszahl der Orte und beschreibt somit die für das modifizierte SIR-Modell wichtige Bevölkerungskapazität in jedem gewählten Knoten aus der Menge $V$. Die Matrix $Matrix_{T}ag_{0}$ gibt uns die Startwerte im Zeitpunkt $t=0$ für unsere Kompartimente. Für die Vereinfachung des Modells haben wir in jedem Knoten einen Infizierten gewählt, der das Infektionsgeschehen in jedem Ort in Gang bringt. Die weiteren wichtigen Variablen für das modifizierte SIR-Modell haben wir anhand unserer Modellierung im Tutorium 3 gewählt.

vec_u0 = [100000 53148 46677 40403 805 472];  
Matrix_Tag_0 = [100000 1 0 0 ; 53148 1 0 0; 46677 1 0 0; 40403 1 0 0; 805 1 0 0 ; 472 1 0 0];
Anzahl_Orte = length(vec_u0)
r = 1 ;
c = 0.3435;
d = 0.03;
T = 100;
w = 1/14;
B = [100000 53148 46677 40403 805 472];

Im nächsten Schritt haben wir in einer Funktion das modifizierte SIR-Modell implementiert.

function wert = SIR_Vec_mod (S, I, R, T,w,c,d,r,B)
 Rneu = w*I - d*I +R  ;
 Ineu =(c/B)*S*I+I - w*I;
 Sneu = S-c*S*I/(r*B) ;
 Tneu = d*I+T ;
 wert=[Sneu Ineu Rneu Tneu];
endfunction

Ergebnismatrix=Matrix_Tag_0;
SIRmatrix = ones(length(vec_u0),4);

Transportmatrix = [0.6 0.3 0.3 0.3 0 0; 0.2 0.5 0.2 0 0 0; 0.1 0.2 0.3 0.2 0.2 0; 0.1 0 0.18 0.498 0 0.15; 0 0 0.02 0 0.8 0; 0 0 0 0.002 0 0.85 ];

Mittels der for-Schleife wird das SIR-Modell für jeden Tag einmal durchgeführt. Die zweite for-Schleife gibt uns die Berechnung der Ergebnisse für jeden Ort.

for i = 1:T
   for j = 1:Anzahl_Orte
     %Epidemiologischer Prozess
      Ergebnismatrix(j,:) = SIR_Vec_mod (Ergebnismatrix(j,1),Ergebnismatrix(j,2),Ergebnismatrix(j,3),Ergebnismatrix(j,4),w,c,d,r,B(1,j));%Räumlicher Prozess
       for k = 1:4  
           if k == 1
             Ergebnismatrix(:,k) = Transportmatrix*Ergebnismatrix(:,k);
           else if k == 2
             Ergebnismatrix(:,k) = Transportmatrix*Ergebnismatrix(:,k);
           else if k == 3 
             Ergebnismatrix(:,k) = Transportmatrix*Ergebnismatrix(:,k);
           else if k == 4
             Ergebnismatrix(:,k) = Ergebnismatrix(:,k); 
           endif
           endif
           endif
           endif 
       endfor
   endfor 
Ergebnismatrix
endfor

Das sich die Transportmatrix lediglich in den ersten paar Zeitschritten bemerkbar macht, liegt daran, dass sich das System durch die Transportmatrix sehr schnell stabilisiert und sich, nach kurzer Zeit, keine weiteren Veränderungen mehr ergeben.

Fazit und weitere Optimierungsmöglichkeiten

Weitere Optimierungsmöglichkeiten wäre die Darstellung in Säulen in den jeweiligen Knoten, um den Verlauf der Pandemie auf unserer diskretisierten Fläche von Rheinland-Pfalz in den Knoten deutlicher darstellen zu können.

Zyklus 4: SIR-Modell Darstellung in 3D

Im ersten Schritt des vierten Zyklus haben wir mit GeoGebra gearbeitet. Zunächst haben wir unsere Karte mit den markierten Knoten in ein 2D-Koordinatensystem eingefügt. Danach haben wir die Karte bzw. die Grafik der Karte so gelegt, dass sie im ersten Quadranten des Koordinatensystem liegt. Dies erleichtert uns im Endeffekt die Bestimmung der Koordinaten. Dementsprechend wurde die Grafik auf der linken Seite auf die y-Achse gelegt und mit dem unteren Rand auf die x-Achse. Daraufhin haben wir durch einfaches Einfügen von Punkten die Koordinaten unserer Knoten erzeugt. Durch Darstellen dieser Punkte in der 3D-Ansicht ergibt sich zumindest der rechteckige Bereich, auf welchem die gewählten Knoten liegen. Anhand unserer somit gefunden Koordinaten konnten wir nun weiterarbeiten.

Einlesen der Koordinaten in Octave

Zunächst haben wir das Gebiet in x- und y-Richtung definiert:

step=0.1;
X=[0:step:7];
Y=[0:step:7.5];
[x,y]=meshgrid(X,Y);

function wert= hallo_funktion(x,y,z1,z2,z3,z4,z5,z6)
 
  if sqrt((x.-3.9).^2+(y.-7.2).^2)<=0.3 wert=z1;
   elseif sqrt((x.-6.7).^2+(y.-6).^2)<=0.28  wert=z2;
   elseif sqrt((x.-6.8).^2+(y.-4.1).^2)<=0.24 wert=z3;
   elseif sqrt((x.-2.7).^2+(y.-4).^2)<=0.2 wert=z4;
   elseif sqrt((x.-6.8).^2+(y.-3.5).^2)<=0.06 wert=z5;
   elseif sqrt((x.-4.1).^2+(y.-2.6).^2)<=0.1 wert=z6;
  else wert=0;
  endif 
endfunction

Fazit & weitere Optimierungsmöglichkeiten

Weiter Optimierungsmöglichkeit wäre das Hinzufügen weiterer Orte. Auch das Impfen als weiterer Faktor für das SIR-Modell könnte hinzugefügt werden. Nach dem Modellieren der Pandemie in den 6 Knoten ergeben sich viele weitere Modellierungsmöglichkeiten. Interessant wäre die Untersuchung und Modellierung der Lockdowns auf das Infektionsgeschehen, ebenso wie der Einfluss der Impfkampagne (angepasst an den realen Verlauf der Impfkampagne). Auch die Entstehung weiterer Varianten des Corona-Virus haben einen großen Einfluss auf den Fortlauf der Pandemie. Auch der Einfluss des Wetters bzw. der Jahreszeit auf die Verbreitung des Virus könnte in einem Modell erarbeitet werden.

Ausblick: Modellierung mit der Impfkampagne. Durch Erstellen einer Funktion, die die Impfquote zu verschiedenen Zeitpunkten abbildet, könnte die vorhergehende Modellierung in Octave angepasst werden.

Verbindung der Werte mit dem SIR-Modell

for p=1:(length(X))
  for q=1:(length(Y))
    Koordinaten_S(q,p)=1*hallo_funktion(x(q,p),y(q,p),Ergebnismatrix(1,1),Ergebnismatrix(2,1),Ergebnismatrix(3,1),Ergebnismatrix(4,1),Ergebnismatrix(5,1),Ergebnismatrix(6,1));
    Koordinaten_I(q,p)=1*hallo_funktion(x(q,p),y(q,p),Ergebnismatrix(1,2),Ergebnismatrix(2,2),Ergebnismatrix(3,2),Ergebnismatrix(4,2),Ergebnismatrix(5,2),Ergebnismatrix(6,2));
    Koordinaten_R(q,p)=1*hallo_funktion(x(q,p),y(q,p),Ergebnismatrix(1,3),Ergebnismatrix(2,3),Ergebnismatrix(3,3),Ergebnismatrix(4,3),Ergebnismatrix(5,3),Ergebnismatrix(6,3));
  endfor
endfor

Kontinuierliche Modellierung

Tutorium 1 - Kontaktmodell

In der ersten Tutoriumsaufgabe soll ein einfaches Kontaktmodell entwickelt werden. Dazu befindet sich unter einer gewissen Anzahl von Individuen ein infiziertes. Nun bewegen sich die alle Individuen. Kommen sie einem Infiziertem zu nahe, infizieren sie sich selbst und tragen so zu einer Verbreitung der Infektionen bei, denn im nächsten Zeitschritt können auch diese die Infektion weitergeben.

Folgende Aufgaben waren zu bearbeiten:

1. Geschwindigkeit, finale Simulationszeit ändern und analysieren, wie sich das auf die Infektionsverbreitung auswirkt.

2. Die Punkte sollen im Verlauf der Simulation die Richtung mehrfach ändern, der Code soll entsprechend angepasst werden.

Umsetzung der Aufgaben

Als erstes muss ein Mashgrid erzeugt werden, wo sich die Induvidien am Anfang der Simulation befinden. Ein Infizierter wird so gewählt, dass er sich im Zentrum aller Individuen befindet.

 Nx=12;% Anzahl der Vektorkomponenten 
 Ny=18;% Anzahl der Vektorkomponenten
 Zeitschritte=3;% pro T=1
 radiusInf=0.55;% Ansteckungsradius
 T=4 ;% Zeiteinheit 
 g=0.90;% Geschwindigkeit der Teilchen/Personen
 % Abstand am Anfang=1
 %------------------------

 x=[1:Nx];%Zeilenvektor mit Nx - Zeilen
 y=[1:Ny];%Spaltenvektor mit Ny - Spalten

 dt=1/Zeitschritte;% Unterteilung der Zeiteinheit in oben genannte Zeitschritte

[x,y]=meshgrid (x,y);% Aufspannung der Koordinatengitter von x und y 

 neuPosX=x.+g*rand(Ny,Nx)-0.5*g;#x-Koordinaten der Teilchen/Passanten
 neuPosY=y.+g*rand(Ny,Nx)-0.5*g;#y-Koordinaten der Teilchen/Passanten

 %erste Infizierte 
 xInf=x(1, Nx/2);%x-Position des ersten Infizierten im Mittelpunkt 
 yInf=y(Ny/2);%y-Position des ersten Infizierten im Mittelpunkt 
 IndInf2=Nx/2;
 IndInf1=Ny/2;
 IndInf1neu=IndInf1;
 IndInf2neu=IndInf2;

Die Laufrichtung der Induvidien wird in einer For-Schleife mit Zufallszahlen generiert. So bewegen sich alle Punkte im Laufe der Simulation in eine unterschiedliche Richtung.

Kritik: Ein Individuum kann nach einer Infektion nicht genesen, ist es einmal infiziert, bleibt es infiziert.

Tutorium 2 - Fundamentallösungen der Diffusionsgleichung

Physikalische Herleitung

Für die Herleitung relevante Funktionen und Größen

$x\in \mathbb {R} ^{n}$ ||Raumveränderliche

$t\in \mathbb {R} ^{+}$ ||Zeitveränderliche

$u:\mathbb {R} ^{n}\times \mathbb {R} ^{+}\to \mathbb {R}$ ||Dichtefunktion der Konzentration der Corona-Infizierten und nach Definition ist die Abbildung auch ein Vektorfeld auf $\mathbb {R} ^{+}$ ´siehe und wird im Fall der Diffusion als Strömungsfeld interpretiert.

$div$ ||Divergenz: Quellen und Senken eines Vektorfelds oder auch Quelldichte genannt. Besagt grundsätzlich wie sehr Vektoren in einer kleinen Umgebung auseinanderstreben.

$\nabla$ ||Gradient: Richtung des steilsten Anstiegs

$\Delta$ ||Laplace-Operator: Ordnet einem zweimal differenzierbaren Skalarfeld die Divergenz seines Gradienten zu. Ist für skalare Funktionen allgemein definiert als $\delta =div\nabla$

$\Delta u=0$ ||Laplace-Gleichung oder auch homogene Poisson-Gleichung.

$\Delta u=f$ ||Poisson-Gleichung

$q(x)$ ||Diffusionsflussrate oder auch Teilchenstromdichte, die angibt wie viele Teilchen links nach rechts strömen.

$c(x)$ ||Ortsabhängiger Diffusionkoeffizient = Maß der Beweglichkeit von Teilchen

Herleitung stationäre Diffusion in räumlicher Dimension

Für die Herleitung der räumlichen Diffusion wird ein dünner horizontaler Stab der Länge $L$ betrachtet. Dementsprechend wird mathematisch gesehen die Konzentration der Teilchen in einer räumlichen Dimension auf dem Intervall $[0,L]$ untersucht. Die Konzentration in vertikaler Richtung wird im System als homogen betrachtet, wodurch nur die Konzentrationsdichte in horizontaler Richtung von Bedeutung ist. Für die Herleitung der Diffusionsgleichung nutzt man nun zwei physikalische Gesetze. Das Fick´sche Gesetz und das Massenerhaltungsgesetz. Zunächst betrachten wir die Herleitung der Diffusionsgleichung in einem von der Zeit $t$ unabhängigen System. Also ein System, dass sich bereits im Gleichgewicht befindet. Dies nennt man auch stationäre Diffusion. Bei der Diffusion werden Teilchen bewegt, präziser gefasst werden sollen sich Teilchen des einen Stoffes in einem zweiten Stoff verteilen. Diese Bewegung ist nicht ungeordnet, sondern besitzt eine konkrete Richtung. Diese Richtung nennt man Transportrichtung. Die Teilchen strömen nach dem zweiten Fick'schen Gesetz vom Ort mit höherer Konzentration in Orte mit niedriger Konzentration. Daraus ergibt sich ein Fluss bzw. eine Größe für den Fluss der Teilchen. Diese Größe nennt man auch die Teilchenstromdichte $q(x)$. Nach dem ersten Fick'schen Gesetz ist diese Teilchenstromdichte proportional zum Konzentrationsgradienten $\nabla u(x)$, der entgegen der Diffusionsrichtung, also der Teilchentransportrichtung verläuft. Die Teilchenstromdichte setzt sich nach dem ersten Fick'schen Gesetz aus dem Konzentrationsgradienten und dem Diffusionskoeffizienten zusammen:

$q(x)=-c(x)\nabla u(x)=-c(x){\frac {\partial u}{\partial x}}=-c(x)u'(x)$

Das zweite bereits oben angesprochene Gesetz ist das Erhaltungsgesetz bzw. Massenerhaltungsgesetz. Dieses besagt grundsätzlich, dass der Stoff in einem Volumen nicht verschwinden kann, noch das zusätzlich Stoff entstehen kann (sofern keine externe Zugabe oder Abnahme erfolgt). Dieser Stoff bleibt also in diesem Volumen erhalten. Nach dem Massenerhaltungsgesetz kann sich die Teilchenanzahl gegeben durch $\int u(s)ds$ in unserem horizontalen Stab, also zwischen den Rändern des Stabes $x$ und $x+\Delta x$ nicht ändern. Durch die Brownsche Bewegung kommt es allerdings zum Austausch der Teilchen an den Rändern des Stababschnnittes $[x,x+\Delta x]$. Der Zufluss wird durch die Teilchenstromdichte am linken Rand $q(x)$ und der Ausfluss durch die Teilchenstromdichte am rechten Rand $q(x+\Delta x)$ gegeben. Aus dem Erhaltungsgesetz folgt nun, dass sich Zufluss und Ausfluss ausgleichen müssen, d.h.:

$q(x)-q(x+\Delta x)=0$

Durch Bilden des Differentialquotienten durch das Teilen durch $\Delta x$ und den Grenzübergang $\Delta x\to 0$ ergibt sich $-q'(x)=0$. Durch Anwendung des ersten Fick'schen Diffusionsgesetz, erhält man die Diffusionsgleichung $-q'(x)=[c(x)u'(x)]'=0$. Bei konstantem Diffusionskoeffizient erhält man nun die Laplace-Gleichung in einer räumlichen Dimension:

$u''(x)=0$

Verallgemeinerung in $\mathbb {R} ^{n}$

Die Gesammtzahl von Teilchen in einem beliebigen Volumen lässt sich auf folgenderweise berechnen: $V\subset \mathbb {R} ^{n}$: $\int _{V}u\,dx$

Der Massenfluss durch den Rand $\partial V$ ergibt sich für Volumen $V$:$\int _{\partial V}q\cdot {\vec {n}}\,dSx$, wobei ${\vec {n}}$ der äußeren Normalenvektor zum Rand von V und die Verknüpfung '$\cdot$' das Skalarprodukt bezeichnet.

Für das Fick'sche Gesetz in $\mathbb {R} ^{n}$ wird die Ableitung von $u$ nach $x$ wieder durch den Konzentrationsgradienten $\nabla u:=(\partial _{x1}u,...,\partial _{xn}u)$ ersetzt.

Für den Massenfluss ergibt sich aus obiger Rechnung ebenfalls die Erhaltungsgleichung $V$:$\int _{\partial V}q\cdot {\vec {n}}\,dSx=0$

Im $\mathbb {R} ^{n}$ wird nun anstelle der Bildung des Differentialquotients der Satz von Stokes und die Lemma Integralmittelwert angewendet. So ergibt sich folgendes:

$\int _{\partial V}q\cdot {\vec {n}}\,dSx=\int _{V}div\ q\ dx=0$

Durch einsetzen des ersten Fick'schen Gesetz bekommt man dann die stationäre (homogene) Diffusionsgleichung in $\mathbb {R} ^{n}$. Diese lautet

$div(-c(x)\nabla u)=0.$ 

Herleitung zeitabhängige Diffusion in räumlicher Dimension

Für die zeitabhängige Diffusion ist das System, welches betrachtet wird, in keinem Gleichgewichtszustand. Dementsprechend ist die gesuchte Dichtefunktion zeitabhängig. Es ergibt sich für die Dichtefunktion der Konzentration $u(x,t)$. Die Gesamtzahl der Teilchen im Stababschnitt mit der Länge $\Delta x$ wird nun als zeitabhängiges Integral angegeben ${\frac {d}{dt}}\int _{x}^{x+\Delta x}u(s,t)ds$ Das Erhaltungsgesetz besagt auch in diesem Fall mit einer zeitlichen Veränderungen, dass die zeitliche Veränderung der Gesamtanzahl der Teilchen durch den Austausch bzw. Zufluss und Ausfluss ausbalanciert sein muss. ${\frac {d}{dt}}\int _{x}^{x+\Delta x}u(s,t)ds=q(x,t)-q(x+\Delta x,t)$

Nach dem Teilen durch $\Delta x$ und dem Grenzübergang $\Delta x\to 0$ erhält man mithilfe des Lemma Integralmittelwert die Erhaltungsgleichung

$\partial _{t}u(x,t)=-\partial _{x}q(x,t).$

Durch erneute Anwendung des ersten Fick'schen Gesetz lässt sich die Teilchenstromdichte $q(x,t)$ ersetzen und man erhält folgende zeitabhängige (homogene) Diffusionsgleichung für die räumliche Dimension:

$\partial _{t}u(x,t)=\partial _{x}(c(x)\partial _{x}u(x,t)).$

Verallgemeinerung in $\mathbb {R} ^{n}$

Wird für die Ableitung nach $x$ der Gradient gesetzt. Wendet man erneut den Satz von Stokes bzw. den Gaußschen Intergalsatz erneut auf den zeitabhängigen Massenfluss an ergibt sich insgesamt folgende (homogene) Diffusionsgleichung:

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

Herleitung Diffusion mit Quellterm

Bei der Diffusion mit Quellterm wird dem System eine externe oder interne Teilchenquelle hinzugefügt. Somit entsteht im System neue Masse. Diese Massenquelle trägt ebenfalls zum Massenerhaltungsgesetz bei und soll mathematisch durch die Dichtefunktion der Materialflussrate $f(x,t)$ beschrieben werden. Diese kennzeichnet die hinzugefügte oder entstandene Masse pro Volumen. Der Quellterm also $f(x,t)$ trägt ebenfalls zur zeitlichen Veränderung der Gesamtzahl von Teilchen im Intervallabschnitt bei, ebenso wie es auch der Austausch an den Rändern tut. Es ergibt sich also folgende Gleichung:

${\frac {d}{dt}}\int _{x}^{x+\Delta x}u(s,t)ds=q(x,t)-q(x+\Delta x,t)+\int _{x}^{x+\Delta x}f(s,t)ds.$

Diese lässt sich nach bekanntem Schema der zeitabhängigen Diffusion über das Teilen durch $\Delta$, den Grenzübergang von $\delta x\to 0$ in folgende (inhomogene) Diffusionsgleichung in der räumlichen Dimension umformen:

$\partial _{t}u(x,t)-\partial _{x}(c(x)\partial _{x}u(x,t))=f(x,t),\ x\in [0,L],t\in \mathbb {R} ^{+}.$

Verallgemeinerung in $\mathbb {R} ^{n}$

Und durch das verallgemeinern in $\mathbb {R} ^{n}$ ergibt sich über das Ersetzen der Ableitung nach $x$ von $u(x,t)$ und die Anwendung des Satz von Stokes folgende inhomogene Diffusionsgleichung:

$\partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=f(x,t),\ x\in \ \mathbb {R} ^{n},\ t\in \mathbb {R} ^{+}.$

Herleitung Diffusion mit Reaktionsterm

Die Diffusion mit Reaktionsterm oder auch die Reaktionsdiffusion wird durch die Reaktionsdiffusionsgleichung beschrieben. Diese ist eine inhomogene zeitabhängige Diffusionsgleichung mit einer Quelldichtefunktion auf der rechten Seite, welche von der Stoffkonzentration abhängig ist. Ähnlich wie bei der Diffusion mit Quellterm wird auch hier auf der rechten Seite der Gleichung eine Quelldichte $f$ hinzugefügt. Der Unterschied liegt jedoch darin, dass der Quellterm von der Dichtefunktion der Stoffkonzentration abhängig ist. Dementsprechend wird die Diffusionsgleichung nicht durch $f(x,t)$ sondern durch $f(u)$ oder auch $f(u(x,t))$ ergänzt. Diese partiellen Differentialgleichungen werden zum Beispiel auch in der Populationsdynamik genutzt.

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

Der Reaktionsterm hat eine bzw. die logistische Form und ist dementsprechend beschränkt. In unserem Fall orientiert sich der Reaktionsterm an Fishers-Gleichung der Populationsdynamik. Diese ist eine Sonderform der Kolmogorov-Petrovsky-Piscounov-Gleichung (KPP-Gleichung) ähnlich der logistischen Funktion, welche den Zusammenhang zwischen vertreichender Zeit und Wachstum beschreibt und durch eine obere Schranke beschränkt wird. In unserem Fall wäre dies die Bevölkerungdichte bzw. die Gesamtanzahl von Personen an einem festen Ort $x$. Der Reaktionsterm hat also folgende Form:

$f(u(x,t))=k.u(x,t){\big (}B(x)-u(x,t){\big )}$

Woraus sich auch die folgende Reaktionsdiffusionsgleichung ergibt:

$\partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=k.u(x,t){\big (}B(x)-u(x,t){\big )},\ x\in \ \mathbb {R} ^{n},\ t\in \mathbb {R} ^{+}.$

Analytische Lösung

Stationäre Diffusion: Laplace Gleichung

Die sogennanten Fundamentallösung der Laplace Gleichung $\Delta u(x)=0$ auf $\mathbb {R} ^{n}$ (Cauchy-Problem) lässt sich anhand der Rotationsinvarianz des Laplace-Operators finden, da sich der selbst durch das Tauschen der Variablen x nicht verändert. Somit ergibt sich die rotationssymmetrische Lösung der Laplace-Gleichung. Formel der Fundamentallösung für $\mathbb {R} ^{n},n=2,\ 3:$

 $u(x)=\Phi (x):=\left\{{\begin{array}{l}-{\frac {1}{2\pi }}\ln {\|x\|}\ ,&n=2\\{\frac {1}{\omega _{n}}}{\frac {1}{\|x\|^{n-2}}}={\frac {3}{4\pi }}{\frac {1}{\|x\|^{n-2}}}\ ,&n\geq 3\\\end{array}}\right.$

wobei $\omega _{n}$ den Volumen der Einheitskugel in $\mathbb {R} ^{n}$ beschreibt und $\|\cdot \|$ die Vektornorm in $\mathbb {R} ^{n}$ ist, vergleiche .

Für n=3 ist $\omega _{3}={\frac {4}{3}}\pi$.

Sowohl für $n=2$ aber auch $n\geq 3$ gilt, dass die Fundamentallösung singulär ist.

 $\lim \limits _{x\to 0}$$\Phi (x)=-\infty$

Stationäre Diffusion: Poisson Gleichung

Ist die Fundamentallösung für einen Differentialoperator bereits gegeben, so erhält man die Lösung der Gleichung $u(x)$ bzw. der der Poissongleichung $-\Delta u(x)=f(x)$ auf $\mathbb {R} ^{n}$ $\forall {x}\in \mathbb {R^{n}}$ durch Faltung der Fundamentallösung mit der Funktion der rechten Seite. (Poissonformel) :

$u(x)=(\Phi *f)(x)=\int _{\mathbb {R} ^{n}}\Phi (x-y)f(y)dy.$ 

Für eine Funktion $f\in C^{2}(\mathbb {R} ^{n})$ mit kompaktem Träger, siehe Definition Kompaktheit gilt,

$\lim _{\|x\|\to \infty }u(x)\to 0.$ 

(für Beweis siehe )

Zeitabhängige Diffusion

Für die homogene (instationäre) Diffusionsgeichung existieren spezielle Lösungsformel für konstanten Diffusionskoeffizient $c(x)\equiv a\in \mathbb {R}$.

Die Fundamentallösung für $\partial _{t}u-a\Delta u=0$, $x\in {\mathbb {R} }^{n},\ t\in {\mathbb {R} }^{+}$ lautet

$u({x},t)=\Psi ({x},t)={\frac {1}{(4\pi at)^{n/2}}}\exp \left(-{\frac {\|{x}\|^{2}}{4at}}\right),$

wobei $\textstyle \|{x}\|^{2}=\sum _{k=1}^{n}x_{k}^{2}$ das Quadrat der euklidischen Norm von ${x}$ ist, siehe .

Die Diffusiongleichung ergänzt durch eine Anfangsbedingung :

$u({x},t=0)=u_{0}({x})\ \forall \ {x}\in \mathbb {R} ^{n}:$

die im Anfangszeitpunkt $t=0$ die Konzentrationsdichte /Dichte der Infizierten beschreibt heißt Anfangswertproblem.

Die Lösung $u({x},t)$ des (homogenen) Anfangswertproblem für $t>0$ ist durch die Faltung der Fundamentallösung $\Psi$ mit der Anfangsfunktion gegeben $u_{0}$:

$u({x},t)=(\Psi *u_{0})({x},t)=\int _{\mathbb {R} ^{n}}\Psi ({x}-{y},t)u_{0}({y})\,d{y}$

siehe . Siehe auch Lösungsformel für das inhomogene Cauchyproblem.

$\lim _{({x},t)\to (0,0)}\Psi (x)=\infty$

Tutoriumsaufgaben

Implementierung der Funktion Φ in Octave

$u(x)=\Phi (x):=\left\{{\begin{array}{l}-{\frac {1}{2\pi }}\ln {\|x\|}\ ,&n=2\\{\frac {1}{\omega _{n}}}{\frac {1}{\|x\|^{n-2}}}={\frac {3}{4\pi }}{\frac {1}{\|x\|^{n-2}}}\ ,&n\geq 3\\\end{array}}\right.$

Grafische Darstellung der Funktion Ψ für n=2

Grafische Darstellung der Funktion Ψ für n≥3

for n=2:3
    if n=2
      x=[-400:1:400];
      phi=(-1/(2*pi))*log(sqrt((x).^2))
      figure(1)
      plot(x, phi,'-')
      axis=('auto')
      grid on 
      xlabel('x-Achse')
      ylabel('y-Achse')
      title("Phi(x)für n=2")
      saveas( 1, "Phix.jpg")
    endif
    if n=3
     for x=-10:0.1:10;
       phi=(3/(4*pi))*(1/(sqrt((x)^2))^(n-2))
       figure (2)
       hold on 
       plot (x, phi,'*b')
       grid on 
       xlabel('x-Achse')
       ylabel('y-Achse')
       title("Phi(x) für n>=2")
       saveas(2, "Phi2x.jpg")
     endfor
    endif
 endfor

Implementierung der Funktion Ψ in Octave

$u(x,t)=\Psi (x,t)={\frac {1}{(4\pi at)^{\frac {n}{2}}}}\exp \left(-{\frac {\|x\|^{2}}{4at}}\right)$

Grafische Darstellung der Funktion Φ

 for t = 1:5;
   x=-10:10;
   psi=(1/(4*0.02*t))*exp(-sqrt((x).^2)/(4*0.02*t))
   hf= figure(t) 
   plot(x,psi)
   title('Grafische Darstellung Psi(x,t)')
   xlabel('x-Achse')
   ylabel('y-Achse')
   axis([-10 10 0 13])
   test=["Figure_", num2str(t),".jpg"]
   saveas(t, test);
 endfor

Zeitabhängige Diffusion als Anfangswertproblem

Bei der zeitabhängigen Diffusion als Anfangswertproblem wird die Diffusionsgleichung durch eine Anfangsbedingung ergänzt, die die Konzentrationsdichte/ Dichte der Infizierten zum Anfangszeitpunkt $t=0$ beschreibt.

Anfangswertfunktion

Anfangswertfunktion u0

Zum Beginn haben wir eine abwechslungsreiche Angfangswertfunktion implementiert

function wert=Anfangswertfunktion(x,y)
  if sqrt((x.-0.0).^2+(y.-2.5).^2)<=0.5 wert=0.5;
  elseif sqrt((x.-2).^2+(y.-0.0).^2)<=0.5 wert=1;
  %elseif sqrt((x.-0.8).^2+(y.-1).^2)<=0.15 wert=3;
  else wert=0;
  endif
endfunction

Implementierung der Formel Psi

Diffusion der Anfangswertfunktion zum Zeitpunkt t = 0.2

Animation der Diffusion über das Zeitintervall [0.2, 1] mit $\Delta t$ = 0.05

   for i=1:(length(X)) 
       for j=1:(length(Y))
       xstar=x(j,i)*ones(length(Y),length(X));
       ystar=y(j,i)*ones(length(Y),length(X));
       % Psi
       psi=(1/(4*pi*a*(t0+k*(1/Zeitschritte))))*(exp(-((xstar.-x).^2+(ystar.-y).^2))/(4*a*(t0+k*(1/Zeitschritte))));
       endfor
   endfor

Integration

       %numerische Integration, alles funktioniert
       I(j,i) = trapz(Y,trapz(X,psi.*u0),2);

Nach dem schon die Integration vonstatten ging, wird nun mit Hilfe der folgenden Schleife die Grafik ausgegeben

%Graphische Ausgaben  
Schritte = k*Zeitschritte+1;
figure (k+1)
meshc(x,y,I)
grid on
title (['Loesung von Poissongleichung zum Zeitpunkt t = ', num2str(k*(1/Zeitschritte))])
xlabel('x - Achse')
ylabel ('y - Achse')
zlabel('z - Achse')
zlim([0 9])

Zum Zeitpunkt t = 0 wird eine Singularität erreicht, wodurch das erste Bild in unserer Animation keine Bedeutung zuzuweisen ist, da die Werte viel zu groß werden.

Tutorium 3 - Modifiziertes SIR Modell

SIR Modell

Für die Bearbeitung der dritten Aufgabe im Rahmen der Tutorien rückten die Modelle Populationsdynamik, präziser gefasst das modifzierte SIR Modell in den Fokus der Vorlesungsinhalte. Ähnlich wie das logistische (SI) Modell gehört auch das SIR Modell zu den Kompartimentmodellen für die Infektionsverbreitung. Das Modell selbst teilt die Gesamtpopulation in drei Gruppen bzw. Kompartimente ein: Susceptibles (Gruppe der Empfänglichen), Infected (Gruppe der Infizierten) und Removed (Gruppe der Genesenen) so wird die Diffusion als Austausch zwischen den Kompartimenten betrachtet.

Die Änderungsraten der jeweiligen Kompartimente werden mittels eines Systems von drei gekoppelten Differentialgleichungen berechnet:

$S'(t)=-kI(t)S(t),$

$I'(t)=kI(t)S(t)-wI(t),$

$R'(t)=wI(t),$

wobei die $w$ die Regenerationsrate bzw. Wechsel- oder Genesensrate und $k$ die Infektionsrate ist. Die Wechselrate selbst ergibt sich durch den Umkehrwert der Genesungszeit $T$. Daraus folgt letztendlich: $w=1/T$. Für COVID 19 ergibt sich somit eine Wechselrate von $w=1/14$, da laut Studien die Inkubationszeit zwischen 2 und 14 Tagen beträgt und 95% aller Befragten auch in diesem Zeitraum Symptome der Viruserkrankung entwickelt haben. Beachtet man nun die Tatsache, dass die Anzahl der Personen in den jeweiligen Gruppen die Gesamtpopulation nicht übersteigen kann, ergibt sich folgende Erhaltungsgleichung:

${\Big (}S+{I}+{R}{\Big )}'=0,$

die im grundlegenden nur besagt, dass die Gesamtpopulation $S+{I}+{R}$ konstant erhalten bleibt, dementsprechend findet kein räumlicher Austausch statt, weswegen das Modell auf eine feste Ortskoordinate $x$ angewendet wird und es sich somit um ein rein dynamisches, also zeitabhängiges Modell handelt.

Modifiziertes SIR Modell

Das modifizierte SIR Modell wird angewendet, wenn die Toten nicht im Kompartiment der Genesenen erfasst werden sollen. Hinzukommend geht man im modifizierten SIR Modell davon aus, dass nur ein gewisser Prozentsatz der Infizierten auch Symptome entwickelt und somit von den Gesundheitsämtern erfasst wird. Unter Betrachtung dieser Modifikationen muss das Modell angepasst werden. So wird das Modell zusätzlich um eine konstante Sterberate $d$, ein konstanter Prozentsatz der durch Tests erfassten Infizierten $r$ und durch eine von $k$ und $B$ abhängige Wachstumsrate $c$ mit $c=kB$ ergänzt. Somit ergibt sich folgendes System von ebenfalls drei gekoppelten Differentialgleichungen:

${\textstyle S'(t)=-{\frac {c}{rB}}S(t)I(t)}$

${\textstyle I'(t)={\frac {c}{B}}S(t)I(T)-wI(t)}$

${\textstyle R'(t)=wI(t)-dI(t).}$

Hinzukommt die Anfangsbedingung der jeweiligen Differentialgleichung:

$\quad S(t_{0})=B-I_{0},I(t_{0})=I_{0},R(t_{0})=0,$

Durch die oben getroffene Anpassung des Modells verändert sich auch die Erhaltungsgleichung aus dem SIR Modell folgendermaßen:

${\Big (}S+{\frac {I}{r}}+{\frac {R}{r}}{\Big )}'=-dI,$

Diese besagt letztendlich nur, dass die Gesamtpopulation ${\Big (}S+{\frac {I}{r}}+{\frac {R}{r}}{\Big )}$ um die verstorbenen Infizierten verringert.

Implementierung des modifizierten SIR-Modells in Octave

 B=55000; % Kapazitätsgrenze, 2/3 der deutschen Bevölkerung (In T) 
 c=0.3545; %Basisinfektionsrate- die Rate des exponentiellen Wachstums am Anfang der Pandemie
 w=1/14; % Wechselrate zu den Genesenen
 d=0.003; % Todesrate
 r=0.1;%Erfassungsfaktor, der besagt, dass 10 Prozent aller Infektionen entdeckt werden
 yo=[B-16/1000;16/1000;0]%Anfangsbedingungen

 %SIR-Modifiziert ohne Slowdown 
 f=@(y,x) [-c*y(1)*y(2)/(r*B);c*y(1)*y(2)/B-w*y(2);w*y(2)-d*y(2)];System von 3 gekoppelten Differentialgleichungen 
 y = lsode (f, yo, timesWHO) %y(1)=S(t) & y(2)=I(t)

 figure(13)%Grafische Darstellung des modifizierten SIR Modells 
 plot (timesWHO, y(:,1), timesWHO, y (:,2), '-.', timesWHO, y (:,3),'-.', timesWHO, y (:,3)+y(:,2))%Darstellung der Lösungen des DGL 
 hold on
 title("SIR modifiziert ohne SLOWDOWN")
 legend ("Susceptible", "Infected", "Removed", " Infected+Removed")%,"German data I+R" ,"German data I ","German data R" ,"location", "west")
 axis([0,n, 0 ,4000])

Aufgabe 1

Erweiterung der Funktion "coronaData"

Die erste Aufgabe war es nun die jeweiligen Zahlen des Infektionsgeschehens in die bereits von Frau Hundertmark vorgefertigte Funktion "coronaData" zu implementieren. Dazu sollten wir die jeweiligen Infektionszahlen in dieser Funktion ergänzen. Die jeweiligen Werte konnte man sehr gut im Internet wiederfinden und so werden bereits auf Wikipedia Tabellen und Datenbanken gepflegt, in welche die neuesten Zahlen direkt eingetragen und aktualisiert werden. Siehe [[1]]

 %Implementierung der gesammelten Corona-Daten 
 A=coronaData_Mai2021();
 inf_falleWHO=A(1,:);%Kumulierten Infektionen 
 inf_falleWHOrecovered=A(3,:);%Kumulierten Genesenen 
 inf_falleWHOtoten=A(2,:);%Kumulierte Gestorbene
 inf_falleWHOaktuell=inf_falleWHO-inf_falleWHOrecovered;%Tägliche Neuinfektionen 
 n=length(inf_falleWHO)
 timesWHO=[0:1:n-1];

 figure(1)
 plot(A(1,:),'-r',A(2,:),'-b', A(3,:),'-g' )
 axis([0 456 0 4000000])
 xlabel ('Tage')
 ylabel ('Zahlen')
 legend('Inizierte','Tote','Genesene')
 saveas(figure(1),"Corona_Daten.jpg")

Anpassung der Infektionsrate mittels einer Slowdown-Funktion

 function val=SIR_slowdown_verbessert(x)
 % Ab Tag t (Argument der Funktion) wird die Infektionsrate in mehrere Schritten sinken:
 if x < 20 
   val=1;
   else if x < 43
     val=0.6;
   else if x < 100%24.05.2020 Infektionsgeschehen hat sich wieder beruhigt
     val=0.095;
   else if x < 220%30.09.2020 die zweite Welle beginnt 
     val=0.0025*(x-100)+0.095;
   else if x < 276
     val=0.27;
   else if x < 310
     val=0.25;
   else if x < 355  
     val=0.15;
   else if x < 385
     val=0.23;
   else if x < 420
     val=0.26;
   else if x < 455 
     val=0.19;
 else  
   val=0.75;
   endif
   endif
   endif 
   endif 
   endif
   endif
   endif
   endif
   endif
 endif
 endfunction

Nachdem wir gemerkt haben, dass die Anpassung der Slowdown-Funktion an das damalige Infektionsgeschehen und nicht nur an die Infektionsrate erfolgen muss, haben wir uns anhand der gesammelten Daten eine neue Slowdown-Funktion gebastelt. Die Ergebnisse sehen sie auf der rechten Seite.

Implementierung der Slowdown-Funktion in Octave

 %Slowdown-Funktion
 for i=1:(length(timesWHO))
   c_funktion(i)=SIR_slowdown_verbessert (timesWHO (i), 0);
 endfor

 T=1;
 c_real(1)=c_funktion(1);
 for i=T+1:n
   c_real(i-T)=(inf_falleWHO(i)-inf_falleWHO(i-T))/(inf_falleWHO(i-T));
 endfor
 c_real_durchschnitt=c_real;
 for i=4:n-4
   c_real_durchschnitt(i)=(c_real(i-3)+ c_real(i-2)+c_real(i-1)+c_real(i)+ c_real(i+1) + c_real(i+2)+c_real(i+3) )/7;
 endfor 

 figure(111)
 plot (timesWHO(1: n-1),c_funktion(1: n-1), timesWHO(1:n-1),c_real_durchschnitt(1:n-1))%,timesWHO(1: n-1), c_real_durchschnitt),timesWHO(1: n-1),c_beschraenkt_durchschnitt)
 grid on
 title ('Tatsächtliche Infektionsrate als 7-Tage Duchschnitt  und slow-down-Funktion')
 legend ('slow-down Funktion',  '7-Tage Durchschnitt', '7-Tage Durchschnitt beschr. Modell')

Bestimmung des relativen Fehlers

 %Berechnung der relativen Fehler und absoluten Fehler 
 aF1=y(:,2)'-(inf_falleWHOaktuell/1000)%Absoluter Fehler der aktuell Infizierten
 aF2=y (:,3)'-(inf_falleWHOrecovered/1000)%Absoluter Fehler der Genesenen 
 aF3=(y (:,3)'+y(:,2)')-(inf_falleWHO./1000)%Absoluter Fehler der kumulierten Infizierten + Genesenen 
 rF1=aF1./(inf_falleWHOaktuell/1000)%Relativer Fehler der aktuell Infizierten
 rF2=aF2./(inf_falleWHOrecovered/1000)%Relativer Fehler der Genesenen 
 rF3=aF3./(inf_falleWHO/1000)%Relativer Fehler der kumulierten Infizierten + Genesenen 

 figure(16)
 plot(timesWHO, rF1*100,'.-.r',timesWHO, rF2*100,'.-.b',timesWHO,rF3*100,'.-.g')
 title('Relativer Fehler Modell und Daten')
 xlabel('Tage ab 27. januar')
 ylabel('Abweichung in Prozent')
 legend('Aktuelle Corona-Infektionen','Genesene','Kumulierte Infektionen und Genesene')

 figure(17) 
 plot(timesWHO, aF1,'.-.r',timesWHO,aF2,'.-.b',timesWHO,aF3,'.-.g')
 title('Absoluter Fehler Modell und Daten')
 xlabel('Tage ab 27. januar')
 ylabel('Abweichung in absoluten Zahlen')
 legend('Aktuelle Corona-Infektionen','Genesene','Kumulierte Infektionen und Genesene')

Anpassung des Erfassungsfaktors durch eine Erfassungsfunktion

Zur Anpassung Erfassungsfunktion haben wir uns zunächst die Zahl der jeweiligen Testungen zur Hand genommen. Dazu haben wir uns auf die Zahlen des RKIs berufen. Diese geben für jede Kalenderwoche die Zahl der jeweils durchgeführten Tests aus und bieten somit die Möglichkeit die Testungen und die Testverfahren für die Bevölkerung transparenter zu gestalten. Grundannahme unserer Anpassung war, dass mit steigenden Testzahlen auch mehr Corona-Infektionen erkannt werden. Dabei haben wir angenommen, dass der Anteil der Infektionen, die vom Gesundheitssystem registriert werden, linear zur Anzahl der Tests steigt. Die Daten des RKIs zeigen, dass die Zahlen bis in Kalenderwoche 17/2020 auf einen durchschnittlichen Wert von ungefähr 300.000 Tests pro Woche angestiegen sind. Wir haben angenommen, dass bei 300.000 Testungen pro Woche 10% der wirklichen Infektionen auch durch das Gesundheitssystem erfasst werden. Aufbauend auf dieser Zahl haben wir unsere Funktion in Orientierung an den RKI-Zahlen angepasst. Dazu haben wir die Tage zwischen den Kalenderwochen gezählt und über den Differenzenquotient ${\frac {y_{1}-y_{0}}{x_{1}-x_{0}}}$, wobei $x$ also den Tagen ab dem 27. Januar jeweils genau ein $y-Wert$ (Erfasste Infektionen in Prozent) zugeordnet wird.

function val=SIR_erfassungsfunktion(x,t)

 % t gibt die Tage an, an denen die Erfassung der Infizierten besser wird.
 % Durch zum Beispiel Tests
 % Annahme bei 300.000 Tests Erfassung von 10% der Fälle 
   if x < t + 147%Signifikant steigende Zahlen der landesweiten Testungen in KW27/2020
     val=0.1 ;
     else if x < t + 196%Erstmalige Zahlen über 1 Mio. Testungen pro Woche 
       val=0.0014*(x-147)+0.1;
     else if x < t + 322%Maximale Testungen in KW 50/2020  
       val= 0.0018*(x-196)+0.3;
     else if x < t + 336%Testungen sinken über Weihnachten und Silvester
       val= -0.02*(x-322)+0.53;
     else if x < t + 455%Testungen steigen auf Mittel von 1.2 Mio pro Woche also ungefähr 40% 
       val= 0.0013*(x-336)+0.25;
   else 
     val= 0.4;
     endif 
     endif 
     endif 
     endif 
   endif 
 endfunction

Implementierung der Erfassungsfunktion

 %Erfassungs-Funktion
 for i=1:length(timesWHO)
   rr(i)=SIR_erfassungsfunktion (timesWHO (i),0);% r 
 endfor

 %Modell mit Erfassungsfunktion
 r=SIR_erfassungsfunktion(timesWHO,0);
 g=@(y,x) [-c*SIR_slowdown_verbessert(x)*y(1)*y(2)/(r*B);c*SIR_slowdown_verbessert(x)*y(1)*y(2)/B-w*y(2);w*y(2)-d*y(2)];
 y = lsode (g, yo, timesWHO);%y(1)= S(t) y(2) = I(t)

 figure(18)
 plot (timesWHO, y (:,2), '-.r', timesWHO, y (:,3),'-.b', timesWHO, y (:,3)+y(:,2),'-.g')
 hold on
 plot (timesWHO, inf_falleWHO./1000, 'sg', "linewidth",4,timesWHO,inf_falleWHOaktuell/1000, '*r' , timesWHO,inf_falleWHOrecovered/1000, '*b' )
 legend ( "Infected", "Removed", " Infected+Removed","German data I+R" ,"German data I ","German data R" ,"location", "west")
 axis([0,n, 0 ,3600])
 title('Vergleich Modell mit Erfassungsfunktion und CoronaDaten ')
 hold off

Bestimmung des relativen Fehlers

 %Berechnung relativer Fehler und absoluter Fehler 
 AF1=y(:,2)'-(inf_falleWHOaktuell/1000)
 AF2=y (:,3)'-(inf_falleWHOrecovered/1000)
 AF3=(y (:,3)'+y(:,2)')-(inf_falleWHO/1000)
 RF1=AF1./(inf_falleWHOaktuell/1000)
 RF2=AF2./(inf_falleWHOrecovered/1000)
 RF3=AF3./(inf_falleWHO/1000)
 
 %Grafische Darstellung 
 figure(19)
 plot(timesWHO, RF1*100,'.-.r',timesWHO, RF2*100,'.-.b',timesWHO,RF3*100,'.-.g')
 title('Relativer Fehler Modell mit Erfassungsfunktion und Daten')
 xlabel('Tage ab 27. januar')
 ylabel('Abweichung in Prozent')
 legend('Aktuelle Corona-Infektionen','Genesene','Kumulierte Infektionen und Genesene')

 figure(20) 
 plot(timesWHO, AF1,'.-.r',timesWHO,AF2,'.-.b',timesWHO,AF3,'.-.g')
 title('Absoluter Fehler Modell mit Erfassungsfunktion und Daten')
 xlabel('Tage ab 27. januar')
 ylabel('Abweichung in absoluten Zahlen')
 legend('Aktuelle Corona-Infektionen','Genesene','Kumulierte Infektionen und Genesene')

Nach der Anpassung der Erfassungsfunktion folgte natürlich auch das Implementieren dieser in unser Modell des modifizierten SIR Modells. Dazu wurden erstmal die Werte der Funktion implementiert und das System der Differentialgleichungen dementsprechend verändert.

Tutorium 4: Finite-Differenzen-Methode (Dirichlet- und Neumann-Randbedingung)

Randwertprobleme

Die genutzten elliptischen Differentialoperatoren wie der Laplace-Operator führen immer auf Randwertprobleme. Dementsprechend werden diese im weiteren Verfahren gelöst. Das Randwertproblem für Diffusionsgleichungen beschreibt ein Diffusionsproblem auf einem beschränkten Gebiet D. Über eine Randwertbedingung wird das Verhalten der approximierten Funktion $u(x,t)$ am Rand des Gebietes festgelegt. Nun gibt es verschiedene Typen von Randwertproblemen für partielle Differentialgleichungen. Für unsere Programmierung in Octave werden das Dirichlet- und das Neumann-Problem betrachtet.

Dirichlet-Problem

Beim Dirichlet-Problem werden die Randwertbedingungen als Funktionswerte vorgegeben.

$-\Delta u(x)=f(x)\qquad x\in D$
$\displaystyle \quad u(x)=g(x)\quad \ \ x\in \partial D.$

Neumann-Problem

Beim Neumann-Problem werden anstatt Funktionswerten Ableitungswerte vorgeschrieben. Präziser gefasst, erfüllt die Funktion $u(x,t)$ die Neumann-Randbedingungen wenn die Richtungsableitungen von $u$ in der Richtung des äußeren Normalenvektors gegeben sind.

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

Finite Differenzen: Grundidee

Aufgabe im vierten Tutorium der Veranstaltung war die Implementierung der Finiten-Differenzen-Methode. Einer der wohl bekanntesten numerischen Diskretisierungsverfahren für partielle Differentialgleichungen. Dabei handelt es sich um ein gitter-basiertes Verfahren, das eine unbekannte Funktion $u(x,t)$ auf endlich vielen Gitterpunkten $P_{ij}=(x_{i},y_{i}),i=1,...,N_{x},j=1,...N_{y}$ approximiert. Dazu werden die Differentialoperatoren der Differentialgleichung durch Punktauswertung des Differenzenquotienten approximiert, wodurch aus (nicht-)linearen Differentialgleichungen Systeme (nicht-) linearer Gleichungen werden. Die Differenzenquotienten ergeben sich als Approximation der Ableitungen von $f$ nach dem Abbruch der Taylorreihe bzw. aus der Taylorformel mit Lagrangeschem Restglied

Taylorreihe

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

Taylorformel mit lagrangeschem Restglied

Sei Funktion $f$ hinreichend oft differenzierbar.

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

Erste Differenzenquotienten

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

Alle drei, also der vorwärts, rückwärts und zentrale Differenzenquotient 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]$.

Zweite Differenzenquotienten

Der zweite Differenzenquotient approximiert 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!}}.$

Beispiel eines Dirichlet-Problems in einer Raumdimension

Grundidee ist also die Differentialoperatoren durch finite Differenzen zu approximieren. Dazu suchen wir eine Funktion $u:[a,b]\to \mathbb {R}$, die für einen Quellterm $f:(a,b)\to \mathbb {R}$ sowohl die Differentialgleichung $u''=f$ (auf dem gesamten Rechengebiet), als auch die Randwertbedingungen $u(a)=u(b)=u_{0}$ erfüllt. Gegeben ist also folgendes Randwertproblem:

$-u''(x)=f(x)$ für $x\in [a,b]$

$u(a)=u(b)=u_{0}$

Zur Lösung mit der Finiten Differenzen Methode wird das Intervall $[a,b]$ durch die Gitterpunkte $x_{}=i\cdot h,i=0,...,N+1$ mit der Gitterweite $h={\frac {b-a}{N+1}}$ diskretisiert. Die zweite Ableitung wird nun durch den oben gezeigten zentralen Differenzenquotient der zweiten Ableitung 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})$.

Dies ergibt an den inneren Gitterpunkten die Differenzengleichungen

$-{\frac {1}{h^{2}}}(-u_{i-1}+2{u_{i}}-u_{i+1})=f$ für $i=1,...,N$

für die numerischen Näherungswerte $u_{i}$ der Lösungswerte $u(x_{i})$. Unter Verwendung der gegebenen Randwerte $u(a)=u(b)=u_{0}$ entsteht ein lineares Gleichungssystem mit $N$ Gleichungen für die $N$ Unbekannten $u_{1},...,u_{N}$.

$A_{h}\cdot {\vec {u}}=f_{h}$ mit der tridiagonalen Systemmatrix $A_{h}$

In Matrixform ergibt sich also folgendes zu lösendes System:

$-{\frac {1}{h^{2}}}{\begin{pmatrix}2&-1&0&\dots &0\\-1&2&-1&\ddots &\vdots \\0&-1&\ddots &\ddots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\dots &0&-1&2\end{pmatrix}}{\begin{pmatrix}u_{1}\\u_{2}\\u_{3}\\\vdots \\u_{n}\end{pmatrix}}={\begin{pmatrix}f+{\frac {u_{0}}{h^{2}}}\\f\\f\\\vdots \\f+{\frac {u_{0}}{h^{2}}}\end{pmatrix}}$

Beispiel eines Neumann-Randwertproblem in einer Raumdimension

Wie oben bereits beschrieben, werden beim Neumann-Randwertproblem die Richtungsableitungen der Funktion $u$ in Richtung der äußeren Normalenvektoren gegeben. Seien die Randbedingungen in den Normalrichtungen $-1,1$ gegeben:

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

Um die Ableitungen in den Randpunkten zu approximieren führen wir neue Unbekannte Funktionswerte in den Randpunkten $a$ und $b$ mit $u_{0}:=u(a),\ u_{n+1}:=u(b)$ ein. Das Gitter geht ebenfalls wie bereits im LGS beim Dirichlet-Randwertproblem von $u_{1},...,u_{n}$. Nun bieten sich zwei Approximationsmöglichkeiten. Die Approximation erster Ordnung über die einseitigen Differenzenquotienten, also den Vorwärts- und Rückwärtsdifferenzenquotient oder die Approximation zweiter Ordnung über den ersten zentralen Differenzenquotient. Gerade die zweite Möglichkeit bietet wohl genauere Lösungen.

Approximation erster Ordnung

Ersetzt man die Ableitung der Funktion $u$ am Punkt $a$ durch den Vorwärtsdifferenzenquotient, erlauben wir uns einen Fehler ${\mathcal {O}}(h)$, der sich durch eine Konstante mal den Wert h nach oben abschätzen lässt.

* den vorwärts -Differenzenquotient $D_{h}^{+}u(a)={\frac {u_{1}-u_{0}}{h}}=u'(a)+{\mathcal {O}}(h)$

Ersetzt man die Ableitung der Funktion $u$ am Punkt $b$ durch den Rückwärtsdifferenzenquotient, erlauben wir uns einen Fehler ${\mathcal {O}}(h)$, der sich durch eine Konstante mal den Wert h nach oben abschätzen lässt.

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

So erhalten wir letztendlich zwei neue Gleichungen. Die Gleichung für $u_{0}$ und $u_{n+1}$, also die 0-te und $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})$

Systemmatrix

So ergibt sich in Matrixform folgende Systemmatrix $\quad A_{h_{1}}\in {\mathbb {R} }^{(n+2)\times (n+2)}$ mit $A_{h_{1}}=-{\frac {1}{h^{2}}}{\begin{pmatrix}h&-h&0&\ldots &0\\-1&2&-1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&-1&2&-1\\0&\ldots &0&-h&h\end{pmatrix}}$ und Vektor u ${\vec {u}}={\begin{pmatrix}u_{0}\\u_{1}\\u_{2}\\\vdots \\u_{n+1}\end{pmatrix}}$.

Lineares Gleichungssystem für Neumann-Problem mit Approximation erster Ordnung

$-{\frac {1}{h^{2}}}{\begin{pmatrix}h&-h&0&\ldots &0\\-1&2&-1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&-1&2&-1\\0&\ldots &0&-h&h\end{pmatrix}}{\begin{pmatrix}u_{0}\\u_{1}\\u_{2}\\\vdots \\u_{n+1}\end{pmatrix}}={\begin{pmatrix}f_{0}-\alpha \\f_{1}\\f_{2}\\\vdots \\f_{n+1}+\beta \end{pmatrix}}$

Approximation zweiter Ordnung

Ersetzt man die Ableitung in den Randpunkten $a,b$ durch den ersten zentralen Differenzenquotienten, approximiert man die Ableitung $u$ in $a$ bzw. $u'(a)$ und $u$ in $b$ bzw. $u'(b)$ durch:

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

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

Dazu werden erneute neue unbekannte Funktionswerte $u_{-1}$ und $u_{n+2}$ festgelegt. Dabei liegt $u_{-1}$ links vom Randpunkt $u_{0}$ und $u_{n+2}$ rechts vom Randpunkt $u_{n+1}$. Durch Umformen erhalten wir nun erneut zwei neue Gleichungen des LGS. DIe $-1$-te Gleichung und die $n+2$-te Gleichung mit

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

Nun möchten wir aber unser System nicht wieder um zwei weitere neue unbekannte ergänzen, weswegen wir auf einen anderen Ansatz zurückgreifen. Dementsprechend ersetzen wir die neuen Unbekannten in den ersten zentralen Differenzenquotienten durch die Approximation der zweiten Ableitung in den Randpunkten, also durch den zweiten Differenzenquotienten an den Randpunkten $u_{o},u_{n+1}$:

$D_{h}^{2}u(a):={\frac {u_{1}-2u_{0}+u_{-1}}{h^{2}}}=u''(a)$ und  $D_{h}^{2}u(b):={\frac {u_{n}-2u_{n+1}+u_{n+2}}{h^{2}}}=u''(b)$

$u_{-1}=u_{1}+2h\alpha +{\mathcal {O}}(h^{3})$ und $u_{n+2}=u_{n}+2h\beta +{\mathcal {O}}(h^{3})$ 

Nun setzen wir in die Approximation der zweiten Ableitung ein:

$D_{h}^{2}u(a):={\frac {u_{1}-2u_{0}+u_{1}+2h\alpha +{\mathcal {O}}(h^{3})}{h^{2}}}=u''(a)$ und $D_{h}^{2}u(b):={\frac {u_{n}+2h\beta +{\mathcal {O}}(h^{3})-2u_{n+1}+u_{n}}{h^{2}}}=u''(b)$

Durch Umformen gelangt man letztendlich zu $u''(a)={\frac {2u_{1}-2u_{0}}{h^{2}}}+{\frac {2\alpha }{h^{2}}}+{\mathcal {O}}(h^{2})$ bzw. ${\frac {2u_{1}-2u_{0}}{h^{2}}}=-2\alpha /h+{\mathcal {O}}(h)$ und $u''(b)={\frac {2u_{n}-2u_{n+1}}{h^{2}}}+{\frac {2\beta }{h^{2}}}+{\mathcal {O}}(h^{2})$ bzw. ${\frac {2u_{n}-2u_{n+1}}{h^{2}}}=-2\beta /h+{\mathcal {O}}(h^{2}).$

Systemmatrix

So ergibt sich in Matrixform folgende Systemmatrix $\quad A_{h_{2}}\in {\mathbb {R} }^{(n+2)\times (n+2)}$ mit $A_{h_{2}}=-{\frac {1}{h^{2}}}{\begin{pmatrix}2&-2&0&\ldots &0\\-1&2&-1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&-1&2&-1\\0&\ldots &0&-2&2\end{pmatrix}}$ und Vektor u ${\vec {u}}{\begin{pmatrix}u_{0}\\u_{1}\\u_{2}\\\vdots \\u_{n+1}\end{pmatrix}}$

Lineares Gleichungssystem für Neumann-Problem für Approximation zweiter Ordnung

$-{\frac {1}{h^{2}}}{\begin{pmatrix}2&-2&0&\ldots &0\\-1&2&-1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&-1&2&-1\\0&\ldots &0&-2&2\end{pmatrix}}{\begin{pmatrix}u_{0}\\u_{1}\\u_{2}\\\vdots \\u_{n+1}\end{pmatrix}}={\begin{pmatrix}f_{0}-{\frac {2}{h}}\alpha \\f_{1}\\f_{2}\\\vdots \\f_{n+1}-{\frac {2}{h}}\beta \end{pmatrix}}$

Regularität der Systemmatrizen

Das somit entstehende LGS wäre nun nur eindeutig lösbar, wenn die Systemmatrizen regulär wären. Das heißt:

Eine $(n\times n)$-Matrix $A$ mit Einträgen aus einem Körper $K$, zum Beispiel die reellen oder komplexen Zahlen, ist genau dann invertierbar, wenn eine der folgenden äquivalenten Bedingungen erfüllt ist:

• Es gibt eine Matrix $B$ mit $AB=I=BA$.
• Die Determinante von $A$ ist ungleich null.
• Die Eigenwerte von $A$ sind alle ungleich null.
• Für alle $b\in K^{n}$ existiert mindestens eine Lösung $x\in K^{n}$ des linearen Gleichungssystems $Ax=b$.
• Für alle $b\in K^{n}$ existiert höchstens eine Lösung $x\in K^{n}$ des linearen Gleichungssystems $Ax=b$.
• Das lineare Gleichungssystem $Ax=0$ besitzt nur die triviale Lösung $x=0$.
• Die Zeilenvektoren sind linear unabhängig.
• Die Zeilenvektoren erzeugen $K^{n}$.
• Die Spaltenvektoren sind linear unabhängig.
• Die Spaltenvektoren erzeugen $K^{n}$.
• Der Rang der Matrix $A$ ist gleich $n$.

Berechnet man nun die Determinante der Matrizen $A_{h_{1}}$ und $A_{h_{2}}$ so ist diese sowohl für die Systemmatrix der Approximation erster Ordnung, als auch für die Approximation zweiter Ordnung $0$. Ebenso gut könnte man allerdings auch zum Beispiel den Rang der Matrix über das Gauß Verfahren bestimmen. Siehe Beispiel $A_{h_{1}}$ für $n=3$. $Rang(A_{h_{1}})=\left\{{\begin{array}{l}3,&h=0\\4,&h\not =0\\\end{array}}\right.$ Des Weiteren sind die Vektoren linear abhängig. Multipliziert man die Systemmatrizen von rechts mit dem Einsvektor erhält man den Nullvektor und somit lässt sich der Nullvektor nicht nur durch die triviale Lösung erzeugen. Daraus folgt aus oben gegebenen äquivalenten Bedingungen, dass die Matrix nicht regulär ist und somit existiert keine eindeutige Lösung für das LGS.

Dirichlet-Randwertproblem in zwei Raumdimensionen

Homogene Randbedingungen

Gegeben ist nun folgendes Dirichlet-Randwertproblem mit homogenen Randbedingungen auf dem rechteckigen Gebiet $D=[a,b]\times [c,d]$ mit $D\subset \mathbb {R^{2}}$:

$-\Delta u(x)=f(x)$ in $D\subset \mathbb {R^{2}}$
      $u=0$ auf $\partial D$

 %----------------------------------
 %Implementierung der Quellfunktion 
 %----------------------------------
 function wert=Quellfunktion_FDM(x,y)
   if sqrt((x.-3.25).^2+(y.-3).^2)<=0.35 wert=1;
     elseif
       sqrt((x.-6).^2+(y.-1).^2)<=0.45 wert=1.5;
  else wert=0;
  endif
 endfunction

%------------------------------------------------
 %Eingabe der homogenen Dirichlet-Randbedingungen
 %------------------------------------------------
 funktion_unten  = @(x) [0];
 funktion_oben   = @(x) [0];
 funktion_links  = @(x) [0];
 funktion_rechts = @(x) [0];

Gesucht wird nun die Funktion $u:\mathbb {R^{2}} \to \mathbb {R}$ mit $(x,y)\in D$

Nun überziehen wir unser rechteckiges Gebiet $D$ mit einem äquidistanten Netz aus Gitterpunkten. Da wir uns im zweidimensionalen Bereich befinden benötigen wir Gitterpunkte in $x-$ und $y-$Richtung mit konstanter Gittergröße $h$.

$X_{i,j}=\{(x_{i},y_{j}),i=1,...,n,j=1,...,m\}$
$h=x_{i}-x_{i-1}=y_{j}-y_{j-1}={\frac {b-a}{n+1}}={\frac {d-c}{m+1}}$

 %--------------------------------------
 %Aufspannen des rechteckigen Gebiets D
 %--------------------------------------
 %definiere Gebiet D
 xend=8;
 yend=5.0; 
 %--------------------------------------
 %Implementierung des Gitters in Octave
 %--------------------------------------
 %Anzahl von X-Punkten, verkleinern Sie diese um die Struktur der Matrix in weiteren Zellen besser zu erkennen.
 N=60
 %Berechnung der konstanten Gittergröße
 hx=xend/(N+1)
 hy=hx;
 M=floor((yend/hy))-1;
 %Koordinatenmatrix der Gitterpunkte
 [x,y]=meshgrid(hx:hx:xend-hx,hy:hy:yend-hy);

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

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

Die Differenzenquotienten mit welchen die Funktionswerte in den Gitterpunkten approximiert werden, müssen nun ebenfalls wie unsere Gitterpunkte an die zweite Raumdimension angepasst werden. Dabei werden die zweiten Differenzenquotienten $D_{h,x}^{2}u(x_{i},y_{j}),\ D_{h,y}^{2}u(x_{i},y_{j})$ genutzt, um die partiellen Ableitungen von der gesuchten Funktion $u$ nach x und y zu approximieren:

${\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}}}$

Ersetzt man nun die zweite Ableitung durch dieses Schema in jedem Gitterpunkt $(x_{i},y_{j})$ ersetzen, bekommen wir unser Approximationsschema. So erhalten wir als Approximation des Laplace-Operators, die Summe dieser zwei Differenzenquotienten:

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

Diese Approximation führen wir für jeden Gitterpunkt durch, wodurch wir $n\times m$- Gleichungen, die durch die Nachbarpunkte des Punktes $u_{i,j}$ verbunden sind. Dieses Schema nennt man auch das 5-Punkt-Sternschema. Dieses ergibt sich nach dem Approximationsschema, indem man 4-mal den zentralen Punkt $u_{i,j}$ nimmt und jeweils einmal den oberen, unteren, rechten und linken Nachbar hinzu addiert. Die Randwerte werden nun approximiert, indem man die Nachbarpunkte im Rand durch die gegebenen Randwertbedingungen ersetzt. 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$.

Nun werden die unbekannten Approximationen der Lösung in der Gitterpunkten in einem langen Vektor eingeordnet. Dabei werden diese Werte in einer bestimmten Reihenfolge in den Vektor eingeordnet. Siehe . Diese Assemblierung findet ebenfalls für den Vektor der rechten Seite statt. Nach der oben beschriebenen numerischen Diskretisierung und Assemblierung des unbekannten Vektors ergibt sich nun ein LGS der Form $A_{h}{\vec {u_{h}}}=f_{h}$ mit der Blocktridiagonalmatrix $A_{h}\in \mathbb {R} ^{(n\cdot m)\times (n\cdot m)}$

size (f);
 %reshape funktioniert Spaltenmässig: es wird Spalte1, dann Spalte2, usw eingeordnet, wobei wir die Matrix f nach Reihen/Zeilen einordnen.
 %Daher soll die Matrix f' mit N(Zeilen)x M(Spalten) in ein langern Vektor eingeordnet.
 b=-1*reshape(f',(N*M),1);

 %----------------------------------------------------------------------------
 %Aufstellen der Systemmatrix durch Zusammensetzen aus verschiedenen Blöcken
 %----------------------------------------------------------------------------
 %Funktion zur Anpassung des Zwischenblocks der Matrix 
 %-----------------------------------------------------
 function wert = Zwischenblock(N,B1,B,C1,zeilennummer)
   wert = [repmat(C1,[1,zeilennummer-2]),B1,B,B1,repmat(C1,[1,N-(zeilennummer+1)])];
 endfunction
 %-----------------------------
 %Erstellen der Blockmatrix
 %-----------------------------
 function wert = Blockmatrix(N,M,h)
   %---------------------------------
   %Erstellen der einzelnen Matrizen
   %---------------------------------
   B = diag(-4*ones(N,1),0)+diag(ones(N-1,1),-1)+diag(ones(N-1,1),1);%Matrix B
   I=diag(ones(N,1),0);%Einheitsmatrix
   C1=diag(zeros(N,1),0);%Nullmatrix
   %---------------------------------
   %Anfang der Blockmatrix 
   %--------------------------------- 
   Ah=[B,I,repmat(C1,[1,M-2]); %1.Zeile
       I,B,I,repmat(C1,[1,M-3])]; %2.Zeile 
   %----------------------------------------------------------
   %Zwischenzeilen der Blockmatrix aus oben genannter Funktion
   %----------------------------------------------------------  
   for i = 1:M-4
     Ah=[Ah;Zwischenblock(M,I,B,C1,i+2)];
   endfor
   %----------------------------------
   %Abschließende Zeilen 
   Ah= [Ah;repmat(C1,[1,M-3]),I,B,I;
        repmat(C1,[1,M-2]),I,B]; 
   %-------------------------------------------------  
   %Definition der oben eingeführten Variable 'Wert'
   %--------------------------------------------------
   wert = Ah;
 endfunction
 
 B=Blockmatrix(N,M,hx);
 B=B*a/hx^2;

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

 %------------------------------------------------------------------
 %Anpassung der rechten Seite auf dem Rand durch die Randbedingungen 
 %------------------------------------------------------------------
 for i=1+1:N  
 for j=1+1:M
 f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i));
 %Dirichlet Randbedingung oben und unten: 
   if j==1 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_unten(i*hx)*a/hx^2;  %(eine Konstante oder eine Funktion von x-also Vektor1(i)); 
   endif
   if j==M 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_oben(i*hx)*a/hx^2;  %(eine Konstanteoder eine Funktion von x -also Vektor2(i)); 
   endif
 %Dirichlet Randbedingung links und rechts:
   if (i==1)
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_links(j*hx)*a/hx^2;    %(eine Konstante oder eine Funktion von y); endif
   endif
   if (i==N )
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_rechts(j*hx)*a/hx^2;    %(eine Konstante oder eine Funktion von y); endif
   endif
 endfor
 endfor

 % Erneute Definition der unteren Reihe, weil es  irgendwie mit der funktion_unten, funktion_links nicht funktioniert...
 f(1,:)= 0*a/hx^2;
 f(:,1)= 0*a/hx^2;
 
 %Eckbeiträge
 f(1,1)=0;
 f(1,N)=0;
 f(M,1)=0;
 f(M,N)=0;

 % Wichtig!: Ecken-Beiträge: da in den Ecken jeweils der Differenzenquotioent in x (Wert Links oder rechts) in der Matrix schon festgelegt ist, muss Null entsprechend auf der rechten Sete stehen.

 %-------------------------------------------------
 % LOESUNGSSCHRITT
 %-------------------------------------------------
 sol=B\b;
 sol_matrix=reshape(sol,N,M);% Matrix mit N(Zeilen)x M(Spalten)
 sol_matrix=sol_matrix'; % Transponiert, wegen reshape

$A_{h}{\vec {u_{h}}}=f_{h}$

$-{\frac {1}{h^{2}}}\left[{\begin{array}{c|c}{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &\cdots &\cdots &0\\\hline {\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &\cdots &0\\\hline 0&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &0\\\hline \vdots &\vdots &\ddots &\ddots &\ddots &\cdots &0\\\vdots &\vdots &\vdots &\ddots &\ddots &\ddots &0\\\hline 0&\cdots &\cdots &\cdots &\cdots &{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}\end{array}}\right]{\begin{pmatrix}u_{1,1}\\u_{1,2}\\\vdots \\u_{1,M}\\\hline u_{2,1}\\u_{2,2}\\\vdots \\u_{2,M}\\\hline u_{3,1}\\u_{3,2}\\\vdots \\u_{3,M}\\\hline \vdots \\\vdots \\\vdots \\\hline u_{N,1}\\u_{N,2}\\\vdots \\u_{N,M}\end{pmatrix}}={\begin{pmatrix}f_{1,1}\\f_{1,2}\\\vdots \\f_{1,M}\\\hline f_{2,1}\\f_{2,2}\\\vdots \\f_{2,M}\\\hline f_{3,1}\\f_{3,2}\\\vdots \\f_{3,M}\\\hline \vdots \\\vdots \\\vdots \\\hline f_{N,1}\\f_{N,2}\\\vdots \\f_{N,M}\end{pmatrix}}$

Inhomogene Randbedingungen

Sind die Randwertbedingungen nicht homogen, also unterschiedlich von Null, dann tragen die bekannten Randwerte der gesuchten Funktion zum Vektor der rechten Seite ${\vec {f_{h}}}$.

 %Eingabe der inhomogenen Dirichlet-Ranbedingungen
 funktion_oben   = @(x) [0.01*x];
 funktion_unten  = @(x) [-0.25];
 funktion_rechts = @(x) [0.25];
 funktion_links  = @(x) [-0.2];

 %------------------------------------------------------------------
 %Anpassung der rechten Seite auf dem Rand durch die Randbedingungen 
 %------------------------------------------------------------------
 for i=1+1:N  
 for j=1+1:M
 f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i));
 %Dirichlet Randbedingung oben und unten: 
   if j==1 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_unten(i*hx)*a/hx^2;  
   endif
   if j==M 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_oben(i*hx)*a/hx^2;  
   endif
 %Dirichlet Randbedingung links und rechts:
   if (i==1)
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_links(j*hx)*a/hx^2;    
   endif
   if (i==N )
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_rechts(j*hx)*a/hx^2;    
   endif
 endfor
 endfor

 %Eckbeitraege!
 f(1,1)=1*Quellfunktion_FDM(x(1,1),y(1,1))+ (funktion_links(1*hx) + funktion_unten(1*hx))*a/hx^2;
 f(1,N)=1*Quellfunktion_FDM(x(1,N),y(1,N))+ (funktion_rechts(1*hx)+ funktion_unten(N*hx))*a/hx^2; %f(1,1);
 f(M,1)=1*Quellfunktion_FDM(x(M,1),y(M,1))+ (funktion_links(M*hx)+ funktion_oben(1*hx))*a/hx^2;
 f(M,N)=1*Quellfunktion_FDM(x(M,N),y(M,N))+ (funktion_rechts(M*hx)+ funktion_oben(N*hx))*a/hx^2; %f(M,1);

$A_{h}{\vec {u_{h}}}=f_{h}$

$-{\frac {1}{h^{2}}}\left[{\begin{array}{c|c}{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &\cdots &\cdots &0\\\hline {\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &\cdots &0\\\hline 0&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}&\cdots &\cdots &0\\\hline \vdots &\vdots &\ddots &\ddots &\ddots &\cdots &0\\\vdots &\vdots &\vdots &\ddots &\ddots &\ddots &0\\\hline 0&\cdots &\cdots &\cdots &\cdots &{\begin{matrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots &1\\0&\ldots &0&1&-4\end{matrix}}&{\begin{matrix}1&0&0&\ldots &0\\0&1&0&\ldots &0\\\vdots &\ddots &\ddots &\ddots &0\\0&\ldots &0&0&1\end{matrix}}\end{array}}\right]{\begin{pmatrix}u_{1,1}\\u_{2,1}\\\vdots \\u_{N,1}\\\hline u_{1,2}\\u_{2,2}\\\vdots \\u_{N,2}\\\hline u_{1,3}\\u_{2,3}\\\vdots \\u_{N,3}\\\hline \vdots \\\vdots \\\vdots \\\hline u_{1,M}\\u_{2,M}\\\vdots \\u_{N,M}\end{pmatrix}}={\begin{pmatrix}f_{1,1}+{\frac {1}{h^{2}}}(lR+uR)\\f_{2,1}+{\frac {1}{h^{2}}}lR\\\vdots \\f_{N,1}+{\frac {1}{h^{2}}}(lR+oR)\\\hline f_{1,2}+{\frac {1}{h^{2}}}uR\\f_{2,2}+0\\\vdots \\f_{N,2}+{\frac {1}{h^{2}}}oR\\\hline f_{1,3}+{\frac {1}{h^{2}}}uR\\f_{2,3}+0\\\vdots \\f_{N,3}+{\frac {1}{h^{2}}}oR\\\hline \vdots \\\vdots \\\vdots \\\hline f_{1,M}+{\frac {1}{h^{2}}}(rR+uR)\\f_{2,M}+{\frac {1}{h^{2}}}rR\\\vdots \\f_{N,M}+{\frac {1}{h^{2}}}(rR+oR)\end{pmatrix}}$

Neumann-Problem in zwei Raumdimensionen

Ähnlich wie beim Dirichlet-Problem wird die Anwendung der zweidimensionalen Differenzenquotienten in ein 5Punkt-Stern-Approximationsschema resultieren, bei dem ebenfalls in jedem Gitterpunkt der zentrale Punkte $u_{i,j}$ viermal genommen wird und man jeweils den oberen, unteren, rechten und linken Nachbar hinzu addiert. Durch Behandlung der Randbedingungen werden zur Änderung einiger Blöcke der tridiagonalen Blockmatrix führen. Analog zum Neumann-Problem in einer Raumdimension wird die Matrix lediglich um die die Unbekannten $u_{0,j},u_{n+1,j},u_{i,0},u_{i,m+1}$ erweitert. Diese werden letztlich wieder zur Bildung der Differenzenquotienten verwendet.

Zunächst werden wieder die Randbedingungen gegeben. Das heißt, die Richtungsableitungen in die Normalenrichtung der unbekannten Funktion in den Randpunkten des Rechtecks werden 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

Bei der Anwendung der einseitigen Differenzenquotienten, analog wie 1-dimensionalem Fall muss der Vektor der 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.)

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

So ergibt sich in Matrixform folgender Block B mit $B=-{\frac {1}{h^{2}}}{\begin{pmatrix}h&-h&0&\ldots &0\\1&4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&1&4&1\\0&\ldots &0&h&-h\end{pmatrix}}$ und Vektor u ${\vec {u}}={\begin{pmatrix}u_{0}\\u_{1}\\u_{2}\\\vdots \\u_{n+1}\end{pmatrix}}$.

Gemischte Randbedingungen

Aufgabe war es nun ein Programm zu schreiben, dass die Unbekannte Funktion $u(x,y)$ in den Gitterpunkten mit gemischten Randbedingungen approximiert. Linkes und rechts sollten homogene Neumann-Randbedingungen behandelten und oben und unten sollten inhomogene Dirichlet-Randbedingungen behandelt werden. Zunächst haben wir dazu unser rechteckiges Gebiet aufgespannt. Und danach das Gebiet mit dem erweiterten Gitter des Neumann-Problems überzogen. Dementsprechend mussten wir auch unsere Dimensionen der Matrizen N und M jeweils um eine Zeile oben und unten und eine Spalte rechts und links erweitern, da sonst das LGS zum Ende des Programms nicht aufgegangen wäre.

 %definiere Gebiet D
 xend=8;
 yend=5.0;
 %Anzahl von X-Punkten, verkleinern Sie diese um die Struktur der Matrix in weiteren Zellen besser zu erkennen.
 N=60
 %--------------------
 hx=xend/(N+1)
 hy=hx;
 M=floor((yend/hy))-1 
 %Koordinatenmatrix der Gitterpunkte
 [x,y]=meshgrid(0:hx:xend,0:hy:yend); 
 %wegen Erweiterung des Gitters bei Neumann-Problem
 N=N+2;
 M=M+2;

 %konstanter Diffusionskoeffizient
 a=0.1;

 %Dirichlet RB:
 funktion_unten  = @(x) [1];% [0.15*x];
 funktion_oben   = @(x) [-1];%[0.1*x];

Ebenfalls wie beim Dirichlet Problem haben wir dann unsere Funktion der rechten Seite implementiert.

 %------------------------------
 %Erstellen der Quellfunktion
 %------------------------------ 
 function wert=Quellfunktion_FDM(x,y)
   if sqrt((x.-3.25).^2+(y.-3).^2)<=0.35 wert=1;
     elseif
     sqrt((x.-6).^2+(y.-1).^2)<=0.45 wert=1.5;
   else wert=0;
   endif
 endfunction

Danach folgte der wohl wichtigste Schritt, nämlich die Anpassung der Blockmatrix, da sich die Blöcke durch die Erweiterung nach Neumann verändern. So haben wir Matrix $B$ um die Unbekannten, siehe Neumann Problem in einer Raumdimension erweitert. Auch die Einheitsmatrix $I$ wurde um die 0-te und $n+1$-Zeile erweitert, die jedoch mit Nullen gefüllt sind.

 %-------------------------------
 %Erstellen der Blockmatrix 
 %-------------------------------
 function wert = Blockmatrix(N,M,h)
   B = diag(-4*ones(N,1),0)+diag(ones(N-1,1),-1)+diag(ones(N-1,1),1);
     %----------------------------------------------
     %Erweiterung des B-Blocks durch die Unbekannten 
     %----------------------------------------------
     B(1,1)=-h;
     B(1,2)=+h;
     B(N,(N-1))=h;
     B(N,N)=-h;
   %---------------------------------
   %Erstellen der einzelnen Matrizen
   %---------------------------------
   I=diag(ones(N,1),0);
     %---------------------------------
     %Anpassen von 0-ter und n+1-Zeile 
     %---------------------------------
     I(1,1)=0;
     I(N,N)=0;
   C1=diag(zeros(N,1),0);
   %--------------------------
   %Anfang der Blockmatrix  
   %--------------------------
   Ah=[B,I,repmat(C1,[1,M-2]); 
       I,B,I,repmat(C1,[1,M-3])]; 
   %-------------------------------   
   %Zwischenzeilen der Blockmatrix  
   %-------------------------------
   for i = 1:M-4
     Ah=[Ah;Zwischenblock(M,I,B,C1,i+2)];
   endfor
    %-----------------------------------
    %Letzen zwei Zeilen der Blockmatrix
    %-----------------------------------
    Ah=[Ah;repmat(C1,[1,M-3]),I,B,I;
        repmat(C1,[1,M-2]),I,B]; 
   
   wert = Ah;
 endfunction

Aufgrund der inhomogenen Dirichlet-Randbedingungen muss der Vektor der rechten Seite angepasst werden. Dazu werden an den oberen und unteren Ränder die Bedingungen behandelt. Da jedoch Neumann für links und rechts gilt und diese homogen sind, verändern diese nichts am Vektor der rechten Seite. %rechte Seite des Systems

 for i=1+1:N  
 for j=1+1:M
 f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i));
   if j==1 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_unten(i*hx)*a/hx^2;  %(eine Konstante oder eine Funktion von x-also Vektor1(i)); 
   endif
   if j==M 
     f(j,i)=1*Quellfunktion_FDM(x(j,i),y(j,i))+ funktion_oben(i*hx)*a/hx^2;  %(eine Konstanteoder eine Funktion von x -also Vektor2(i)); 
   endif
 endfor
 endfor

Da an den Ecken jeweils die homogenen Randbedingungen gelten, müssen die Eckbeiträge nicht angepasst werden, wie es im reinen Dirichlet-Problem mit inhomogenen Randbedingungen ist. Aufgrund der

 %Eckbeiräge
 f(1,1)=0;
 f(1,N)=0;
 f(M,1)=0;
 f(M,N)=0;
 % Wichtig!: Ecken-Beiträge: da in den Ecken jeweils der Differenzenquotion in x (Wert Links oder rechts) in der Matrix schon festgelegt ist, muss Null entsprechend auf der rechten Sete stehen

Nun werden die unbekannten Approximationen der Lösung in der Gitterpunkten in einem langen Vektor eingeordnet. Dabei werden diese Werte in einer bestimmten Reihenfolge in den Vektor eingeordnet. Siehe [5]. Diese Assemblierung findet ebenfalls für den Vektor der rechten Seite statt. Nach der oben beschriebenen numerischen Diskretisierung und Assemblierung des unbekannten Vektors ergibt sich nun ein LGS der Form ${\displaystyle A_{h}{\vec {u_{h}}}=f_{h}}$

 size (f);
 %reshape funktioniert Spaltenmässig: es wird Spalte1, dann Spalte2, usw eingeordnet, wobei wir die Matrix f nach Reihen/Zeilen einordnen.
 %Daher soll die Matrix f' mit N(Zeilen)x M(Spalten) in ein langern Vektor eingeordnet.
 b=-1*reshape(f',(N*M),1);
 
 %-------------------------------------------------
 % LOESUNGSSCHRITT
 %-------------------------------------------------
 sol=B\b;
 sol_matrix=reshape(sol,N,M);% Matrix mit N(Zeilen)x M(Spalten)
 sol_matrix=sol_matrix'; % Transponiert, wegen reshape

Tutorium 5: Abschlussprojekt

Theoretischer Input

Für den Abschluss der letzten Tutoriumsaufgaben wird letztendlich sämtliche Theorie in einer letzten Modellierung zusammengebracht. Beginnend mit der Reaktionsdiffusionsgleichung, fortlaufend mit den Kompartimentmodellen und abschließend mit dem FDM-Verfahren sowie den Randwertproblemen.

Numerische Diskretisierung der Reaktionsdiffusionsgleichung

Die instationäre, inhomogene partielle Differentialgleichung der Reaktionsdiffusionsgleichung ist bereits aus Tutorium 2 und 3 bekannt, bietet sie uns die Möglichkeit über den Reaktionsterm der rechten Seite, die in Tutorium 3 behandelten Modelle, mit der Diffusion zu verbinden. Somit stellt der Reaktionsterm den Zusammenhang zwischen der räumlichen Infektionsverbreitung und den dynamischen Kompartimentmodellen her.

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

Die Funktion u sei dabei die Populationsdichte der (aktuell oder kummulierten) Infizierten Individuen, wobei der Reaktionsterm $f=f(u)$ die zeitliche Dynamik der aktuell oder der kummulierten Infizierten steuert.

Räumliche Semidiskretisierung für die Diffusiongleichung

(Genommen aus Skript Frau Prof. Dr. Hundertmark)

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.

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.

Ergänzungen Octave-Tutorial

Diese Ergänzungen wurden bei dem Octave-Tutorial vorgenommen, um die Implementation der Gruppe bzgl. der verwendeten Octave-Befehle nachvollziehbar zu machen.

• Matrix-Multiplikation in Octave - notwendig für die Modellierung diskreter Transportprozesse.

Skripte in Octave

Referenzieren Sie hier die von Ihnen verwendeten Skripte und Bibliotheken in Octave oder R/Studio

Literatur

Notieren Sie hier die von Ihnen verwendete Literatur

• Quelle 1