Kurs:Räumliche Modellbildung/Gruppe 11

Gruppenseite - BSM

Diese Seite enthält die gesammelten Materialien der Gruppe 11 - BSM für die Portfolioprüfung. (Bitte noch nicht löschen. Diese Seite soll als Grundlage für die Modulprüfung am 09.06.21 genutzt werden. Vielen Dank)

COVID-19/Mathematische Modellbildung

Teilnehmer-innen

Böhmer, Hannah
Müller, Kevin
Schneider, Sebastian

1. Problemstellung

Ausgangslage unserer Modellierung ist die COVID-19-Pandemie. Seit März beschäftigt diese Infektion der Atemwege weltweit Menschen mit der Eindämmung des Virus. Vor allem sogenannte Risikogruppen sind stark von dem Virus betroffen. Zu den Risikogruppen zählen (Stand RKI 21.08.2020) ältere Personen, Raucher, stark adipöse Menschen und Personen mit gewissen Vorerkrankungen. Der Fall-Verstorbenen-Anteil liegt in Deutschland bei 4% und die Basisreproduktionszahl liegt zwischen 2 und 3,8. Da die Symptome (Husten, Fieber, Schnupfen, Störung des Geruchs- und Geschmackssinns und Pneumonie) der Infektion häufig einer gewöhnlichen Erkältung ähneln, führen viele Personen ihr Leben wie gewohnt weiter, wodurch eine schnellere Verbreitung der Infektion stattfindet. Die folgende Modellierung soll die Ausbreitung der Krankheit veranschaulichen. Zudem soll der Einfluss der Schutzmaßnahmen (vor allem Schließung der Schulen und Gaststätten, sowie Einschränkung der Freizeitaktivitäten) dargestellt werden.

2. Theoretische Grundlagen

Für die Modellierung der Pandemie wird zunächste eine Dichteverteilungsfunktion infizierter Personen gesucht. Integriert man über diese Funktion, erhält man die absolute Anzahl der Infizierten.
Der Infektionsprozess kann mathematisch beschrieben werden durch:

$\partial _{t}u(x,t)-div(c(x)\nabla u(x,t))=0$
Dies ist eine homogene partielle Differentialgleichung.
Hierbei gibt u die gesuchte Dichteverteilungsfunktion in Abhängigkeit von Ort und Zeit und c den Infektionskoeffizienten an.

$\partial _{t}u(x,t)$ stellt die Ableitung der orts- und zeitabhängigen Dichtefunktion nach der Zeit dar.
Die Divergenz beschreibt das Skalarprodukt des Nabla-Operators $\nabla$ (den Gradienten) mit einem Vektorfeld.
Der Gradient $\nabla$ ist definiert durch:

{\vec {\nabla }}=\left({\frac {\partial }{\partial x_{1}}},\ldots ,{\frac {\partial }{\partial x_{n}}}\right)

und gibt somit die Ableitung nach jeder Raumrichtung an.

Für den Spezialfall einer gleichmäßigen Infektion $c(x)=const$ erhalten wir für den stationären (zeitunabhängigen) Fall die Laplacegleichung $\Delta u(x)=0$ , da gilt:

$\Delta :=div\nabla$ .

Betrachtet man einen Infektionsprozess mit Quellterm, erhält man folgende inhomogene Differentialgleichung:

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

$f$ gibt dabei die Funktion der Infektionsverbreitung an.

Zur Lösung der gesuchten Infektionsdichtefunktion gibt es mehrere Ansätze:

Stationäre Lösung:

1. Laplace Gleichung in $n$ Dimensionen
Voraussetzung: Zeitunabhängigkeit, keinen Quellterm (Homogenität der Differentialgleichung) und gleichmäßige Infektion (siehe Spezialfall von oben)
Fundamentallösung:

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

Für $n=3$ gilt:
$\omega _{3}={\frac {4}{3}}\pi$

2. Poisson Gleichung
Voraussetzung: Zeitunabhängigkeit
Mit der Poisson Gleichung ist somit die Lösung folgender inhomogene Differentialgleichung $-\Delta u(x)=f(x)$ mittels Faltung möglich:
$u(x)=(\Phi *f)(x)=\int _{\mathbb {R} ^{n}}\Phi (x-y)f(y)dy.$
$\Phi$ entspricht der Fundamentallösung der Laplace Gleichung.
Wenn f ein kompakter Träger ist, gilt: $\lim _{\|x\|\to \infty }u(x)\to 0.$

Instationäre Lösung

1. Zeitabhängige Infektion
Voraussetzung: kein Quellterm (Homogenität der Differentialgleichung), gleichmäßige Infektion

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

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

