Diese Seite ist noch in Bearbeitung, und wird noch komplett geändert

Gruppenseite - VT

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

COVID-19/Mathematische Modellbildung

Teilnehmer-innen

Toker, Merve
Vatter, Melanie

Ziel der Veranstaltung

Das Modul M11: Entwicklung der Mathematik in Längs- und Querschnitten unterteilt in zwei Bereiche:

diskrete Prozesse (Herr Niehaus)
Kontinuierlichen Diffusionsprozesse (Frau Hundertmark)

Beide dienen der Modellierung der COVID-19-Pandemie.

In beiden Bereichen wurden zeitliche und räumliche Betrachtungen vorgenommen.

Für beide Teilbereiche sollte während eines Lehrveranstaltung Aufgaben bearbeitet werden.

Die Aufgaben zu den diskreten Prozessen sind hier zu finden^[1].

Die detaillierten Aufgabenstellungen zu den kontinuierlichen Diffusionsprozessen sind hier einzusehen ^[2].

In einigen Beispielen werden Beispielrechnungen für die Stadt Mainz gemacht. Hierbei ist zu beachten, dass wir nur das "Zentrum" der Stadt betrachten, also nur die Stadtteile: Altstadt, Neustadt, Oberstadt und Hartenberg-Münchfeld. Wenn also im Folgenden Einwohnerzahlen, Infektionsraten oder Bevölkerungsdichten genannt werden, beziehen sich diese nur auf die soeben genannten Stadtteile.

Diskrete Diffusionsprozesse (Niehaus)

SIR Modell

Die erste Aufgabe bestand darin, einen Zusammenhang zwischen S, I und R (Susceptible, Infected und Recovered) darzustellen. Hierfür haben wir das SIR Modell gewählt. Um die Entwicklung der Infektionszahlen am Beispiel der Stadt Mainz darzustellen wird das SIR Modell über Funktionen angewandt.

Dafür wurden zuerst Anfangswerte generiert. Dabei wurde Bevölkerung B, die Basisinfektionsrate beta, die Genesungsrate gamma und die Todesrate d definiert. Dabei besteht die Bevölkerung B nur aus zwei Drittel der gesamten Mainzer Bevölkerung dar. Denn bei einer Epidemie geht man meistens davon aus, dass sich nur zwei Drittel aller Menschen infizieren (können).

Des Weiteren werden die Anzahl der Tage, für die das Modell laufen soll generiert.

Anschließend wurden die Anfangswerte der Funktionen und die Funktionen für Folgetage für S (Susceptible), I (Infected), R (Removed) und kI (kummulierte Infizierte) generiert.

Zum Schluss wurden die Funktionen geplottet und wir erhielten das Bild, was unter den Implementierungen von Octave zu sehen sind.

  B = 141600;    #Kapazitaetsgrenze, 2/3 der Mainzer Bevölkerung (212400)
  beta = 0.29;   #Basisinfektionsrate (die Rate des exponentiellen Wachstums am Anfang der Pandemie)
  gamma = 0.1;   #Genesungsrate (Wechselrate der Infizierten in die Gruppe der Genesenen)
  d = 0.003;     #Todesrate

  N_t=125;
  t=linspace(0,N_t,N_t+1);

  S(1)=B-I(1);
  I(1)=50;       #Am Anfang der Modellierung gibt es 50 Infizierte in Mainz 
  R(1)=0;        #Am Anfang der Pandemie ist noch niemand Genesen

  for n=1:N_t;
  S(n+1)=S(n)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)=R(n)+gamma*I(n)-d*I(n);
  D(n+1)=d*I(n);
  kI(n+1)=I(n+1)+R(n+1)+D(n+1);
  endfor

  hold on
  plot(t,S,'b');
  plot(t,I,'m');
  plot(t,R,'g');
  plot(t,D,'c');
  plot(t,kI,'r');
  axis([0 130 0 142000]);
  legend("S","I","R","D" "kI","location","east")
  hold off

Beim Plot der Inzizierten entstand eine Glockenkurve mit dem Höhepunkt nach ca. 48 Tagen. Das heißt, dass nach 48 Tagen die Infektionen zurück gehen und sich immer weniger Menschen infizieren. Die Infektion wäre nach etwa 120 Tagen überstanden.

Das Modell geht davon aus, dass die Menschen nach einer überstandenen Infektion immun gegenüber des Virus sind.

SIR modifiziert

Dies ist aber in der Realität nicht der Fall. Es sind über die Corona-Infektion noch keine Langzeitstudien bekannt, allerdings sagt man pauschal, dass man wenn man eine Infizierung überstanden hat, man etwa ein halbes Jahr (ca. 183 Tage) immun ist. Daher haben wir das SIR nochmal neu modelliert, mit der Annahme, dass die Genesenen nach ca. 183 Tagen wieder zurück in die Gruppe der Anfälligen wechseln. Dabei wird unsere neue Kurve für R nicht mehr die kummulierten Infizierten anzeigen, da für die Berechnung der S, immer die R der 183 Tage vorher abgezogen werden müssen.

Der Code dafür sieht nun folgendermaßen aus:

  for n=1:N_t;
  S(n+1)=S(n)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)=R(n)+gamma*I(n)-d*I(n);
  kI(n+1)=I(n+1)+R(n+1);

   if n>183 && S(n)+R(n-183)<=B && R(n)-R(n-183)>=0; #Bedingung benötigt da sonst Modell negative Werte bekommt 
   S(n+1)=S(n)+R(n-183)-(beta/B)*S(n)*I(n);
   I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
   R(n+1)= R(n)-R(n-183)+gamma*I(n)-d*I(n);
   endif
 
  endfor

Nahm man, wie bei dem ursprünglichen SIR-Modell an, dass nach einer Infektion eine vollständige Immunisierung stattfindet, war die Infektion nach etwa 120 Tagen überstanden. Wenn man nun beachtet, dass man nach einem halben Jahr erneut wieder infiziert werden kann, dauert die Pandemie laut diesem Modell nach etwa 6000 Tagen erst überstanden (also erst nach etwa 26 Jahren und 4 Monaten, Ende Juli 2036).

Hoffen wir mal, dass die Modellierung nicht wahr wird.

SIRV-Modell

Des Weiteren werden seit Dezember 2020 Menschen gegen das Corona-Virus geimpft. Da man davon ausgehen kann, dass die Impfungen die Pandemie schwächen, haben wir zusätzlich eine Gruppe V für Vaccinated (geimpft) angelegt. Um diese zu implementierten waren einige Vorüberlegungen notwendig

Startdatum Impfungen + 14 Tage (Wirkung)
Impfquote

Die Impfungen starteten am Tag 305 unseres SIR-Modells, zum ersten Mal waren am 15.01.21 Menschen vollständig geimpft, mit der 14-Tage Wirkungszeit ist also unser erster Tag der Auswirkungen der Impfung der 29.01.21 (Tag 339)

Wie berechnet man die Impfquote?

Die Zahlen der Impfungen in Deutschland wurden vom Robert-Koch-Institut zur Verfügung gestellt. Diese haben wir uns als Graph ausgeben lassen:

Anhand der ersten und letzten bereitgestellten Zahl des RKI haben wir eine Exponentialfunktion berechnet (siehe blauer Graph), die den Daten ähneln sollte. Anfang und Endpunkt waren zwar gleich, die anderen Punkte aber eher weniger. Daher haben wir uns dazu entschieden die Funktion der Impfungen zu splitten. Wie auf dem obigen Bild gut zu erkennen ist, ähnelt der Anfang des Graphen eher einer linearen Funktion und wird erst zum Ende hin ähnlich einer Exponentialfunktion.

So haben wir die Werte aufgeteilt in die ersten 103 Tage der Imfpungen (lineare Funktion) und die danach. Dafür haben wir die Werte aber erstmal für die Stadt Mainz anteilsmäßig runtergerechnet:

Anzahl_Impfungen_Mainz = Anzahl_Impfungen Deutschland : 83 Mio. x 212400

Für die ersten 103 Tage der Impfung ergab sich folgendes Diagramm, für welches wir Excel die Lineare Funktion haben ausrechnen lassen:

Anschließend haben wir für die restliche Funktion die Wachstumsrate selbst bestimmt durch den Tag 104 und den letzten gegeben Tag. Wir kamen auf eine Wachstumsrate von delta=0.027

Unsere Implementierung sah dann folgendermaßen aus:

  B = 141600;    #Kapazitaetsgrenze, 2/3 der Mainzer Bevoelkerung (212400)
  beta = 0.29;          #Basisinfektionsrate- die Rate des exponenziellen Wachstums am Anfang der Pandemie
  gamma = 0.1;         #Genesungsrate
  d = 0.003;             #Todesrate
  delta=0.027;
  N_t=700;
  t=linspace(0,N_t,N_t+1);

  S(1)=B-I(1);
  I(1)=5; 
  R(1)=0;
  V(1)=0;
  V(339)=1.56;

  for n=1:N_t;
  S(n+1)=S(n)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)=R(n)+gamma*I(n)-d*I(n);
  D(n+1)=d*I(n);
  V(n+1)=V(n);

  if n>183 && n<=339 && S(n)+R(n-183)<=B && R(n)-R(n-183)>=0;
  S(n+1)=S(n)+R(n-183)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)= R(n)-R(n-183)+gamma*I(n)-d*I(n);
  D(n+1)=d*I(n);
  V(n+1)=V(n);
  endif

  if n>339 && n<=443 && S(n)+R(n-183)<=B && R(n)-R(n-183)>=0;
  S(n+1)=S(n)+R(n-183)-V(n)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)=R(n)-R(n-183)+gamma*I(n)-d*I(n);
  D(n+1)=d*I(n);
  V(n+1)=V(n)+141.89;
  endif

  if n>443 && n<=N_t && S(n)+R(n-183)<=B && R(n)-R(n-183)>=0;
  S(n+1)=S(n)+R(n-183)-V(n)-(beta/B)*S(n)*I(n);
  I(n+1)=I(n)+(beta/B)*S(n)*I(n)-gamma*I(n);
  R(n+1)=R(n)-R(n-183)+gamma*I(n)-d*I(n);
  D(n+1)=d*I(n);
  V(n+1)=V(n)+delta*(B-V(n));

Graphische Darstellung

Man kann sehen, dass die Impfung in unserem Modell positive Auswirkungen hat. Denn die Pandemie dauert "nur" noch 500 Tage. Das heißt die Pandemie wäre am 9.Juli 2021 vorbei. Hoffen können wir ja mal....

Verbesserungsideen

Wäre die Infektionsrate kleiner gewählt worden, wäre der Effekt der Immunisierungsaufhebung deutlicher geworden
Man könnte zusätzlich noch eine Slowdownfunktion (Kontakbeschränkungen, Maskenpflicht, Aufhebungen etc.) einbringen

SIR-Modell mit Bewegung

Diese Modellierung basiert auf den Implementierungen der Gruppe 6.

Das soeben vorgeführte SIR-Modell ist statisch. Das heißt, es findet keine Bewegung statt und es interessiert auch nicht, wo die Menschen sind. 1 Infizierter steckt im Schnitt während seiner Krankheit 2,9 weitere Menschen an (beta:gamma=2,9) egal wo sich diese befinden.

