Kurs:Räumliche Modellbildung/Gruppe 8

Gruppenseite - ST

Diese Seite enthält die gesammelten Materialien der Gruppe 8 (Spaniol, Traub) für die Portfolioprüfung.

COVID-19/Mathematische Modellbildung

Teilnehmer-innen

Teilnehmer-innen treffen sich standardmäßig in Videokonferenz-Breakoutroom 8

(S) Spaniol, Michelle
(T) Traub, Sebastian

Diskrete Modellierung von Diffusionsprozessen

Einfaches SIR-Modell

Bei einem SIR Modell geht man von drei Gruppen aus:

Susceptible (engl. Empfängliche)
Infected (engl. Infizierte)
Recovered (engl. Genesene)

Dabei können Susceptible zu den Infected und Infected zu den Recovered wechseln. Wenn man nun Raten definiert, mit denen diese Wechselprozesse stattfinden, kann man Krankheitsverläufe rudimentär modellieren.

Anfangszustand in Deutschland

In unserem Fall sind wir für die Susceptible von der Gesamtbevölkerung von Deutschland ausgegangen sowie einem Infizierten

 Infected=1;            
 Susceptible=83000000-Infected;
 Recovered=0;

Festlegung der Parameter

Im nächsten Schritt wird die Basisreproduktionszahl (BRZ) sowie die Anzahl der Zeitschritte festgelegt. Die Basisreproduktionsrate gibt an, wie viele Individuen ein Infizierter in einem Zeitschritt ansteckt. Ein Zeitschritt entspricht hier 14 Tagen, da ein Infizierter in etwa so lange seine Krankheit weitergeben kann.

 BRZ=1.5;
 Zeitschritte=70;
 t=linspace(0, Zeitschritte, Zeitschritte+1);

Der Befehl linspace erzeugt einen Vektor mit den Einträgen aller Zeitschritte in einer Zeile.

Aufstellen von Vektoren für S, I und R

Uns interessiert die Entwicklung der Zahlen für S, I und R in allen Zeitschritten. Diese Zahlen sollen jeweils in einem Vektor gespeichert werden. Demnach werden drei Vektoren benötigt, die die gleiche Länge des Zeitvektors t haben. Zu Anfang haben diese Vektoren nur Nullen als Einträge.

 S=(zeros(Zeitschritte+1, 1));      
 I=(zeros(Zeitschritte+1, 1));
 R=(zeros(Zeitschritte+1, 1));

Der erste Eintrag soll die Anfangszustände beinhalten.

 S(1)=Susceptible;
 I(1)=Infected;
 R(1)=Recovered;

Epidemiologischer Prozess

Im nächsten Schritt folgt die Berechnung des Krankheitsverlaufs auf Grundlage der Anfangsbedingungen. Grundgerüst ist eine Schleife, die für jeden nächsten Zeitschritt die Werte für S, I und R berechnet und in die Vektoren speichert.

 for n=1:Zeitschritte
     I(n+1)=I(n)*BRZ*(S(n)/S(1));
     S(n+1)=S(n)-I(n+1);
     R(n+1)=R(n)+I(n);
 end

Die Infizierten des nächsten Zeitschrittes berechnen sich durch die jetzigen Infizierten mal der Basisreproduktionszahl. Da im Laufe der Zeit die Zahl der Susceptible sinkt, kann ein Individuum weniger andere anstecken. Dieser Umstand wird durch den Faktor (S(n)/S(1)) berücksichtigt.
Die Susceptible des nächsten Zeitschrittes berechnen sich durch die jetzigen Susceptibles abzüglich der neuen Infizierten
Da nach einem Zeitschritt die Infizierten die Krankheit überwunden haben, berechnen sich die nächsten Recovered durch die Summe aus aktuellen Recovered und Infizierten

Skript in Cocalc

Einfaches SIR-Modell

Einfaches SIR-Modell mit Transportprozessen

Bei dem vorherigen Modell sind wir von einer konstanten Bevölkerungszahl an einem Ort ausgegangen. Die Erweiterung besteht jetzt darin, dass sich diese konstante Bevölkerungszahl auf 5 verschiedene Orte verteilt. Zusätzlich zum epidemiologischen Prozess werden jetzt noch zusätzlich Transportprozesse zwischen den 5 Orten implementiert.

Festlegung der Parameter

In unserem Fall wurde eine 5x5-Matrix erstellt, die den Transportprozess beschreibt. Hierbei geben konstante Raten den Wechsel der Personen aus einem Kompartiment in ein anderes an. Außerdem haben wir die Parameter Basisreproduktionszahl (BRZ), Zeitschritte und Dimension definiert. Mit der Festlegung der Dimension auf den Wert 5 gehen wir im Modell von ebenso vielen Orten aus.

 BRZ=1.5;   
 Zeitschritte=100; 
 Dimension=5;
 TM=[0.89, 0.03, 0.01, 0.04, 0.03;    
   0.01, 0.91, 0.05, 0.01, 0.02;
   0.03, 0.02, 0.90, 0.04, 0.01;
   0.05, 0.02, 0.03, 0.88, 0.02;
   0.02, 0.02, 0.01, 0.03, 0.92];
 t= linspace(0, Zeitschritte, Zeitschritte+1);

Anfangszustände von S, I und R

Im nächsten Schritt werden zuerst Nullmatrizen für die Werte der Susceptible, Infected und Recovered aufgestellt, die dann folgend mit Werten gefüllt werden. Des Weiteren definieren wir die Anfangszustände, wobei für S, I und R jeweils eine 5x1 Matrix aufgestellt werden muss. In besagter Matrix werden für die 5 Orte jeweils Anfangszustände festgelegt, die den Personenstand in dem jeweiligen Kompartiment (Ort) beschreiben. Hierbei ist der Anschauungsraum Deutschland, weshalb wir wieder von der entsprechenden Bevölkerungszahl ausgehen. Die Susceptible belaufen sich dahingehend zusammen auf 83 Mio (-1) und in einem einzigen Ort befindet sich zu Beginn ein Infected.

 Snull=zeros(Zeitschritte-1, Dimension);   
 Inull=zeros(Zeitschritte-1, Dimension);
 Rnull=zeros(Zeitschritte-1, Dimension);
 AS=[20000000, 13500000, 25000000, 8500000, 15999999];
 AI=[0,0,0,0,1];  
 AR=[0,0,0,0,0];

Im weiteren Verlauf setzen wir die Anfangszustände und die Nullmatrizen zusammen und definieren diese neuen Matrizen als unsere S-, I- und R-Matrizen für die kommenden Berechnungen.

 S=[AS; Snull];   
 I=[AI; Inull];
 R=[AR; Rnull];

Ablauf der Simulation

Ausgehend von den Anfangszuständen folgt im nächsten Abschnitt zuerst die Berechnung des Krankheitsverlaufs. Hier werden Susceptible, Infected und Recovered in einer Schleife für jeden Zeitschritt berechnet und entsprechend in den Matrizen festgehalten. Der epidemiologische Prozess wird hier durch ausgelagerte Funktionen beschrieben, die zu einem späteren Zeitpunkt näher erläutert werden. Damit ist die erste Schleife abgeschlossen und wir begeben uns in die zweite Schleife und damit in den Transportprozess. Hier kommt die Transportmatrix ins Spiel, die den Austausch der Personen zwischen den Kompartimenten beschreibt und damit ein Voranschreiten der Infektion begünstigt.

 for n=1:Zeitschritte-1
     for k=1:Dimension
          I=InachZeitschritt(I, S, BRZ, n, k);
          S=SnachZeitschritt(S, I, n, k);
          R=RnachZeitschritt(R, I, n, k);
     end
     S=Transportprozess(S, n, Dimension, TM);
     I=Transportprozess(I, n, Dimension, TM);
     R=Transportprozess(R, n, Dimension, TM);
 end

Epidemiologischer Prozess I

Nun werden die definierten Funktionen für den epidemiologischen Prozess näher erläutert. Die erste Funktion beschreibt die Berechnung der Infected im nächsten Zeitschritt. Die Vorgehensweise unterscheidet sich dabei nur hinsichtlich der hinzugekommenen Dimension von dem ersten SIR-Modell. Allerdings befinden sich hier die Einträge des nächsten Zeitschrittes nicht in der nächsten Spalte, sondern in der nächsten Zeile. Auch hier berechnen sich die Infected des nächsten Zeitschrittes, aus dem Produkt der aktuellen Infected, der Basisreproduktionszahl und einem weiteren Faktor. Dieser Faktor (S(n)/S(1)) berücksichtigt, wie im vorherigen Modell bereits erläutert, die sinkende Anzahl der Susceptible im Vergleich zum Anfangszustand.

 function[A]=InachZeitschritt(A, S, BRZ, n, k);
     A(n+1,k)=A(n,k)*BRZ*(S(n,k)/S(1,k));
 end