2. Zeitabhängige Infektion mittels AWP
Voraussetzung: kein Quellterm (Homogenität der Differentialgleichung, Anfangsbedingung zum Zeitpunkt $t=0$

Lösungsformel:

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

3. Kontaktmodell

Einführung

Das Kontaktmodell ist ein einfaches Modell zur Veranschaulichung der Verbreitung des Corona-Virus.
Hierbei wird davon ausgegangen, dass sich die Personen zufällig im Raum mit konstanter Geschwindigkeit bewegen. Sobald eine Person zu nah an einen Infizierten herankommt, wird auch diese infiziert.

Octave Implementierung

Das Modell wurde folgendermaßen mittels Octave modelliert:

%Verbreitung der Infektion aufgrund zufälliger Bewegungen


close all 

clear all 

tic %Starte Zeitmessung 

%Anfangsparameter 

nx=10; %Anz. Breite 

ny=20; %Anz. Höhe 

nt=4; %Anz. Zeitschritte pro Zeiteinheit 

dt=1/nt; %Zeitschritt 

r_Inf=0.8; %Radius, bei der Infizierte Gesunde anstecken 

T=10; %Länge der Simulation in Zeiteinheiten 

v=3; %Geschwindigkeit, mit der sich die Fußgänger bewegen 

 

%Erstellen der Vektoren/Matrizen 

x=[1:nx]; 

y=[1:ny]; 

[x,y]=meshgrid(x,y); 



%Setze ersten Infizierten 

xInf=x(1,nx/2); 

yInf=y(ny/2,1); 



%Übergeben der KO des ersten Infizierten 

IndInf2=xInf; 

IndInf1=yInf; 



%Grafische Ausgabe Anfangssituation 

figure (1) 

hold on 

plot(x,y, '*b') 

plot(xInf,yInf, '*r') 

axis([-5 15 -5 25]) 

set(gca,'Visible','off') 

title (['Positionen nach 0s']) 

hold off 



%______________________________________________________________________________% 



%Simulation zufälliger Bewegungen und Ansteckung bei r<r_Inf 

%Erstelle Geschwindigkeitsmatrizen 



M=rand(ny,nx).-0.5; 

Mx=rand(ny,nx).-0.5; 

My=rand(ny,nx).-0.5; 

x_neu=x; 

y_neu=y; 



for i=1:T*nt 



  %Bestimmen neuer Richtungen nach jeder ganzen Zeiteinheit 

  if mod(i,4)==0 

    Mx=rand(ny,nx).-0.5; 

    My=rand(ny,nx).-0.5; 

  endif 



 %____________________________________________________________________________% 

 %Vektorpfeile anzeigen lassen - Teil 1 

  %x_alt=x_neu 

  %y_alt=y_neu   

 %____________________________________________________________________________% 

 

 %Berechne neue Pos für jeden Punkt 

  x_neu=x_neu.+(v*Mx)*dt; 

  y_neu=y_neu.+(v*My)*dt; 
   
 

  %Speichere den letzten Plot 

  test=["Fig_", num2str(i)]; 

  saveas(i, test, 'jpg') 

 

  figure (i+1) 

  title (['Positionen nach ' num2str(i*dt) 's']) 

  axis([-5 15 -5 25]) 

  hold on 

  plot(x_neu,y_neu, '*b') 



 %____________________________________________________________________________% 

  %Vektorpfeile anzeigen lassen - Teil 2 

   %quiver(x_neu,y_neu,x_neu-x_alt,y_neu-y_alt); 

  %____________________________________________________________________________% 

 

  %Überprüfen ob Infiziert und wenn, dann markiert 

  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) 

      r=norm([x_neu(l,j)-x_neu(IndInf1(k1),IndInf2(k1)), y_neu(l,j)-y_neu(IndInf1(k1),IndInf2(k1))]); 

     

        if r<r_Inf 

         %wenn der Abstand zum Infizierten unter dem Ansteckunsradius liegt, wird  der l,j-te  Punkt infiziert (mit rotem Stern gekenzeichnet) 

            

            plot(x_neu(l,j),y_neu(l,j), '*r') ; 

            set(gca,'Visible','off') %Achsanzeige deaktiviert 

            %-------------Speicheroptimierung-------------- 

            k2=1; 

            while k2 <=length(IndInf1)&&(l~=IndInf1neu(k2) || j~=IndInf2neu(k2)) 

             k2=k2+1; 

            endwhile 

            if k2 >length(IndInf1) 

              IndInf1=[IndInf1 l];

              IndInf2=[IndInf2 j]; 

            endif 

        endif 

     endfor 

  endfor 

 

 

  endfor 


  hold off 

 