Wir wollen uns das Ganze jetzt einmal mit möglichen Bewegungen ansehen. Damit die Zahlen nicht so groß werden, teilen wir die Bevölkerungszahl durch tausend. Somit müssen die Ergebnisse der Modellierung am Ende mit Tausend multipliziert werden um realistische Werte zu betrachten. Bei diesem Modell, haben wir zur Vereinfachung die Todesrate nicht beachtet.

  B = 2/3*212.4;    #Kapazitaetsgrenze, 2/3 der Mainzer Bevölkerung
  beta = 0.29;      #Basisinfektionsrate- die Rate des exponenziellen Wachstums am Anfang der Pandemie
  gamma = 0.1;      #Genesungsrate
  T = 130;          #Anzahl der Tage
  dt = [1:T];

Erstellung Schachbrettmuster

Um eine räumliche Modellierung erstellen zu können, erstellen wir zuerst ein Schachbrettmuster. Also ein Gitter, auf dem Menschen verteilt sind.

  xend    = 12; %in km    #Seitenlänge Schachbrett 11km
  yend    = 10; %in km    #Seitenlänge Schachbrett 9km da erst ab 1 begonnen wird
  X       = 12 ; 
  hx      = 1;
  hy      = hx;
  Y       = 10;
  h       = hx;
  [x,y]   = meshgrid(1:hx:12, 1:hy:10);

  q = 3/4;                #Anteil der pro Zeitschritt von einer Zelle in benachbarte Zellen der Moore-Matrix wandert

Zusätzlich werden Matrizen erstellt, in denen später die Werte für S, I und R gespeichert werden können. Da wir noch keine Funktionen erstellt haben, sind diese zur Zeit alle noch mit $0$ -en gefüllt:

  S = zeros(Y,X,T);         #Verteilung der Gesamtbevölkerung Bd auf dem Schachbrett (Susceptible)
  I = zeros(Y,X,T);         #Infected
  R = zeros(Y,X,T);         #Recovered

Bevölkerungsmatrix

Des Weiteren wird eine Bevölkerungsmatrix erstellt, die anschließend auf dem Schachbrett verteilt wird. Die Bevölkerungsmatrix ist entstanden aus der Bevölkerungsdichte pro Quadratmeter in den einzelnen Stadteilen. Auch bei diesen Zahlen wird in Tausendsteln gearbeitet, da die exakten Zahlen von Octave aufgrund der Größe nicht berechnet werden konnten.

Die Bevölkerungsmatrix ist folgendermaßen entstanden:

Die Bevölkerungsmatrix in Octave sah dann wie folgt aus:

  function val = Bd(x,y);
  # Leute pro qkm * 1000
  density = [  2 2	8 8	7 7	4 3.5 3.5 1 1;
               2 2	8 8	7 4 4 3.5 3.5 1 1;
               2 2	2 5	4 4	4 1.1 1.1 1 1;
               3 3	5 5	5 2	2 1.1 1.1 1 1;
               3 3	3 3	2 2	2 1.1 1.1 1.1 1;
               3 3	3 2	2 2	2 1.1 1.1 1.1 1.1;
              1.5 1.5	1.5	1 1	1.5	2.0	2.0 1.1 1.1 1.1;
              1.5 1.5 1.5 1.0 2.5 1.5 0.6 0.6 0.6 0.6 0.6;
              1.5 1.5 1.5 1.5 2.5 1.5 0.6 0.6 0.6 0.6 0.6;
              ];
   val = density(x,y);
   endfunction

Bestimmt man anschließend, dass alle sich in der Matrix Bd(x,y) befindlichen Menschen Susceptible sind, kann man sich das Graphisch wie folgt ausgeben lassen:

  for i=1:X;
  for j=1:Y;
   S(j,i,1) = Bd(j,i);
  endfor
  endfor
  S(:,:,1)=flipud(S(:,:,1));

  #plot Succeptible Startkonfiguration
  figure(1)
  surfc(S(:,:,1))
  xlabel('x-Achse')
  ylabel('y-Achse')
  colormap ("jet")
  colorbar
  axis ([1 xend 1 yend 0 max(max(S(:,:,1)))])
  view(0,90)
  title(["Bevoelkerungsmatrix"])

Entweder 1-Dimensional: oder 3-Dimensional:

Die Colorbar an den Seiten gibt an wie viele Einwohner pro Quadratkilometer sind. Die Zahl muss man für die Realität noch mit 1000 multiplizieren.

Infizierte

Anschließend platzieren wir die ersten Infizierten. Hierfür arbeiten wir wieder in Tausendstel und setzen an drei Orten Infizierte die wir auch von den Susceptible/Anfälligen abziehen

   I(5,8,1) = 5/1000;  #5 Infizierte in Hartenberg (Bsp. Schule)
   S(5,8,1) = S(5,8,1)-I(5,8,1);
   I(6,9,1) = 11/1000; #11 in Altstadt (Bsp. Einkaufszentrum)
   S(6,9,1) = S(6,9,1)-I(6,9,1);
   I(8,4,1) = 9/1000; # 9 Infizierte in Bretzenheim (Bsp. Bushaltestelle)
   S(8,4,1) = S(8,4,1) - I(8,4,1);

Transportprozess

Moore-Matrix

Als nächstes werden anhand der Moore-Nachbarschaft eine Matrix erstellen, welche die Anzahl der nächsten Nachbarn enthält.

"Die Moore-Nachbarschaft ist eine Nachbarschaftsbeziehung in einem quadratischen Raster. Alle Flächen, welche mindestens eine Ecke mit der Basisfläche gemeinsam haben, gelten als Nachbarn. Sie ist nach Edward F. Moore benannt und wird auch als 8er-Nachbarschaft bezeichnet. Sie entspricht den Zugmöglichkeiten des Königs beim Schach"^[3].

Das bedeutet für unsere Moore-Matrix, dass alle Punkte die im inneren der Matrix liegen, jeweils $8$ Nachbarn haben. Die Punkte, die auf den Rändern liegen, haben nur $5$ Nachbarn. Die Punkte in den Ecken haben nur $3$ Nachbarn. Somit sieht unser Moore-Matrix folgendermaßen aus:

$MM={\begin{pmatrix}3&5&5&5&5&5&5&5&5&5&3\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\5&8&8&8&8&8&8&8&8&8&5\\3&5&5&5&5&5&5&5&5&5&3\\\end{pmatrix}}$

Die dazugehörige Octave-Implikation sieht folgendermaßen aus:

  MM = ones(Y,X)*8; #Zuerst werden alle Einträge der Matrix auf 8 gesetzt              
  #Randpunkte haben 5 Nachbarn, daher werden alle diese auf 5 gesetzt.
  for i=1:X                      
  MM(1,i)=5;
  MM(Y,i)=5;
  endfor
  for j=1:Y                      
  MM(j,1)=5;
  MM(j,X)=5;
  endfor
  #Eckpunkte haben 3 Nachbarn, somit werden diese auf 3 festgelegt
  MM(1,1)=3;                     
  MM(1,X)=3;
  MM(Y,1)=3;
  MM(Y,X)=3;

Matrizen zum Zwischenspeichern

Anschließend werden 3 Matrizen zum zwischenspeichern der Werte erstellt, welche bei Bewegung entstehen. Auch diese sind zu Anfang mit $0$ -en gefüllt.

  ZS = zeros(Y,X);        #Zwischenspeicher Succeptible
  ZI = zeros(Y,X);        #Zwischenspeicher Infected
  ZR = zeros(Y,X);        #Zwischenspeicher Recovered

Funktionen zur Bewegung

Am Anfang der Implementierung haben wir festgelegt, dass sich $3/4$ der Menschen in den jeweiligen Matrizen bewegen. Nun haben wir je eine Funktion erstellt, die angibt, wie viele Menschen sich wohin bewegen, und wie viele Menschen bleiben.

Diese sehen folgendermaßen aus:

Bewegungsfunktion:

  function wert=anteil_Bewegung(l,k,j,i,q,t,MM,Z,SIR)
  wert=Z(l,k)+q/MM(j,i)*SIR(j,i,t-1);
  endfunction

keine_Bewegung-Funktion:

  function wert=anteil_keine_Bewegung(j,i,q,t,Z,SIR)
  wert=Z(j,i)+(1-q)*SIR(j,i,t-1);
  endfunction

Nun haben wir alle erforderlichen Schritte durchgeführt, um einen Transportprozess durchlaufen zu lassen.

Es werden nun für alle Werte, die innerhalb des Schachbrettes liegen, auf einen Teil die Bewegungs- auf den anderen die keine_Bewegung-Funktion angewandt. Dabei wird zuerst festgelegt, welcher Anteil bleibt und erst dann auf die nicht Bleibenden die Bewegungsfunktion angewandt. Die Funktionen laufen sowohl auf den Matrizen S, I und R. Die neu berechneten Werte werden dann in den Zwischenspeichermatrizen gespeichert.

  #Alle Zeitschritte werden durchlaufen
  for t=2:T;
  #ERSTENS: MENSCHENWANDERUNG
  #Pro Zeitschritt werden alle Felder durchlaufen
  for i=1:X
   for j=1:Y
    #Pro Feld geht der Anteil q in die MM
     for k=(i-1):(i+1)
       for l=(j-1):(j+1)

         #Liegt k oder l außerhalb des Definitionsbereiches, soll nichts getan werden
         if (((k==0 || k==X+1) || l==0) || l==Y+1) 
           a=0; #do nothing
                   
         #Der Anteil (1-q) bleibt in der Zelle selbst
         elseif (j==l && i==k)
           ZS(j,i)=anteil_keine_Bewegung(j,i,q,t,ZS,S);
           ZI(j,i)=anteil_keine_Bewegung(j,i,q,t,ZI,I);
           ZR(j,i)=anteil_keine_Bewegung(j,i,q,t,ZR,R);
        
                              
         #Auf alle anderen MM wird der Anteil q aufgeteilt
         else
           ZS(l,k)=anteil_Bewegung(l,k,j,i,q,t,MM,ZS,S);
           ZI(l,k)=anteil_Bewegung(l,k,j,i,q,t,MM,ZI,I);
           ZR(l,k)=anteil_Bewegung(l,k,j,i,q,t,MM,ZR,R);
         endif
       endfor
     endfor
   endfor
  endfor

Funktionen für S, I und R

Nun haben wir neue Werte für die Anzahl der Menschen in den jeweiligen Matrizen für S, I und R erhalten. Diese sind noch in den Zwischenspeichermatrizen enthalten. Nun müssen wir auf die Werte noch die Funktionen für S, I und R anwenden um das SIR-Modell zu erhalten. Dafür werden folgende Funktionen benötigt:

Susceptible:

  function wert=S_funktion(j,i,t,gamma,S,ZS,ZI,ZR)
  wert = ZS(j,i)-(0.29*ZS(j,i)*ZI(j,i));
  endfunction

Infected:

  function wert=I_funktion(j,i,t,gamma,S,ZS,ZI,ZR)
  wert = ZI(j,i)+0.29*ZS(j,i)*ZI(j,i)-gamma*ZI(j,i);
  endfunction

Recovered:

  function wert=R_funktion(j,i,t,gamma,B,ZS,ZI,ZR)
  wert = ZR(j,i)+gamma*ZI(j,i);
  endfunction

Diese Funktionen werden nun auf die Zwischenspeichermatrizen angewandt. Die Ergebnisse davon werden dann in die Ursprungsmatrizen von S, I und R gespeichert und die Zwischenspeicher wieder auf $0$ gesetzt.

  for i=1:X
   for j=1:Y
       S(j,i,t) = S_funktion(j,i,t,gamma,S,ZS,ZI,ZR);
       I(j,i,t) = I_funktion(j,i,t,gamma,S,ZS,ZI,ZR);
       R(j,i,t) = R_funktion(j,i,t,gamma,B,ZS,ZI,ZR);
    endfor 
   endfor
 
   #Zwischenspeichermatrix gleich Null setzen
   ZS = zeros(Y,X);
   ZI = zeros(Y,X);
   ZR = zeros(Y,X);
   endfor