Epidemiologischer Prozess S

Bei der zweiten Funktion wird die Berechnung der Susceptible ausgeführt. Hier bestand zuerst das Problem, dass die Differenz aus jetzigen Susceptible und Infected des nächsten Zeitschrittes teilweise negative Werte lieferte. Dies wird durch den Transportprozess hervorgerufen, wodurch die Infected des nächsten Zeitschrittes teilweise größer sind als die jetzigen Susceptible. Dieses Problem haben wir mithilfe einer if-Abfrage behoben. Solange die Infected Zahl die Susceptibles nicht überschreitet, wird wie im vorherigen Modell die Differenz berechnet. Sind in einem Kompartiment mehr oder genauso viele Personen Infected wie Susceptible wird keine Differenz berechnet, sondern die Susceptible des nächsten Zeitschrittes werden auf 0 gesetzt. Somit wird das Problem mit den negativen Personenzahlen korrigiert.

 function[A]=SnachZeitschritt(A, I, n, k)
      if I(n+1,k)<A(n,k)
              A(n+1,k)=A(n,k)-I(n+1,k);
              else A(n+1,k)=0;
      endif
 end

Epidemiologischer Prozess R

Die dritte Funktion führt die Berechnung der Recovered aus. Dabei wird wie im Modell zuvor die Summe aus den jetzigen Recovered und den jetzigen Infected gebildet, um die Recovered des nächsten Zeitschrittes zu bestimmen.

 function[A]=RnachZeitschritt(A, I, n, k);
     A(n+1,k)=A(n,k)+I(n,k);
 endfunction

Transportprozess

Die letzte Funktion beschreibt den Transportprozess, der zwischen den Kompartimenten stattfindet. Hier werden die zuvor berechneten neuen Infected, Susceptible und Recovered jeweils mit der Transportmatrix multipliziert. Das Matrixprodukt erzeugt somit den Wechsel von Personen in andere Orte und treibt die Verbreitung des Virus weiter an.

 function[A]=Transportprozess(A, n, Dimension, TM);
     A(n+1,1:Dimension)=A(n+1,1:Dimension)*TM;
 end

Epidemiologischer Prozess und Transportprozess in den einzelnen Kompartimenten

Skript in Cocalc

Einfaches SIR-Modell mit Transportprozessen

Kontinuierliche Modellierung von Diffusionsprozessen

Einfaches Kontaktmodell

Eine der größten Ausbreitungsmöglichkeiten für Krankheiten besteht im direkten Kontakt eines infizierten Individuums mit einem anderen Menschen. Das folgende Modell probiert diese Ausbreitungsmöglichkeit darzustellen. Ausgangssituation ist ein Koordinatensystem, in dem sich viele in einem Rechteck homogen verteilte Punkte befinden, die jeder einen Menschen darstellen. Am Anfang soll ein Individuum infiziert sein. Dies wird durch einen roten Kasten um die jeweilige Position des Infizierten gekennzeichnet. Wenn durch zufällige Bewegung der Abstand eines Infizierten und eines Gesunden unter einen festen Wert fällt, wird letzterer infiziert und kann andere anstecken.

Festlegen der Koordinaten für alle Individuen

Zuerst müssen die Koordinaten aller Menschen zu Anfang definiert werden. Dazu werden Vektoren x und y erstellt, die in Einer-Schritten Einträge von 1 bis zur Anzahl der Personen pro Raumdimension enthalten. Der Befehl meshgrid erstellt aus diesen Vektoren zwei Matrizen x und y der Dimension Nx x Ny, die an festen Indizes die Koordinaten einer Person beinhalten. So sind die Koordinaten der Person an Ausgangsposition (3|5) immer in x(3|5) und y(3|5) zu finden.

 Nx=10;  %Anzahl Individuen in x-Richtung
 Ny=20;  %Anzahl Individuen in y-Richtung
 x=(1:Nx);
 y=(1:Ny);
 [x, y]=meshgrid(x, y);

Dokumentation der anfänglich Infizierten

Als nächstes muss die Anfangsposition des ersten Infizierten festgelegt und gespeichert werden. Die Festlegung erfolgt in der Mitte aller anderen Personen. zur Speicherung wird eine dritte Matrix Inf der selben Dimension wie x und y erstellt. Diese wirkt quasi wie eine dritte Koordinate zu jeder Person. Ist der Eintrag an einer Stelle der Matrix Inf 0, ist die Person dieser Stelle gesund und ist er 1, ist die Person infiziert.

 xInf=x(1, Nx/2);  %Anfangsposition des 1. Infizierten
 yInf=y(Ny/2, 1);
 Inf=zeros(Ny, Nx); %Erstellung Matrix, die anzeigt an welchen Koordinaten Infizierte Personen sind 
 Inf(yInf, xInf)=1;

Die Ausgangssituation sieht folgendermaßen aus:

Allgemeine Parameter

Da davon ausgegangen wird, dass eine Infektionsübertragung bei einer Abstandsunterschreitung stattfindet, muss dieser Mindestabstand auch definiert werden. Ebenso wird die Bewegungsgeschwindigkeit, die Anzahl der Zeitschritte und zwei Zähler festgelegt. Die Bedeutung von letzterem wird später weiter erläutert. Die beiden Matrizen Rx und Ry geben die Richtung der Bewegung in x- bzw. y Richtung an und enthält zufallsgenerierte Zahlen zwischen 0 und 1.

 g=0.5;       %Variable für Geschwindigkeit
 radInf=0.5;  %Distanz unter derer eine Infektionsübertragung stattfindet
 p=3;         %Zähler für Zeitschritte wegen Richtungsänderung
 q=3;         %Zähler für Zeitschritte wegen Richtungsänderung
 Rx=rand(Ny, Nx);
 Ry=rand(Ny, Nx);
 
 Zeitschritte=30;

Ablauf der Simulation

Die eigentliche Simulation des Modells enthält 3 Schritte, welche pro Zeitschritt wiederholt werden

 for i=1:Zeitschritte
     [Rx, p]=richtungsaenderung(p, Nx, Ny, Rx);
     [Ry, q]=richtungsaenderung(q, Nx, Ny, Ry);
     p=p+1;
     q=q+1;
 
     x=NeuePosition(x, Rx, g);
     y=NeuePosition(y, Ry, g);
 
     Inf=abstandsabfrage(Nx, Ny, Inf, x, y, radInf);
 end

Richtungsänderung

Im ersten Schritt wird die Bewegungsrichtung geändert. Dies funktioniert über die Funktion "richtungsaenderung".

 function[R, n]=richtungsaenderung(n, Nx, Ny, Rjetzt);
     if n==3          %Änderung der Richtung nach 3 Schritten
     R=rand(Ny, Nx);  %Richtungsänderung in x Richtung
     n=0;
     else R=Rjetzt;
     endif
 end

Hier spielen die oben genannten Zähler p und q eine Rolle (Hier werden sie als Variable n in der Schleife genutzt). Diese starten bei Null und werden im Hauptprogramm in jeden Durchgang der Schleife um 1 erhöht. Solange der Zähler kleiner als 3 ist, wird die else-Bedingung erfüllt und es ändert sich nichts. Ist der Zähler allerdings 3, werden die Richtungsmatrizen Rx und Ry mit neuen Zufallszahlen gefüllt. Zudem wird der Zähler zurück auf Null gesetzt. Dies resultiert in einer Änderung der Bewegungsrichtung alle drei Zeitschritte.

Berechnung der neuen Position

Im zweiten Schritt wird auf Grundlage der aktuellen Position, der Richtungsmatrizen und der Geschwindigkeit die Neue Position bestimmt. Dies passiert für die x- und die y-Komponente separat.

 function[Position]=NeuePosition(Position, R, g)
     Position=Position.+(R.-0.5)*g;
 end

Es werden schlicht komponentenweise auf die aktuelle Position die Einträge der Richtungsmatrix addiert. Diese enthält Zufallswerte zwischen 0 und 1. Da aber auch eine Bewegung in negativer Richtung möglich sein soll, wird von allen Zufallszahlen 0,5 abgezogen. g stellt die Geschwindigkeit dar. Mit ihr werden die Werte der Richtungsmatrix multipliziert, was bei höherer Geschwindigkeit zu einer größeren Änderung der Position führt.

Abstandsabfrage