endfor 


test=["Fig_", num2str(i+1)]; %Speichere letztes Bild 

saveas(i+1, test, 'jpg') 

toc %Beende Zeitmessung

Darstellung Kontaktmodell

Darstellung des Kontaktmodells:

Darstellung des Kontaktmodells mit erhöhter Geschwindigkeit:

In dieser Darstellung wurde die Geschwindigkeit erhöht. Hierbei ist eine schnellere Infektion zu beobachten, wodurch die Einschränkung der Freizeitaktivitäten zu legitimieren ist.

Darstellung des Kontaktmodells mit verringertem Infektionsradius:

In dieser Darstellung wurde der Infektionsradius verringert. Hierbei ist eine verlangsamte Ausbreitung des Virus zu beobachten. Diese Modellierung soll die Übertragung während der Mundschutzpflicht darstellen. Die Viren werden überwiegen über Nase und Mund aufgenommen und abgegeben. Durch das Bedecken des Mund und Nasenraums verringert sich der Infektionsradius.

Problem

In diesem Modell wird die Genesung nach 14 Tagen nicht berücksichtigt. Es wird davon ausgegangen, dass die Personen weiterhin bei zu geringem Abstand weitere Personen infizieren. Zudem betrachtet dieses Modell weder Sterbe- noch Geburtenrate.

4. SI Modell

Einführung

Eine geeignete Modellierung bietet das SI Modell (Kompartimentmodell).
Wir betrachten hier nur das begrenzte oder auch logistische Wachstumsmodell.

In diesem Modell wird zunächst die Gesamtpopulation B betrachtet. Wir wählen eine geschätze Zahl von 2/3 der Gesamtbevölkerung.

Somit gilt: B= 55.000.000.
Diese Gesamtpopulation wird eingeteilt in die Infektionsanfälligen S (susceptible) und die kummulierten Infizierten I (infected).

Sobald sich eine Person aus der Gruppe der Infektionsanfälligen infiziert hat, wechselt er irreversibel in die Gruppe der Infizierten.
Zu Beginn der Modellierung $t_{0}$ ist eine bestimmte Anzahl an Personen $I_{0}$ infiziert. Somit gilt:
$S(t_{0})=B-I_{0}$
Zudem gilt:

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

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

und
$I(t_{0})=I_{0}.$

Octave Implementierung

Folgende Modellierung wurde mittels Octave gemacht:

In dieser Modellierung wird mit 16 Infizierten gestartet.

clear all ; 



%Hole RKI-Daten von coronaData: ; 

%RKI-Daten: ; 

RKI=coronaData ; 

RKI=[0:size(RKI,2)-1;RKI]; %[Tage;Infizierte;Tote;Geheilte];


% Festlegung der Parameter:

dt=0.1;                     %Zeitschritt in d 

I0=16;                      %Anfangsinfizierte 

B=(2/3)*280;                %Kapazität in Tausend 

k=0.22;                     %Infektionsrate 



%Zeitarray: 

tmax=180; 

t=0:dt:tmax; 



%Anfangswerte:

y0=[B-I0/1000; I0/1000]; 



%logistisches Modell (DGL-System): 

f1=@(y) [-k*y(1)*y(2)/B ; k*y(2)*y(1)/B]; %y(1)=S, y(2)=I 



%Erzeuge numerische Lösung:

y=lsode(f1,y0,t); 



%Graphische Ausgabe: 

figure() 

plot(t,y(:,1),t,y(:,2),RKI(1,:),RKI(2,:)/1000,'.',"linestyle", "none") 

hold on 

xlabel('Zeit [d]') 

axis([0 tmax 0 B])

ylabel('Bevölkerung (in Tsd.)')

legend(['Infektionsanfaellige';'Infizierte (kummuliert)';'RKI-Daten'],"location","east") 

title(['SI-Modell in',num2str(tmax),' Tagen mit ',num2str(I0),...

' Anfangsinfizierten und Infektionsrate k=',num2str(k)])

hold off

Darstellung SI Modell

In dieser Darstellung sind neben den modellierten Zahlen für Infektionsanfällige und Infizierte die Daten des RKI aufgetragen. Hierbei ist eine gute Übereinstimmung nur zu beobachten, wenn man die Zahl der Gesamtpopulation (B) reduziert. Deshalb wurde bei der Modellierung eine Gesamtbevölkerung von 280.000 angenommen.

Problem

Das Problem dieses Modells ist, dass weder Tote noch Geburten beachtet werden. Die gestrichelte Linie zeigt die RKI Daten. In weitester Betrachtung stimmen die Funktionen zwar überein, jedoch sollte eine Verbesserung des Modells vorgenommen werden.