SIR-Modell geplottet

Nun konnten wir uns die Werte zu bestimmten Zeiten in unserem Gitter durch Plotten anzeigen lassen. Wir haben sie für jeden 5. Zeitschritt geplottet. Es folgt ein Beispiel der Implementierung anhand der kummulierten Infizierten. Die Implementierung der anderen Ausgaben ist ähnlich, nur dass $I+R$ ersetzt werden muss.

  for i=1:T/5;
  t=5*i;
  figure(t)
  Z=(I(:,:,t)+R(:,:,t));
  surf((Z), "EdgeColor", "none")
  colormap ("jet")
  colorbar
  caxis ([0 max(max(max(I+R)))])
  axis([1 12 1 10 -0.1 max(max(max(I+R)))])
  title(["kummulierte Infizierte t=",num2str(t)])
  xlabel("x")
  ylabel("y")
  test=["Fig_", num2str(t),".jpg"]
  saveas(t, test)
  endfor

Kummulierte Infizierte:

Infizierte: Anfällige:

Genesene:

SIRV-Modell geplottet

Das Gleiche haben wir mit unserem SIRV Modell versucht.

Dazu haben wir die Implementierungen für die SIR Funktionen geändert und die V-Funktion hinzugefügt:

S-Funktionen

  function wert=S_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZS(j,i)-(0.29*ZS(j,i)*ZI(j,i));
  endfunction

  function wert=S_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZS(j,i)+R(j,i,t-183)-(0.29*ZS(j,i)*ZI(j,i));
  endfunction

  function wert=S_funktion3(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZS(j,i)+R(j,i,t-183)-ZV(j,i)-(0.29*ZS(j,i)*ZI(j,i));
  endfunction

I-Funktion

  function wert=I_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZI(j,i)+0.29*ZS(j,i)*ZI(j,i)-gamma*ZI(j,i);
  endfunction

R-Funktionen

  function wert=R_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZR(j,i)+gamma*ZI(j,i);
  endfunction

  function wert=R_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZR(j,i)-R(j,i,t-183)+gamma*ZI(j,i);
  endfunction

V-Funktionen

  function wert=V_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZV(j,i);
  endfunction

  function wert=V_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZV(j,i)+141.89;
  endfunction

  function wert=V_funktion3(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV)
  wert = ZV(j,i)+0.027*(Bd(j,i)-ZV(j,i));
  endfunction

Ebenfalls mussten wir Matrizen für V und ZV erstellen, die am Anfang mit Nullen belegt wurden.

So sah dann unsere Modellierung aus:

  for i=1:X
   for j=1:Y
       S(j,i,t) = S_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
       I(j,i,t) = I_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
       R(j,i,t) = R_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
       V(j,i,t) = V_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
  
  if   t>183   && t<=339 && S(t)+R(t-183)<=B && R(t)-R(t-183)>=0
    S(j,i,t) = S_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
    I(j,i,t) = I_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
    R(j,i,t) = R_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
    V(j,i,t) = V_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
  endif

  if t >339 && t<=T && S(t)+R(t-183)<=B && R(t)-R(t-183)>=0
   S(j,i,t)=S_funktion3(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   I(j,i,t) = I_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   R(j,i,t) = R_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   V(j,i,t) = V_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
  endif 
 
   if t>443 && t<=T && S(t)+R(t-183)<=B && R(t)-R(t-183)>=0;
   S(j,i,t) = S_funktion3(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   I(j,i,t) = I_funktion(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   R(j,i,t) = R_funktion2(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   V(j,i,t) = V_funktion3(j,i,t,T,gamma,S,R,B,ZS,ZI,ZR,ZV);
   endif 
   endfor
   endfor

Kontinuierliche Diffusionsprozesse (Hundertmark)

Wie am Anfang dieser Seite bereits erwähnt, wurden uns Aufgaben zur Bearbeitung gestellt. Dabei gab es 5 Übungsaufgaben, die hier chronologisch aufgezeigt werden. Die Abkürzung FDM, die in einigen Überschriften vorkommt bezeichnet hier die Finite-Differenzen-Methode. Die Übungsaufgaben sollen Versuche sein, durch Modelle realistische Aussagen über die Pandemie-Lage treffen zu können. Für alle Übungsaufgaben wurden uns von Frau Hundertmark Ansätze und Beispiel gezeigt, sodass wir bei einigen nicht mehr viel ändern mussten. Bei Aufgabe 4 und 5 haben wir uns größtenteils an Gruppe 1 orientiert. ^[4]

1. Übungsaufgabe SW1: Kontaktmodell

Die Arbeitsaufträge der ersten Übungsaufgabe waren:

einfaches Kontaktmodell erstellen
Simulation für verschiedene Ansteckungsradien
unterschiedliche Position des Infizierten
verschiedene Bewegungsgeschwindigkeiten
verschiedene Bewegungsrichtungen

Einfaches Kontaktmodell

Was ist ein Kontaktmodell? Ein Kontaktmodell ist ein Modell, bei dem sich Punkte bewegen, die sich unter einer gewissen Abstandsschranke mit einem rotem Punkt ebenfalls rot einfärben. Die Punkte sollen selbstverständlich Menschen darstellen. Die roten Punkte sind infizierte Menschen, die andere ebenfalls unter einem gewissen Abstand infizieren.

Schachbrett und Anfangswerte

Zuerst erstellen wir eine 20x20 Matrix bzw. ein Schachbrett mit den längen 20x20 und legen unsere Anfangswerte fest:

  Nx=20;             %Anzahl von Punkten in x und y Richtung
  Ny=20;
  Zeitschritte=3;    %pro eine Zeiteinheit T=1, pro Zeiteinheit 3 Schritte
  radiusInf=0.50;    %Radius für die Verbreitung der Infektion bei Kontakt (Ansteckungsradius)
  T=3;               %Die Berechnung wird für 3 Zeiteinheiten laufen
  g=0.50;            %Geschwindigkeit mit der sich die Fußgänger bewegen
  x=[1:Nx];          % Vektoren der x und y Koordinaten
  y=[1:Ny];
  dt=1/Zeitschritte;
  [x,y]=meshgrid (x,y);      % Erstellung des Schachbrettes. Ein Punkt steht für eine Person

Position Infizierte Person

Anschließend setzen wir den ersten Infizierten in die Mitte des Schachbrettes und lassen uns dies Graphisch ausgeben:

  xInf=x(1, Nx/2);
  yInf=y(Ny/2);
  IndInf2=Nx/2;
  IndInf1=Ny/2;
  IndInf1neu=IndInf1;
  IndInf2neu=IndInf2;
  close all;
  figure (1)
  hold on;
  plot(x,y, '*b') ;
  plot(xInf,yInf, 'sr') ;
  axis([-5 25 -5 25]);
  hold off;

Schachbrett mit einem Infizierten in der Mitte

Bewegungsgenerierung

Danach werden zufällige Bewegungen generiert. Dafür wird zuerst durch eine Bewegungsgleichung die Zielposition der Punkte bestimmt.

  neuPosX=x.+g.*rand(Ny,Nx)-0.5*g;  %Neue zufällige Zielposition der sich bewegenden Punkte (Personen)
  neuPosY=y.+g.*rand(Ny,Nx)-0.5*g;

Anschließend werden die neuen Positionen für jeden Zeitschritt berechnet, mit den Infiziertenpositionen abgeglichen und ausgegeben:

  %Zeitschleife für 3 Zeitschritte pro eine Zeiteinheit T, zur Überprüfung ob neue Infizierte pro 
  Zeitschritt dazukommen
  for i=1:T*Zeitschritte
  neuX=x.+1*(neuPosX.-x)*i*dt;
  neuY=y.+1*(neuPosY.-y)*i*dt;
  figure (i+20)
  title(['t=' num2str(i*dt)])
  axis ([-5 25 -5 25]);
  hold on;
  plot(neuX,neuY,'*b');
  for k1=1:length(IndInf1)
  for j=1:Nx
  for l=1:Ny

  %die aktuellen Positionen der Fußgänger-Punkte werden mit den Positionen der Infizierten abgeglichen 
  (Abstand wird berechnet)
     abstand=norm ( [neuX(l,j)-neuX(IndInf1(k1),IndInf2(k1)), neuY(l,j)-neuY(IndInf1(k1),IndInf2(k1))]);
       if abstand<radiusInf
        %wenn der Abstand zum Infizierten unter dem Ansteckunsradius liegt, wird  der l,j-te  Punkt 
        infiziert (mit rotem Quadrat gekenzeichnet)
        abstand;    
           plot(neuX(l,j),neuY(l,j), 'sr'); 
           k2=1;

           % der Vektor der Infizierten wird untersucht, falls die l,j-Indizes der aktuellen Punkte nicht unter den infizierten Indizes auftauchen, werden  die l,j-Indizes des aktuell Infizierten zu den neuen Vektoren der Infizierten-Indizes
           % IndInf1neu, IndInf2neu zugefügt.(Grund: Speicheroptimierung)
           while k2 <=length(IndInf1)&&(l!=IndInf1neu(k2) || j!=IndInf2neu(k2)) 
            k2=k2+1;
           endwhile
           if k2 >length(IndInf1)
             IndInf1neu=[IndInf1neu l];
             IndInf2neu=[IndInf2neu j];
           endif
       endif
    endfor
  endfor
  endfor
  IndInf1=IndInf1neu;
  IndInf2=IndInf2neu;
  hold off
  lengt=length(IndInf1);
  endfor
  %Zeitschleife abgeschlossen

Wir erhalten folgende grafisches Ausgabe:

Da wir den Abstand, ab dem sich Menschen infizieren ziemlich klein gewählt haben und auch die Geschwindigkeit der Individuen relativ gering war, hat sich in 3 Zeiteinheiten nur eine weitere Person angesteckt. Im Folgenden werden wir noch sehen, was eine Geschwindigkeitsänderung ausmachen kann.

Versechsfachung der Geschwindigkeit

Versechsfachen wir die Geschwindigkeit der Individuen wie folgt, stecken sich wesentlich mehr Personen an, als vorher.

Die Zeitschleife bleibt fast unverändert. Lediglich die Berechnung der neuen Positionen verändert sich zu:

  neuX=x.+Param*(neuPosX.-x)*i*dt;
  neuY=y.+Param*(neuPosY.-y)*i*dt;

Bei einer langsameren Geschwindigkeit würde sich niemand neues anstecken, weswegen wir diesen Plot auslassen.

Änderung Infiziertenradius

Da kaum neue Infektionen entstanden sind, haben wir das gleiche nochmal mit einem größeren Ansteckungsradius von 0.9 durchgeführt. Einmal mit normaler Geschwindigkeit und einmal Versechsfacht.

  radiusInf=0.9;

Ursprungsgeschwindigkeit:

Hier sieht man eine riesen Veränderung zu vorher. Ein größerer Ansteckungsradius lässt viel mehr Infizierte entstehen.

Versechsfacht:

Auch hier sind viel mehr neue Infizierte entstanden.

In den weiteren Beispielen arbeiten wir wieder mit dem kleineren Infektionsradius.

Änderung Anfangsposition

Rand

Wir wollen nun schauen was passiert, wenn die Infizierte Person nicht in der Mitte des Gitters, sondern am Rand des Gitters ist. Dafür ändern wir im Skript die Befehle für die Infiziertenposition wie folgt ab:

  %erste Infizierte 
  xInf=x(1, Nx/2);
  yInf=y(1);
  IndInf2=Nx/2;
  IndInf1=1;
  IndInf1neu=IndInf1;
  IndInf2neu=IndInf2;

Lässt man dieses Modell genauso durchlaufen wie das erste, erhalten wir folgende graphische Ausgabe:

Auch hier infizieren sich nur wenige.

Was passiert dieses Mal bei einer Versechsfachung der Geschwindigkeit?

Auch hier sieht man, dass eine schnellere Bewegung, zu mehr Infektionen führt, allerdings sind es weniger Infektionen wie bei dem Infizierten in der Mitte.

Ecke

   %erste Infizierte 
   xInf=x(1, 1);
   yInf=y(1);
   IndInf2=1;
   IndInf1=1;
   IndInf1neu=IndInf1;
   IndInf2neu=IndInf2;

Versechsfachung:

Änderung Bewegungsrichtung

Wir ändern nun die Richtung der Bewegung in jedem Zeitschritt. Durchgeführt wird dies, mit der Abstandsschranke von 0.9, da sich bei kleinerem Ansteckungsradius niemand angesteckt hatte.

  for i=1:T*Zeitschritte
  neuX=x.+1*(neuPosX.-x)*1*dt;
  neuY=y.+1*(neuPosY.-y)*1*dt;
  x=neuX;
  y=neuY;
  neuPosX=neuX.+g.*rand(Ny,Nx)-0.5*g;
  neuPosY=neuY.+g.*rand(Ny,Nx)-0.5*g;
  figure (i+20)
  title(['t=' num2str(i*dt)])
  plot(neuX,neuY, '*b') ;
  hold on;
  axis([-5 25 -5 25]);
  quiver(neuX,neuY,neuPosX.-x,neuPosY.-y,'m')
  title(['t=' num2str(i*dt)])
  for k1=1:length(IndInf1)
  for j=1:Nx
   for l=1:Ny
     abstand=norm ( [neuX(l,j)-neuX(IndInf1(k1),IndInf2(k1)), neuY(l,j)-neuY(IndInf1(k1),IndInf2(k1))]);
       if abstand<radiusInf
            abstand;    
           plot(neuX(l,j),neuY(l,j), 'sr') ;
           %---------------------------
           k2=1;
           % der Vektor der Infizierten wird untersucht, falls die l,j-Indexen des aktuellen Punkten nicht unter den
           %infizierten Indexen auftauchen, werden  die l,j-Indexen des aktuell Infizierten zu den neuen Vektoren der Infizierten-Indexen
           % IndInf1neu, IndInf2neu zugefügt.
           %(Grund: Speicheroptimierung)
           while k2 <=length(IndInf1)&&(l!=IndInf1neu(k2) || j!=IndInf2neu(k2)) 
            k2=k2+1;
           endwhile
           if k2 >length(IndInf1)
             IndInf1neu=[IndInf1neu l];
             IndInf2neu=[IndInf2neu j];
           endif
       endif
    endfor
  endfor
  endfor
  IndInf1=IndInf1neu;
  IndInf2=IndInf2neu;
  hold off
  lengt=length(IndInf1);
  endfor

2. Übungsaufgabe SW2: Fundamentallösungen

Bei der zweiten Aufgabe sollten wir Diffusionsgleichungen durch Fundamentallösungen lösen. Sowohl für stationäre als auch für instationäre Fälle mit inhomogenen Differentialgleichungen.

Stationärer Fall

Zuerst betrachten wie den stationären Fall. Hierfür legen wir zuerst ein quadratisches Gebiet in $\mathbb {R^{2}}$ fest, über welches die Implementierung laufen soll. Die Schrittweite des Gitters legen wir auf 0,03 fest.

  step=0.03
  X=[0:step:2]; 
  Y=[0:step:2];
  [x,y]=meshgrid(X,Y);

Des weiteren soll nun eine Infektionsquelle auf dem Gebiet festgelegt werden. Diese haben wir definiert mit:

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

Um die Funktionswerte der Funktion $f$ in jedem Gitterpunkt $f(x(j,i),y(j,i))$ der Matrix speichern zu können benötigen wir folgende Implementierung:

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

Unsere Infektionsquelle auf dem Gebiet sieht somit folgendermaßen aus:

Die Lösung der inhomogenen Laplace-Gleichug, der Poissongleichung $-\Delta u(x)=f(x)$ auf ${R}^{n}$ lässt sich mithilfe der Faltung der Fundamentallösung und der Funktion der rechten Seite berechnen.

Fundamentallösung: $u(x)=(\Phi *f)(x)=\int _{\mathbb {R} ^{n}}\Phi (x-y)f(y)dy$ ,

wobei die Formel der Fundamentallösung für $\mathbb {R} ^{n},n=2,\ 3:$ für $Phi$ ist:

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

$\omega _{n}$ = Volumen der Einheitskugel in $\mathbb {R} ^{n}$ , $\|\cdot \|$ = Vektornorm in $\mathbb {R} ^{n}$
Für n=3 ist $\omega _{3}={\frac {4}{3}}\pi$ .

Dabei ist die Faltung von $Phi$ und $f$ die doppelte Integration von $Phi$ und $f$

Da wir unser Gebiet auf $\mathbb {R} ^{2}$ definiert haben, nutzen wir die obere Gleichung für $Phi$ .

Das doppelte Integral kann dann mit der Trapezregel berechnet werden.

Octave-Code:

  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=-log(sqrt((xstar.-X).^2+(ystar.-Y).^2))/(2*pi);
  phi(j,i)=0;
  I(j,i) = trapz(Y,trapz(X,phi.*f,2));
  endfor  
  endfor

Als Lösung der Poissongleichung erhalten wir nun:

Bei der Quellfunktion erhielten wir eine Funktion mit 4 Trägern. Diese könnten in der Corona-Pandemie Hotspots darstellen. Bei der Lösung der Poissongleichung kann man erkennen, wie sich die Infektion von den Hotspots auf verteilt bis das ganze Gebiet Infiziert ist. Dabei sind zu den Ecken hin. mehr Infizierte als in der Mitte des Gebietes, da in den Ecken die Quellen waren.

Instationärer Fall

Die homogene Differentialgleichung $\partial _{t}u-a\Delta u=0,x\in \mathbb {R} ^{n},t\in \mathbb {R} ^{+}$ mit dem konstanten Diffusionskoeffizienten $c(x)\equiv a\in \mathbb {R}$ wird mit folgender Fundamentallösung gelöst: $u({x},t)=(\Psi *u_{0})({x},t)=\int _{\mathbb {R} ^{n}}\Psi ({x}-{y},t)u_{0}({y})\,d{y}$

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

$\textstyle \|{x}\|^{2}=\sum _{k=1}^{n}x_{k}^{2}$ = Quadrat der Euklidischen Norm von ${x}$ .

1. Anfangslösung

Den Instationären Fall haben wir mehrmals durchgeführt. Einmal mit der Quellfunktion als Anfangslösung aus dem stationären Fall (mit verschiedenen Konstanten für $a$ ) und einmal mit einer anderen Anfangslösung.

Lassen wir oben genannte Quellfunktion als Anfangslösung u_0 über 3 Zeitschritte laufen mit folgenden Befehlen:

Lassen wir die Funktion über 3 Zeiteinheiten mit je 3 Zeitschritten laufen:

  T=3
  Zeitschritte=3
  dt=1/Zeitschritte;
  t0=0.001;
  a=10
  for k=0:T*Zeitschritte
   t=t0+k*dt
   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(-((xstar.-x).^2+(ystar.-y).^2))/(4*a*t);
   #Faltung von Anfangslosung und psi
   J(j,i) = trapz(Y,trapz(X,psi.*u0,2));
   endfor
   endfor
   Losung(:,:,k+1)=J; #Damit nach der Zeitschleife auf verschiedene Lösungen
   #zugegriffen werden kann, werden sie hier gespeichert
   endfor

Die Lösung nach den jeweiligen Zeitschritten wird durch folgende Implementierung graphisch ausgegeben:

  for k=0:T*Zeitschritte
  figure(3+k+1)
  meshc(x,y,Losung(:,:,k+1))
  grid on
  title(['Loesung der Differentialgleichung in t=',num2str(t0+k*dt)])
  xlabel('x')
  ylabel('y')
  test=["Fig_", num2str(k+1),".jpg"]
  saveas(k+1, test)
  endfor

Wir erhalten folgende grafische Darstellungen:

Im abgebildeten Octave Plot wählten wir den Diffusionskoeffizienten $a=10$ :

Ändert man diese auf $a=5$ erhalten wir:

$a=1$ :

$a=0.5$ :

Hieran kann man erkennen, dass je größer kleiner der Diffusionkoeffizient $a$ ist, desto größer werden die Werte der Lösung der Differentialgleichung.

2. Anfangslösung

Nun nutzen wir anstatt der Quellfunktion die Anfangslösung, also eine Funktion, die zum Zeitpunkt $u_{0}$ in das Gebiet eintritt. Diese sieht folgendermaßen aus:

  function wert=Anfangslosung(x,y)
  if sqrt((x.-0.2).^2+(y.-0.4).^2)<=0.05 wert=1;
  else  if sqrt((x.-1.2).^2+(y.-1.4).^2)<=0.04 wert =0.5;
  else  wert=0;
  endif
  endif
  endfunction

Um die Funktionswerte der Funktion $u_{0}$ in jedem Gitterpunkt $u_{0}(x(j,i),y(j,i))$ der Matrix speichern zu können benötigen wir folgende Implementierung:

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

Wir erhalten:

und die dazugehörige lösung der Differentialgleichung

Hier sieht man, dass die Verteilung der Lösung Zeitabhängig ist. Desto länger die Lösung im System ist, desto mehr verteilt sie sich. Da die Lösung aber relativ wenig für das große Gebiet war, wird die Funktion am Ende so flach, dass sie kaum noch zu sehen ist.

3. Übungsaufgabe SW4: Kompartimentmodelle

Aufgaben waren dieses Mal Kompartimentmodelle für Deutschland zu erstellen:

Erstellen SI-Modell
- verschiedene Infektionsraten
Erstellen SIR-Modell
- Vergleich zu den realen Daten des RKI -exponentielles Wachstum
Modifiziertes SIR-Modell
- Faktor r
- Slowdownfunktionen einbauen

SI-Modell

Zuerst legen wir die Anfangswerte fest. Es sollte angenommen werden, dass sich nur 2 Drittel der deutschen Bevölkerung infizieren kann. Das entspricht 55 Millionen Menschen.

  B = 55000;              %Kapazität in Tausend
  beta=0.29;              %Infektionsrate 
  gamma= 1/14;            %Genesungsrate
  times=(0:0.1:180);      %Zeitspanne
  I0=16/1000;             %Anzahl Infizierte, zum Zeitpunkt t0
  y0=[B-I0;I0;];          %Anfangswerte Spaltenvektor

Dann implementieren wir die Funktion der rechten Seite mit der Octave arbeiten soll. Schaut man sich die Funktionen in der eckigen Klammer bei $f$ an, erkennt man die Gleichungen des SI-Modells wobei $y(1)=S$ , $y(2)=I$ und $k=beta$ sind:

{\textstyle S'(t)=-kS(t)I(t),}

{\textstyle I'(t)=kI(t)S(t),\quad S(t_{0})=B,\ I(t_{0})=I_{0}}

  %Funktion der rechten Seite einbauen
  %Definiere die Funktion f der rechten Seite des DGL Systems y'=f(y,t), Vektorwertige Funktion f
  f=@(y,x)[-beta*y(1)*y(2)/B;beta*y(1)*y(2)/B];

Dies soll nun für jeden Zeitschritt ausgerechnet und geplottet werden:

  %numerische Lösung 
  y=lsode(f,y0,times);
  figure 1
  plot (times, y (:,1),'.',times, y (:,2), '.');
  legend("Succeptible","Infected","location","east")
  title ('SI Modell')

Wir erhalten:

verschiedene Infektionsraten

Wir haben nun verschiedene Infektionsraten gewählt und plotten lassen:

  %Betrachtung verschiedene Infektionsraten:
  beta0=0; 
  beta1=0.25;
  beta2=0.5;
  beta3=0.75;
  beta4=1;
  %Vektorwertige Funktionen für vers. Infektionsraten
  t=@(y,x)[-beta0*y(1)*y(2)/B;beta0*y(1)*y(2)/B];
  g=@(y,x)[-beta1*y(1)*y(2)/B;beta1*y(1)*y(2)/B];
  h=@(y,x)[-beta2*y(1)*y(2)/B;beta2*y(1)*y(2)/B];
  m=@(y,x)[-beta3*y(1)*y(2)/B;beta3*y(1)*y(2)/B];
  n=@(y,x)[-beta4*y(1)*y(2)/B;beta4*y(1)*y(2)/B];

Zur besseren Anschauung haben wir die ausgegebenen Bilder als Animation dargestellt und auch der beta-Wert des ersten Modells eingebaut:

Man sieht, je größer beta, desto schneller breitet sich die Infektion aus. Je kleiner beta, desto länger dauert die Infektion. Allerdings sind Werte von $beta$ ≥ $0.5$ eher unrealistisch.

SIR-Modell

Die Anfangswerte vom SIR-Modell entsprechen denen des SI-Modells. Lediglich die Anfangswerte des Spaltenvektors und die Funktion der rechten Seite ändern sich. Zusätzlich kommt noch der Anfangswert der Genesenen hinzu.

Auch hier entspricht die Funktion der rechten Seite denen des SIR-Modells wobei $beta=k$ , $y(1)=S$ , $y(2)=I$ und $gamma=w$ :

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

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

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

  R0=0;
  y0=[B-I0; I0; R0];
  %Funktion der rechten Seite 
  %Definiere die Funktion f der rechten Seite des DGL Systems y'=f(y,t), Vektorwertige Funktion f
  f=@(y,x)[-beta*y(1)*y(2)/B;beta*y(1)*y(2)/B-gamma*y(2);gamma*y(2)];
  %numerische Lösung 
  y=lsode(f,y0,times);
  plot (times, y (:,1),'.',times, y (:,2), '.', times, y (:,3),'.',times, y (:,3)+y (:,2),'.')
  legend("Succeptible","Infected","Removed","cummulative Infected","location","east")
  title ('SIR Modell, beta=0.29')

Dazugehörige Graphik:

Änderung Infektionsraten

Auch hier haben wir die Infektionsraten geändert und ein GIF erstellt. Die Infektionsraten entsprechen den gleichen wie denen des SI-Modells:

Modifiziertes SIR-Modell

Bei dem modifizierten SIR-Modell sollten wir die tatsächlichen Werte für Genesene, Infizierte, Tode und Geheilte des RKI betrachten.

Bisher haben wir für unsere SI und SIR-Modelle die Ableitungen von S, I und R genutzt. Möchte man die Gesamtanzahl der Infizierten berechnen benötigt man die I-Funktion. Welche wie folgt lautet:

 $I(t)=I_{0}e^{k(t-t_{0})}.$

Nutzen wir die I-Funktion mit unserem beta=0.29 anstelle von k, plotten diese für die ersten 30 Tage (da sie nur da einer Exponentialfunktion ähnelt und vergleichen dies mit den Daten des RKI erhalten wir folgendes:

Man sieht, dass die Daten des RKI mit den ersten 20 Tagen der Exponentialfunktion übereinstimmen, danach aber auseinander gehen.

Bestimmung Beta

Wir suchen nun ein beta, für das die Daten sich noch mehr gleichen. Dafür bestimmen wir ein Residuum: "Als Residuum bezeichnet man in der numerischen Mathematik die Abweichung vom gewünschten Ergebnis, welche entsteht, wenn in eine Gleichung Näherungslösungen eingesetzt werden."^[5]

Wir berechnen nun zuerst das Residuum von unserer Funktion mit beta=0.29 zu den Daten des RKI's:

  %Bestimmung des Residuums (Abweichung der exakten Lösung)
  for i=1:30       %Berechnung für jeden der 30 tage
  Infizierte=A(1,i);
  I(i)=I_0*exp(beta*(i-t_0));
  R(i)=(Infizierte.-I(i)).^2;
  endfor
  j=[1:30];
  R(j)
  Summe=(sum(R(j)));    %Gesamtabweichung aller Tage

Summe=97157

Anschließend haben wir uns von Cocalc ermitteln lassen, bei welchem Wert das Residuum den kleinsten Wert hat, bei welchem $beta$ also der kleinste Fehler für die ersten 30 Tage entsteht. Dafür wird die Residuenquadratsumme gebildet: "Bildet man die Summe der quadrierten Residuen für alle Beobachtungen, so erhält man die Residuenquadratsumme:

SQR:=\sum _{i=1}^{n}{\hat {\varepsilon }}_{i}^{2}=\sum _{i=1}^{n}(y_{i}-{\hat {y}}_{i})^{2}

Diese spezielle Abweichungsquadratsumme taucht in vielen statistischen Maßen, wie z. B. dem Bestimmtheitsmaß, der F-Statistik und diversen Standardfehlern, wie dem Standardfehler der Regression auf. Die Minimierung der Residuenquadratsumme führt zum Kleinste-Quadrate-Schätzer."^[6]:

  Summe=0; %Die Summe aller Residuen (also aller Abweichungen) muss 0 ergeben
  for beta=[0.2:0.002:0.4]
  for i=1:30
  Infizierte=A(1,i);
  I(i)=I_0*exp(beta*(i-t_0));
  R(i)=sum((Infizierte.-I(i)).^2);
  endfor
  j=[1:30];
  R(j);
  Summe=[Summe , (sum(R(j)))];
  endfor
  Summe/1000;

Ermittlung beta:

  MIN=min (Summe(2:length(Summe)))
  index=find(Summe==MIN)
  beta=[0.2:0.002:0.4];
  plot (beta,Summe(2:length(Summe)), '*')
  title ('Residuum')
  optim_Infrate=beta(index)

Dabei erhielten wir folgende Werte:

  MIN =  6839.1
  index =  28
  optim_Infrate =  0.25400

Diese bedeuten, dass das kleinste Residuum 6839.1 ist, wenn $beta=0.254$

Vergleichen wir nun das SIR-Modell mit $beta=0.254$ mit den Daten des RKI durch nachstehende Implementierung erhalten wir folgende Abbildung:

  B=55000;        % Kapazitätsgrenze, 2/3 der deutschen Bevölkerung
  beta=0.254;     % Basisinfektionsrate- die Rate des exponentiellen Wachstums 
  gamma=1/14;     % Wechselrate zu den Genesenen
  d=0.003;        % Todesrate
  r=1;            % Anteil der erfassten Infizierten
  I0=16/1000;     % Anzahl der Infizierten in T0
  R0=0;           % Anzahl der Genesenen in T0
  times=(0:0.1:30);
  plot (times, y (:,2), 'r.', times, y (:,3),'b.',times, y (:,3)+y (:,2),'k .')
  plot (timesWHO, inf_falleWHO./1000, 'k *', timesWHO,inf_falleWHOaktuell/1000, 'r*' , 
  timesWHO,inf_falleWHOrecovered/1000, 'b*' )
  legend ("Infected", "Removed", "cummulative Infected " ,"German data I +R" ,"German data I ","German 
  data R" ,"location", "west")
  hold off

Man kann sehen, dass die Graphen nicht übereinstimmen. Das liegt daran, dass wir durch das Residuum eine Infektionsrate bestimmt haben, die am wenigsten von den Daten des RKI abweicht, nicht aber eine die gar nicht abweicht. Des Weiteren ist zu erkennen, dass die Kurven erst ab Tag 10 auseinandergehen und die Infektionsrate des RKI anscheinend höher ist als unsere. Die Kurven für R und I+R haben wir weggelassen, da die laut den Daten es in den ersten 30 Tagen noch keine Genesenen gab, somit $I+R=I$ und $R=0$ ist.

Einflussfaktor r

Oben in der Implementierung ist dem ein oder anderen eventuell schon das $r$ aufgefallen, welches in vorherigen Implementierungen nicht zu sehen war. R ist der Einflussfaktor der gemeldeten Daten. Ist $r=1$ sind alle Infizierten auch beim RKI gemeldet. Ist $r$ < $1$ , so wurden nicht alle Infizierten erfasst.

Wir wollen nun die Auswirkungen von r auf das SIR Modell betrachten:

Gewählt haben wir 5 Werte für r:

   r=1
   r1=0.75;
   r2=0.5;
   r3=0.25;
   r4=0.1;
    % Funktion der rechten Seite angepasst:  
   f=@(y,x) [-beta*y(1)*y(2)/(r*B);beta*y(1)*y(2)/B-gamma*y(2);gamma*y(2)-d*y(2)];
   y = lsode (f, yo, times);

Wir erhalten folgende Bilder:

Slowdownfunktion

Als nächstes wollen wir eine Slowdownfunktion erstellen die sowohl zu den Daten des RKI passt als auch die Maßnahmen und Lockerungen seit Beginn der Pandemie darstellen.

  function val=slowdown(x,t)
  % Ab Tag t (Argument der Funktion: 25) wir die Infektionsrate in mehreren Schritten sinken und steigen:
  if x<t val=1 ; 
            else if x<t +8 val=t./(x.^1.1);     %Schulschließung
            else if x<t +16 val=t./(x.^1.2);    %Kontaktbeschränkung
            else if x<t +51 val=t./(x.^1.3);    %Langzeitauswirkungen
            else if x<t +56 val=t./(x.^1.2);    %Lockerungen
            else if x<t +63 val=t./(x.^1.1);    %Demo Querdenker1
            else if x<t +88 val=t./(x.^1);      %Demo Querdenker2
            else if x<t +140 val=t./(x.^0.9);   %Urlaubssaison
            else if x<t +221 val=t./(x.^0.8);   %zweite Welle
            else if x<t +235 val=t./(x.^0.9);   %Verschärfung der Maßnahmen
            else if x<t +281 val=t./(x.^1);     %Lockdown light
            else if x<t +292 val=t./(x.^1.1);   %Schulschließung
            else if x<t +295 val=t./(x.^1);     %Weihnachten
            else if x<t +318 val=t./(x.^1.01);  %Impfungen (nur sehr wenig Auswirkungen)
            else if x<t +328 val=t./(x.^1.11);  %FFP2-Maskenpflicht
            else if x<t +366 val=t./(x.^1.12);  %zweiter Impfstoff
            else if x<t +373 val=t./(x.^0.92);  %Schulöffnungen + mehr Testungen
            else if x<t +377 val=t./(x.^0.91);  %Aussetzung Astra Zeneca
            else if x<t +390 val=t./(x.^0.92);  %Astra wieder aktiv
            else if x<t +393 val=t./(x.^0.93);  %Impfung Hausarzt
            else val=t./(x.^0.83) ;             %Ostern / schönes Wetter Menschen gehen mehr raus
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
  endif
  if x> 418 val=0.2; endif
  endfunction

Wir haben die Funktion mit den Daten des RKI abgeglichen und auch hier eine geeignete Infektionsrate durch die Residuenquadratsumme bestimmt, denn die zuvor bestimmte, passte nicht mehr so sehr zu der Slowdownfunktion. Für $beta$ ergab sich nun der Wert $0.34$ .

Geplottet auf die ersten 100 Tage und verglichen mit den Daten des RKI kann man erkennen, dass die tatsächlichen Werte sich den jeweiligen Funktionen sehr ähneln.

Implementiert man die Slowdownfunktion komplett, erhält man folgende Graphik:

Man sieht, dass diese das Infektionsgeschehen zwar verlängert, aber dennoch nicht die Realtität abbildet, denn sonst hätten wir die Pandemie schon Ende letzten Jahres überstanden.

4. Übungsaufgabe SW7: Randwertaufgabe für Poissongleichung mit FDM

In dieser Aufgabe sollten wir die stationäre Diffusion mithilfe der Poisson-Gleichung auf einem rechteckigem Gebiet mit der Finiten-Differenzen Methode implementieren. Die Methode ist ein Gitter-basiertes Verfahren, welches die unbekannte Funktion $u(x,t)$ auf endlich vielen Gitterpunkten $P_{i,j}=(x_{i},y_{j}),\ i=1,\ldots N_{x},\ j=1,\ldots N_{y}$ approximiert.

$-a\Delta u=f~in~D=(0,X_{end})x(0,Y_{end}),$

$u(x)=g(x)~auf~\partial D$

Zunächst müssen einige Eigenschaften wie die Quelldichtefunktion, das Gebiet etc. festgelegt werden:

Die "rechte Seite" der Poisson-Gleichung

Die "rechte Seite" der Poissongleichung besteht aus einer Quelldichtefunktion $f$ welche die Flussgeschwindigkeit der Quelle in das System hinein beschreibt. Bezogen auf die Corona-Pandemie, also die Geschwindigkeit in der Neuinfizierungen kontinuierlich in das System/in einem Gebiet hinzukommen.

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

Gebiet der Modellierung

Um ein Gebiet D festzulegen, müssen die Kantenlängen des Rechtecks angegeben werden. Des Weiteren benötigt man die Anzahl an Gitterunkten in X-Richtung:

  xend=9
  yend=6
  N=30  %Anzahl von X-Punkten (ohne Randpunkte)

Um das äquidistante Gitter erstellen zu können, muss der Abstand zwischen den einzelnen Gitterpunkten festgelegt werden. Dafür muss die $x$ -Kantenlänge durch $N+1$ geteilt werden, da der $0$ -te und $N$ -te Index den Rändern entspricht und nicht zu den Gitterpunkten zählen:

  hx=xend/(N+1)
  hy=hx;
  M=floor((yend/hy))-1
  [x,y]=meshgrid(hx:hx:xend-hx,hy:hy:yend-hy);

Zusätzlich benötigen wir noch einen konstanten Diffusionskoeffizienten (Verbreitungsgeschwindigkeit). Je größer dieser ist, desto schneller erfolgt die Verbreitung.

  a=1.5

Randbedingungen

Dirichlet-Randwertproblem

Bei dem Dirichlet Randwertproblem ist die gesuchte Funktion $u$ am Rand $\partial D$ durch eine Funktion $g$ gegeben.

-\Delta u(x)=f(x)~~~~{\mbox{in}}~~D\subset \mathbb {R} ^{2},

u(x)=g(x)~~~~~~~~~~{\mbox{auf}}~~\partial D

Wir wählen für unsere Funktion $g$ für die verschiedenen Ränder einen konstanten Wert:

  %Dirichlet-RB:
  wert_unten=-0.035
  wert_oben=0.035
  wert_links=0.035
  wert_rechts=-0.035

Systemmatrix

Um die gesuchte Funktion $u$ zu berechnen, muss $\Delta u(x)$ approximiert werden. Eindimensional kann dabei $u''$ betrachtet werden. Die zweite Ableitung kann dabei mit Hilfe des zentralen Differenzenquotienten diskretisiert werden:

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

Dieser ergibt sich aus der Differenz des Rückwärts- und des Vorwärtsdifferentenquotienten $D_{h}^{+}u(x_{i})-D_{h}^{-}u(x_{i})$ .

Bei einer zweidimensionalen Betrachtung wie bei unserem Fall (rechteckiges Gebiet $D$ ), müssen zwei Raumrichtungen x und y berücksichtigt werden und damit zwei Differenzenquotienten $D_{h,x}^{2}u(x_{i},y_{j}),\ D_{h,y}^{2}u(x_{i},y_{j})$ aufgestellt werden.

Diese werden durch die partiellen Ableitungen im Punkt $(x_{i},y_{j})$ approximiert:

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

Summiert man diese erhält man: $\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.$

Um die Systemmatrix aufzustellen zu können, müssen die unbekannten Funktionswerte $u(x_{i},y_{j})$ aller Gitterpunkte, sowie die Funktionswerte der rechten Seite $f(x_{i},y_{j})$ in je einen langen Vektor eingeordnet werden. Diese sollen wie folgt aufgestellt werden:

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

{\vec {f}}_{h}=(f_{1,1},f_{2,1},\ldots f_{n,1}|\ldots |f_{1,m},\ldots f_{n,m})^{T}

Daraus lässt sich ein lineares Gleichungssystem $A_{h}{\vec {u}}_{h}={\vec {f}}_{h}$ aufstellen, bei welchem $A_{h}\in {\mathbb {R} }^{n.m\times n.m}$ die blockdiagonale Matrix ist.

Es folgt:

A_{h}={\frac {1}{h^{2}}}{\begin{pmatrix}B&I&&\ldots &0\\I&B&I&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\&&I&B&I\\0&\ldots &0&I&B\end{pmatrix}}

mit Diagonalblock

B={\begin{pmatrix}-4&1&0&\ldots &0\\1&-4&1&\ldots &0\\\vdots &\ddots &\ddots &\ddots \\&&1&-4&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.

Praxis

Um die oben genannte Theorie zu implementieren wird zuerst die Systemmatrix folgendermaßen aufgestellt:

  Vec1=ones(N*M,1)  %erzeug Vektor der Länge N*M
  BB=diag(-4.*Vec1,0)+diag(ones(N*M-1,1),1)+diag(ones(N*M-1,1)-1)+diag(ones(N*(M-1),1),N)+diag(ones(N*(M-1),1),-N); 
  %Systemmatrix N^2xM^2 mit -4 auf der Diagonale, 1 auf den Nebendiagonalen und 1 auf den N und -N-ter Nebendiagonalen

Es muss zusätzlich beachtet werden, dass beim Übergang in die nächste Blockmatrix B die Nebendiagonalen aus Einsen unterbrochen werden und dort je ein Nulleintrag vorliegt:

${\begin{pmatrix}-4&1&0~&1&0&0\\1&-4&1~&0&1&0\\0&1&-4~&0&0&1\\\\1&0&0~&-4&1&0\\0&1&0~&1&-4&1\\0&0&1~&0&1&-4\end{pmatrix}}$

Diese Eigenschaft wird durch folgenden Befehl implementiert:

  %Korrektur der Matrix (Blockdiagonalität)
  for i=1:(M-1)
  BB(i*N+1,i*N)=0;BB(i*N,i*N+1)=0;
  endfor

Anschließend wird die Matrix mit dem Diffusionskoeffizient/h² multipliziert:

  BB=BB*a/hx^2

Matrixdarstellung: (s.Bild: Sytemmatrix)

  figure(1);
  spy(BB);
  xlabel("Spalten")
  ylabel("Zeilen")
  zlabel("Matrix B")
  title ("Systemmatrix B");

Rechte Seite der Poissongleichung

Um die rechte Seite des Systems aufzustellen, wird die oben definierte Quellfunktion für jeden Gitterpunkt gespeichert. Danach werden die Dirichlet-Randbedingungen hinzugefügt.

  clear f;
  for i=1:N
  for j=1:M
  f(j,i)=1*Quellfunktion(x(j,i),y(j,i));  %*hx^2;
  %Dirichlet Randbedingung unten und oben:
  if j==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_unten*a/hx^2; endif #(wegen Poisson muss NRB 
  1/h^2)
  if j==M f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_oben*a/hx^2; endif
  %Dirichlet Randbedingung links und rechts:
  if i==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+ wert_links*a/hx^2; endif
  if i==N f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+ wert_rechts*a/hx^2; endif
  endfor
  endfor

Ecken

Die Ecken müssen anders betrachtet werden, da hier zwei Randbedingungen aufeinandertreffen:

  f(1,1)=1*Quellfunktion(x(j,i),y(j,i))+(wert_links+wert_unten)*a/hx^2; 
  f(1,N)=1*Quellfunktion(x(j,i),y(j,i))+(wert_rechts+wert_unten)*a/hx^2;
  f(M,1)=1*Quellfunktion(x(j,i),y(j,i))+(wert_links+wert_oben)*a/hx^2;
  f(M,N)=1*Quellfunktion(x(j,i),y(j,i))+(wert_rechts+wert_oben)*a/hx^2;

Nun kann die rechte Seite in den Vektor $b$ umgeformt werden. Des Weiteren fließt das $-1$ aus der Poissungleichung $-\Delta u(x)=f(x)$ ein.

  size (f) #rechte Seite kann in Vektor b umgeformt werden. Negativ wegen PG
  b=-1*reshape(f',N*M,1); %Quellfunktion in Spaltenvektor um Systemmatrix drauf zu werfen. (wird 
  transponiert, damit die Zeilen untereinander aufgestellt werden)
  %reshape funktioniert Spaltenmässig: es wird Sp1, dann Sp2, usw eingeordnet, wobei wir die Matrix f 
  nach Reihen/Zeilen einordnen.
  %Daher soll die Matrix f' mit N(Zeilen)x M(Spalten) in einen langen Vektor eingeordnet werden

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

Graphische Ausgaben Im Folgenden sind die Bilder der Lösung und der Quellfunktion zu sehen.

  %Öffne ein Bildfenster und plotte die Lösung auf dem realen Koordinatenfeld
  figure(2);
  surfc(x,y,sol_matrix);
  title ("Loesung");
  ylabel("y")
  xlabel("x")
  
  %weiteres Bild mit Quellfunktion- rechte Seite des Systems
  figure(3);
  surfc(x,y,f);
  title ("Quellfunktion");
  ylabel("y")
  xlabel("x")

Implementierung der Neumann-Randbedingungen

Nun werden die Dirichlet-Randbedingungen mit den Neumann-Randbedingungen kombiniert. Für die Neumannn-Randbedingungen müssen die Richtungsableitungen der Funktion $u$ in Richtung des äußeren Normalenvektors gegeben sein:

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

Da wir auf einem rechteckigen Gebiet arbeiten gilt:

Es wird nun die Finite-Differenzen-Methode (FDM) für die Neumann-Randbedingung auf dem Rechteck $D$ mit homogenen Randbedingungen mit $\alpha =\beta =\gamma =\delta =0$ implementiert. Unser betrachtetes Gitter muss somit um die Randwerte erweitert werden:

  [x,y]=meshgrid(0:hx:xend,0:hy:yend);
  N=N+2;
  M=M+2;
  %konstanter Diffusionskoeffizient
  a=1.0;

Systemmatrix

Dadurch dass das Gitter erweitert wurde, gibt es nun weitere unbekannte Funktionswerte $u(x_{i},y_{j})$ . Der lange Vektor wird dann wie folgt aussehen:

{\vec {u}}_{h}=(u_{0,0},u_{1,0},\ldots u_{n,0},u_{n+1,0}|u_{0,1},u_{1,1},\ldots u_{n,1},u_{n+1,1}|\ldots |u_{0,m},u_{1,m},\ldots u_{n+1,m}|u_{0,m+1},u_{1,m+1},\ldots u_{n+1,m+1})^{T}

Dabei entspricht die erste Spalte dem unterem Gitterrand, die letzte dem oberen. Außerdem wurden die Spaltenabschnitte um je zwei Punkte erweitert, wobei der erste auf dem linken Rand, der letzte auf dem rechten Rand liegt.

Die Dirichlet-Randbedingungen sollen nur an bestimmten Rändern durch Neumann-Randbedingungen ergänzt werden, jedoch nicht vollständig ersetzt, da sonst kein brauchbares Ergebnis herauskommt. Dies geschieht hier auf dem linken und rechten Rand. Das Gitter wird jedoch nicht darauf angepasst, wodurch ein kleiner Fehler entsteht.

Um das lineare Gleichungssystem $A_{h}{\vec {u}}_{h}={\vec {f}}_{h}$ aufzustellen muss unsere blockdiagonale Matrix $A_{h}$ um Einträge erweitert werden, damit $A_{h}\in {\mathbb {R} }^{(n+2)(m+2)\times (n+2)(m+2)}$ . Dabei wird der Diagonalblock $B$ jeweils um eine Zeile oben und unten erweitert. Diese enstprechen den Werten des linken und rechten Randes, weshalb hier ein anderer Eintrag erfolgen muss. Dieser leitet sich aus dem Differenzenquotienten her.

  ${\frac {\partial u}{\partial {\vec {n}}}}=\pm {\frac {\partial u}{\partial x}}$  oder  $\pm {\frac {\partial u}{\partial y}}$

- 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})=-\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})=\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})=-\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})=\delta ,\ i=1,\ldots n,$

siehe ^[7]

Die Systemmatrix sieht dann wie folgt aus:

A_{h}={\frac {-1}{h^{2}}}{\begin{pmatrix}B&I^{*}&&\ldots &0\\I^{*}&B&I^{*}&&\vdots \\\vdots &\ddots &\ddots &\ddots &\\&&I^{*}&B&I^{*}\\0&\ldots &0&I^{*}&B\end{pmatrix}}

mit Diagonalblock

B={\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}},\quad B\in {\mathbb {R} }^{n+2\times n+2}

