Kurs:Räumliche Modellbildung/Aufgaben Tutorium kontinuerliche räumliche Modellierung

Übungsaufgaben zur kontinuerlichen Modellierungsansätzen

Die Aufgaben sollen in 3er Gruppen mithilfe von Octave bzw. als Jupyter Notebook mit Octave Kernel bearbeitet werden. In der kollaborativen Umgebung CoCalc via Kurs:räumliche_Modellierung wird der Austausch der bearbeiteten und korrigierten Skripte stattfinden. Alternativ können diese über Teilnehmerordner SOSE 21 des OLAT-Kurses Diffusionsprozesse hochgeladen werden.

1. Übungsaufgabe SW1: Kontakt- Modell

Besprechung im SW3

Animation: Räumliche Ausbreitung einer Infektion durch direkten Kontakt von sich bewegenden Fußgängern - (2020) Anna Hundertmark

Animation: Räumliche Ausbreitung einer Infektion durch direkten Kontakt von sich bewegenden Individuen mit wechselnden Richtungen - (2021) Anna Hundertmark

Machen Sie sich vetraut mit dem Software Octave, gehen Sie dabei nach dem internativen Jupyter Notebook Einfuhrung in Octave.ipynb im Cocalc Kurs:räumliche_Modellierung. Diese Datei finden Sie auch im OLAT Materialordner SOSE 21/Tutorium. Alternativ können Sie direkt nach dem dort enthaltenem Skript zur Einführung in Octave Octave_intro.pdf vorgehen.
Testen Sie das Octave Skript zum Kurs:Räumliche_Modellbildung/Deterministisches_Kontaktmodell für die räumliche Verbreitung der Infektion durch zufällige Bewegung der Individuen und Kontakt und passen Sie diese an. Erstellen Sie dabei Simulationen für verschiedene Simulationszeiten, Ansteckungsradien, Bewegungsgeschwindigkeiten und Bewegungsrichtungen. Als Tutorium zur Bearbeitung dieser Aufgabe dient das interaktive Jupyter Notebook Kontakt Modell.ipynb in unserem CoCalc-Kurs und die dazugehörige statische Datei Kontakt_Modell_Ausarbeitung.pdf aus

OLAT Materialordner SOSE 21/Tutorium. Folgende Befehle und Funktionen sind nutzlich:

Legen Sie die Anzahl von sich bewegenden Punkten durch die Anzahl der x- und,y-Koordinaten der Punkten fest:

Nx=...;
Ny=..;

Erstellen Sie die Koordinaten-Vektoren der Punkten:

x=[1:Nx];
y=[1:Ny];

Durch diesen Befehl wird der Punktegitter erstellt; x und y Koordinaten der Punkte sind Matrizen

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

Dieser Befehl generiert die Matrix der Zufallszahlen der Dimension Nx mal Ny:

rand(Ny,Nx)

Der Abstand von einem infizierten Individuum(Punkt) mit Koordinaten (Xinf, Yinf) zum (i,j)-ten Punkt mit Koordinaten (x(i,j),y(i,j)) wird Mithilfe der Funktion für Euklidische Norm berechnet:

Abstand=[x(i,j)-Xinf y(i,j)-Yinf];                                                                                      
norm(Abstand,2)

2. Übungsaufgabe SW2 - Fundamentallösungen

Besprechung SW4
a) Implementieren Sie die Formel für die Fundamentallösungen $\Phi ,\Psi$ der Laplace und der instationären Diffusionsgleichung in Octave. Stellen Sie diese graphisch als Funktionen auf $\mathbb {R} ^{2}$ bzw. $\mathbb {R} \times \mathbb {R} ^{+}$ dar. Achten Sie dabei darauf, dass Sie die Singularitätspunkte aus dem Definitionsbereich (aus dem Punktegitter) ausschneiden, d.h. dass die Singularitätspunkte nicht im Gitter enthalten sind.

Bei der Imlementierung in Octave können Sie den Gitter herstellen:

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

Definition der Funktion des natürlichen Logarithmus und der Exponenzialfunktion $\log(\|x\|),\ e^{x/t}$

f=@(X,Y) log((X).^2+(Y).^2))
g=@(X,t) exp(X./t)