5. SIR Modell

Einführung

Das SIR Modell ist ebenso wie das SI Modell den Kompartimentmodellen zuzuordnen.
Zusätzlich zur Gruppe der Infektionsanfälligen (S) werden hier auch die aktuell Infizierten (I) betrachtet und eine Gruppe der Genesenen inklusive der Verstorbenen (R) ("removed") berücksichtigt.

Zu Beginn der Modellierung $t_{0}$ ist eine bestimmte Anzahl an Personen $I_{0}$ aktuell infiziert. Die Anzahl der Genesenen bzw. Verstorbenen wird zunächst null gesetzt. Somit gilt:
$S(t_{0})=B-I_{0},$

$R(t_{0})=0$

Zudem gilt:

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

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

$R'(t)=wI(t),$
wobei $w=1/T$ der Kehrwert der Genesungszeit T ist.

Weiterhin gilt:
$I(t_{0})=I_{0}.$

Octave Implementierung

Folgende Modellierung wurde mit Octave erstellt:

%Hole RKI-Daten von coronaData: 

%RKI-Daten: 

RKI=coronaData; 

RKI=[0:size(RKI,2)-1;RKI]; %[Tage;Infizierte(kum);Tote;Geheilte] 

RKI=[RKI;RKI(2,:)-(RKI(3,:)+RKI(4,:))] %Erzeuge Aktuell Infizierte 



%Parameter: 

dt=0.1;                          %Zeitschritt in d 

I0=RKI(2,1);                     %Anfangsinfizierte 

k=0.3;                          %Infektionsrate 

R0=0;                            %Anfangsgenesene bzw. -tote 

B=(2/3)*250;                   %Kapazität in Tausend 

T=14;                            %Genesungsdauer in Tagen 

w=1/T;                           %Genesungsrate 



%Zeitarray:

tmax=180; 

t=0:dt:tmax; 





%Anfangswerte: 

y0=[B-(I0+R0)/1000,I0/1000,R0/1000]; 


 
f2=@(y,x) [-k*y(2)*y(1)/B ; k*y(2)*y(1)/B-w*y(2) ; w*y(2)];  %y(1)=S, y(2)=I, y(3)=R 



%Erzeuge numerische Lösung 

y=lsode(f2,y0,t); 

 

%Graphische Ausgabe: 

figure() 

plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,2)+y(:,3),RKI(1,:),RKI(2,:)/1000,'.',... 

RKI(1,:),(RKI(3,:)+RKI(4,:))/1000,'.',RKI(1,:),RKI(5,:)/1000,'.')

hold on 

xlabel('Zeit [d]') 

axis([0 tmax 0 B]) 

ylabel('Bevölkerung (in Tsd.)') 

legend(['Infektionsanfaellige';'Infizierte';... 

'Genesene bzw. Verstorbene';'Infizierte (kummuliert)';... 

'RKI Infizierte';'RKI Genesene/Verstorbene';'RKI Infiziert (kum.)'],"location","east") 

title(['SIR-Modell in',num2str(tmax),' Tagen mit ',num2str(I0),... 

' Anfangsinfizierten und Infektionsrate k=',num2str(k)]) 

hold off

Darstellung SIR Modell

Problem

Ebenso wie im SI Modell stimmen die Zahlen des RKI nur mit der Kurve überein, wenn man die Gesamtbevölkerung geringer wählt. Dies ist auf eine hohe Dunkelziffer zurückzuführen.

6. Modifiziertes SIR Modell

Einführung

Das modifizierte SIR Modell berücksichtigt einerseits, dass die gestorbenen Infizierten nicht mehr in der Gruppe der Genesenen erfasst werden. Da nicht alle Infizieten durch Tests identifiziert werden, wird zusätzlich die Dunkelziffer einbezogen. Somit gilt:

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

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

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

$S(t_{0})=B-I_{0}$

$I(t_{0})=I_{0}$

$R(t_{0})=0$

Hierbei ist $c=kB$ , wobei k und B wie im SI bzw. SIR Modell gewählt werden.
d bezeichnent die Sterbrate der Infizierten und r den Anteil der durch Tests identifizierten Infizierten.

Wählt man r=1 und d=0 erhält man das SIR Modell.

Octave Implementierung

Folgendes Modell wurde mit Octave erstellt:

%Implementierung SIR-Modell modifiziert:



clear all 



%Hole RKI-Daten von coronaData: 

RKI=coronaData; 

RKI=[0:size(RKI,2)-1;RKI]; %[Tage;Infizierte(kum);Tote;Geheilte] 