und null-erweiterter Einheitsmatrix

I^{*}={\begin{pmatrix}0&&&\ldots &0\\&1&&\ldots &\\\vdots &\ddots &\ddots &\ddots \\&&&1&\\0&\ldots &&&0\end{pmatrix}},\quad I^{*}\in {\mathbb {R} }^{n+2\times n+2}

.

Zugehörige Octave-Implementierung

  Vec1=ones(N*M,1); % erzeugt Vektor der Länge N*M
  BB=diag(-4.*Vec1,0)+diag(ones(N*M-1,1),1)+ diag(ones(N*M-1,1),-1)+diag(ones(N*(M
  -1),1),N)+diag(ones(N*(M-1),1),-N);
  %Systemmatrix N^2xM^2 mit -4 auf der Diagonale und 1 auf den Nebendiagonalen und 1 auf den N und 
  -N-ter Nebendiagonale
  %Korrektur der Matrix (Blockdiagonalität)
  for i=1:(M-1)
  BB(i*N+1,i*N)=0;BB(i*N,i*N+1)=0;
  endfor
  % Matrix mit Diffkoeff/h^2 multiplizieren
  BB=BB*a/hx^2;

Um die obere und untere Zeile der Blockmatrizen zu ändern folgt:

  % Neumann RB für y: c* partial_y u=0 - einzelne Blöcke der Matrix ersetzen
  h=hx;
  block=diag((a/h)*ones(N,1));
  % Neumann RB für x: c* partial_x u=0 - einzelne Zeilen der Matrix ersetzen
  for i=1:(M-2)               #NRB linker und rechter Rand - einzelne Zeilen der
                              #Matrix ersetzen. oberste und unterste Zeile der Blockmatrix
  vector_up=zeros(1,N*M);     #Müssen geändert werden. Es werden die RB für rechts und
  vector_up(N*i+1)=a/h;       #links in die Systemmatrix eingepflegt. 
                              #a/h, denn a/h² *h (Matrix wurde schon mit a/h² multipliziert
  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 ;    #Ersetzt entsprechende Zeile mit Vektor
  BB(i*N+N,:)=vector_down ;
  endfor