Als Letztes folgt die Überprüfung, ob sich zwei Individuen näher als der definierte Mindestabstand gekommen sind. Falls dies zutrifft, werden beide als infiziert gekennzeichnet.

 function[Inf]=abstandsabfrage(Nx, Ny, Inf, xNeu, yNeu, radInf);
  (1)   for k=1:Nx   
                for l=1:Ny
  (2)              if Inf(l, k)==1  %Falls ein ein Infizierter gefunden wird...
  (3)                   for m=1:Nx   % Abstandsberechnung aller anderen Objekte zum Infizierten
                            for n=1:Ny
  (4)                           abstand=norm([x(l,k)-x(n,m), y(l,k)-y(n,m)]);
  (5)                           if abstand<radInf
  (6)                               Inf(n, m)=1;
                                end
                            end
                        end
                    end
                 end
        end
 end

Die Informationen, an welcher Position sich ein Infizierter befindet, sind in der Matrix Inf gespeichert. Diese wird Zeilen- und Spaltenweise durchsucht (1). Wenn ein Infizierter gefunden wird (2), werden wieder alle Individuen durchgegangen (3) und jeweils der Abstand zwischen ihnen und der in (2) gefundenen infizierten Person ermittelt (4). In einer letzten if-Bedingung findet die Überprüfung statt, ob zwei Personen den Mindestabstand zueinander unterschritten haben (5). Wenn dies der Fall ist, wird die andere Person ebenfalls als infiziert markiert und kann im nächsten Zeitschritt ebenfalls die Krankheit übertragen.

Ergebnis

In jedem Zeitschritt werden Bilder der aktuellen Situation gespeichert. Zusammengefügt ergeben sich folgende Resultate:

Skript in Cocalc

Kontaktmodell

Fundamentallösung der Poisson- und Diffusionsgleichung

stationäre Fundamentallösung für Laplace

In diesem Abschnitt wird die Fundamentallösung $\Phi$ für die homogene Laplace-Gleichung untersucht und graphisch als Funktionen auf $\mathbb {R} ^{2}$ dargestellt.

Erzeugung eines Gitters

 x = linspace(0,10,100);
 y = linspace(0,10,100);
 [x,y] = meshgrid(x,y);

Um die stationäre Fundamentallösung für Laplace zu implementieren, müssen wir zuerst ein Gitter definieren, welches das Gebiet der Fundamentallösung beschreibt.

Implementierung der stationären Diffusion

 u=@(x,y) -1/(2*pi)*log(sqrt(x.^2+y.^2));
 meshc(x,y,u(x,y))

Die Formel der Fundamentallösung für $\mathbb {R} ^{2}$ lautet daher: $u(x)=\Phi (x)=-{\frac {1}{2\pi }}\ln \|x\|$

Der Befehl meshc liefert dann die graphische Ausgabe, wobei x und y die Koordinatenmatrizen sind und u(x,y) die Wertematrix darstellt. Folglich entsteht das Ergebnis der stationären Laplace-Gleichung:

Skript in Cocalc

Laplace stationär

instationäre Fundamentallösung für Laplace

Nachfolgend wird die Fundamentallösung $\Psi$ der Laplace-Gleichung untersucht und wie zuvor graphisch als Funktion auf $\mathbb {R} ^{2}$ dargestellt. Dabei unterscheidet sich die instationäre Diffusion von der stationären durch die Hinzunahme der Zeit t.

Erzeugung eines Gitters

 x = linspace(-10,10,100);
 y = linspace(-10,10,100);
 [x,y] = meshgrid(x,y);
 a=0.05;  %konstanter Diffusionskoeffizient
 T=20;

Wie zuvor wird zuerst ein Gitter erzeugt, um das Gebiet zu definieren. Allerdings werden hier auch noch die Parameter a und T festgelegt, wobei a den konstanten Diffusionskoeffizienten darstellt und T die Anzahl der Schritte im Prozess.

Implementierung der instationären Diffusion

for t=1:T
   u=@(x,y,t) 1/(4*pi*a*t)*exp(-sqrt(x.^2+y.^2)/(4*a*t));
   figure(t)
   meshc(x,y,u(x,y,t))
 endfor

Da die instationäre Diffusion zeitabhängig ist, muss in der Schleife zuerst der Parameter t definiert werden, der die Zeit für jeden Schritt angibt. Für die homogene Diffusionsgleichung mit konstantem Diffusionskoeffizienten existiert eine Lösungsformel, die wie folgt lautet: $u(x,t)=\Psi (x,t)={\frac {1}{4\pi at}}\exp \left(-{\frac {{\|x\|}^{2}}{4at}}\right)$

Skript in Cocalc

Laplace instationär

stationäre Fundamentallösung der Poissongleichung

Die Poissongleichung $-\Delta u(x)=f(x)$ auf $\mathbb {R} ^{n}$ ist die Lösung der inhomogenen Laplace Gleichung. Die Poissongleichung lässt sich anhand der Faltung der Fundamentallösung $\Phi$ und der Funktion der rechten Seite bestimmen. Die Poissonformel lautet: $u(x)=(\Phi *f)(x)=\int _{\mathbb {R} ^{n}}\Phi (x-y)f(y)dy$ Wichtig hierbei ist, dass es sich um eine Funktion mit kompaktem Träger handeln muss, also eine Funktion, die außerhalb eines Kompaktums nur den Wert 0 annimmt. Für diese Funktion gilt dann $\lim _{\|x\|\to \infty }u(x)\to 0$

Definition des Gebietes

 s=0.1;  
 X=0:s:10;
 Y=0:s:10;
 [x,y]=meshgrid(X,Y);

Zu Beginn muss ein rechteckiges Gebiet in $\mathbb {R} ^{2}$ festgelegt werden, mit einem Gitter, das sich aus den Vektoren X und Y, sowie der entsprechenden Schrittweite ergibt.

Funktion der rechten Seite

Um die Funktion der rechten Seite darzustellen, müssen wir eine Quellfunktion aufstellen. Diese Funktion ist in unserem Modell wie folgt definiert:

 function wert=Quellfunktion(x,y)
 if sqrt((x.-4).^2+(y.-4).^2)<=1 wert=1;
         else wert=0;
 endif
 endfunction

Die Funktionswerte f müssen dann in jedem Gitterpunkt $f(x(j,i),y(j,i))$ in eine Matrix gespeichert werden. Das geschieht in der nachfolgenden Schleife.

 for i=1:(length(X))
  for j=1:(length(Y))
    f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
  endfor
 endfor

Die folgende Grafik zeigt die Darstellung der oben definierten Funktion der rechten Seite.

Implementierung der Poissonformel

 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));
 Phi mit verschobenem Argument:
     phi=-log((sqrt((xstar.-x).^2+(ystar.-y).^2)))/(2*pi);
     phi(j,i)=0;
 numerische Integration
     I(j,i) = trapz(Y,trapz(X,phi.*f,2));
     endfor
 endfor
 close all

Nun müssen wir feste Matrizen $x^{*}$ und $y^{*}$ mit konstanten Einträgen festlegen. Danach werden dann die Werte der Fundamentalfunktion $\phi$ um ein festes $x^{*}$ und $y^{*}$ verschoben und in der Matrix $phi$ gespeichert. Am Schluss wird die numerische Integration mit Hilfe der Trapezregel durchgeführt.

Lösung der Poissongleichung

Als Ergebnis erhalten wir somit die Lösung der Poissongleichung, die in der nächsten Abbildung graphisch dargestellt ist.

Skript in Cocalc

Poisson stationär

instationäre Fundamentallösung der Poissongleichung

Nachfolgend wird die Fundamentallösung $\Psi$ der Poissongleichung untersucht und wie zuvor graphisch als Funktion auf $\mathbb {R} ^{2}$ dargestellt. Dabei unterscheidet sich die instationäre Diffusion von der stationären durch die Hinzunahme der Zeit t.

Definition des Gebietes

 s=0.5; 
 X=0:s:10;
 Y=0:s:10;
 [x,y]=meshgrid(X,Y);
 a=0.1;       
 T=10;     
 dt=0.5;

Zu Beginn muss wie im Modell zuvor ein rechteckiges Gebiet in $\mathbb {R} ^{2}$ festgelegt werden, mit einem Gitter, das sich aus den Vektoren X und Y, sowie der entsprechenden Schrittweite ergibt. Außerdem muss ein konstanter Diffusionskoeffizient festgelegt werden, sowie die Grundlagen für die zeitliche Abhängigkeit. Die Anzahl der Zeitschritte wird durch T gegeben und die Zeitschritte durch dt.

Definition der Anfangsfunktion

Um die Funktion der rechten Seite darzustellen, müssen wir eine Anfangsfunktion aufstellen. Diese Funktion ist in unserem Modell wie folgt definiert:

 function wert=Anfangsfunktion(x,y)
   if 0.3*x+0.4*y<=1 wert=1;
           else wert=0;
   endif
 endfunction