RKI=[RKI;RKI(2,:)-(RKI(3,:)+RKI(4,:))]; %Erzeuge Zeile aktuell Infizierte 



%Parameter:
dt=0.1;                          %Zeitschritt in d 

tmax=180;                        %Betrachtungsdauer in d 

I0=RKI(2,1);                     %Anfangsinfizierte 



d=0.003;                         %Sterberate 

R0=0;                            %Anfangsgenesene bzw. -tote 

B=(2/3)*83000;                   %Kapazität in Tausend 

T=14;                            %Genesungsdauer in Tagen 

w=1/T;                           %Genesungsrate 

r=0.9                            %rel. Anteil durch Tests erfasste Infizierte 



k=[0.3;0.1;0.6];                %Infektionsraten (wurde durch Ausprobieren festgelegt) 

it1=56;                          %Ende 1./Start 2. Infektionsperiode (56=22.03.) 

it2=122;                         %Ende 2./Start 3. Infektionsperiode (122=27.05.) 



y0=[B-(I0+R0)/1000,I0/1000,R0/1000];  %Anfangswerte [S,I,R] 

 

%________________________% 

%Zeitarrays: 

t1=0:dt:it1-dt; 

t2=it1:dt:it2-dt; 

t3=it2:dt:tmax; 

t=0:dt:tmax;        
         
 

%Modifiziertes SIR Modell (DGL-System): 

f=@(y,x) [-k(1)*y(2)*y(1)/(r*B) ; k(1)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)];  %y(1)=S, y(2)=I, y(3)=R 

 

%Erzeuge numerische Lösungen: 

y1=lsode(f,y0,t1);  %Berechne 1. Phase 

f=@(y,x) [-k(2)*y(2)*y(1)/(r*B) ; k(2)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)]; %Ändern der Infektionsrate 

y2=lsode(f,y1(length(t1),:),t2);  %Berechne 2. Phase 

f=@(y,x) [-k(3)*y(2)*y(1)/(r*B) ; k(3)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)]; %Ändern der Infektionsrate 

y3=lsode(f,y2(length(t2),:),t3);  %Berechne 3. Phase 

y=[y1;y2;y3]; 

%yy=ones(length(y(1)),1)*B; 

%yy=yy-(y(:,1)+y(:,2)+y(:,3)); %Berechnen der Verstorbenen (?) 

%y=[y,yy]; 

 

%Graphische Ausgabe 

figure() 

plot(t,y(:,1),"b",t,y(:,2),"r",t,y(:,3),'color','g',t,y(:,2)+y(:,3),'color','c',... 

RKI(1,:),RKI(5,:)/1000,'.','color','r',RKI(1,:),(RKI(3,:)+RKI(4,:))/1000,'.','color','g',... 

RKI(1,:),RKI(2,:)/1000,'.','color','c') 

xlabel('Zeit [d]') 

axis([0 tmax 0 B]) 

ylabel('Bevölkerung (in Tsd.)') 

legend(['Infektionsanfaellige';'Infizierte';... 

'Genesene bzw. Verstorbene';'Infizierte  (kummuliert)';'Verstorbene';... 

'RKI Infizierte';'RKI Genesene/Verstorbene';'RKI Infizierte (kum.)'],"location","east") 

title(['Modifiziertes SIR-Modell in ',num2str(tmax),' Tagen mit ',num2str(I0),... 

' Anfangsinfizierten und Kapazität bei ',num2str(round(B*1000)),'.']) 

%_______________________% 

%Da r den relativen Anteil der von den Tests erfassten Infizierten darstellt, gilt 0<r<1. 

%In dieser Darstellung ist r mit 0,9 recht groß gewählt. 

%Dies bedeutet, dass der Test fast alle Infizierte erfassen würde. 

%Hierdurch ist die Konvergenz am rechten Rand deutlich höher, wie wenn r=0,1 o.ï¿½. gewählt wird. 

%Gibt man für r=0,1 ein, sieht man eine vergleichsweise geringe Zahl der Infizierten. 

%Jedoch würde man dann eine große Dunkelziffer haben. 



%______________________________________________________________________________% 

%Anpassung der Kapazität und Infektionsrate an RKI-Daten: 

%Implementierung SIR-Modell modifiziert: 

clear all  







%Hole RKI-Daten von coronaData: 

RKI=coronaData; 

RKI=[0:size(RKI,2)-1;RKI]; %[Tage;Infizierte(kum);Tote;Geheilte] 

RKI=[RKI;RKI(2,:)-(RKI(3,:)+RKI(4,:))]; %Erzeuge Zeile Aktuell Infizierte 



%Parameter: 

dt=0.1;                          %Zeitschritt in d 