Matrixdarstellung

Stellt man das ganze als Matrix da, sieht es wie folgt aus: s. Bild Systemmatrix2

  figure (4);
  %surface (BB);
  spy (BB);
  xlabel("Spalten")
  ylabel("Zeilen")
  zlabel("Matrix B")
  title ("Systematrix B");

rechte Seite des Systems

  for i=1:N
  for j=1:M
  f(j,i)=1*Quellfunktion(x(j,i),y(j,i));%*hx^2;
  %Neumann Randbedingung links und rechts:
  if i==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i)); 
  endif
  if i==N f(j,i)=1*Quellfunktion(x(j,i),y(j,i)); 
  endif
  %Dirichlet Randbedingung oben und unten:
  if j==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_unten*a/hx^2; 
  endif
  if j==M f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_oben*a/hx^2; 
  endif
  endfor
  endfor

Ecken

Die Eckfunktionen werden folgendermaßen implementiert:

  f(1,1)=1*Quellfunktion(x(j,i),y(j,i));
  f(1,N)=1*Quellfunktion(x(j,i),y(j,i));
  f(M,1)=1*Quellfunktion(x(j,i),y(j,i));
  f(M,N)=1*Quellfunktion(x(j,i),y(j,i));
  size (f)
  b=-1*reshape(f',N*M,1);
  %reshape funktioniert Spaltenmässig: es wird Sp1, dann Sp2, usw eingeordnet, wobei wir die Matrix f 
  nach 
  Reihen/Zeilen einordnen.
  %Daher soll die Matrix f' mit N(Zeilen)x M(Spalten) in ein langen Vektor eingeordnet werden.
  % LOESUNGSSCHRITT:
  sol=BB\b;
  sol_matrix=reshape(sol,N,M);% Matrix mit N(Zeilen)x M(Spalten)
  sol_matrix=sol_matrix'; % Transponiert, wegen reshape

Graphik Graphisch sieht das ganze folgendermaßen aus: Zuerst ist wieder die Lösung und danach die Quellfunktion zu sehen.

  figure(5); 
  surfc(x,y,sol_matrix);
  title ("Loesung");
  ylabel("y")
  xlabel("x")

  %weiteres Bild mit Quellfunktion- rechte Seite des Systems
  figure(6);
  surfc(x,y,f);
  title ("Quellfunktion");
  ylabel("y")
  xlabel("x")

5. Übungsaufgabe SW10: Epidemiologische Modellierung als Reaktionsdiffusionsprozess mit FDM

Arbeitsaufträge:

Implementierung FDM-Methode der Reaktionsdiffusionsgleichung für die zeitlich-räumliche Modellierung der Epidemie auf einem Rechteckgebiet $D$ mit homogenen Neumannrandbedingungen ${\frac {\partial u}{\partial {\vec {n}}}}|_{\partial D}=0$ .
Ergänzung um geographisch differenzierte Bevölkerungsdichte. Funktion von Raumvariablen.
zusätzliche nicht konstanter Diffusionskoeffizient.

Modellierung Corona-Ausbreitung der Stadt Mainz

Hierfür benötigen wir wie in der vorigen Aufgabe die Finite-Differenzen-Methode zur Lösung der Reaktionsdiffusionsgleichung $\partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=f(u(x,t)).$ .

Wir betrachten in der Modellierung die Bevölkerungsdichte der Stadt Mainz Zentrum (Altstadt, Neustadt, Oberstadt und Hartenberg-Münchfeld).

Fläche: 98km² (wir haben diese auf 100km² aufgerundet um eine quadratische Fläche von 10kmx10km erstellen zu können). Bevölkerungsdichte: 6111,25 Einwohner/km²

Datei:Bild von Mainz.png

Bild von Mainz

Auf der Abbildung der Stadt Mainz kann man sehen, dass man die Stadt, wie oben schon erwähnt, grob als Quadrat ansehen.

  xend=10.0; %in km
  yend=10.0; %in km
  N=80; % Anzahl von Gitterpunkten in x-Richtung(0-te und N-te Index 
  entsprechen den Rändern, keine Gitterpunkte am Rand)
  a=0.01; %konstanter Diffusionskoeffizient
  T=450; %Zeitintervall in Tagen
  delta_t=0.03; %Zeitschritt

Um das Gitter zu konstruieren und um die $N$ Gitterpunkte in gleichgroßem Abstand auf der Kante zu verteilen, bzw. den Abstand der Gitterpunkte zu berechnen muss die Kantenlänge des Quadrates durch $N+1$ geteilt werden. Da ein Quadrat vorliegt, ist die Anzahl der Gitterpunkte in x- und y-Richtung gleich groß.

  hx=xend/(N+1) %Abstand Gitterpunkte/um Gitterpunkte im gleichen Abstand auf der Karte zu verteilen
  hy=hx; % setze so, damit sich Quadrate bilden 
  %M=floor((yend/hy))-1 -> Anzahl der Gitterpunkte in y-Richtung
  M=N; %Im Falle eines Quadrates, gleiche Anzahl an Gitterpunkten 
  Nr_timesteps=floor (T/delta_t) %Anzahl der Zeitschritte

Anfangswerte

  B=6111.25; %Bevölkerungsdichte Einwohner pro Quadratkilometer in Mainz Zentrum (Altstadt, Neustadt, 
  Oberstadt, Hartenberg-Münchfeld),Durchschnittlich: 2236
  c=0.34; %Infektionsrate aus SIR-Abgabe3
  k=c./B; %Infektionsrate k.u.(B-u)
  w=1/14;%Wechselrate

Anfangsfunktion

    function wert = initialfunction (x, y) #Quellfunktion
    wert = 0;
    r = 0.31;  # Fläche des Kreises: 0.3km^2 entspricht ca. 1833 Einwohnern 
    if sqrt ((x - 6.5) .^ 2 + (y - 6.5) .^ 2) <= r
   wert = 1 / (pi * r ^ 2); # Infektionen im Zentrum von Mainz (ein Infizierter)
  endif
  endfunction

Systemmatrix

Das erstellte Gitter muss (für die NRB) um die Randpunkte erweitert werden:

  [x,y]=meshgrid(0:hx:xend,0:hy:yend);
  N=N+2;
  M=M+2;

Die Systemmatrix wird nach dem gleichen Schema wie in Übungsaufgabe 4 erstellt:

  Vec1=ones(N*M,1); % erzeugt Vektor der Länge N*M
  BB=diag(-4.*Vec1,0)+diag(ones(N*M-1,1),1)+ diag(ones(N*M
  -1,1),-1)+diag(ones(N*(M-1),1),N)+diag(ones(N*(M-1),1),-N);
  % Systemmatrix N^2xN^2 mit -4 auf der Diagonale und 1 auf den 
  Nebendiagonalen und 1 auf den N und -N-ter Nebendiagonale
  %Korrektur der Matrix (Blockdiagonalität)
  for i=1:(M-1)
  BB(i*N+1,i*N)=0;BB(i*N,i*N+1)=0;
  endfor
  % Matrix mit Diffkoeff/h^2 multiplizieren
  BB=BB*a/hx^2;

homogene Neumann-Randbedingungen integrieren

Bei den homogenen Neumann-Randbedingungen ist die Ableitung über den Rand ${\frac {\partial u}{\partial {\vec {n}}}}|_{\partial D}=0$ . Damit wird ein freier Fluss der Infektion in die angrenzenden Gebiete beschrieben. Für diese Randbedingungen werden jeweils die erste und letzte Blockzeile der Systemmatrix, welche dem oberen und unteren Rand entsprechen, angepasst. Auf dem linken und rechten Rand müssen die oberste und unterste Zeile der Diagonalblöcke der Systemmatrix ersetzt werden:

  % Neumann RB für y: a* partial_y u=0 - einzelne Blöcke der Matrix ersetzen
  block=diag((a/hx)*ones(N,1));
  %unterer Rand
  BB(1:N,1:N)=-block;
  BB(1:N,N+1:2*N)=block;
  %oberer Rand
  BB(N*(M-1)+1:N*M,N*(M-1)+1:N*M)=-block;
  BB(N*(M-1)+1:N*M,N*(M-2)+1:N*(M-1))=block;

  % Neumann RB f¸r x: a* 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/hx;
  vector_up(N*i+2)=-a/hx;
  vector_down=zeros(1,N*M);
  vector_down(N*i+N)=a/hx;
  vector_down(N*i+N-1)=-a/hx;
  BB(i*N+1,:) =-vector_up ; #Vektor der ersten Zeile
  BB(i*N+N,:) =-vector_down ; #Vektor der letzten Zeile
  endfor

Hieraus ergibt sich die folgende Systemmatrix (s.rechts):

Anfangslösung in Matrix

Die oben definierte Anfangslösung muss nun als Matrix aufgestellt werden und Infizierte müssen gesetzt werden.

  I0=1 % Anzahl der Anfangsinfizierungen pro gegebene Fläche pi*r^2
  for i=1:N
  for j=1:M
  initial(j,i)=I0*initialfunction(x(j,i),y(j,i));
  endfor
  endfor
  size (initial)
  sol_old=1*reshape(initial',N*M,1);
  %reshape funktioniert Spaltenmäßig: Sp1,Sp2, wir wollen Zeilenmäßig
  %Matrix f' ist N(Zeilen)x M(Spalten)

Grafik zur Einwohnerzahl implementieren

In diesem Teil wird eine Grafik zur Einwohnerzahl

als Bevölkerungsmatrix eingelesen. Wir haben uns für einen geeigneten Bildausschnitt für die Stadt Mainz entschieden:

  %Graphik einlesen und
  %Geeigneter Ausschnitt wählen für Mainz, Matrix A ist eine 3D Matrix
  A=imread ('Monitor einwohnerzahl zoom1.png');
  size(A); % Zeigt die Dimension der Bildmatrix
  S=A(83:92,105:114);
  imshow (S); % Bild wird angezeigt
  S=double(S);
  [s1 s2]=size(S);
  Max=max(max(S))
  S=(Max-S); % Farbumkehrung/ Anpassung

Anschließend werden die Werte der oben aufgestellten Matrix (10x10) an das am Anfang erstellte Gitter angepasst, denn die Matrizen sind unterschiedlich eingeteilt. Die Schritte des Gitters sind kleiner, als die der Matrix, weswegen die Matrix interpoliert werden muss:

  Sint=interp2(S,xend*(1*x.+2)/(xend+2),yend*(1*y.+2)/(yend+2), 
  'cubic'); 
  %y-Koordinate im Image läuft von oben nach unten, Bild flippen oben 
  <-> unten
  Sint=flipud (Sint);
  size (Sint)
  figure(67);
  surface(x,y, (B/Max)*Sint)
  title ("Bevoelkerungsdichte");
  ylabel("y")
  xlabel("x")
  colorbar
  %Skalierung auf max. Wert B-> Anpassung der Bevölkerungsdichte auf 
  die Karte
  %Korrektur der Bevölkerungsdichte (keine Null-Einträge) +0.1*Max 
  -> minimalen Wert der Bevölkerungsdichte setzen wir auf 10 % des 
  Maximum
  %Assemblierung der Matrix in ein Vektor
  B=reshape((B/Max)*Sint',N*M,1).+0.1*B; %Angepasste Bevölkerungsdichte
  size (B)

Dabei erhalten wir nun das Bild Bevölkerungsdichte Mainz. Vergleicht man den Bildausschnitt von Mainz mit dem Bild der gewählten Bevölkerungsdichte kann man schließen, dass die beiden Bilder annähernd übereinstimmen.

Reaktionsterm für kumulierte Infizierte

In der nächsten Teilaufgabe sollten wir die Bevölkerungsdichte im Reaktionsterm geographisch differenziert als Funktion von Raumvariablen betrachten. Für die Berechnung der Dichte der kumulierten Infizierten werden folgende Funktionen benötigt: $f(u(x,t))={\frac {c}{B(x)}}u(x,t){\big (}B(x)-u(x,t){\big )}$ (Rekationsterm)

  F=@(u) (c./B).*u.*(B.-u);

Zeitschleife

   for i=1:Nr_timesteps
  #Term F wird nun zu unserem System vom GDGL des räumlich diskretisierten NR-Problem für die Diffusion hinzugefügt.
  sol= sol_old+ BB*sol_old*delta_t + 1*F(sol_old)*delta_t;
  sol_old=sol;
  matrix_solution(:,i)=sol;
  endfor
  % Auswahl an Bildern drucken, jedes fünfzehnte Bild
  fig_index=floor([0.15:0.15:1]*Nr_timesteps)
  j=0
  for i=fig_index
  sol_matrix=reshape(matrix_solution(:,i),N,M);% Matrix mit N(Zeilen)x M(Spalten)
  sol_matrix=sol_matrix';
  disp(['Figure ',num2str(round(i/fig_index(1)))]);
  j=j+1;
  figure(j);
  % Alternativ surfc(x,y,sol_matrix*hx^2, "EdgeColor", "none") -> 3D
  surface(x,y,sol_matrix*hx^2, "EdgeColor", "none") % multipliziert mit hx^2, damit Anzahl der Infizierten pro 
  Rasterdargestellt wird (und nicht die Dichte)
  colormap ("jet")
  axis([0 xend 0 yend 0 1.1*max(max(B))*hx^2])
  colorbar
  title (["Loesung in t=", num2str(round(delta_t*i))]);
  ylabel("y")
  xlabel("x")
  test=["Fig_", num2str(j),".jpg"]
  endfor

Anfangslösung

  sol_matrix=reshape(matrix_solution(:,1),N,M);% Matrix mit N(Zeilen)x M(Spalten)
  sol_matrix=sol_matrix';
  disp(['Figure ',num2str(0)]);
  j=j+1;
  figure(j);
  %%%Alternativ surfc(x,y,sol_matrix*hx^2, "EdgeColor", "none")
  surface(x,y,sol_matrix*hx^2, "EdgeColor", "white")
  colormap ("jet")
  axis([0 xend 0 yend 0 1.1*max(max(B))*hx^2])
  colorbar
  title (['Loesung in t=', num2str(round(delta_t))]);
  ylabel("y")
  xlabel("x")

Die Abbildung zeigt, dass wir zu Beginn der Corona-Pandemie ca. einen Infizierten in Mainz haben:

Fazit

Die Legende gibt die Anzahl der kumulierten Infizierten pro Gitterpunkt an. Auf einem Gitterpunkt(Pixel) sind durchschnittlich 61 Einwohner. Wenn davon am Ende 45 Einwohner infiziert sind, sind das 2/3 aller Einwohner. Ab Tag 134 passiert nichts mehr. Hier ist schon der Höhepunkt erreicht und es infiziert sich niemand neues mehr.

Reaktionsterm für aktuell Infizierte

f_{I}(u(x,t))={\frac {c}{B(x)}}u_{I}(x,t)u_{s}(x,t)-wu_{I}(x,t)

Dafür wird die alte Lösung genutzt

  sol_old=1*reshape(initial',N*M,1); % Aktuell Infizierte
  sol_old_susp=B.-sol_old; % Susceptibles-Infektions-Anfällige

Unsere Slowdownfunktion aus Übungsaufgabe 3 wird ebenfalls benötigt:

  function val=slowdown4(x,t)
  % Ab Tag t (Argument der Funktion) wird die Infektionsrate in mehreren Schritten sinken:
  if x<t val=1 ; 
            else if x<t +8 val=t./(x.^1.1); 
            else if x<t +16 val=t./(x.^1.2); 
            else if x<t +51 val=t./(x.^1.3); 
            else if x<t +56 val=t./(x.^1.2);
            else if x<t +63 val=t./(x.^1.1);
            else if x<t +88 val=t./(x.^1);
            else if x<t +140 val=t./(x.^0.9);
            else if x<t +221 val=t./(x.^0.8);
            else if x<t +235 val=t./(x.^0.9);
            else if x<t +281 val=t./(x.^1);
            else if x<t +292 val=t./(x.^1.1);
            else if x<t +295 val=t./(x.^1);
            else if x<t +318 val=t./(x.^1.01);
            else if x<t +328 val=t./(x.^1.11);
            else if x<t +366 val=t./(x.^1.12);
            else if x<t +373 val=t./(x.^0.92); %Schulöffnungen + mehr Testungen
            else if x<t +377 val=t./(x.^0.91);
            else if x<t +390 val=t./(x.^0.92);
            else if x<t +393 val=t./(x.^0.93);
            else val=t./(x.^0.83) ;
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
                  endif
  endif
  if x> 418 val=0.2; endif

Reaktionsterm

  F=@(u,v,t) (c*slowdown4(25,t)./B).*u.*v;

Zeitschleife

  for i=1:Nr_timesteps #Die Zeitschleife erfolgt erneut mit Explizitem Eulerverfahren
  %zuerts nichtlin. F-tion des Reaktionsterms aus alten Werten ausrechnen, wir weiter benutzt:
  Reakt=F(sol_old,sol_old_susp, i*delta_t);
  %Erste Gleichung - für Susceptibles
  sol_susp= sol_old_susp+ BB*sol_old_susp*delta_t - Reakt*delta_t;
  %Zweite Gleichung - für Infected (aktuell), ( um Term -w*I erweitert)
  sol= sol_old+ BB*sol_old*delta_t+ (Reakt-w.*sol_old)*delta_t;
  %Speicherung für den nächsten Zeitschritt
  sol_old=sol;
  sol_old_susp=sol_susp;
  % Speicherung der aktuell Infizierten in einer Matrix, Speicherung der Anfälligen (S) nicht
  implementiert
  matrix_solution1(:,i)=sol;
  endfor

Animation

Anschließend erhalten wir mit dem gleichen Octave-Befehl wie im ersten Lösungsschritt unsere Animation für die aktuell Infizierten.

Zum Vergleich erfolgt hier noch die Graphik für die aktuell Infizierten ohne Slowdown-Funktion.

[1] ttps://de.wikiversity.org/wiki/Kurs:R%C3%A4umliche_Modellbildung

[2] ttps://de.wikiversity.org/wiki/Kurs:R%C3%A4umliche_Modellbildung/Aufgaben_Tutorium_SS2020

[3] ttps://de.wikipedia.org/w/index.php?title=Moore-Nachbarschaft&oldid=207626306

[4] ttps://de.wikiversity.org/wiki/Kurs:R%C3%A4umliche_Modellbildung/Gruppe_1

[5] ttps://de.wikipedia.org/wiki/Residuum_(numerische_Mathematik)

[6] ttps://de.wikipedia.org/wiki/St%C3%B6rgr%C3%B6%C3%9Fe_und_Residuum

[7] ttps://de.wikiversity.org/wiki/Numerische_Modellierung_der_Diffusion#Inhomogene_Randbedingungen

[1]

[2]

[3]

[4]

[5]

[6]

[7]