Kurs:Räumliche Modellbildung/Aufgaben Tutorium SS2020

Übungsaufgaben zur kontinuerlichen Modellierungsansätzen (AH)

Die Aufgaben sollen in 2er oder 3er Gruppen bearbeitet werden und im Teilnehmerordner von OLAT-Kurs Diffusionsprozesse hochgeladen werden.

1. Übungsaufgabe SW1: Kontaktmodell

Besprechung im SW3

Machen Sie sich vetraut mit dem Software Octave, gehen Sie dabei nach dem Skript Einführung in Octave Octave_intro.pdf und octave_first.m aus OLAT-Materialordner/Octave. Alternativ können Sie mit dem internativen Jupyter Notebook Einfuhrung in Octave.ipynb auf CoCalc arbeiten. Diese Datei finden Sie im OLAT-Materialordner/Tutoriumaufgaben.
Schreiben Sie selbs ein Octave Skript zum Einfachen Kontaktmodell für die räumliche Verbreitung der Infektion durch zufällige Bewegung der Individuen und Kontakt. Erstellen Sie Simulationen für verschiedene Ansteckungsradien, Bewegungsgeschwindigkeiten und Bewegungsrichtungen. Als Tutorium zur Bearbeitung dieser Aufgabe dient das interaktive Jupyter Notebook Kontakt Modell.ipynb und die dazugehörige statische Datei Kontakt_Modell_Ausarbeitung.pdf aus OLAT-Materialordner/Tutoriumaufgaben.

Folgende Befehle und Funktionen können nutzlich werden:

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) Untersuchen Sie die Fundamentallösungen $\Phi ,\Psi$ für die Laplace und für die instationäre Diffusionsgleichung. Stellen Sie diese graphisch als Funktionen auf $\mathbb {R} ^{2}$ bzw. $\mathbb {R} \times \mathbb {R} ^{+}$ in Octave dar. Achten Sie dabei darauf, dass Sie die Singularitätspunkte aus dem Definitionsbereich (aus dem Punktegitter) ausschneiden. Bei der Imlementierung in Octave können Sie verwenden:

[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 Koordintenmatrizen, 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 $\int \int z(x,y)dx,dy$ mithilfe der Trapezregel:

trapz(y,trapz(x,z,2))

Hinweis
Im OLAT-Materialordner/Octave/Tutoriumaufgabe 2 finden Sie ein Octave- Skript zur Implementierung der Poissonformel, alternativ können Sie mit dem Jupyter Notebook Losungsformel_Poisson.ipynb arbeiten, diesen finden Sie im OLAT-Materialordner/Tutoriumaufgaben. Schreiben Sie dieses Skript für die Lösungsformel $(\Psi *u_{0})(x,t)$ des Anfangswertproblems um.

3. Übungsaufgabe SW4: Kompartimentmodelle

Besprechung: SW6.

a) Implementieren Sie in Octave die numerische Lösung der Differentialgleichung des logistischen Wachstumsmodells (oder auch SI Modell bekannt) und des Systems von Differenzialgleichungen des SIR-Modells für die Infektionsverbreitung. Benutzen Sie realistische Anfangswerte (die vom Robert Koch Institut gemeldeten Daten am Anfang der Datenerfassung). Für die freie Kapazität (gleichzeitig den Anfangswert für die Gruppe S des SIR-Modells) nehmen Sie an, dass sich etwa 2/3 der Gesammtbefölkerung infiziert.

Experimentieren Sie mit verschiedenen Infektionsraten. Vergleichen Sie die Ergebnisse des logistischen (SI) und des SIR-Modells. (Achten Sie auf den Unterschied der Zahl der aktuell Infizierten und der kummulativen Infiziertenzahl.)

Numerische Lösung des Systems von Differenzialgleichungen in Octave mit 'lsode':

Parameter:

B=... %Kapazität: 2/3 der deutschen Befölkerung (in Tausend)
c=... %Infektionsrate
w=... ; % Genesungsrate, hier T-Genesungszeit ist 14 Tage
%Zeiten in den die Lösungs berechnet werden soll
times=(0:0.1:180);

%Anfangswerte als Spaltenvektor (zahlen in Tausend): SIR Modell:
yo=[...; ... ; 0];
%logist. Modell: yo=...;

Definiere die Funktion $f$ der rechten Seite des DGL-Systems $y'=f(y,t)$

%SIR Modell: Vektorwertige Funktion f
f=@(y,x) [-c*y(1)*y(2)/B;c*y(1)*y(2)/B-w*y(2);w*y(2)];
%logistisches Modell (skalarvertige Funktion f):
%f=@(y,x) c*y(1)*y(2)/B;

%numerische Lösung
y = lsode (f, yo, times);

Graphische Darstellung der Lösung des SIR Modells (drei Lösungsvektoren):