tmax=180;                        %Betrachtungsdauer in d 

I0=RKI(2,1);                     %Anfangsinfizierte 



d=0.003;                         %Sterberate 

R0=0;                            %Anfangsgenesene bzw. -tote 

B=(2/3)*300;                     %Kapazität in Tausend 

T=14;                            %Genesungsdauer in Tagen 

w=1/T;                           %Genesungsrate 

r=0.9;                            %rel. Anteil durch Tests erfasste Infizierte 



k=[0.33;0.15;0.45];             %Infektionsraten (wurde durch Ausprobieren festgelegt) (3. Wert für Rate nach Lockerungen) 

it1=56-17;                      %Ende 1./Start 2. Infektionsperiode (56=22.03.) 

it2=122-45;                     %Ende 2./Start 3. Infektionsperiode (122=27.05.) 



%Durch Rückdatierung Anpassung an RKI-Daten 



y0=[B-(I0+R0)/1000,I0/1000,R0/1000];  %Anfangswerte [S,I,R] 



%________________________% 

%Zeitarrays: 

t1=0:dt:it1-dt; 

t2=it1:dt:it2-dt; 

t3=it2:dt:tmax; 

t=0:dt:tmax;        
         


%Modifiziertes SIR Modell (DGL-System): 

f=@(y,x) [-k(1)*y(2)*y(1)/(r*B) ; k(1)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)];  %y(1)=S, y(2)=I, y(3)=R 

 

%Erzeuge numerische Lösungen: 

y1=lsode(f,y0,t1);  %Berechne 1. Phase 

f=@(y,x) [-k(2)*y(2)*y(1)/(r*B) ; k(2)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)]; %Ändern der Infektionsrate 

y2=lsode(f,y1(length(t1),:),t2);  %Berechne 2. Phase 

f=@(y,x) [-k(3)*y(2)*y(1)/(r*B) ; k(3)*y(2)*y(1)/B-w*y(2) ; w*y(2)-d*y(2)]; %Ändern der Infektionsrate 

y3=lsode(f,y2(length(t2),:),t3);  %Berechne 3. Phase 

y=[y1;y2;y3]; 

%yy=ones(length(y(1)),1)*B; 

%yy=yy-(y(:,1)+y(:,2)+y(:,3)); %Berechnen der Verstorbenen 

%y=[y,yy]; 



%Graphische Ausgabe 

figure() 

plot(t,y(:,1),"b",t,y(:,2),"r",t,y(:,3),'color','g',t,y(:,2)+y(:,3),'color','c',... 

RKI(1,:),RKI(5,:)/1000,'.','color','r',RKI(1,:),(RKI(3,:)+RKI(4,:))/1000,'.','color','g',... 

RKI(1,:),RKI(2,:)/1000,'.','color','c') 

xlabel('Zeit [d]') 

axis([0 tmax 0 B]) 

ylabel('Bevölkerung (in Tsd.)') 

legend(['Infektionsanfaellige';'Infizierte';... 

'Genesene bzw. Verstorbene';'Infizierte  (kummuliert)';... 

'RKI Infizierte';'RKI Genesene/Verstorbene';'RKI Infizierte (kum.)'],"location","east") 

title(['Modifiziertes SIR-Modell in ',num2str(tmax),' Tagen mit ',num2str(I0),... 

' Anfangsinfizierten und Kapazität bei ',num2str(round(B*1000)),'.'])

Darstellung modifiziertes SIR Modell

Problem

Bei diesem Modell wird die Geburtenrate und die allgemeine Sterberate nicht berücksichtigt. Zudem wird die räumliche Infektionsverbreitung sowie die unterschiedliche räumliche Bevölkerungsdichte vernachlässigt.

7. Modell des Reaktionsdiffusionsprozess

Einführung

Das Modell der Reaktionsdiffusion wird auch als populationsbasiertes räumlich-kontinuierliches Modell bezeichnet. Bei dem Reaktionsdiffusionsprozess wird die Reaktionsdiffusionsgleichung herangezogen. Dies ist eine inhomogene Gleichung mit Quelldichtefunktion (vgl. Abschnitt 2.). Diese ist folgendermaßen definiert:

$\partial _{t}u(x,t)-div\ (c(x)\nabla u(x,t))=f(u(x,t))$
f wird auch als Reaktionsterm bezeichnet.

Im Gegensatz zu den rein dynamischen Kompartimentmodellen wird in diesem Modell der räumliche Austausch betrachtet. Dies liegt daran, dass die Kompartimentmodelle auf der Massenerhaltung beruhen.
Somit beschreibt dieses Modell nicht nur die zeitliche, sondern auch die räumliche Ausbreitung.