Graphische Befehle, hier x,y sind die Koordinatenmatrizen, z ist die Wertematrix.

meshc(x,y,z)
surface (x,y,z)

b) Implementieren Sie die Lösungsdarstellunsformel für Poissongleichung $(\Phi *f)(x)$ und für den Anfangswertproblem der instationären Diffusionsgleichung $(\Psi *u_{0})(x,t)$ mithilfe numerischen Integration (z.B. Trapezregel). Wählen Sie beliebige beschränkte Quellfunktion $f$ bzw. Anfangsbedingung $u_{o}$ mit kompakten Träger.

Numerische Integration über ein Rechteck $\int \int u(x,y)dxdy$ mithilfe der Trapezregel:

trapz(y,trapz(x,u,2))

Hinweis

Im OLAT Materialordner SOSE 21/Tutorium/Fundamenallösunen finden Sie den Jupyter Notebook Losungsformel_Poisson.ipynb mit der Implementierung der Poissonformel $(\Phi *f)(x)$ , den Sie für die Implementierung der Fundamentrallösung der Diffusionsgleichung $(\Psi *u_{0})(x,t)$ (in einer Zeitschleife) umschreiben können.

3. Übungsaufgabe SW5: Kompartimentmodelle

Besprechung/Abgabe: SW7.

Verwenden Sie das modifizierte SIR Modell zur Modellierung der Infektionsverbreitung in Deutschland unter der Annahme, dass die gemeldeten Infizierten-Daten nur ein Bruchteil $0<r<1$ der tatsächtlichen Anzahl der Infizierten erfassen (den Rest bilden die nicht-erkannte / asymptomatische / nicht gemeldete Infizierten.) Kalibrieren Sie dieses Modell, sodass die zugehörige numerische Lösung die vom Robert Koch Institut gemeldeten Daten zur Corona-Infizierten approximert und treffen Sie Zukunftsvorhersagen in dem Sie auch Impfstrategie im Modell abbilden:

Als Anfangsbedingungen verwenden Sie die vom Robert Koch Institut gemeldeten Daten am Anfang der Datenerfassung.
Wählen Sie den Faktor $r$ - der Datenerfassung, dieser kann sich im Laufe der Pandemie ändern, da im späterem Velauf der Pandemie mehr getestet wurde.
Verwenden Sie eine zeitabhängige Infektionsrate $c.s(t)$ , wobei $c$ die Basisinfektionsrate, die das exponenziellen Wachstum am Anfang der Pandemie modelliert und $s(t)$ den Einfluß der Kontaktreduzierenden Maßnahmen (Schulschließungen, Ausgangsperren usw.), oder des warmes Wetters an die Infektionsverbreitung.
Implementieren Sie die numerische Lösung des modifizierten SIR-Modell in Octave zum Beispiel mit Hilfe von Octave-Function lsode() und wählen Sie die Infektionsrate $c.s(t)$ und Erfassungsfaktor $r$ bzw. $r(t)$ passend, sodass die Kurven der numerischen Lösung den RKI Daten mindestens bis zum Ende 2020 näherunsweise enstprechen (betrachten Sie hierzu den relativen Fehler).
Ergänzen Sie das modifizierte Modell um Terme, die zusätzlich die Wirkung der Impfstrategie sei Anfang 2021 berücksichtigen und treffen Sie Zukunftsvorhersagen für unterschiedliche Impfraten.

Hinweis

Im OLAT Materialordner SOSE 21/Tutorium/SIR_Modell finden Sie den Jupyter Notebook SIR_Modell.ipynb zur numerischen Lösung des SIR Modells und die Octave Funktion coronaData.m zu den RKI Daten vom Anfang der Pandemie. Diese Funktion sollen Sie um Daten aus dem Sommer, Herbst, Winter 202o und Frühjahr 2021 aktualisieren

4. Übungsaufgabe SW7: Poissongleichung mit FDM

Besprechung/Abgabe: SW 9

Hinweis
Im OLAT Materialordner SOSE 21/Tutorium/FDM finden Sie den Jupyter Notebook FDM_Poissongleichung_Dirichlet.ipynb zur numerischen Diskretisierung der stationärer Poissongleichung mit homogenen Dirichlet- Randbedingungen. Diese können Sie als Grundlage für die Bearbeitung der folgenden Teilaufgaben a) und b) verwenden und adaptieren.