Die Funktionswerte f müssen dann wie zuvor in jedem Gitterpunkt $f(x(j,i),y(j,i))$ in eine Matrix gespeichert werden. Das geschieht in der nachfolgenden Schleife.

 for i=1:(length(X))
     for j=1:(length(Y))
     u0(j,i)=1*Anfangsfunktion(x(j,i),y(j,i));
     endfor
 endfor

Die Funktion der rechten Seite ist in der folgenden Grafik dargestellt.

Implementierung der Poissongleichung

 for t=1:T
    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=(1/(4*pi*a*t))*exp(-((0.3*(xstar.-x)+0.4*(ystar.-y)))/(4*a*t));
         %numerische Integration
         I(j,i,t+1) = trapz(Y,trapz(X,psi.*u0,2));
         endfor
     endfor
     figure (t+1)
     meshc(x,y,I(:,:,t+1))
     title(['Losung von Diffusionsgleichung in t=', num2str(t)])
     axis([0 10 0 10 0 1.2])
 endfor

Aufgrund der Zeitabhängigkeit der instationären Diffusionsgleichung muss zu Beginn der Schleife der Parameter t definiert werden, der die Zeit für jeden Schritt angibt. Für die homogene Diffusionsgleichung mit konstantem Diffusionskoeffizienten existiert eine Lösungsformel, die wie folgt lautet: $u(x,t)=\Psi (x,t)={\frac {1}{4\pi at}}\exp \left(-{\frac {{\|x\|}^{2}}{4at}}\right)$ Nun müssen wir feste Matrizen $x^{*}$ und $y^{*}$ mit konstanten Einträgen festlegen. Danach werden dann die Werte der Fundamentalfunktion $\Phi$ um ein festes $x^{*}$ und $y^{*}$ verschoben und in der Matrix $phi$ gespeichert. Am Schluss wird die numerische Integration mit Hilfe der Trapezregel durchgeführt. Die graphische Lösung der instationären Diffusionsgleichung sieht wie folgt aus:

Skript in Cocalc

Poisson instationär

Kompartimentmodelle für die dynamische Beschreibung der Infektionsverbreitung auf Grundlage der aktuellen Datenlage

Nachfolgend werden verschiedene Modelle vorgestellt, die zum ersten mal auf authentischen Daten beruhen. In Aufgabenteil 1 wird die Infektionsverbreitung nach dem logistischen Wachstumsmodell, einamel als SI- und einmal als SIR-Modell, nachgebildet. Aus den realen Anfangswerten wird eine Infektionsrate bestimmt. Zudem werden die Ergebnisse der Modelle untereinander sowie mit den tatsächlichen Daten verglichen. In Aufgabenteil 2 soll die bisher konstant gehaltene Infektionsrate zeitlich variabel gemacht werden und damit eine möglichst genaue Näherungsfunktion für die vom RKI bereitgestellten Daten erstellt werden. Maßnahmen wie Kontaktbeschränkungen, Maskenpflicht oder auch wieder die Lockerung dieser sollen dafür als zeitlicher Indikator für eine Änderung der Infektionsrate dienen.

Teil 1 - SI-Modell

Das logistische Modell gehört zu den Kompartimentmodellen für die Infektionsverbreitung. Dabei wird die Gesamtpopulation in zwei Gruppen aufgeteilt, nämlich in die Susceptible (Infektionsanfällige) und die Infected (kummulierte Infizierte). Hier wird vorausgesetzt, dass die infizierten Individuen nach der Genesung nicht wieder infektionsanfällig werden.

Festlegung der Parameter

 Gesamtbev=83000;
 B=(2/3)*Gesamtbev;
 I_0=16/Gesamtbev;
 c=0.3219;
 t0=0;
 dt=1;
 times=(0:0.1:100);
 y0=[B;I_0];

Durch die der WHO entnommenen Anfangswerte $I_{0}={\frac {16}{83Mio}}$ erhält man die Konstante $\beta =0,3219$ . Die Infektionsrate c wird also durch $\beta =0,3219$ beschrieben. Als Kapazität $B$ nehmen wir ${\frac {2}{3}}$ der Gesamtbevölkerung Deutschlands an, wobei die Gesamtbevölkerung in Tausend angegeben ist. Die Gesamtpopulation ist somit stets konstant. Der Anfangszustand der Susceptible und der Infected ist durch $S(t0)=B$ und $I(t0)=I_{0}$ gegeben und wird in der Matrix $y0$ gespeichert.

Numerische Lösung

 f=@(y,x) [-c*y(1)*y(2)/B; c*y(1)*y(2)/B];
 %numerische Lösung
 y=lsode(f, y0, times);

Um die Veränderung der Susceptible und der Infected zu bestimmen, werden die jeweiligen Gleichungen benötigt. Diese sind wie folgt gegeben:

$S'(t)=-kS(t)I(t)$

$I'(t)=kS(t)I(t)$

Hier wird deutlich, dass die Gesamtpopulation konstant ist, da die Erhaltungsgleichung $S'(t)+I'(t)=0$ gilt. Somit ist die Summe aus Susceptible und Infectet stets gleich der Konstanten $B$ . Die beiden Differentialgleichungen werden in einer Funktion f gespeichert. Die numerische Lösung des Systems von Differentialgleichungen folgt in Octave über den Befehl "lsode", der Anhand der Ausgangswerte der Susceptible und Infected die Differentialgleichungen löst.

Graphische Darstellung des SI-Modells

Teil 2 - SIR-Modell

Das SIR-Modell betrachtet nun eine zusätzliche Gruppe, nämlich die der Removed. Diese stellen die genesenen Individuen dar, die nach der Genesung von der Gruppe der Infizierten in die Gruppe der Genesenen (und Toten) $R$ wandert. Dieser Welchsel wird mit der Genesungsrate $g$ beschrieben. Insgesamt sind somit die kummulierten Infizierten als Summe der aktuellen Infizierten und der Genesenen darstellbar.

Festlegung der Parameter

 Gesamtbev=83000;
 B=(2/3)*Gesamtbev;
 I_0=16/Gesamtbev;
 R_0=0;
 c=0.3219;   %Infektionsrate
 g=0.07143;  %Genesungsrate
 t_0=0;
 dt=1;
 times=(0:0.1:180);
 y0=[B;I_0;R_0];

Auch hier müssen zuerst einige Parameter festgelegt werden. Diese überschneiden sich größtenteils mit den zuvor festgelegten Parametern des SI-Modells. Als neuer Parameter wird der Anfangszustand der Recovered benötigt, der zu Beginn 0 beträgt. Außerdem werden die Parameter um eine Genesungsrate erweitert. Die Matrix der Anfangszustände enthält nun drei Werte, da auch hier um den Anfangswert der Removed erweitert wurde.

Numerische Lösung

 f=@(y,x) [-c*y(1)*y(2)/B;c*y(1)*y(2)/B-g*y(2);g*y(2)];
 y=lsode(f, y0, times);

Ähnlich wie im vorherigen Modell wird die numerische Lösung des Systems von Differentialgleichungen mit dem Befehl "lsode" durchgeführt, der Anhand der Ausgangswerte der Susceptible, Infected und Recovered die Differentialgleichungen löst. Mithilfe der Differentialgleichungen wird die Funktion $f$ der rechten Seite des Differentialgleichungssystems definiert. Dazu werden die Differentialgleichungen wie folgt definiert:

$S'(t)=-kS(t)I(t)$

$I'(t)=kS(t)I(t)-gI(t)$

$R'(t)=gI(t)$

Auch hier gilt wie zuvor die Erhaltungsgleichung $(S+I+R)'=0$ die besagt, dass die Gesamtpopulation konstant ist.

Graphische Darstellung des SIR-Modells

Teil 3 - Modifiziertes SIR-Modell

In diesem Modell wird davon ausgegangen, dass die Zahl der erfassten Infizierten unvollständig ist und dass die Infektionsrate aufgrund mehrerer Kontaktbeschränkungen zur Verlangsamung der Infektionsausbreitung führt. Dafür werden die vom RKI ermittelten Daten bezüglich der Infektionszahlen benutzt. Das modifizierte SIR-Modell wird dann bestmöglich an die tatsächlichen Daten des RKI angepasst.

Festlegung der Parameter

 NM=55000; % Kapazitätsgrenze, 2/3 der deutschen Bevölkerung
 c=0.3219; %Basisinfektionsrate- die Rate des exponenziellen Wachstums am Anfang der Pandemie
 w=1/14; % Wechselrate zu den Genesenen
 d=0.003;
 %Anteil der erfassten Infizierten
 r=0.1;
 e1=1;
 e2=1.32;
 e3=1.26;
 e4=0.96;
 e5=0.92;