Wir betrachten auch in diesem Modell wieder nur das beschränkte Wachstum. Somit gilt für den Reaktionsterm:
$f(u(x,t))=k.u(x,t){\big (}B(x)-u(x,t){\big )},$
wobei u die Dichtefunktion der kumulierten Infizierten, k die logistische Infektionsrate und B die Dichtefunktion der Gesamtpopulation angibt.

Teil man die kumulierten Infizierten in aktuell Infizierte und Genesene bzw. Tote (vgl. SIR Modell), erhält man für die Dichtefunktionen der Infektionsanfälligen $u_{S}$ und der aktuell Infizierten $u_{I}$ folgende Differentialgleichungen:

$\partial _{t}u_{S}(x,t)-div(c(x)\nabla u_{S}(x,t))=f_{S}(u_{S}(x,t),u_{I}(x,t))$
$\partial _{t}u_{I}(x,t)-div(c(x)\nabla u_{I}(x,t))=f_{I}(u_{S}(x,t),u_{I}(x,t))$

Die Reaktionsterme sind definiert durch:

$f_{S}(u_{S}(x,t),u_{I}(x,t))=-k.u_{I}(x,t)u_{S}(x,t)$ und
$f_{I}(u_{S}(x,t),u_{I}(x,t))=k.u_{I}(x,t)u_{S}(x,t)-wu_{I}(x,t)$
Hierbei gibt w die Wechselrate zwischen I und R an.

Randnotizen:

Differenzenquotienten:

erster zenraler Differenzenquotient: $D_{h}f(x)=f'(x)+{\cal {O(h^{2}):D_{h}f(x):={\frac {f(x+h)-f(x-h)}{2h}}}}$

zweiter zentraler Differenzenquotient: $D_{h}^{2}f(x)=f''(x)+{\cal {O(h^{2}):D_{h}^{2}f(x):={\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}}}$

vorwärts(+) / rückwärts(-) Differentienquotient: $D_{h}^{\pm }f(x)=f'(x)+{\cal {O(h):D_{h}^{+}f(x):={\frac {f(x+h)-f(x)}{h}}bzw.D_{h}^{-}f(x):={\frac {f(x)-f(x-h)}{h}}}}$

Octave Implementierung

%Implementierung Diffusion-Reaktion 



close all 

clear all 



%________________________________Startparameter________________________________% 

%Definiere Groesse Anfangsgebiet und Schrittweiten 

xend=5; 

yend=10; 

dx=0.1; 

dy=dx; 

[x,y]=meshgrid(0:dx:xend,0:dy:yend); 

N=size(x,2); %Anzahl x-Schritte 

M=size(y,1); %Anzahl y-Schritte 





%Definiere Dauer, Zeitschrittlaenge und Anz. Zeitschritte 

T=200; 

dt=0.1; 

nt=T/dt; 

 

%Definiere (konst.) Diffusionskoeffizienten 

D=0.01; 

c=0.3; %Infektionsrate 

w=1/14; %Wechselrate 

t_lock=28; %Tag des Lockdowns 

 

%___________Einlesen, interpolieren und anpassen Bevoelkerungsdichte___________% 
 
pic=imread('abc.png'); 

pic=pic(:,:,2); 

pic=double(pic); 

[s1 s2]=size(pic); 

Max=max(max(pic)); 

pic_int=interp2(pic,xend*(1*x.+2)/(xend+2),yend*(1*y.+2)/(yend+2), 'cubic')/Max+0.1; 

pic_int=flipud (pic_int); 