a)

Implementieren Sie in Octave die Finite-Differenzen-Methode (FDM) für Dirichlet-Randwertaufgabe mit inhomogenen Randbedingungen auf einem rechteckigem Gebiet $D=(0,X_{end})\times (0,Y_{end})$ .

Benutzen Sie äquidistantes Gitter der inneren Punkte $\{(x_{i},y_{j}),\ i=1,\ldots N,\ j=1,\ldots M\}$ , wobei $x_{i}=i.h,\ x_{N+1}=X_{end},\ y_{j}=j.h,\ y_{M+1}{\stackrel {\approx }{=}}Y_{end}$ .

Gehen Sie wie folgt vor

1. Ordnen Sie gedanklich die Matrix der unbekannten Werte $u_{i,j},\ i=1,\ldots N,\ j=1,\ldots M$ in ein langer Vektor ${\vec {u}}_{h}\in {\mathbb {R} }^{NM}$ zeilenweise, siehe Assemblierung.

2. Implementieren Sie die entsprechende Systemmatrix $A_{h}\in {\mathbb {R} }^{NM\times NM}$ , beispielsweise indem Sie den entsprechenden Vektor auf die Diagonale und die Nebendiagonalen der Matrix setzen:

Ah=(1/{h^2})*diag(vector1,0)+diag(vector2,1)+ diag(vector3,-1)+diag(vector4,N)+diag(vector5,-N);

Bemerkung: Unterbrechen Sie danach noch Nebendiagonalen mit 0 auf entsprechenden Stellen, um die Blocktridiagonalität zu erzeugen.

3. Erweitern Sie die Matrix der Werte der Funktion der rechten Seite $f_{i,j}:=f(x_{i},y_{j}),\ i=1,\ldots N,\ j=1,\ldots M$ um die Beiträge aus den Dirichlet-Randbedingungen. Ordnen/assemblieren Sie $f_{i,j}$ in ein langer Spalten-Vektor ${\vec {f}}_{h}\in {\mathbb {R} }^{NM}$ (zeilenweise) mit dem Befehl 'reshape(Matrix,Anzahl_Zeilen,Anzahl_Spalten)'.

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

Bemerkung: 'reshape' ordnet eine Matrix spaltenmäßig um , d.h. zuerst Spalte 1, danach Spalte 2 von A usw. in ein langen Vektor. Probieren Sie aus:

 A=[1 2 3 ; 4 5 6; 7 8 9]
 reshape  (A,9,1)

Da Spaltenindexen $i=1,\ldots N,$ für x-Koordinaten stehen und in Assemblierung zeilenmäßig vorgegangen ist, muss hier die Matrix f der rechten Seite transponiert werden.

4. Lösen Sie das Gleichungssystem $-{c}A_{h}{\vec {u}}_{h}={\vec {f}}_{h}\in {\mathbb {R} }^{NM},\ c-$ konstanter Diffusinskoeffizient, ordnen Sie die Lösung wieder in eine Matrix um und stellen Sie die Lösung graphisch dar:

sol=Ah\fh;
   
 sol_matrix=reshape(sol,N,M);% Matrix mit N(Zeilen)x M(Spalten)
 sol_matrix=sol_matrix';
 ...

b)

Implementieren Sie in Octave die Finite-Differenzen-Methode (FDM) für Neumann-Randwertaufgabe auf dem Recheck D mit homogenen Randbedingungen mit $\alpha =\beta =\gamma =\delta =0$ . Benutzen Sie die Diskretisierung der Randbedingung von der ersten Ordnung.
Gehen Sie wie folg vor:

1. Erweitern Sie den äquidistantes Gitter um die Randpunkte, $\{(x_{i},y_{j}),\ i=0,\ldots N+1,\ j=0,\ldots M+1\}$ .

