Kurs:Räumliche Modellbildung/Gruppe Nr
Gruppenseite - (GBN)
BearbeitenDiese Seite enthält die gesammelten Materialien der Gruppe 17 für die Portfolioprüfung.
Teilnehmer-innen
Bearbeiten- Sophia Gazelkowski
- Marie Baudy
- Anna-Lena Neumeyer
Diskrete Modellierung
BearbeitenZiel unserer diskreten Modellierung ist es, das Infektionsgeschehen auf kleinskaligem Niveau (beispielsweise innerhalb eines Klassen- oder Büroraumes) mit Hilfe des SIR-Modells abzubilden. Dabei sollen Transportprozess und epidemiologischer Prozess abwechselnd ausgeführt werden.
Implementierung in Excel
BearbeitenDie Implementierung in Excel erfolgte nach der Ausarbeitung der Ausgangsidee über zwei Zyklen.
Unsere Ausgangsidee: In einem 10x10 Gebiet bewegen sich drei Personen. Dabei wird die x-Koordinate als zeitlicher Verlauf angesehen und nur mit der y-Koordinate wird die räumliche Ausbreitung bestimmt. Daher wird die Ausgangsidee für eine räumliche Ausbreitung sowohl in x- als auch in y-Richtung im 1. Zyklus überarbeitet.
In unserem ersten Modellierungs-Zyklus gehen wir vereinfacht von 5 Personen innerhalb eines 20x20-Gebietes aus. Jeder dieser Personen wird pro Zeitschritt mit dem Befehl =ZUFALLSBEREICH() eine zufällige x- und y-Koordinate zugeordnet. Über den Infektionsstatus kann festgelegt werden, ob die Personen zu Beginn, im Zeitpunkt t=0, infiziert (1) oder anfällig (0) sind. Zu Beginn setzen wir nur eine Person als infiziert fest. In jedem Zeitschritt kommt es zu einer Positionsänderung der Personen und einer anschließenden Überprüfung, ob eine Infizierung stattfinden kann.
Eine Infizierung findet unter folgenden Bedingungen statt:
1) Eine der beiden Personen muss infiziert sein.
2) Kritischer Abstand von 0.2 unterschritten
3) Ansteckungswahrscheinlichkeit erfüllt
Im zweiten Modellierungszyklus wird der unnatürliche, sprunghafte Bewegungsablauf zu einem natürlicheren, kontinuierlichen Bewegungsablauf optimiert. Die Position zu dem Zeitpunkt t=0 wird weiterhin über den Zufallsbereich bestimmt. Für jeden darauffolgenden Zeitpunkt wird über eine Zufallszahl bestimmt, ob sich die Koordinaten der jeweiligen Person um einen Wert von 1 erhöhen, verringern oder gleich bleiben. Läuft eine Person am Rand des Gebiets heraus, läuft sie auf der anderen Seite wieder herein. Um dem SIR-Modell zu genügen, wird eingeführt, dass Infizierte nach einer Genesungsdauer von 14 Tagen automatisch in die Gruppe der Genesenen wechseln. Dabei wird ihnen ein Infektionsstatus von 2 zugeordnet.
Implementieren der Personengruppen S,I und R
BearbeitenZu Beginn wird zunächst eine bestimmte Anzahl von Personen festgelegt, die sich insgesamt in dem Raum befinden sollen. Dabei geht man von 6 Personen gemäß dem SIR-Modell aus. Es wird für jede einzelne Gruppe eine Matrix erstellt, in der die x- und y- Koordinate der Personen gespeichert werden.
I_alt=[5,10;2,10;1,5;-1,-1;-1,-1;-1,-1] %Matrix der Infizierten
I_neu=[5,10;2,10;1,5;-1,-1;-1,-1;-1,-1]
S_alt=[-1,-1;-1,-1;-1,-1;1,5;7,8;1,4] %Matrix der Anfälligen
S_neu=[-1,-1;-1,-1;-1,-1;1,5;7,8;1,4]
R_alt=[-1,-1;-1,-1;-1,-1;-1,-1;-1,-1;-1,-1] %Matrix der Genesenen
R_neu=[-1,-1;-1,-1;-1,-1;-1,-1;-1,-1;-1,-1]
Transport der einzelnen Personengruppen
BearbeitenBeim Transportprozess gibt es grundsätzlich zwei Möglichkeiten, damit die im Raum befindlichen Personen nicht aus dem Gitter verschwinden können:
- Die Person kommt am anderen Rand wieder in das Gitter hinein. (zyklische Randbedingung)
- Die Person kann garnicht erst über den Rand hinauslaufen, sondern behält ihre Randposition, wenn die neue Position größer wäre als die Länge oder Breite des Gitters.
In der vorliegenden Modellierung haben wir uns dazu entschieden, die Personen nicht über den Rand hinauslaufen zu lassen. Da wir uns beispielhaft in einem Klassenraum befinden, spiegelt diese Variante das Unterrichtsgeschehen besser wider. Stellt man sich den Klassenraum als unser gewähltes Gitter vor, erscheint es unrealistisch, dass sich Personen auf der einen Seite heraus bewegen und dann wieder auf der anderen Seite hereinkommen. Die Schülerinnen und Schüler verweilen in unserem Modell vereinfacht im vorgegebenen Gitter.
function P_neu=transport(P_alt) %Funktionsname: transport, Einsetzen von P_alt, um P_neu zu erhalten
Breite=10
Hoehe=20 %des Gitters
uG=-1/3; %untere Grenze für die Zufallszahl
oG=1/3; %obere Grenze für die Zufallszahl
Gesamtzahl=6; %6 Personen im Gitter
P_neu=P_alt;
for i=1:Gesamtzahl
Zufall1=(rand(1)*2-1); %Zufallszahl zwischen -1 und 1 für jedes i, wirkt auf die erste Komponente (x-Koordinate)
Zufall2=(rand(1)*2-1); %Zufallszahl zwischen -1 und 1 für jedes i, wirkt auf die zweite Komponente (y-Koordinate)
if (P_alt(i,1)>-1 && P_alt(i,2)>-1) %Person i darf bewegt werden, beide Einträge in Zeile i größer als -1, d.h. Person fällt in die entsprechende Gruppe des SIR-Modells
%Zufall in erster Komponente
if (Zufall1<uG)
P_neu(i,1)=P_alt(i,1)-1
if (P_neu(i,1)<0)
P_neu(i,1)=P_alt(i,1)
endif
elseif (Zufall1>oG)
P_neu(i,1)=P_alt(i,1)+1; %x-Koordinate wird um 1 größer = Bewegung nach rechts
if (P_neu(i,1)>Breite)
P_neu(i,1)=P_alt(i,1) %kann nicht über den Rand hinauslaufen, z.B. Klassenraum --> Person bleibt am Rand stehen, wenn die neue Position größer wäre als die Breite des Gitters
endif
else
P_neu(i,1)=P_alt(i,1); %Person bleibt zeilenmäßig vor Ort
endif
%Zufall in zweiter Komponente
if (Zufall2<uG)
P_neu(i,2)=P_alt(i,2)-1; %y-Koordinate wird um 1 kleiner = Bewegung nach unten
if (P_neu(i,2)<0)
P_neu(i,2)=P_alt(i,2)
endif
elseif (Zufall2>oG)
P_neu(i,2)=P_alt(i,2)+1; %y-Koordinate wird um 1 größer = Bewegung nach oben
if (P_neu(i,1)>Hoehe)
P_neu(i,1)=P_alt(i,1) %kann nicht über den Rand hinauslaufen, z.B. Klassenraum --> Person bleibt am Rand stehen, wenn die neue Position größer wäre als die Breite des Gitters
endif
else
P_neu(i,2)=P_alt(i,2); %Person bleibt zeilenmäßig vor Ort
endif
endif
endfor
endfunction
Überschreiben der alten Positionen:
I_neu=transport(I_alt)
S_neu=transport(S_alt)
R_neu=transport(R_alt)
S_alt=S_neu
I_alt=I_neu
R_alt=R_neu
Epidemiologischer Prozess
Bearbeiten- Ansteckung I: Anfällige werden infiziert
- Ansteckung R: Infizierte werden genesen
- Ansteckung S: Anfällige wechseln durch Impfung in die Gruppe der Genesenen
Ansteckung I: Infektion, Wechseln von S zu I
BearbeitenZunächst wird der Abstand zwischen den Infizierten und anfälligen Personen berechnet. Im nächsten Schritt wird die Ansteckung beschrieben, bei der eine infizierte Person eine anfällige Person ansteckt
function d=abstand_Personen( Person_I, Person_S)
d=sqrt( (Person_S(1)-Person_I(1))^2 + (Person_S(2)-Person_I(2))^2 ) %euklidische Norm
endfunction
function SI_neu=AnsteckungI(S_alt,I_alt,R_alt) %Funktionsname: AnsteckungI, Einsetzen der alten Infizierte, Anfälligen und Genesenen um die neuen Infizierten zu erhalten
d=4 %Abstand
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=1:Gesamtzahl
if (I_alt(j,1)>-1 && I_alt(j,2)>-1) %1. Komponente ist x-Koordinate, 2. Komponente ist y-Koordinate der Infizierten
%hier ist klar, dass die i-te Person infektiös ist
Person_I = [I_alt(j,1),I_alt(j,2)]
if (S_alt(s,1)>-1 && S_alt(s,2)>-1) %Index an die Gruppe angepasst
%hier ist klar, dass die Person anfällig ist
Person_S = [S_alt(s,1),S_alt(s,2)] %euklidische Norm ausrechnen zwischen S und I
Zufall=rand(1)
if ( (abstand_Personen( Person_I, Person_S)<d) && Zufall>0.7)
%s-te Person verschwindet aus S
S_neu(s,1)=-1
S_neu(s,2)=-1
%s-te Person wird infektiös
I_alt(s,1)=S_alt(s,1)
I_alt(s,2)=S_alt(s,2)
else I_neu=I_alt; S_neu=S_alt;
endif
endif
endif
endfor
endfor
SI_neu=[S_neu, I_neu]
endfunction
SI_neu=AnsteckungI(S_alt,I_alt,R_alt)
Da die Infizierten und Anfälligen in einer großen Matrix (6x4) gespeichert sind, müssen diese im folgenden wieder in 2 getrennte Matrizen der Dimension (6x2) getrennt werden, dies geschieht mit der Extrahier-Funktion
Separierung der Personengruppe I:
%Personengruppen müssen aus der Matrix separiert werden
function I_neu=SI_extrahiere_I(SI_neu)
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=3:4
I_neu(j,s-2)=SI_neu(j,s)
endfor
endfor
endfunction
I_neu=SI_extrahiere_I(SI_neu)
I_alt=I_neu
Separierung der Personengruppe S:
%Personengruppen müssen aus der Matrix separiert werden
function S_neu=SI_extrahiere_S(SI_neu)
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=1:2
S_neu(j,s)=SI_neu(j,s)
endfor
endfor
endfunction
S_neu=SI_extrahiere_S(SI_neu)
S_alt=S_neu
Ansteckung R: Genesung, Wechsel von I zur R
BearbeitenHier wird eine bestimmte Wahrscheinlichkeit angenommen, mit der ein Infizierte in die Gruppe der Genesenen wechselt.
function IR_neu=AnsteckungR(S_alt,I_alt,R_alt) %Funktionsname: AnsteckungI, Einsetzen der alten Infizierte, Anfälligen und Genesenen um die neuen Infizierten zu erhalten
G=0.86
Gesamtzahl=6;
for j=1:Gesamtzahl
Zufall1=rand(1)
if (Zufall1>G && (I_alt(j,1)>-1 && I_alt(j,2)>-1))
%j-te Person verschwindet aus I
I_neu(j,1)=-1
I_neu(j,2)=-1
%j-te Person wird recovered
R_neu(j,1)=I_alt(j,1)
R_neu(j,2)=I_alt(j,2)
else
I_neu=I_alt
R_neu=R_alt
endif
endfor
%wir geben hier ein 2-dimensonales Vektorelement für S und R zurückgegeben (4x6 Matrix)
IR_neu=[I_neu,R_neu]
endfunction
IR_neu=AnsteckungR(S_alt,I_alt,R_alt)
Separierung der Personengruppe R:
%Personengruppen müssen aus der Matrix separiert werden
function R_neu=IR_extrahiere_R(IR_neu)
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=3:4
R_neu(j,s-2)=IR_neu(j,s)
endfor
endfor
endfunction
R_neu=IR_extrahiere_R(IR_neu)
R_alt=R_neu
Separierung der Personengruppe I:
%Personengruppen müssen aus der Matrix separiert werden
function I_neu=IR_extrahiere_I(IR_neu)
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=1:2
I_neu(j,s)=IR_neu(j,s)
endfor
endfor
endfunction
I_neu=IR_extrahiere_I(IR_neu)
I_alt=I_neu
Ansteckung S: Impfung, Wechsel von S zu R
BearbeitenHier wechseln Anfällige in die Gruppe der Genesenen mit einer bestimmten Impfrate, die abhängig vom Tag ist.
function SR_neu=AnsteckungS(S_alt,I_alt,R_alt, t ) %Impfung
%Zeitschritte=5
Impfrate=0.6 %Personen pro Zeitschritt
Impfung=floor(Impfrate*t)
%m=0:1:t
Gesamtzahl=6;
m=0
for s=1:Gesamtzahl
if (S_alt(s,1)>-1 && S_alt(s,2)>-1)
%s-te Person verschwindet aus S
S_neu(s,1)=-1
S_neu(s,2)=-1
%s-te Person wird recovered
R_neu(s,1)=S_alt(s,1)
R_neu(s,2)=S_alt(s,2)
m=m+1
else
S_neu=S_alt
R_neu=R_alt
endif
if m==Impfung
break
endif
endfor
%hier wird ein 2-dimensonales Vektorelement für S und R zurückgegeben (4x6 Matrix)
SR_neu=[S_neu,R_neu] %S und R sollen zurückgegeben werden
endfunction
SR_neu=AnsteckungS(S_alt,I_alt,R_alt,4)
Separierung der Personengruppe R:
function R_neu=SR_extrahiere_R(SR_neu) %extrahiere R aus SR_neu
Gesamtzahl=6;
for j=1:Gesamtzahl %Zeilen
for s=3:4 %Spalten
R_neu(j,s-2)=SR_neu(j,s) %die erste und zweite Spalte von R_neu soll der 3. und vierten Spalte von SR_neu entsprechen
endfor
endfor
endfunction
R_neu=SR_extrahiere_R(SR_neu)
R_alt=R_neu
Separierung der Personengruppe S:
%Personengruppen müssen aus der Matrix separiert werden
function S_neu=SR_extrahiere_S(SR_neu)
Gesamtzahl=6;
for j=1:Gesamtzahl
for s=1:2
S_neu(j,s)=SR_neu(j,s)
endfor
endfor
endfunction
S_neu=SR_extrahiere_S(SR_neu)
S_alt=S_neu
Alternative mit Zufallszahl für die Impfrate
function S_neu=AnsteckungS(S_alt,I_alt,R_alt) %Funktionsname: AnsteckungI, Einsetzen der alten Infizierte, Anfälligen und Genesenen um die neuen Infizierten zu erhalten
G=0.5
Gesamtzahl=6;
for s=1:Gesamtzahl
Zufall1=rand(1)
if (Zufall1>G && (S_alt(s,1)>-1 && S_alt(s,2)>-1))
%s-te Person verschwindet aus S
S_neu(s,1)=-1
S_neu(s,2)=-1
%s-te Person wird recovered
R_neu(s,1)=S_alt(s,1)
R_neu(s,2)=S_alt(s,2)
else
S_neu=S_alt
R_neu=R_alt
endif
endfor
endfunction
Kontinuierliche Modellierung
BearbeitenTutorium 1 - Kontaktmodell
BearbeitenEinführung zum Kontaktmodell
Bearbeiten
Mögliche Szenarien
BearbeitenSzenarien:
BearbeitenSzenario 1: Hohe Anzahl der Infizierten: Bei hohem Ansteckungsradius (z.B. Britische Mutante bei 2m) und gleichzeitig hohe Mobilität der Bevölkerung (Geschwindigkeit und Bewegungsrichtung) steigen die Infektionszahlen im zeitlichen Verlauf am stärksten.
Szenario 2: Niedrige Anzahl der Infizierten: Bei niedrigem Ansteckungsradius (z.B. herbeigeführt durch das Tragen medizinischer Masken) und geringer Mobilität der Bevölkerung (durch Ausgangsbeschränkungen und Kontaktminimierung) bleibt die Anzahl der Infizierten im zeitlichen Verlauf auf niedrigem Niveau
Vergleich Szenario 1 und Szenario 2
BearbeitenDas Szenario 1 soll ein Infektionsgeschehen mit hohen Infektionszahlen simulieren. Dafür wurde der Ansteckungsradius auf 1.5 erhöht, sowie Bewegungsgeschwindigkeit moderat auf 2. Die Position des Ausgangsinfizierten wird mittig gesetzt und die einzelnen Personen laufen pro Durchlauf der Zeitschleife in zufällige verschiedene Richtungen. Das Szenario 2 soll-gegensätzlich zum ersten Szenario- ein Infektionsgeschehen mit möglichst geringer Anzahl an Infizierten simulieren. Daher wird der Ansteckungsradius auf 0.3 verringert, sowie die Bewegungsgeschwindigkeit auf 0.4. Die Position des Ausgangsinfizierten wird an den Rand gesetzt. Die einzelnen Personen laufen pro Durchlauf der Zeitschliefe in dieselbe Richtung, da der Befehl
neuPosX=x.+g.*rand(Ny,Nx)-0.5*g;
neuPosY=y.+g*rand(Ny,Nx)-0.5*g;
vor die Zeitschleife gesetzt wurde, so dass die Schleife wiederholt auf dieselbe Bewegungsrichtung zugreifen kann. In den beiden GIFs kann man das Infektionsgeschehen der beiden Szenarien gut miteinander vergleichen.
Szenario 3: Vergleich unterschiedlicher Positionen des Patient-0 (Erstinfizierte*r)
Befindet sich die Infizierte Person am Rand der Menschenmenge, so werden unter Umständen bei der Bewegung der Personen weniger anliegende Personen infiziert. Befindet sich die Person jedoch im Zentrum der Menschenmenge, so werden mit großer Wahrscheinlichkeit eine höhere Anzahl von Personen infiziert.
Ergebnisse
BearbeitenDurch mehrfache Änderung der Parameter sind wir auf folgende Ergebnisse gekommen:
Ansteckungsradien: Je größer der Ansteckungsra
dius gewählt wird, desto schneller steigt die Anzahl der Infizierten, da sich die Personen nicht zu nahe kommen müssen, um sich zu infizieren. Wenn der Infektionsradius größer als 1 gewählt wird, werden im nächsten Schritt alle Nachbarn der Infizierten angesteckt, da die Personen einen Abstand von 1 zueinander haben. Die Infektion verbreitet sich in diesem Fall ohne Bewegung, wobei das Einbauen von Bewegung sogar zu einer Verlangsamung der Infektionsentwicklung führen könnte.
Bewegungsgeschwindigkeit: Je höher die Bewegungsgeschwindigkeit, desto schneller entsteht eine hohe Anzahl an Infizierten. Allerdings sind die Personen nach gewisser Zeit so weit voneinander entfernt, dass die Infektionsentwicklung plateauförmig verläuft. Wählt man die Bewegungsgeschwindigkeit zu hoch, so tritt die Plateauwirkung schon früher bei noch geringerer Infiziertenanzahl ein.
Zeiteinheit: Eine Erhöhung der Zeiteinheit T bedeutet, dass mehr Infizierte entstehen, da die zufällige Bewegung über einen längeren Zeitraum (mehrere Stunden oder Tage) läuft. Erhöhung der Zeitschritte liefert exaktere Zwischenergebnisse, da innerhalb einer Zeiteinheit T, häufiger die Anzahl der Infizierten abgerufen wird.
Änderung der Bewegungsrichtung: Je öfter die zufällige Bewegungsrichtung einer Person geändert wird, desto mehr Infizierte können unter Umständen entstehen. Jedoch ist die Änderung der Bewegungsrichtung nicht der alleinige Faktor für hohe Infektionszahlen. Es bedarf einer Kombination aus zentraler Position, mittlerer Bewegungsgeschwindigkeit und mittlerem bis großem Infektionsradius.
Eine Veränderung der Ausgangsposition des Infizierten beeinflusst die weitere Infektionsausbreitung. Je zentraler sich der Infizierte ausgangs in der Menschenmenge befindet, desto mehr Menschen kann er innerhalb der nächsten Zeitschritte anstecken. Befindet er sich dagegen am Rand der Menschenmenge, kommt er in der gleichen Zeit mit deutlich weniger Personen in Kontakt.
Anmerkung: Durch den "rand"-Befehl werden bei jeder Ausführung des Algorithmus neue Zufallszahlen generiert, sodass sich die Ergebnisse der Ausführungen mit gleichen Parametern unterscheiden.
Ausgangscode
BearbeitenAnmerkung: Dieser Code wurde uns über CoCalc von Frau Prof. Dr. Hundertmark zur Verfügung gestellt.
%Anzahl von Punkten in x und y Richtung
Nx=10;
Ny=20;
Zeitschritte=3; %pro eine Zeiteinheit T
radiusInf=0.80;
T=6;%Zeiteinheiten
g=0.50;%Bewegungsgeschwindigkeit
%------------------------
x=[1:Nx];
y=[1:Ny];
dt=1/Zeitschritte;
[x,y]=meshgrid (x,y);
%erste Infizierte
xInf=x(1, Nx/2);
yInf=y(Ny/2);
IndInf2=Nx/2; %Spalte
IndInf1=Ny/2; %Zeile
IndInf1neu=IndInf1;
IndInf2neu=IndInf2;
close all;
figure (1)
hold on;
plot(x,y, '*b') ;
plot(xInf,yInf, 'sr') ;
axis([-5 15 -5 25]);
test=["Fig_", num2str(1),".jpg"]
saveas(1, test)
hold off;
neuPosX=x.+g.*rand(Ny,Nx)-0.5*g;
neuPosY=y.+g*rand(Ny,Nx)-0.5*g;
Anzahl=0; % Zur zeitlichen Darstellung der Infiziertenzahl
%Zeitschleife für 3 Zeitschritte pro eine Zeiteinheit T:
%---------------------------
for i=1:T*Zeitschritte
neuX=x.+1*(neuPosX.-x)*i*dt;
neuY=y.+1*(neuPosY.-y)*i*dt;
figure (i+1)
title (['t=' num2str(i*dt)])
axis([-5 15 -5 25]);
hold on
plot(neuX,neuY, '*b') ;
%neue Infizierungen:
%In Vektoren IndInf1, IndInf2 sind die i- und j-Indexen der (x,y)-Koordinaten der bereits Infizierten gespeichert.
%---------------------------
zahler=1;
for k1=1:length(IndInf1)
for j=1:Nx
for l=1:Ny
%die aktuellen Positionen der Füß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-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)
Anzahl= [Anzahl lengt];
endfor
%Zeitschleife abgeschlossen
%---------------------------
hold off;
% ---------------------------
Infiz1=[IndInf1; IndInf2] %Ausabe Inizierten Indexen
hold off;
Anzahl
plot((1:T*Zeitschritte)*dt, Anzahl(2:length(Anzahl)), 'sr-')
xlabel('Zeit')
ylabel('Anzahl von Infizierten')
grid on
Änderung der Bewegungsrichtung
Bearbeiten%Anzahl von Punkten in x und y Richtung
Nx=10;
Ny=20;
g=1; %Geschwindigkeit mit der sich die Personen bewegen
radiusInf=0.80;
T=6;
Zeitschritte=3;
dt=1/Zeitschritte;
x=[1:Nx];
y=[1:Ny];
[x,y]=meshgrid (x,y);
%erste infizierte Person
xInf=x(1,Nx/2); %Position des 1. Infizierten mit konkreten Koordinaten
yInf=y(Ny/2);
IndInf1=Ny/2;
IndInf2=Nx/2;
IndInf1neu=IndInf1;
IndInf2neu=IndInf2;
figure (1)
hold on;
plot(x,y, '*b');
plot(xInf, yInf, 'sr');
axis([-5 15 -5 25]);
%Speichern der Bilder
test=["Fig_", num2str(1),".jpg"]
saveas(1, test)
hold off;
Anzahl=0; % Zur zeitlichen Darstellung der Infiziertenzahl
%Zeitschleife für 3 Zeitschritte pro eine Zeiteinheit T: 18 Durchläufe
for i=1:T*Zeitschritte
%Bestimmen der neuen Position pro Durchlauf
neuPosX=x.+1*g.*(rand(Ny,Nx)-0.5); %Zufallszahl zwischen -0,5 und 0,5
neuPosY=y.+1*g.*(rand(Ny,Nx)-0.5);
neuX=x.+(neuPosX.-x)*1*dt;
neuY=y.+(neuPosY.-y)*1*dt;
x=neuX;
y=neuY;
figure (i+1)
title (['t=' num2str(i*dt)])
axis([-5 15 -5 25]);
hold on
plot(neuX, neuY, '*b') ;
quiver(neuX, neuY, neuPosX.-x, neuPosY.-y) %Richtungspfeile
%neue Infizierungen:
%In Vektoren IndInf1, IndInf2 sind die i- und j-Indexen der (x,y)-Koordinaten der bereits Infizierten gespeichert.
for k1=1:length(IndInf1) %Länge von Zeile (Ny/2) im ersten Durchlauf
for j=1:Nx
for l=1:Ny
%Abgleichen der aktuellen Positionen der Personen mit denen der Infizierten durch Berechnung des Abstands
abstand=norm ( [neuX(l,j)-neuX(IndInf1(k1),IndInf2(k1)), neuY(l,j)-neuY(IndInf1(k1),IndInf2(k1))]);
if abstand<radiusInf %es kommt zu einer Infizierung am l,j-ten Punkt (rotes Quadrat)
plot(neuX(l,j),neuY(l,j),'sr') ;
k2=1;
%Untersuchung des Vektors der Infizierten: l,j-Indexen der aktuellen Personen sollen nicht unter den infizierten Indexen auftauchen
while k2 <=length(IndInf1)&&(l!=IndInf1neu(k2) || j!=IndInf2neu(k2))
k2=k2+1;
endwhile
%Speicherung der l,j-Indexen der neuen infizierten Person im neuen Vektor der Infizierten-Indexen zur Speicheroptimierung
if k2 >length(IndInf1)
IndInf1neu=[IndInf1neu l];
IndInf2neu=[IndInf2neu j];
endif
endif
IndInf1=IndInf1neu;
IndInf2=IndInf2neu;
endfor
endfor
endfor
hold off
lengt=length(IndInf1)
Anzahl= [Anzahl lengt];
%Speichern der Bilder
test=["Fig_", num2str(i+1),".jpg"]
saveas(i+1, test)
endfor
Infiz1=[IndInf1; IndInf2] %Ausgabe der Infizierten-Indexen
hold off;
Tutorium 2 - Fundamentallösungen
BearbeitenIm ersten Teil der Aufgabe geht es darum, die Lösungsformel für die LaPlace-Gleichung und die Diffusionsgleichung zu implementieren.
Fundamentallösung für die LaPlace-Gleichung in zweidimensionalen Fall
LaPlace-Gleichung: Die Fundamentallösung der LaPlace-Gleichung lautet folgendermaßen: Die || . || beschreibt die euklidische Norm und wn das Volumen der Einheitskugel, welches sich jedoch in Abhängigkeit der Dimension verändert. In der obigen Formel ist das Volumen in dreidimensionalen Fall gegeben.
clear all
phi=@(x,y) -log(sqrt(x.^2+y.^2))/(2*pi); %Lösungsformel für die Fundamentallösung der LaPlace-Gleichung (vgl. Skript)
step=0.25;
X=[-5:step:5];
Y=[-5:step:5];
[x,y]=meshgrid(X,Y);
figure 3
meshc(x ,y, phi(x,y))
axis([-5 5 -5 5 -0.4 0.2])
title('Fundamentalloesung Laplace')
test=["Fig_.jpg"]
saveas(3, test)
Lösungsformel für die Diffusionsgleichung
Für die Diffusionsgleichung gibt es eine Lösungsformel, welche die Gleichung löst. Jedoch muss als Voraussetzung gelten, dass der Diffusionskoeffizient konstant ist.
Die Fundamentallösung für die Gleichung lautet:
clear all
%Diffusionskoeffizient c(x)=a= 1 fest, kann jedoch beliebig variert werden
% t= Zeitschritte welche durchlaufen werden
t0=0;
dt=0.1;
for k=1:10
t=t0+k*dt %t>0, d.h. kritische Punkte kommen nicht vor, folglich ist die Singularität ausgeschlossen
%(Wenn beide Paramter gegen Null streben, so geht psi gegen unendlich)
psi=@(x,y,t) 1/(4*pi*t)*exp(-(x.^2+y.^2)/(4*t));
step=0.1;
X=[-4:step:4];
Y=[-4:step:4];
[x,y]=meshgrid(X,Y);
figure (k)
surf(x,y,psi(x,y,t))
axis([-3 3 -3 3 0 0.8])
grid on
title (['t=' num2str(t)])
xlabel('x')
xlabel ('y')
test=["Fig_",num2str(t),".jpg"]
saveas(k, test)
endfor
Implementierung der Poissonformel
Die Poissongleichug ist eine inhomogene LaPlace-Gleichung:
Die Lösung dieser Gleichung ergibt sich aus einer Faltung mit der Fundamentallösung der LaPlace-Gleichung, welche oben dargestellt wurde, und der Funktion der rechten Seite.
Zunächst wird die rechte Seite der Funktion definiert:
function wert=Quellfunktion(x,y)
if abs(x.-1)+abs(y.-1)<=0.5 wert=1; %Betragsnorm
else wert=0;
endif
endfunction
%Definition des Gebietes
step=0.03;
X=[0:step:2];
Y=[0:step:2];
[x,y]=meshgrid(X,Y);
%Die Funktionswerte werden in jedem Gitterpunkt in einer Matrix gespeichert
for i=1:(length(X))
for j=1:(length(Y))
f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
endfor
endfor
for i=1:(length(X))
for j=1:(length(Y))
xstar=x(j,i)*ones(length(Y),length(X));
ystar=y(j,i)*ones(length(Y),length(X));
% Phi mit verschobenem Argument:
phi=-log(sqrt((xstar.-x).^2+(ystar.-y).^2))/(2*pi);
phi(j,i)=0;
%numerische Integration
I(j,i) = trapz(Y,trapz(X,phi.*f,2));
endfor
endfor
figure 1
meshc(x,y,I)
grid on
title ('Losung von Poissongleichung')
xlabel('x')
xlabel ('y')
test=["Fig_1_2",".jpg"]
saveas(2, test)
figure 2
surfc(x,y,f)
title ('Funktion der rechten Seite')
test=["Fig_2_2",".jpg"]
saveas(3, test)
Änderung der rechten Seite:
clear all;
function wert=Quellfunktion(x,y)
if nthroot((x).^3+(y).^3,3)<=0.5 wert=1; %p-Norm
else wert=0;
endif
endfunction
Lösung der Poissiongleichung für diese rechte Seite:
Implementierung der zeitabhängigen Diffusionsgleichung
Die Lösung dieser Gleichung ergibt sich aus einer Faltung mit der Fundamentallösung der Diffusionsgleichung, welche oben dargestellt wurde, und der Anfangsfunktion.
Zunächst wird die Anfangsfunktion definiert:
clear all;
function wert=Quellfunktion(x,y)
if sqrt((x.-1.0).^2+(y.-1).^2)<=0.1 wert=1; %euklidische Norm
elseif sqrt((x.-0.8).^2+(y.-1.8).^2)<=0.2 wert=1;
else wert=0;
endif
endfunction
a=1;
n=2;
step=0.02;
X=[0:step:2];
Y=[0:step:2];
[x,y]=meshgrid(X,Y);
for i=1:(length(X))
for j=1:(length(Y))
f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
endfor
endfor
t0=0;
dt=0.001;
for k=1:20
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 mit verschobenem Argument:
psi=(1/((4*a*pi*t)^(n/2)))*exp(-((xstar.-x).^2+(ystar.-y).^2)/(4*t*a));
%psi(j,i)=0; %hiermit wird ananlog zu LaPlace und Diffusionsgleichug die Sigularität ausgeschlossen
%numerische Integration
I(j,i) = trapz(Y,trapz(X,psi.*f,2));
endfor
endfor
figure (k)
surf(x,y,I)
hold on
axis([0 2 0 2 0 1]);%Normierung der Achsen
grid on
title (['t=' num2str(t)])
xlabel('x')
xlabel ('y')
hold off
test=["Fig_", num2str(t),".jpg"]
saveas(k, test)
endfor
Tutorium 3 - SIR-Modell
BearbeitenAufgabe 1: Beginn der Pandemie
x=1:67;
x=1:30;
y=[ 16.0 18.0 21.0 26.0 53.0 66.0 117.0 150.0 188.0 240.0 400.0 639.0 795.0 902.0 1139.0 1296.0 1567.0 2369.0 3062.0 3795.0 4838.0 6012.0 7156.0 8198.0 10999.0 13957.0 16662.0 18610.0 22672.0 27436.0 31554.0 36508.0 42288.0 48582.0 52547.0 57298.0 61913.0 67366.0 73522.0 79696.0 85778.0 91714.0 95391.0 99225.0 103228.0 108202.0 113525.0 117658.0 120479.0 123016.0 125098.0 127584.0 130450.0 133830.0 137439.0 139897.0 141672.0 143457.0 145694.0 148046.0 150383.0 152438.0 154175.0 155193.0 156337.0 157641.0 159119.0];
y=y(1:30);
f = @(b,x) b(1).*exp(b(2).*x); % Objective Function
B = fminsearch(@(b) norm(y - f(b,x)),[16;0.1]) % geschätzte Parameter; 16 wird durch die Ausgangslage vorgegeben
figure
plot(x, y, 's')
hold on
plot(x, f(B,x))
hold off
grid
xlabel('x')
ylabel('f(x)')
text(50, 150, sprintf('f(x) = %.1f\\cdote^{%.3f\\cdotx}%+.1f', B))
B(1)
Aufgabe 2: Slowdown
NM=55000 % Kapazitätsgrenze, 2/3 der deutschen Bevölkerung in 1000 umgerechnet
c=0.32450
w=1/13.90; % Wechselrate zu den Genesenen
d=0.003 % Todesrate
%Anteil der erfassten Infizierten
r=0.1 % 10%
%%aktuelle Daten: (passen nicht mehr zur exponenziellen Funktion)
A=coronaDataneu();
inf_falleWHO=A(1,:);
inf_falleWHOrecovered=A(3,:);
inf_falleWHOtoten=A(2,:);
inf_falleWHOaktuell=inf_falleWHO-inf_falleWHOrecovered; % Infizierten - Genesene ergeben die aktuellen Infizierten
n=length(inf_falleWHO) %wichtig, Vektoren benötigen die gleiche Länge
timesWHO=[0:1:n-1];
figure(1)
title ('Corona- infizierten Detuschland, Quelle RKI')
plot (timesWHO, inf_falleWHO./1000, 'markersize', 1, 's', "linewidth",4, timesWHO,inf_falleWHOaktuell/1000,'markersize', 1, '*' , timesWHO,inf_falleWHOrecovered/1000,'markersize', 1,'s', timesWHO, inf_falleWHOtoten./1000 ,'markersize', 1, '*' ) %durch 1000 teilen, da die Daten absolute Zahlen sind
xlabel ('Tage von 26. Februar')
ylabel ('Zahlen in Tausend')
legend ( "German data I +R" ,"German data I ","German data R", "German data deaths" ,"location", "north")
grid on
xticks([50,100,150,200,250,300,350,400,450,500])
axis([0,n, 0 ,4000])
c=0.32450; %von unrsprünglichem Notebook entnommen
function val=slowdown(x)
% Ab Tag t (Argument der Funktion) wir die Infektionsrate in mehreren Schritten sinken:
if x<29 val=1.1 ;
else if x<36 val=val= (-0.0106*x+ 0.4675)/0.32450;
else if x<51 val=(-0.0054*x+0.2831)/0.32450;
else if x<81 val=(-0.0006*x+0.0516 +0.05)/0.32450 ;
else if x<201 val=(0.00003*x - 0.00091 +0.06)/0.32450;
else if x<251 val=(0.0006 *x - 0.1244+0.12)/0.32450;
else if x<281 val=(-0.000561 *x + 0.1724+0.1)/0.32450;
else if x<381 val=(-0.0002*x+ 0.0807+0.1)/0.32450;
else if x<455 val=(0.0003*x +0.15 )/0.32450;
else val=1.3;
endif
endif
endif
endif
endif
endif
endif
endif
endif
endfunction
%Slowdown-Funktion
for i=1:length(timesWHO)
cc(i)=slowdown (timesWHO (i) );
endfor
plot (timesWHO,cc)
hold on
%infektionsrate aus den RKI-Daten berechnet:
T=1
ccc(1)=cc(1);
for i=T+1:n
ccc(i-T)=(inf_falleWHO(i)-inf_falleWHO(i-T))/(inf_falleWHO(i-T));
% INfektionsrate aus dem beschränktem Modell
ccc2(i-T)=NM*(inf_falleWHO(i)-inf_falleWHO(i-T))/(inf_falleWHO(i-T)*(NM-0.001*inf_falleWHO(i-T)));
ccc3(i-T)=NM*(0.001*inf_falleWHO(i)-0.001*inf_falleWHO(i-T) +w*0.001*inf_falleWHO(i-T) )/(0.001*inf_falleWHO(i-T)*(NM-0.001*inf_falleWHO(i-T)-0.001*inf_falleWHOrecovered(i-T)));
endfor
% Plot
x=0:1:500;
for i=1:1:length(x)
vec2(i)=slowdown(x(i));
endfor
plot(x,vec2)
hold on
plot (ccc2, 'markersize', 0.5, '*', ccc3,'markersize', 0.5, 's')
Aufgabe 3: Lösen der Differentialgleichung mit Anfangsbedingungen
%Lösung der Differetialgleichung mit angepasster slowdown-Funktion und in Abhängigkeit von den Testzahlen angepasstem Erfassungsfaktor r
times=(0:0.1:454);
c=0.3245;
%Anfangsbedingungen
yo=[NM-16/1000;16/1000;0];%Anfällige im Zeitpunkt t0, Startwert darf nicht 0 sein, hier -16 ; Infizierte ; 0 Genesene
f=@(y,x) [-c*y(1)*y(2)/(1*NM);c*y(1)*y(2)/NM-w*y(2);w*y(2)-0*y(2)];
%f=@(y,x) [-slowdown(x,20)*y(1)*y(2)/(0.1*NM);slowdown(x,20)*y(1)*y(2)/NM-w*y(2);w*y(2)-0*y(2)];
yy = lsode (f, yo, times);%Lösung:
figure(1)
plot ( times, yy (:,2), '.', times, yy (:,3),'.',times, yy (:,3)+yy (:,2),'.',times, yy (:,1), '*' )
% Er hat vorher die alten berechneten Erbenisse geplottet
hold on
plot (timesWHO, inf_falleWHO./1000, 'markersize', 1,'s', "linewidth",4, timesWHO,inf_falleWHOaktuell/1000, 'markersize', 1, '*' , timesWHO,inf_falleWHOrecovered/1000,'markersize', 1.0, '*' ) %times arbeitet mit 1/10 eines Tages als Schrittweite
legend ( "Infected", "Removed", " Infected+Removed","German data I +R" ,"German data I ","German data R" ,"location", "north")
axis([0,n, 0 ,5000])
hold off
figure (2)
plot ( times, yy (:,2), '.', times, yy (:,3),'.',times, yy (:,3)+yy (:,2),'.',times, yy (:,1), '*' )
figure (3)
plot (timesWHO, inf_falleWHO./1000, 'markersize', 1,'s', "linewidth",4, timesWHO,inf_falleWHOaktuell/1000, 'markersize', 1, '*' , timesWHO,inf_falleWHOrecovered/1000,'markersize', 1.0, '*' )
Aufgabe 4: Einfluss des Erfassungsfaktors r
Je größer r, desto genauer entsprechen die erhobenen Daten der Realität und die Dunkelziffer der unentdeckten Fälle verringert sich.
r ist von mehreren Faktoren abhängig. Zunächst beeinflussen die Eigenschaften des Virus den Erfassungsfaktor. Würden alle Erkrankten eindeutige Symptome aufweisen, würde man auch alle erfassen können. Da sich Covid-19 allerdings dadurch auszeichnet, dass ein Teil der Erkrankungen asymptomatisch verläuft, spielt r eine wichtige Rolle in der Modellierung des Infektionsgeschehens. Darüber hinaus hat die Anzahl der Testungen einen wesentlichen Einfluss auf den r-Wert: Je mehr Menschen getestet werden, desto mehr Asymptomatisch-Erkrankte können entdeckt werden und desto niedriger die Dunkelziffer. Wählt man den Erfassungsfaktor r=1, so würden alle Infizierten auch erfasst werden. Allerdings ist dieser Fall unrealistisch, da selbst bei regelmäßigen Testungen nicht alle infizierten Personen erfasst werden können, da die Tests nur in einer ganz gewissen Zeitspanne mit hoher Virenlast überhaupt erst anschlagen.
Im Folgenden werden drei unterschiedliche Szenarien mit den Erfassungsfaktoren von 0.1, 0.5 und 0.8 entworfen und vergleichend gegenübergestellt. Dazu wird die Bevölkerung in der Differentialgleichung der Anfälligen mit dem jeweiligen Faktor multipliziert.
Anschließend wurde eine zeitlich-abhängige Funktion ähnlich der Slowdown-Funktion für den Erfassungsfaktor r aufgestellt, welche die Testzahlen in Deutschland im Verlauf des Pandemiegeschehens abbilden soll. Dabei wurden die Testzahlen des RKI als Orientierungshilfe genommen, um eine stückweise konstante Erfassungsfaktor-Funktion zu generieren.
Berücksichtigt werden hierbei besonders:
- kostenlose Bereitstellung von Tests für Hotelliere und Gastronomen (ab 2. November 2020/Kalenderwoche 45, Tag 251), gleichzeitig vermehrte Testungen durch das Aufkommen der mutierten Varianten
- Bereitstellung von einem kostenlosen Test pro Woche und Bürger (Kalenderwoche 9 2021, Tag 370)
- Einführung von verpflichtenden Selbsttests in Grundschulen und weiterführenden Schulen nach den Osterferien 2021 (etwa ab Kalenderwoche 14, Tag 405)
- dann: Wirkung der Impfungen, dementsprechend gehen die Testzahlen wieder etwas zurück (Zweifach Geimpfte und bereits Genesene müssen keinen Negativtest erbringen, um Gastronomie oder Geschäfte zu nutzen.) --> Erfassungsfaktor bzw. Dunkelziffer bleibt gleich!
function val=Erfassungsfaktor_r(x)
% t=13, bis einschließlich Kalenderwoche 10 wurden die Testzahlen in einem gemeinsamen Wert dargestellt
if x<13 val=0.1 ; %geringe Testzahlen zu Beginn der Pandemie
else if x<13 +111 val=0.2; % ab Tag 13 bis Tag 123 (Kalenderwoche 11 bis )
else if x<13 +140 val=0.25; % ab Tag 124 bis Tag 152 (Kalenderwoche bis)
else if x<13+238 val=0.5; % ab Tag 153 bis Tag 250 (Kalenderwoche bis )
else if x<13+280 val=0.8; % ab Tag 251 bis Tag 292 (Kalenderwoche 45 bis 51)
else if x<405 val=0.45; % ab Tag 293 bis Tag 404
else if x<431 val=0.7; % ab Tag 405 bis Tag 431 (Kalenderwoche 14 bis 17)
else if x<455 val=0.7; % ab Tag 432 bis Tag 454
endif
endif
endif
endif
endif
endif
endif
endif
endfunction
Das Ergebnis der Differentialgleichungen mit zeitlich-abhängigem Erfassungsfaktor ähnelt dem Szenario 2 mit konstantem mittlerem Erfassungsfaktor.
Aufgabe 5: Einbauen der Impfungen
Es wurden zunächst Daten des RKI zur Anzahl an täglich Geimpften in Deutschland recherchiert. Da man erst nach der Zweitimpfung optimal geschützt ist und im SIR-Modell erst dann in die Gruppe der Genesenen wechselt, haben wir lediglich die Daten zu den Zweitimpfungen betrachtet. Mit dem Befehl "movmean" wurden die Daten zunächst geglättet, um Schwankungen auszugleichen, die aufgrund von geringeren Impfzahlen am Wochenende zustande kommen.
x=325:454;
v1=[846 1021 17787 15843 36442 54580 41311 33571 46112 21833 38761 46899 58531 45873 57774 53643 31112 61561 67211 97846 72151 77913 48345 25046 56057 73818 77130 68512 75574 45456 25711 55781 56403 56302 52994 55581 36268 28378 57730 60472 59149 58958 66891 39892 24064 49866 52105 70506 62832 68539 48656 33483 51592 54801 67041 60896 70860 50277 35680 59226 68170 80927 75592 84239 53270 38190 76462 84056 89606 81159 88540 62335 53264 88613 92602 101517 90710 67980 64374 50368 70215 94704 96368 87674 88923 65056 50320 76376 73697 83212 70304 74194 56616 39006 61415 62771 73137 67195 68931 58948 49048 73911 91455 126260 132603 136224 73675 65574 114770 158780 227837 219092 234216 152356 101493 189805 284220 360967 158205 222738 161250 117245 196259 346892 521981 481187 433280 244691 151487 176811];
v=v1.*0.001;
plot(x,v)
hold on
M=movmean(v,20)
plot(x,M)
hold off
Schließlich erfolgte eine lineare Approximation der geglätteten Daten, um eine Impffunktion für die täglich Geimpften aufzustellen.
%Impffunktion für die täglich Geimpften
function val=Impffunktion(x)
if x<325 val=0;
else if x<345 val= 1.545*x-475.22;
else if x<400 val= 0.342*x-60.127;
else if x<420 val= -0.305*x+198.6;
else if x<454 val=6.088*x-2486.459;
else val=0;
endif
endif
endif
endif
endif
endfunction
%Durch die beiden Daten wird eine lineare Funktion gelegt.
% f(x)=1.545*x+-475.225 % 1. Abschnitt (Tag 325 bis 345)
% f(x)=0.342*x+-60.127 % 2. Abschnitt (Tag 345 bis 400)
%f(x)=-0.305*x+198.6 % 3. Abschnitt (Tag 400 bis 420)
%f(x)=6.088*x+-2486.459 % 4. Abschnitt (Tag 420 bis 454)
%Plot der Impffunktion
x=0:1:500;
for i=1:1:length(x)
vac(i)=Impffunktion(x(i));
endfor
plot(x,vac)
Zur besseren Interpretation werden nachfolgend die Ergebnisse der Differentialgleichung mit und ohne Impfungen gegenübergestellt. Im Szenario mit Impfungen wurde die Impffunktion in der Gruppe der Anfälligen abgezogen und zur Gruppe der Genesenen addiert. Zusätzlich haben wir die Impffunktion mit dem Faktor 1,65 multipliziert und somit etwas hochgesetzt, da durch Anwendung des Moving Average auf die realen Daten vor allem gegen Ende des betrachteten Zeitraums deutlich niedrigere Werte entstanden.
Tutorium 4 - Finite-Differenzen-Verfahren
Bearbeiten1.Teil: Implementierung der Dirichlet-Randbedingung
BearbeitenZunächst wurden die Dirichlet Randbedingungen auf allen Rändern gleich gewählt.
Zu Beginn wurde eine Quellfunktion für die Funktion der rechten Seite implmentiert:
clear all
function wert=Quellfunktion(x,y)
if sqrt((x.-3.25).^2+(y.-3).^2)<=0.35 wert=1;
else wert=0;
endif
endfunction
Anschließend wurde das Gebiet und die Gittergröße festgelegt:
%definiere Gebiet D
xend=8;
yend=5.0;
%Anzahl von X-Punkten, verkleinern Sie diese um die Struktur der Matrix in weiteren Zellen besser zu erkennen.
N=10
%--------------------
hx=xend/(N+1)
hy=hx; %damit ein äquidistantes Gitter entsteht
M=floor((yend/hy))-1
%Koordinatenmatrix der Gitterpunkte
[x,y]=meshgrid(hx:hx:xend-hx,hy:hy:yend-hy);
Im Anschluss wurde die Finite Differenzen Methode angewandt:
%konstanter Diffusionskoeffizient
a=1.0; %konstanter Diffusionskoeffizient
%--------------------------
%RANDBEDINGUNGEN Vorbereitung fuer eine Teilaufgabe]
%-------------------------------------------
%Dirichlet RB mit Konstanten und gleichen Rändern, alle Ränder haben den gleichen Wert:
wert_unten=0.15;
wert_oben=0.15;
wert_links=0.15;
wert_rechts=0.15;
%-------------------------------------------
% Systemmatrix
%-------------------------------------------
Vec1=ones(N*M,1); %erzeugt Vektor aus Einsern der Laenge N*M
%Systemmatrix N^2xM^2 mit -4 auf der Diagonale und 1 auf den Nebendiagonalen und 1 auf den N und -N-ter Nebendiagonale
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);
%Korrektur der Matrix (Blocktridiagonalitaet)
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; %die aufgestellte Matrix muss mit h^2 (h=hx=hy; damit ein äqudistantes Gitter entsteht) dividiert werden.hx^2 entsteht, da Laplace u mit Hilfe der 2. Differenzenquotienten approximiert wird. Der Diffusionskoeffizient hat in diesem Fall keine Auswirkungen, da dieser gleich 1 gesetzt wurde
%Matrixdarstellung - optional
%-----------------------------
figure (1);
%surface (BB);
spy (BB);
xlabel("Spalten")
ylabel("Zeilen")
zlabel("Matrix B")
title ("Systematrix B");
Einbezug der rechten Seite des Systems:
%rechte Seite des Systems
for i=1:N
for j=1:M
f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
%Vorgehensweise nach dem Sternschema:
%Dirichlet Randbedingung unten und oben
if j==M f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+ wert_oben*a/hx^2; endif
if j==1 f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+ wert_unten*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
% Wichtig!: Ecken-Beiträge: da in den Ecken jeweils der Differenzenquotion in x (Wert Links oder rechts)
% als auch in y Richtung (oben/unten) approximiert und aufadiert werden muss
%Beispiel: linke untere Ecke wird durch linken und unteren Rand definiert/beeinflusst
%andere Ecken gehen analog
f(1,1)=1*Quellfunktion(x(j,i),y(j,i))+ wert_links*a/hx^2 + wert_unten*a/hx^2; %%die Randwerte müssen mit h^2 dividiert werden, a=1
f(1,N)= 1*Quellfunktion(x(j,i),y(j,i))+ wert_rechts*a/hx^2+ wert_unten*a/hx^2;
f(M,1)=1*Quellfunktion(x(j,i),y(j,i))+ wert_links*a/hx^2+ wert_oben*a/hx^2;
f(M,N)=1*Quellfunktion(x(j,i),y(j,i))+ wert_rechts*a/hx^2 + wert_oben*a/hx^2;
size (f)
%reshape funktioniert Spaltenmässig: es wird Spalte1, dann Spalte2, usw eingeordnet, wobei wir die Matrix f nach Reihen/Zeilen einordnen.
%Daher soll die Matrix f' mit N(Zeilen)x M(Spalten) in ein langern Vektor eingeordnet. Assemblierung
b=-1*reshape(f',N*M,1);
Im letzten Schritt wird das System gelöst:
% LOESUNGSSCHRITT
%-------------------------------------------------
sol=BB\b;
sol_matrix=reshape(sol,N,M);% Matrix mit N(Zeilen)x M(Spalten)
sol_matrix=sol_matrix'; % Transponiert, wegen reshape
%-------------------------------------------------
%Oeffne ein Bildfenster und plotte die Lösung auf dem realen Koordinatenfeld
%-------------------------------------------------
figure(2);
surfc(x,y,sol_matrix);
title ("Loesung");
ylabel("y")
xlabel("x")
colormap("hsv")
colorbar
%weiteres Bild mit Quellfunktion- rechte Seite des Systems
figure(3);
surfc(x,y,f);
title ("Quelfunktion");
ylabel("y")
xlabel("x")
%%colormap("hsv")
Im Anschluss wurden Veränderungen vorgenommen:
Zum einen an der Quellfunktion:
clear all
function wert=Quellfunktion(x,y)
if exp((x.-5).^2+(y.-3).^2)<=0.35 wert=1; %die vorgegebene Quellfunktion wurde in diesem Schritt abgeändert
else wert=0;
endif
endfunction
Zum anderen an den Dirichlet-Randbedingungen. Dies sind zwar weiterhin konstant, aber nicht mehr auf allen Rändern gleich, sondern verschieden:
%Dirichlet RB mit konstanten aber unterschiedlichen Randbedingungen:
wert_unten=0.5;
wert_oben=-0.7;
wert_links=0.15;
wert_rechts=0.8;
Die nachfolgenden Schritte verlaufen analog.
Fazit: Die Quellfunktion an sich hat nicht so einen starken Einfluss auf die Lösung des Systems, wie die Dirichlet-Randbedinungen
2. Teil: Gemischte Ränder: Implementierung sowohl der Neumann, als auch Dirichlet Randbedingungen
Zunächst wird auch hier wieder eine Quellfunktion und ein Gebiet definiert. Dies ist analog zu der Implementation des 1. Teils.
Anschließend werden die Dirichlet-Bedingungen festgelegt:
wert_links = 0.15;
wert_rechts = 0.15;
Implementierung der Neumann-Bedingungen:
% Hilfestellung von der Wikiversity-Seite
% Neumann Randbedingung:
N=N+2; %da das Gitter für das Neumann-Randwertproblem erweitert wurde, um den zentralen Differenzenquotient anzuwenden
M=M+2;
%Systemmatrix aufstellen
Vec1=ones(N*M,1);
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);
%Korrektur der Matrix (Blocktridiagonalitaet)
for i=1:(M-1)
BB(i*N+1,i*N)=0;BB(i*N,i*N+1)=0;
endfor
BB=BB*(a/hx^2);
%Neumann Randbedingung muss nur für den oberen und unteren Rand implementiert werden, da auf den anderen Rändern die Dirichlet Randbedingung gilt
% Neumann RB für y: c* partial_y u=0 - einzelne Blöcke der Matrix ersetzen
%----------------------------------------------------------------------
% da die Neumann Bedingung auf dem oberen und unteren Rand implementiert wurde, müssen demensprechend diese Ränder angepasst werden
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: c* partial_x u=0 - einzelne Zeilen der Matrix ersetzen
%----------------------------------------------------------------------
% for i=1:(M-2)
%vector_up=zeros(1,N*M); %0 Zeile
%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 ;
%BB(i*N+N,:)=vector_down ;
%endfor
Anpassung der rechten Seite:
% im nachfolgenden wird die rechte Seite des Systems betrachtet
for i=1:N
for j=1:M
f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
%Dirichlet Randbedingung auf dem rechten und linken Rand
if (i==N ) f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_rechts*a/hx^2; endif
if (i==1) f(j,i)=1*Quellfunktion(x(j,i),y(j,i))+wert_links*a/hx^2; endif
endfor
endfor
%im nachfolgenden Schritt müssen die Ecken angepasst werden, da in allen 4 Ecken Dirichlet und Neumann aufeinandertreffen kommt es zu keiner Summation wie etwa in Aufgabe a)
%Ecke rechts unten
% AH richtig mit den Eckpunkten, weil durch die Struktur der Matrix realisieren sich in den Eckpunkten die y-Ableitungen Randbedingungen.
f(1,N)= 1*Quellfunktion(x(1,N),y(1,N))+(0*wert_rechts*a/hx^2);
%Ecke rechts oben
f(M,N) = 1*Quellfunktion(x(M,N),y(M,N))+(0*wert_rechts*a/hx^2);
%Ecke links unten
f(1,1) = 1*Quellfunktion(x(1,1),y(1,1))+(0*wert_links*a/hx^2);
%Ecke links oben
f(M,1) = 1*Quellfunktion(x(M,1),y(M,1))+(0*wert_links*a/hx^2);
b=-1*reshape(f',N*M,1);
Abschließend wird das lineare System analog zum 1. Teil gelöst
3. Teil: Implementierung der Neumann-Bedingung
Achtung: Diese Bedingung kann nicht erfüllt werden, da dadurch die Matrix singulär wird. Jedoch gibt Octave eine Lösung aus mit der entsprechenden Warnung.
Definition der Quellfunktion:
function wert=Quellfunktion(x,y)
if sqrt((x.-3.25).^2+(y.-3).^2)<=0.35 wert=1;
else wert=0;
endif
endfunction
Implementierung der Neumann-Bedingung auf allen Rändern:
%Nutzung der Hilfestellung von der Wikiversity-Seite
% Neumann Randbedingung:
N=N+2;
M=M+2;
%Systemmatrix aufstellen
Vec1=ones(N*M,1);
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);
%Korrektur der Matrix (Blocktridiagonalitaet)
for i=1:(M-1)
BB(i*N+1,i*N)=0;BB(i*N,i*N+1)=0;
endfor
BB=BB*(a/hx^2);
% Neumann RB für y: c* partial_y u=0 - einzelne Blöcke der Matrix ersetzen
%----------------------------------------------------------------------
%da die Neumann Bedingung auf dem unteren Rand implementiert wurde, muss auch dementsprechend der untere Rand angepasst werden
block=diag((a/hx)*ones(N,1));
%unterer Rand
% AH am unterem Rand Neumann!
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) %AH hier bis M-2, richtir, weil Blocke 2,3, bis vorletzter Block muessen x-Ableitung approximeiren
vector_up=zeros(1,N*M); %0 Zeile
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 ;
BB(i*N+N,:)=vector_down ;
endfor
Anpassung der rechten Seite:
for i=1:N
for j=1:M
f(j,i)=1*Quellfunktion(x(j,i),y(j,i));
endfor
endfor
b=-1*reshape(f',N*M,1);
Anschließend wird das lineare System mit den bekannte Befehlen gelöst.
Tutorium 5 - Räumliche Epidemieausbreitung mit Reaktionsdiffusionsprozess
BearbeitenAusgangsparameter
Bearbeiten %Gebiet wurde vergrößert
xend=30; %in km
yend=40; %in km
N=60 %Unterteilung der x-Achse
a=0.1; %(konstanter Diffusionskoeffizient)
T=240 %(Zeitintervall in Tagen) so lange werden die Effekte der Maßnahmen durch unsere Slowdown abgebildet, (jedoch keine 454 Tage, da dies zu viel Zeit zum Ausführen benötigt)
delta=0.1
%Für Reaktionsdiffusionsgleichung:
B=0.206; %(Bevoelkerungsdichte: 206 Einwohner pro km^2 in Rheinland-Pfalz)
c=0.32450; %(Infektionsrate, Vorsicht k=c/B in unserem Modell k.u.(B-u)-w.u)
w=1/13.9; %Wechselrate
%--------------
hx=xend/(N+1)
hy=hx;
M=floor((yend/hy))-1
Nr_timesteps=floor (T/delta)
Ausgangsfunktion
Bearbeiten% Anfangsfunktion
% 2 Hotspots in Anlehung an die Bevölkerungsdichte
function wert=initialfunction(x,y)
wert=0;
r1=0.5;
r2=0.3;
if sqrt((x.-6).^2+(y.-8).^2)<=r1
wert=1/(pi*r1^2);
elseif sqrt((x.-16).^2+(y.-4).^2)<=r2 wert=0.5/(pi*r2^2);
endif
endfunction
Ergänzungen Octave-Tutorial
BearbeitenDiese 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
BearbeitenReferenzieren Sie hier die von Ihnen verwendeten Skripte und Bibliotheken in Octave oder R/Studio
Literatur
BearbeitenNotieren Sie hier die von Ihnen verwendete Literatur
- Quelle 1