pic_int=reshape(pic_int',M*N,1); 

pic_int=(pic_int+0.5)/1.5; %Korrektur Dichte 





%________________________sonstige Parameter____________________________________% 



k=1; %Zaehler graphische Ausgaben 

p=10/dt; %Anz. Zeitschritte bevor neues Bild ausgegeben wird 

%______________________________________________________________________________% 




%_______________________Erzeugen der Systemmatrix______________________________% 

%Hier noch ohne Randbedingungen 

B=diag(-4.*ones(N*M,1),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); 

%Korrektur 

for i=1:(M-1)  

    B(i*N+1,i*N)=0; %Korrektur linker Rand 

    B(i*N,i*N+1)=0; %Korrektur rechter Rand 

endfor 



%Gewichtung mit Diff.koeff und Schrittweite 

B=B*D/dx^2; 

 

%______________________________________________________________________________% 

%Neumann'sche Randbedingungen: Setze die Richtungsableitung an allen Raendern 0 

%Neumann-Korrektur: 

for i=1:N 

  B(i,i)=B(i,i)+1; %Oberer Rand 

endfor 



for i=N*(M-1)+1:N*M 

  B(i,i)=B(i,i)+1; %Unterer Rand 

endfor 

 

for i=1:N:N*(M-1)+1 

  B(i,i)=B(i,i)+1; %Linker Rand 

endfor 

 

for i=N:N:N*M 

  B(i,i)=B(i,i)+1; %Rechter Rand 

endfor 

%______________________________________________________________________________% 

  



%Erzeugen der Anfangsloesung 

for  i=1:N 

  for j=1:M 

   start(j,i)=0.005*anfang2(x(j,i),y(j,i)); 

  endfor   

endfor  

 

%Erzeuge 'Vektor' aus Anfangssituation 

inf_alt=1*reshape(start',N*M,1); 

sus_alt=pic_int-inf_alt; 

inf_anz=sum(inf_alt); 

sus_anz=sum(sus_alt); 

t=0:dt:T; 



%Reaktions-Term/Funktion (SIR)

fS=@(uI,uS,x) -c*slowdown(x,t_lock)*uI.*uS./pic_int; 

fI=@(uI,uS,x) c*slowdown(x,t_lock)*uI.*uS./pic_int-w*uI; 





%_________________________Berechnung Diff-Reak_________________________________% 

for i=1:nt 

  sus=sus_alt+B*sus_alt*dt+fS(inf_alt,sus_alt,i*dt)*dt; 

  inf=inf_alt+B*inf_alt*dt+fI(inf_alt,sus_alt,i*dt)*dt; 

  inf_anz=[inf_anz sum(inf)]; 

  sus_anz=[sus_anz sum(sus)]; 

  sus_alt=sus; 

  inf_alt=inf; 

  sus_mat(:,i)=sus; 

 inf_mat(:,i)=inf; 

   

  %Ausgabe jedes p-ten Bildes 

  if mod(i,p)==0 

    matrix_ausgabe=reshape(inf_mat(:,i),N,M)';% Matrix mit N(Spalten)x M(Zeilen) 

    figure(k) 

    surfc(x,y,matrix_ausgabe, "EdgeColor", "none") 

    colormap("jet") 

    colorbar 

    axis([0 xend 0 yend -0.001 1]) 

    %caxis([0 0.3]) 

    title (["Normierte Infiziertendichte nach ", num2str(dt*i), "Tagen"]); 

    ylabel("y in km") 

    xlabel("x in km") 

    view(2) 

    test=["Fig_", num2str(k)]; 

    saveas(k, test, 'jpg') 

    

    k=k+1; 

  endif 

 

endfor 



figure(k) 

plot(t,inf_anz,t,sus_anz) 

xlabel('Zeit in Tagen') 

ylabel('BevÃ¶lkerungsdichte pro Flaecheneinheit') 

legend(['Infiziertendichte';'Anfaelligendichte']) 

title('Darstellung der Dichten in Abhaengigkeit der Zeit') 

test=["Fig_", num2str(k)]; 

saveas(k, test, 'jpg')

Graphische Darstellung

Problem

In dieser Modellierung wurde ein konstanter Diffusionskoeffizient angenommen. Jedoch hängt dieser in der Realität von der Bevölkerungsdichte ab. Dieses Problem kann behoben werden, indem der Diffusionskoeffizient in Abhängigkeit der Populationsdichte modelliert wird.

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^[1] oder R/Studio

Quellenangaben

Notieren Sie hier die von Ihnen verwendete Literatur

Boto von Querenburg (2013) Mengentheoretische Topologie - Springer-Lehrbuch
Wikipedia-Artikel zu Epidemiologie (2020)^[2]

Literatur

↑ Silva, I., & Moody, G. B. (2014). An open-source toolbox for analysing and processing physionet databases in matlab and octave. Journal of open research software, 2(1).
↑ „Epidemiologie“. In: Wikipedia, Die freie Enzyklopädie. Bearbeitungsstand: 25. Juni 2020, 09:27 UTC. URL: https://de.wikipedia.org/w/index.php?title=Epidemiologie&oldid=201288758 (Abgerufen: 25. Juni 2020, 10:24 UTC)

[1] Silva, I., & Moody, G. B. (2014). An open-source toolbox for analysing and processing physionet databases in matlab and octave. Journal of open research software, 2(1).

[2] „Epidemiologie“. In: Wikipedia, Die freie Enzyklopädie. Bearbeitungsstand: 25. Juni 2020, 09:27 UTC. URL: https://de.wikipedia.org/w/index.php?title=Epidemiologie&oldid=201288758 (Abgerufen: 25. Juni 2020, 10:24 UTC)

[1]

[2]