2. Erweitern Sie die Systemmatrix um neue Blockzeilen/Blockspalten, ${\tilde {A}}_{h}\in {\mathbb {R} }^{(N+2)(M+2)\times (N+2)(M+2)}$ , hierbei ist jeder Block dieser Matrix eine Matrix $(N+2)\times (N+2)$ . Hier können Sie wie im
Punkt a) 2. vorgehen, wenn sie nun N mit N+2 und M mit M+2 ersetzen, d.h. auf die Diagonale wird ein Vektor Länge $(N+2)\times (M+2)$ gesetzt werden usw. Anschließend, um die Approximation der Richtungsableitungen an den Rändern zu implementieren sollen die erste zwei Blöcke der ersten Blockzeile und die letzte zwei Blöcke der letzen Blockzeile ersetzt werden, wie auch die erste und letzte Zeile der einzelnen Blöcke etwa so:

% Neumann Randbedingung:
N=N+2;
M=M+2;
 ...
 
 % Neumann RB für y:   c* partial_y u=0 - einzelne Blöcke der Matrix ersetzen
 %----------------------------------------------------------------------
  block=diag((c/h)*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:  c* partial_x u=0 - einzelne Zeilen der Matrix ersetzen
  %----------------------------------------------------------------------
   for i=1:(M-2)
   vector_up=zeros(1,N*M);
   vector_up(N*i+1)=c/h;
   vector_up(N*i+2)=-c/h;
   vector_down=zeros(1,N*M);
   vector_down(N*i+N)=c/h;
   vector_down(N*i+N-1)=-c/h;
   BB(i*N+1,:)=vector_up ;
   BB(i*N+N,:)=vector_down ;
  endfor

Warum erhalten Sie kein gebrauchbares Ergebniss wenn Sie Neumann Randbedingungen an allen Rändern verwenden?

Um dies zu vermeiden verwenden Sie Dirichlet Bedingungen an zumindest einem der Rändern. Hierfür müssen Sie den Punktegitter nicht wieder verkleinern, sondern die Dirichlet-Randwerte gedanklich in die entsprechende äußere Punkte mit Koordinaten $x_{-1},\ x_{N+2},y_{-1},\ y_{M+2}$ übertragen.

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

Besprechung: SW 13

Hinweis

Für die Bearbeitung der Teilaufgabe a) finden Sie im Cocalc Kurs /Tutorium 5 oder im OLAT Materialordner SOSE 21/Tutorium/Abschlussprojekt den Jupyter Notebook Raumliches_Epidemiemodel_mit_RDGl.ipynb. Hier ist die Neumann-Randwertaufgabe für die zeitabhänhige Diffusionsgleichung mit FDM-Verfahren unter Anwendung des expliziten Eulerverfahren für die Zeitsdiskretisierung implementiert.

Teilaufgabe a)

Implemenieren Sie das FDM-Verfahren zur Lösung der Reaktionsdiffusionsgleichung für die zeitlich-räumliche Modellierung einer Epidemie auf einem Rechteckgebiet $D$ mit homogenen Neumann-Randbedingungen ${\frac {\partial u}{\partial {\vec {n}}}}|_{\partial D}=0$ .

Vewenden Sie dabei:

einen konstanten Diffusionskoffizient $c(x)=c$ ,
konstante aber realistische Bevölkerungsdichtefunktion $B(x)\equiv B$ , siehe Bevölkerunsgdichte Deutschland
die in Tutoriumaufgabe 3 anhand der RKI Daten festgestelte Infektionsraten (konstant und zeitabhängig-'slowdown').

Gehen Sie wie folg vor:
Definieren Sie eine Anfangslösung, die Sie als Startvektor verwenden und berechnen Sie die weitere Zeitschritte mit (beispielsweise) 'for' Schleife, ergänzen Sie den Reaktionsterm.

 sol_initial=reshape(initial',N*M,1);
 sol_old=sol_initial;
 
 %------------
 % Implementieren Sie die nichlineare Funktion des Reaktionstermes für beschränktes Wachstum der  kummulativen Infizierten 
 % (ggf. die Reaktionsterme fuer die aktuellen Infizierte u_I und  die anfaellige u_S)
 F=@(u) ....;
 %------------
 
 for i=1:Nr_timesteps
 sol= sol_old+ BB*sol_old*delta_t+ F(sol_old)*delta_t;
 
 sol_old=sol;
 %Spechern der aktuell  berechneten Lösungen als Vektor in eine Matrix, die am Ende Loesunen aus allen Schiweien enthaelt
 matrix_solution(:,i)=sol;
 endfor

Hinweis: Achten Sie darauf, das explizite Eulerverfahren hinreichend kleine Schrittweite $\Delta t$ benötigt.

Teilaufgabe b)