In diesem Modell beschreibt NM die Kapazitätsgrenze der Bevölkerung Deutschlands, die sich tatsächlich ansteckt (in Tausend). Der Faktor c gibt die Basisinfektionsrate an, also die des exponentiellen Wachstums am Anfang der Epidemie. Die tägliche Wechselrate w gibt die Individuen an, die von der Gruppe der Infizierten in die Gruppe der Recovered oder der Toten wechselt. Die Todesrate wird durch d beschrieben. Der durch r festgelegte konstante Faktor gibt den prozentualen Anteil der durch Tests erfassten Infizierten an. Wir gehen also davon aus, dass nur 10% der Erkrankten tatsächlich erfasst werden. Die Parameter e1-e5 werden an späterer Stelle erläutert.

Daten des RKI

Die Funktion 'coronaData_aktuell' enthält die vom RKI ermittelten Daten zu Corona-Infizierten in Deutschland im Zeitraum vom 24.02.2020 - 29.07.2020. Dabei wird unterschieden zwischen kumulierten Infizierten, Recovered und Todesfällen.

 %aktuelle Daten: (passen nicht mehr zur exponentieller Funktion)
  A=coronaData_aktuell();
  IWHO=A(1,:);
  RWHO=A(3,:);
  DWHO=A(2,:);
  IWHOaktuell=IWHO-RWHO-DWHO;  
  n=length(IWHO); %n:=Anzahl der Tage
  timesWHO=[0:1:n-1];

Die aktuellen Infizierten berechnen sich somit mithilfe der vom WHO ermittelten kumulierten Infizierten. Hier müssen die vom WHO erfassten Recovered und Todesfälle von den Infizierten abgezogen werden. Damit erhalten wir die Daten der aktuell von der WHO erfassten Infizierten. Die Daten des RKI werden in der folgenden Grafik veranschaulicht:

Verlangsamung der Infektionsrate

An dieser Stelle wird die zeitabhängige Infektionsrate Infolge der Kontaktbeschränkungen implementiert.

 function val=slowdown(x,t,e1,e2,e3,e4,e5)
 %Tag 1 24.02. ; Tag 21  15.03. Kontaktbeschränkungen; Tag 66 29.04. Maskenpflicht; Tag 78  11.05. Lockerungen der Kontaktbeschränkungen; Tag 
 113 15.06. Öffnung des Schengenraums
 if x<t val=1 ; else  if x<t+21 val=t./(x.^e1); else if x<t+66 val=t./(x.^e2); else if x<t+78 val=t./(x.^e3); else if x<t+113 val=t./(x.^e4); 
 else val=t./(x.^e5);
 endif
 endif
 endif
 endif
 endif
 endfunction

Dafür definieren wir die Slowdownfunktion, die ab Tag t in mehreren Schritten sinken soll. Dafür wurde recherchiert, wann welche Maßnahmen eingetreten sind, die zu einer Verlangsamung der Infektionsrate beigetragen haben. Hier wurde zuerst eine Parametersuche durchgeführt, die passende Parameter für die entsprechenden Abschnitte der Funktion berechnet hat.

Wie genau die Parameter e1 bis e5 ermittelt wurden, folgt später.

Vergleich der Infektionsraten

In diesem Abschnitt werden die Infektionsraten der Slowdownfunktion mit der tatsächlichen Infektionsrate der RKI-Daten verglichen.

 %Slowdown-Funktion
 for i=1:length(timesWHO)
 cc(i)=slowdown2(timesWHO (i),21, e1, e2, e3, e4, e5)*c;   %ab Tag 21 erste Maßnahmen
 endfor 
 plot  (timesWHO,cc)
 hold on
 %infektionsrate aus den RKI-Daten berechnet:
 T=1;
 ccc(1)=cc(1);
 for i=T+1:n
 ccc(i-T)=(IWHO(i)-IWHO(i-T))/(IWHO(i-T));
 endfor

Die tatsächliche Infektionsrate (ccc) muss aus den RKI-Daten berechnet werden. Dazu wird die Differenz der kumulierten Infizierten eines Tages i und der kumulierten Infizierten des vorherigen Tages gebildet. Mit der Division des Ergebnisses durch die kumulierten Infizierten am Vortag ergibt sich die Infektionsrate aus den RKI-Daten für den jeweiligen Vortag. Im Folgenden werden die Raten graphisch dargestellt:

Numerische Lösung des SIR-Systems

In diesem Schritt soll nun das modifizierte SIR-Modell numerisch gelöst werden.

 %Lösung des ODE Systems
 times=(0:0.1:n);
 %Anfangsbedingungen
 yo=[NM-16/1000;16/1000;0];
 f=@(y,x) [-c*slowdown2(x, 21, e1, e2, e3, e4, e5)*y(1)*y(2)/(r*NM);c*slowdown2(x, 21, e1, e2, e3, e4, e5)*y(1)*y(2)/NM-w*y(2);w*y(2)-d*y(2)];
 y = lsode (f, yo, times);

Zuerst müssen die Anfangsbedingungen der Susceptible, Infected und Recovered in einem Vektor gespeichert werden. Wir bereits im vorherigen SIR-Modell wird die numerische Lösung des Systems von Differentialgleichungen mit dem Befehl "lsode" durchgeführt, der Anhand der Ausgangswerte der Susceptible, Infected und Recovered die Differentialgleichungen löst. Die Anfangsbedingungen werden aus den Daten des RKI entnommen. Mithilfe der Differentialgleichungen wird die Funktion $f$ der rechten Seite des Differentialgleichungssystems definiert. Dazu werden die Differentialgleichungen wie folgt definiert:

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

$I'(t)={\frac {c(t)}{NM}}(t)I(t)-wI(t)$

$R'(t)=wI(t)-dI(t)$

Es gilt die Erhaltungsgleichung $(S+{\frac {I}{r}}+{\frac {R}{r}})'=-dI$ die besagt, dass die Gesamtpopulation sich um die verstorbenen Infizierten verringert. Die numerische Lösung wird nun in der Grafik gemeinsam mit den Daten des RKI dargestellt:

Relativer Fehler Infected

Im folgenden wird der relative Fehler der numerischen Lösung des modifizierten SIR-Modells zu den RKI.Daten berechnet. Dabei wird in unserem Modell nur der relative Fehler der Infected betrachtet.

 %Interpolation der Funktionswerte an die RKI-Zeiten (Tage)
 I=interp1 (times, y(:,2), timesWHO+1);
 rel_I=abs(I-IWHOaktuell/1000)./(IWHOaktuell/1000);
 %weiterer Fehlermaß Durchschnitt und Euklid-Norm
 rel_I_mean=mean(rel_I)

Mithilfe dieser Fehlerberechnung wurden die Parameter e1 bis e5 bestimmt. Dabei wurden schlicht alle möglichen Kombinationen von den Parametern zwischen 0,8 und 1,5 ausprobiert und jeweils der relative Fehler berechnet. Daraus ergibt sich dann eine optimale Kombination, die den Fehler minimiert.

 for e1=0.8:0.1:1.5
 for e2=0.8:0.1:1.5
 for e3=0.8:0.1:1.5
 for e4=0.8:0.1:1.5
 for e5=0.8:0.1:1.5
 
   ...........  Lösung des ODE Systems ............
 
 
   ...........  Berechnung des Fehlers ............
 
           if rel_I_mean<rel_I_opt
               rel_I_opt=rel_I_mean
               e1opt=e1;
               e2opt=e2;
               e3opt=e3;
               e4opt=e4;
               e5opt=e5;
           endif
 endfor
 endfor
 endfor
 endfor
 z=z+1
 endfor

Bei dieser Funktion wurden die Parameter in Schritten von 0,1 variiert. Um das Ergebnis weiter zu verfeinern, wurde eine zweite Parametersuche mit kleinerer Schrittweite aber auch kleinerem Bereich (basierend auf den ersten Ergebnissen) durchgeführt:

 for e1=0.9:0.02:1.1
 for e2=1.2:0.02:1.4
 for e3=1.1:0.02:1.3
 for e4=0.9:0.02:1.1
 for e5=0.8:0.02:1
 
 .........
 
 endfor
 endfor
 endfor
 endfor
 z=z+1
 endfor

Der durchschnittliche relative Fehler der Infizierten in unserem Modell beläuft sich mit somit auf 0.13524.

Darstellung Relativer Fehler Infected

Randwertaufgabe für Poissongleichung mit finite Differenzen Methode (FDM)