plot (times, y (:,1),'.',times, y (:,2), '.', times, y (:,3),'.',times, y (:,3)+y (:,2),'.')
legend ("Succeptible","Infected", "Removed", "cummulative  Infected ","location", "east")

b) Verwenden Sie das modifizierte SIR Modell zur Modellierung der Infektionsverbreitung in Deutschland unter der Annahme, dass die gemeldeten Infizierten-Daten nur ein Bruchteil der tatsächtlichen Anzahl der Infizierten erfassen (den Rest bilden die nicht-erkrankten Infizierten.) Kalibrieren Sie dieses Modell, sodass die zugehörige numerische Lösung die vom Robert Koch Institut gemeldeten Daten zur Corona-Infizierten approximert:

Als Anfangsbedingungen verwenden Sie die vom Robert Koch Institut gemeldeten Daten am Anfang der Datenerfassung.
Wählen Sie den Faktor r - der Datenerfassung.
Verwenden Sie eine Basisinfektionsrate $c\equiv k.N$ , die das exponenziellen Wachstum am Anfang der Pandemie modelliert.
Im Monat März wurden Kontaktreduzierende Maßnahmen wie Schulschließung und Kontaktverbort umgesetzt, die darauf folgende Verhaltensänderung hat sich auf die Infektionszahlen ausgewirkt. Modellieren Sie diesen Einfluß der Maßnahmen durch die Reduzierung/Anstieg der Infektionsrate in ihrem Octave-Skript, indem Sie die Infektionsrate als zeitabhängige Funktion definieren,

$c=c(t)$ .
Mithilfe des modifizierten SIR-Modells untersuchen Sie den Einnfluß der Kontaktreduzierenden Maßnahmen und der darauf folgenden Lockerung dieser Maßnahmen auf die Entwicklung der Pandemie. Welchen Einfluß hat Parameter r in dieser Modellierung?

Den Jupyter Notebook SIR_Modell.ipynb zu der Teilaufgabe b) und die Octave Funktion coronaData.m zu den RKI Daten finden Sie im OLAT-Materialordner/Tutoriumaufgaben.

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

Besprechung: SW 9

Als Grundlage für die Bearbeitung der folgenden Teilaufgaben a) und b) können Sie den Jupyter Notebook FDM_Poissongleichung_Dirichlet.ipynb verwenden, den Sie im OLAT-Materialordner/Tutoriumaufgaben finden. Hier ist die homogene Dirichlet-Randwertaufgabe für die Poissongleichung (für Null-Randwerte) implementiert.

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=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((a/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)=a/h;
   vector_up(N*i+2)=-a/h;
   vector_down=zeros(1,N*M);
   vector_down(N*i+N)=a/h;
   vector_down(N*i+N-1)=-a/h;
   BB(i*N+1,:)=vector_up ;
   BB(i*N+N,:)=vector_down ;
  endfor

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üsssen 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 SW10: Epidemiologische Modellierung als Reaktionsdiffusionsprozess mit FDM

Besprechung: SW 11

Für die Bearbeitung der Teilaufgabe a) finden Sie im OLAT-Materialordner/Tutoriumaufgaben den Octave Skript solve_diffusion.m Hier ist die Neumann-Randwertaufgabe für die zeitabhänhige Diffusionsgleichung mit FDM-Verfahren unter Anwendung des expliziten Eulerverfahren für die Zeitsdiskretisierung implementiert.

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

Gehen Sie wie folg vor:

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

Definieren Sie eine Anfangslösung, die Sie als Startvektor verwenden und berechnen Sie die weitere Zeitschritte mit (beispielsweise) 'for' Schleife

 sol_initial=reshape(initial',N*M,1);
 sol_old=sol_initial;
 
 %------------
 % Implementieren Sie die nichlineare Funktion des Reaktionstermes für beschränktes Wachstum der aktuellen oder kummulativen Infizierten
 F=@(u) ....;
 %------------
 
 for i=1:Nr_timesteps
 sol= sol_old+ BB*sol_old*delta_t+ F(sol_old)*delta_t;
 
 sol_old=sol;
 %Specheirn aller Lösungen in Vektorformat
 matrix_solution(:,i)=sol;
 endfor

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

b)
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. Speichern Sie dazu ein Ausschnitt einer für Ihre Modellierung relevanten Graphik, importieren (mit 'imread') und bearbeiten Sie diese Graphik weiter in Octave:

%----------------------------------------------------
%Graphik 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));
 ...
 S=fl
 
%Interpolation der Werte an die Koordinatenmatrix der numerischen Berechnung
 %----------------------------------------------------
 Sint=interp2(S,xend*(1*x.+2)/(xend+2),yend*(1*y.+2)/(yend+2), 'cubic');
 %% alternativ 
 %%Sint=interp2(S,x,y, 'cubic');
 
 %y-Koordinate im Image läuft von oben nach unten, 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 max. 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$ stern positiv ist.

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.