Inhomogene Populationsdichte in unserem Beispiel, angepasst von Monitor Einwohnerzahl 2011

Animation der Epidemieausbreitung mit geographisch differenzierter Populationsdichte, Anzahl der Gesammtinfizierten pro Rasterzelle (220mx220m), konstanter Diffusionskoeffizient c=0.1, Infektionsrate 0.3.

Ergänzen Sie Ihr Programm aus Teilaufgabe a) in dem Sie die Bevölkerunsgdiche im Reaktionsterm geographisch differenziert als Funktion von Raumvariable betrachten,

f(u(x,t))={\frac {c}{B(x)}}u(x,t){\big (}B(x)-u(x,t){\big )}

für Berechnung der Dichte der kummulierten Infizierten, oder

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

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

für Berechnung der Dichte der aktuellen Infizierten nach dem Kompartiment-Prinzip, siehe dazu die Diskretisierung des Reaktionstermes.

Verwenden Sie realistische Bevölkerungsverteilungsfunktion in einem ausgewählten Teil Deutschlands.

Fügen Sie der rechten Seite der Reaktionsdiffusionsgleichung, bzw. des Systems die notwendige Terme um die Wirkung der Impfungen abzubilden, diese wurden in

3. Tutoriumaufgabe- SIR Modell ausgearbeitet. Experimentieren Sie mit unterschiedlicher Impfraten.

Hinweise zu der Umwandlung eines graphischen Formats:

Speichern Sie ein Ausschnitt einer für Ihre Modellierung relevanten Grafik ab, importieren (mit 'imread') und bearbeiten Sie diese Graphik weiter in Octave:

%----------------------------------------------------
%Grafik einlesen und 
%Geeigneter Auschnitt wählen, Matrix A ist eine 3D Matrix
%----------------------------------------------------
 A=imread ('Monitor_Einwohnerzahl_detail.png');
 size(A) % Zeigt die Dimension der Bildmatrix
 S=A(20:24,20:27,2);
 
%Bilddarstellung 
%----------------------------------------------------
imshow(S)

%Inversion von unit8 zu double format
%----------------------------------------------------
 S=double (S);
 [s1 s2]=size(S);
 Max=max(max(S));
 
 
%Interpolation der Werte an die Koordinatenmatrix der numerischen Berechnung
 %----------------------------------------------------
 Sint=interp2(S,x,y, 'cubic');
 
 %y-Koordinate im Image läuft von oben nach unten, bei uns von unen nach oben, daher Bild flippen, oben <-> unten
 Sint=flipud (Sint)
  
 figure(Nr_timesteps+2);
 surface(x,y, (B/Max)*Sint)
 title ("Bevölkerungsdichte");
 ylabel("y")
 xlabel("x")
 colorbar
 
 %----------------------------------------------------
 %Skalierung auf einen gewünschten maxim. Wert B, 
 %Korrektur der Bevölkerungsdichte (keine Null-Einträge) und
 %Assemblierung der Matrix in ein Vektor
 %------------------------------------------------------
 B=reshape((B/Max)*Sint',N*M,1) .+0.1;
 size (B)

Hinweis: Achten Sie darauf, dass die gewählte Bevölkerungsdichtefunktion $B(x)>0$ streng positiv ist.

Teilaufgabe c)

Modellieren Sie den Reaktionsdiffusionsprozess der räumlichen Epidemieausbreitung mit nicht-konstantem Diffusionkoeffizienten $c(x)={\tilde {c}}(B(x))$ (Ausbreitungsgeschwindigkeit abhängig von der Populationsdichte). Verwenden Sie die entsprechende Diffusionsmatrix $A_{h}$ die für die Diskretisierung von $div(c(x,y)\nabla u(x,y))=\partial _{x}(c(x,y)\partial _{x}u(x,y))+\partial _{y}(c(x,y)\partial _{y}u(x,y))$ angepasst werden muss. Das Approximationsschema ergibt sich nach zweifacher Anwendung des zentralen Differenzenquotienten im obigen Ausdruck.