In dieser Aufgabe soll die Lösung der Poissongleichung mit zwei verschiedenen Randwertproblemen mit Hilfe der finiten Differenzen Methode modelliert werden. Die Grundidee dieser ist, die Ableitungen einer Funktion in einer Differentialgleichung nicht genau zu bestimmen, sondern sich nur äquidistante Punkte (in Abstand h) in dieser Funktion anzuschauen und den Differentialquotienten durch einen Differenzenquotienten zwischen diesen Punkten zu nähern. Hier relevant ist vor allem der zentrale Differenzenquotient $D_{h}f(x):={\frac {f(x+h)-f(x-h)}{2h}}$ , welcher den Mittelwert zwischen der Steigung aus einem Gitterpunkt vor (x-h) und einem Gitterpunkt nach (x+h) dem untersuchten x-Wert bildet.

Da in der Poissongleichung mit dem Laplace-Operator die zweite Ableitung eine Rolle spielt, braucht man den zweiten Differenzenquotienten $D_{h}^{2}f(x):={\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}$ . Dieser kann ganz einfach durch zweifache Anwendung des ersten Differenzenquotienten mit halber Schrittweite hergeleitet werden.

Oben wurde gezeigt, wie man den zweiten Differenzenquotienten herleitet. Allerdings geschah dies nur in einer Raumdimension. Jetzt liegt allerdings ein zweidimensionales Problem vor. Es gibt also je einen zweiten Differenzenquotienten in x- und in y-Richtung.
$D_{h,x}^{2}u(x_{i},y_{j}):={\frac {u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}}$ $D_{h,y}^{2}u(x_{i},y_{j}):={\frac {u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}}$
Der zweite Differenzenquotient in zwei Raumdimensionen folgt durch Addition der beiden oberen 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.$

Dirichlet Randwertbedingungen - inhomogen

Die Dirichlet Randbedingungen sagen aus, dass die Werte einer Funktion $u(x,y)$ mit $x,y\in D$ am Rand $x,y\in \partial D$ durch eine andere Funktion gegeben sind. Im jetzige Fall werden die Randwerte an den Rändern als konstant definiert. Zudem wird der Diffusionskoeffizient festgelegt:

 #Diffusionskoeffizient
 a=1;
 
 wert_unten=0.01;   
 wert_oben=-0.01;
 wert_links=-0.01;
 wert_rechts=0.01;

Danach muss ein zweidimensionales Gitter definiert werden:

 #Definiere Endpunkte
 x_end = 8.71;
 y_end = 5.38;
 
 #Definiere Schrittweite
 hx = 0.1;
 hy = hx;
 
 #Definiere Gitter
 x = [0:hx:x_end];
 y = [0:hy:y_end];
 [X,Y] = meshgrid(x,y);
 
 #Bestimme N & M
 N = floor(x_end/hx+1)
 M = floor(y_end/hy+1)

Dabei werden die Endpunkte der Achsen sowie der Abstand der Gitterpunkte bestimmt. Wenn die Division von Endpunkt durch die Schrittweite nicht nicht eine natürliche Zahl ergibt, sind n und m auch keine natürlichen Zahlen, was später Probleme verursacht. Aus diesem Grund werden die Ergebnisse der Division mit dem Befehl "floor" auf eine natürliche Zahl abgerundet.

Da die Poissongleichung inhomogen ist, gibt es wie in vorherigen Aufgaben eine Funktion f der rechten Seite, die auch auf ähnliche Weise definiert wird.

 for i=1:N
 for j=1:M
   f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
   %Dirichlet Randbedingung unten und Oben:
   if j==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))-wert_links*a/hx^2; endif
   if j==M f(j,i)=1*Quellfunktion(x(j,i),y(j,i))-wert_rechts*a/hx^2; endif
   %Dirichlet Randbedingung links und rechts:
   if i==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))-wert_oben*a/hx^2; endif
   if i==N f(j,i)=1*Quellfunktion(x(j,i),y(j,i))-wert_unten*a/hx^2; endif
 endfor
 endfor
 
 % Ecken-Beiträge:
 f(1,1)=1*Quellfunktion(x(j,i),y(j,i))-(wert_links+wert_oben)/hx^2;
 f(1,N)= 1*Quellfunktion(x(j,i),y(j,i))-(wert_links+wert_unten)/hx^2;
 f(M,1)=1*Quellfunktion(x(j,i),y(j,i))-(wert_rechts+wert_oben)/hx^2;
 f(M,N)=1*Quellfunktion(x(j,i),y(j,i))-(wert_rechts+wert_unten)/hx^2;

Neu sind hier die if-Abfragen, die die oben definierten Dirichlet-Randbedingungen in die Funktion eintragen. Die Ecken des Definitionsbereiches liegen an zwei Rändern, weshalb hier die Randbedingungen separat als die Summe der beiden Werte der anliegenden Ränder bestimmt werden. Dies können die if-Bedingungen nicht leisten.

Um die zu lösende Gleichung $\Delta u(x_{i},y_{j})=f(x_{i},y_{j})$ in ein gewöhnliches lineares Gleichungssystem zu schreiben, müssen die Matrizen $u(x_{i},y_{j})$ und $f(x_{i},y_{j})$ in einen Vektor umgeschrieben werden. Dabei werden beide Matrizen zeilenweisen in einen neuen Vektor einsortiert.

 b=-1*reshape(f',N*M,1);

Das zu lösende Gleichungssystem ist somit $BB*{\vec {u}}={\vec {f}}$

Als nächstes folgt die Erstellung des Approximationsschemas BB gemäß der Gleichung $\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.$ , die den Laplace Operator approximiert Diese Matrix hat auf der auf der ersten und -ersten Nebendiagonale die Einheitsmatrix I sowie eine weitere Matrix B auf der Hauptdiagonale. B wiederum hat auf der Hauptdiegoanle den Wert -4 und auf der ersten und -ersten Nebendiagonale den Wert 1. Zuerst werden die Matrizen B und I erstellt. BB wird auch erstellt, ist hier allerdings noch leer.

 I=eye(N,N);
 B=-4*eye(N,N)+diag(ones(N-1,1),1)+diag(ones(N-1,1),-1);
 BB=zeros(N*M,N*M);

Im nächsten Schritt werden Die Matrizen B und I in die Matrix BB eingesetzt sowie mit a/h_x² multipliziert.

 for i=1:M    %Hauptdiagonale besetzen
   BB((i-1)*N+1:i*N,(i-1)*N+1:i*N)=B;
 endfor
 
 for i=1:M-1    %linke Nebendiagonale
   BB(i*N+1:(i+1)*N,(i-1)*N+1:i*N)=I;
 endfor
 
 for i=1:M-1    %rechte Nebendiagonale
   BB((i-1)*N+1:i*N,i*N+1:(i+1)*N)=I;
 endfor
 
 BB=BB*a/hx^2;

Jetzt sind alle Bestandteile des Gleichungssystems definiert, weswegen dieses nur noch gelöst werden muss.

 sol=BB\b;
 sol_matrix=reshape(sol,N,M);
 sol_matrix=sol_matrix';

Bei der Quellfunktion

 if sqrt((x.-3.25).^2+(y.-3).^2)<=0.35 wert=1;
   else wert=0;
 endif

ergibt sich folgende Lösung:

Neumann Randwertbedingungen

Die Modellierung der Neumann Randbedingungen unterscheidet sich nur in manchen Punkten von der der Dirichlet Randbedingungen.

Die Neumann Randbedingungen lautet $\nabla (u)*{\vec {n}}=g(x)$ (Hier: $g(x)=0$ . In der Randbedingung taucht mit dem Nabla-Operator die Ableitung auf, welche mittels des vorwärts- bzw. rückwärts Differenzenquotienten (je nach dem an welchem Rand) approximiert wird.

linker Rand: $D_{h,x}^{+}u(x_{0},y_{j})={\frac {u_{1,j}-u_{0,j}}{h}}=0$
rechter Rand: $D_{h,x}^{-}u(x_{n+1},y_{j})={\frac {u_{n+1,j}-u_{n,j}}{h}}=0$
unterer Rand: $D_{h,y}^{+}u(x_{i},y_{0})={\frac {u_{i,1}-u_{i,0}}{h}}=0$
oberer Rand: $D_{h,y}^{-}u(x_{i},y_{m+1})={\frac {u_{i,m+1}-u_{i,m}}{h}}=0$

Eine Implementierung von Neumann-Randbedingungen an allen Rändern ist jedoch physikalisch unsinnvoll, weswegen diese nur am linken und rechten Rand implementiert werden. Am oberen und unteren Rand bleiben die Dirichlet Randbedingungen stehen. Somit läuft der Index $j$ jetzt nicht mehr nur von $1$ bis $m$ , sondern von $0$ bis $m+1$ . Aus diesem Grund müssen die Matrizen B und I um Spalten erweitert werden.

 #Definiere Endpunkte
 xend = 8.71;
 yend = 5.38;
 
 #Definiere Schrittweite
 hx = 0.1;
 hy = hx;
 h=hx;
 
 #Definiere Gitter
 x = [0:hx:xend];
 y = [0:hy:yend+2*hy];
 [x,y] = meshgrid(x,y);
 
 a=1;
   
 %Dirichlet Randbedingungen
 oberer_rand=0.01;
 unterer_rand=-0.01;
 N = floor(xend/hx+1);
 
 #auf linkem/rechten Rand Neumann Randbedingungen
 M = floor(yend/hy+1)+2;

Die Erstellung der Systemmatrix funktioniert genauso wie bei den Dirichlet-Randbedingungen. Lediglich die erste und letzte Zeile der mittleren Blockmatrizen müssen zur Approximation der Neumann-Randbedingungen am linken und rechten Rand abgeändert werden.

 I=eye(N,N);
 B=-4*eye(N,N)+diag(ones(N-1,1),1)+diag(ones(N-1,1),-1);
 BB=zeros(N*M,N*M);
 for i=1:M    %Hauptdiagonale besetzen
   BB((i-1)*N+1:i*N,(i-1)*N+1:i*N)=B;
 endfor
 
 for i=1:M-1    %linke Nebendiagonale
   BB(i*N+1:(i+1)*N,(i-1)*N+1:i*N)=I;
 endfor
 
 for i=1:M-1    %rechte Nebendiagonale
   BB((i-1)*N+1:i*N,i*N+1:(i+1)*N)=I;
 endfor
 BB=BB*a/hx^2; % Matrix mit Diffkoeff/h^2 multiplizieren
 
   % Neumann RB für x:  c* partial_x u=0 - einzelne Zeilen der Matrix ersetzen
 
  for i=1:(M-2)
  vector_up=zeros(1,N*M);
  vector_up(N*i+1)=a/h;
  vector_up(N*i+2)=-a/h;
  vector_down=zeros(1,N*M);
  vector_down(N*i+N)=a/h;
  vector_down(N*i+N-1)=-a/h;
  BB(i*N+1,:)=vector_up ;
  BB(i*N+N,:)=vector_down ;
 endfor

Da am Rand nur die ersten Differenzenquotienten der Neumann Randbedingungen, die sich wie oben geschrieben nur aus zwei benachbarten Gitterpunkten berechnen, relevant sind, sind pro abgeänderter Zeile nur zwei Einträge zu finden. Im nächsten Schritt wird die Funktion der rechten Seite implementiert. Dort sind auch die Dirichlet Randbedingungen wiederzufinden.

 for i=1:N
   for j=1:M
     f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
     %Dirichlet Randbedingung unten und Oben:
     if j==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+unterer_rand*a/hx^2; endif
     if j==M f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+oberer_rand*a/hx^2; endif
   endfor
 endfor
   
 b=-1*reshape(f',N*M,1);

Der Lösungsschritt des Systems ist der gleiche wie zuvor:

 sol=BB\b;
 sol_matrix=reshape(sol,N,M);
 sol_matrix=sol_matrix';

Epedimiologische Modellierung durch numerische Diskretisierung des Reaktionsdiffusionsprozesses

In diesem Abschnitt wird die Reaktionsdiffusiongleichung mit bevölkerungsdichteabhängigen Diffusionskoeffizienten betrachtet:

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

Zuerst müssen notwendige Parameter definiert werden. Diese beinhalten Werte für das betrachtete Gebiet, Basisinfektionsrate und -diffusionskoeffizient, Zeitparameter zur Berechnung der Lösung sowie die in Aufgabe 3 bestimmten Parameter der slowdown-Funktion.

 % Startparameter Gebiet
 OstWest=10; %in 20km
 NordSud=10; %in 20km
 N=40;
 hx= OstWest/(N+1);
 hy= hx;
 M=floor((NordSud/hy))-1;
 [x,y]=meshgrid(0:hx:OstWest,0:hy:NordSud);
 N= N+2;
 M= M+2;
 
 a_basis =0.1;          % Diffusionskoeffizient
 c_basis=0.3219;        % Basisinfektionsrate
 k=0.5;                 % Proportionalitätskonstante bevölkerungsdichteabhängiger Diffusionskoeffizient   
 
 % Zeitparameter
 T=100;                          %Zeitintervall in Tagen
 delta_t=0.1;                  %Zeitschritt
 Nr_timesteps=floor (T/delta_t);
 
 % Parameter Slowdownfunktion
 e1=1;
 e2=1.32;
 e3=1.26;
 e4=0.96;
 e5=0.92;

Als nächstes wird die Matrix der Bevölkerungsdichte bestimmt.

 for i=1:N;
   for j=1:M;
     B(j,i) =Bevoelkerungsdichte(1+floor(j/M*(NordSud-1)),1+floor(i/N*(OstWest-1)));
   endfor
 endfor
 B=1*reshape(B',N*M,1);

Die Funktion "Bevoelkerungsdichte" sieht folgendermaßen aus und ist auf die Einheit 1000Menschen/km² skaliert:

 function val=Bevoelkerungsdichte(x,y);
   Bd=[91300 834200 94255 90280 62250 120115 97850 117220 78515 47030,  
   276060 510890 508725 378685 324115 81600 97315 69795 93490 30470,
   501140 969720 1084940 793545 205940 120870 53165 40515 31435 52425,
   545350 848445 353345 335845 196005 64635 51000 51860 47800 44360,
   216805 546320 606705 128440 142680 78495 30790 25735 21250 59200,
   176035 1130850 403645 107065 117745 80090 27660 35145 32190 43375,
   138030 406515 786840 111695 103125 180895 32610 39995 65370 50510,
   144880 167845 519990 118455 70720 115955 78165 91930 109360 34930,
   99185 78060 214615 76610 59656 59810 74505 151230 60890 22660,
   32900 30000 77070 137725 66245 90385 77065 159960 68315 23245];
 
   val=Bd(x,y)*(1/400)*(1/1000);  %Bevölkerungsdichte in 1000 Menschen pro km² 
 endfunction

Die Daten der Matrix stellen einen Ausschnitt Deutschlands dar, der wie folgt aussieht:

Um die Einwohnerzahlen der einzelnen Gebiete zu ermitteln, wurden Infos der statistischen Ämter des Bundes und der Länder hinzugezogen, insbesondere die interaktive Karte "Zenusatlas". In dieser interaktiven Karte ist es möglich, ein Gebiet auszuwählen und danach folgen automatisch die jeweiligen Bevölkerungszahlen. Die Seite ist über folgenden Link erreichbar: https://atlas.zensus2011.de/.

Weiterhin muss auch die Startverteilung an Infizierten gegeben sein

 % Anfangskonzentration an Infizierten
 for  i=1:N;
   for j=1:M;
    initial(j,i)=0.01*initialfunction(x(j,i),y(j,i));   
   endfor;
 endfor;

Die Anfangsverteilung der Infizierten wird durch folgende Funktion beschrieben:

function wert=initialfunction(x,y)
 wert=0;
 r1=0.35;
 r2=0.15;
 if sqrt((x.-6.25).^2+(y.-3).^2)<=r1 wert=1/(pi*r1^2);
   elseif  sqrt((x.-2.25).^2+(y.-2).^2)<=r2 wert=0.5/(pi*r2^2);
 endif
endfunction

Sie weist zwei Runden Hotspots mit dem Radius 0,35*20km=7km und 0,15*20km=3km die Infiziertendichte 1 und 0,5 zu. Durch Multiplikation dieser Werte mit 0.01 im vorangegangenen Schritt ergeben sich somit in diesen Hotspots Dichten von 10 bzw 5 Infizierten pro 1000 Einwohnern. Auch diese Funktion muss wieder in einen einzelnen Vektor geschrieben werden.

 % anfängliche Infiziertendichte als Vektor schreiben
 sol_old = 1*reshape(initial',N*M,1);

Die Funktion $f(x)$ der rechten Seite beschreibt, wie viel "Material" dem System zugeführt wird. Bisher war diese Funktion entweder 0 oder konstant. Es ist allerdings sinnvoll anzunehmen, dass es bei einer höheren Dichte an Infizierten mehr Ansteckungen gibt als sonst, weswegen in diesem Schritt die Funktion $f(x)$ von der Infiziertendichte abhängig gemacht wird. Es wird also zu $f(u(x,t))$ . Die genaue Funktionsvorschrift ergibt sich aus dem logistischen Wachstumsmodell und lautet $f({\vec {u}}(t))\approx F({\vec {u}}_{h}(t)):={\frac {c}{\vec {B}}}*{\vec {u}}_{h}(t)*({\vec {B}}-{\vec {u}}_{h}(t))$ .

Implementiert wird dies folgendermaßen:

 F=@(U_old,t) c_basis*slowdown2(x,25,e1,e2,e3,e4,e5)*U_old./B.*(B-U_old);

Die Infektionsrate ist hier durch $c_{basis}*slowdown2(x,25,e1,e2,e3,e4,e5)$ wie in Aufgabe 3 beschrieben. Der Reaktionsterm wird verwendet, um die Lösung des nächsten berechneten Zeitschrittes zu ermitteln. $U_{old}$ ist deswegen die Lösung des Zeitschrittes vor des aktuell zu berechnenden, welche mit ( $sol_{old}$ ) bezeichnet wurde.

Als nächstes gilt es den anfänglichen Diffusionskoeffizienten zu bestimmen. Da dieser nun von den Koordinaten x und y abhängig ist, kann man ihn in der Diffusionsgleichung : $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 )}$ nicht mehr vor die Divergenz ziehen, weswegen das Approximationsschema angepasst werden muss. Nach Ausführung der partiellen Ableitungen, welche durch den zentralen Differenzenquotienten approximiert werden, ergibt sich folgendes Approximationsschema: $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 )}.$ Ein neues Problem besteht darin, dass durch die nur einfache Anwendung der partiellen Ableitung auf a(x,y), welche durch den zentralen Differenzenquotienten approximiert wird, Schritte mit halber Schrittweite h entstehen. Es werden somit Zwischenwerte von a(x,y) benötigt. Um diese Zwischenwerte zu erhalten, wird eine Matrix der doppelten Größe von u(x,y) erstellt, welche an jeder zweiten Stelle in x- und y-Richtung mit neuen von der Bevölkerungsdichte linear abhängigen Werte für a gefüllt wird. Die freien Zwischenwerte werden mit Durchschnittsbildung der umliegenden Werte genähert.

 a = a_function(k,b,a_basis,b_mittel,N,M);

 function val = a_function(k,b,a_basis,b_mittel,N,M)
   c=zeros(2*M,2*N);
   #a auf jeden zweiten Gitterpunkt von c verteilen
   for i=1:N
     for j=1:M
       a(2*j,2*i)=a_basis+k*(b(j,i)-b_mittel);
     endfor
   endfor
 
   #Zwischenwerte bestimmen als Mittelwert der umliegenden
   #Erste Zeile bestimmen
   for i=1:N
     a(1,2*i) = a(2,2*i);
   endfor
   for i=2:N
     a(1,2*i-1) = 1/2*(a(1,2*i-2)+a(1,2*i));
   endfor
 
   #Erste Spalte bestimmen
   for j=1:M
     a(2*j,1) = a(2*j,2);
   endfor
   for j=2:M
     a(2*j-1,1) = 1/2*(a(2*j-2,1)+a(2*j,1));
   endfor
 
   #rechte ober Ecke
   a(1,1)=1/2*(a(1,2)+a(2,1));
 
   #Zeilen auffüllen
   for j=1:M 
    for i=1:N-1
      a(2*j,2*i+1)=1/2*(a(2*j,2*i)+a(2*j,2*i+2));
    endfor
   endfor
 
   for i=2:2*N 
     for j=1:M-1
       a(2*j+1,i)=1/2*(a(2*j,i)+a(2*j+2,i));
     endfor
   endfor
 
   val = a;
 endfunction

Gemäß des obigen Approximationsschemas muss eine neue Systemmatrix $A_{h}$ erstellt werden. Dies wird mit der Funktion BB_function gemacht.

 #Block-Tri-Diagonale Matrix
 function val = BB_function(a,h,N,M)
 BB=zeros(N*M,N*M);
 
 #oben links
 BB(1:N,1:N) = -h*I_Eins_Matrix(j,a,h,N,M);
 BB(1:N,N+1:2*N) = h*I_Eins_Matrix(j,a,h,N,M);
 
 #unten rechts
 BB((M-1)*N+1:N*M,(M-1)*N+1:N*M) = -h*I_M_Matrix(j,a,h,N,M);
 BB((M-1)*N+1:N*M,(M-2)*N+1:N*(M-1)) = h*I_M_Matrix(j,a,h,N,M);
 
 #Hauptdiagonale
 for j=2:M-1
   BB((j-1)*N+1:j*N,(j-1)*N+1:j*N) = B_j_Matrix(j,a,h,N,M);
 endfor
 
 #linke Nebendiagonale
 for j=1:M-2
   BB(j*N+1:(j+1)*N, (j-1)*N+1:j*N) = I_l_j_Matrix(j,a,h,N,M);
 endfor
 
 #rechte Nebendiagonale
 for j=2:M-1
   BB((j-1)*N+1:j*N,j*N+1:(j+1)*N) = I_r_j_Matrix(j,a,h,N,M);
 endfor
 # Matrix mit Diffkoeff/h^2 multiplizieren
 val = BB/h^2;
 endfunction

Dabei werden alle benötigten Blockmatrizen in die Matrix BB eingesetzt, welche wieder durch eigene Unterfunktionen bestimmt werden:

 # I_0 Matrix erstellen
 function y = I_Eins_Matrix(j,c,h,N,M)
   I_Eins = zeros(N,N);
   for i=1:N
     I_Eins(i,i) = c(2*1,2*i);
   endfor
   y = I_Eins;
 endfunction

 # I_m+1 Matrix erstellen
 function y = I_M_Matrix(j,c,h,N,M)
   I_M = zeros(N,N);
   for i=1:N
     I_M(i,i) = c(2*M,2*i);
   endfor
   y = I_M;
 endfunction

 #I_r_j Matrix erstellen
 function y = I_r_j_Matrix(j,c,h,N,M)

I_r_j = zeros(N,N);

   for i=2:N-1
     I_r_j(i,i) = c(2*(j+1/2),2*i);
   endfor
   y = I_r_j;
 endfunction

 # Î_l_j Matrix erstellen
 function y = I_l_j_Matrix(j,c,h,N,M)
   I_l_j = zeros(N,N);
   for i=2:N-1
     I_l_j(i,i) = c(2*(j-1/2),2*i);
   endfor
   y = I_l_j;
 endfunction

 #Gibt uns die Matrix B_j aus
 function y = B_j_Matrix(j,c,h,N,M)
 	
   # B_j Matrix erstellen
   B_j = zeros(N,N);
   #oberer linker Eintrag
   B_j(1,1) = -h*c(2*j,2*1);
   #unterer rechter Eintrag
   B_j(N,N) = -h*c(j,N);
   #untere Nebendiagonale
   #oben rechts
   B_j (N,N-1) = h*c(2*j,2*N);
   for i=1:N-2
     B_j(2+i-1,1+i-1) = c(2*j,2*(i+1/2));
   endfor
 
   #obere Nebendiagonale
   #oben links
   B_j(1,2) = h*c(2*j,2*1);
   for i=2:N-1
     B_j(1+i-1,2+i-1) = c(2*j,2*(i+1/2));
   endfor
 
   #Hauptdiagnonale
   for i=2:N-1
     B_j(i,i) = -(c(2*j,2*(i-1/2))+c(2*j,2*(i+1/2))+c(2*(j+1/2),2*i)+c(2*(j-1/2),2*i));
   endfor
   y = B_j;
   endfunction

Jetzt sind alle Ausgangswerte definiert und das Gleichungssystem kann gelöst werden. Die Lösung des nachfolgenden Zeitschrittes ergibt sich durch die Lösung vorherigen Zeitschrittes plus die Änderung im Zeitschritt:

{\vec {u}}_{h}(t_{k+1})={\vec {u}}_{h}(t_{k})+{\frac {d}{dt}}{\vec {u}}_{h}(t)

Die Diffusionsgleichung:

\partial _{t}u(x,t)=div\ (c(B(x,y)\nabla u(x,t))+f(u(x,t))

wird genähert durch:

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

womit sich die Gleichung für den nächsten Zeitschritt ergibt:

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

 for i=1:Nr_timesteps
   i
   sol= sol_old + BB_function(a,hx,N,M)*sol_old*delta_t + F(sol_old,i*delta_t)*delta_t;
   sol_old=sol;
   matrix_solution(:,i)=sol;
 
   # neue Basisinfektionsmatrix aufstellen für nächsten Zeitschritt
   a = a_function(k,b,a_basis,b_mittel,N,M);
 endfor

Die Lösungen für jeden Zeitschritt werden in der Matrix matrix_solution spaltenweise gespeichert.

Grafische Lösung der Reaktionsdiffusionsgleichung

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 (siehe auch Octave-Tutorial)

Literatur

Notieren Sie hier die von Ihnen verwendete Literatur

Boto von Querenburg (2013) Mengentheoretische Topologie - Springer-Lehrbuch