Kurs:Räumliche Modellbildung/Gruppe Nr

Aus Wikiversity

Gruppenseite - (GBN)[Bearbeiten]

Diese 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[Bearbeiten]

Ziel 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[Bearbeiten]

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

Ausgangsidee


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

Zyklus 1

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.

Zyklus 2

Implementieren der Personengruppen S,I und R[Bearbeiten]

Zu 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[Bearbeiten]

Beim Transportprozess gibt es grundsätzlich zwei Möglichkeiten, damit die im Raum befindlichen Personen nicht aus dem Gitter verschwinden können:

  1. Die Person kommt am anderen Rand wieder in das Gitter hinein. (zyklische Randbedingung)
  2. 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[Bearbeiten]

Zunä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[Bearbeiten]

Hier 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[Bearbeiten]

Hier 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[Bearbeiten]

Tutorium 1 - Kontaktmodell[Bearbeiten]

Einführung zum Kontaktmodell[Bearbeiten]

Mögliche Szenarien[Bearbeiten]

Szenarien:[Bearbeiten]

Szenario 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 1: Kontaktmodell mit hoher Anzahl von Infektionen

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

Szenario 2: Kontaktmodell mit einer geringen Anzahl von Infektionen

Vergleich Szenario 1 und Szenario 2[Bearbeiten]

Das 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)

Szenario 3: Infizierter befindet sich in der Mitte
Szenario 3: Infizierter befindet sich am Rand


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[Bearbeiten]

Durch 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[Bearbeiten]

Anmerkung: 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[Bearbeiten]

Im 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

Fundamentallösung LaPlace-Gleichung 2D

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 der Diffusionsgleichung

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:

Funktion der rechten Seite
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
Lösung Poissiongleichung für spezielle rechte Seite
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)
Funktion der rechten Seite

Ä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:

Lösung der Poissiongleichung für die rechte Seite
Zeitabhängige Diffusionsgleichung

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[Bearbeiten]

Aufgabe 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)
Pandemiebeginn

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])
Daten Deutschland
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')
Gegenüberstellung: Slowdown-Funktion, ccc2 und ccc3


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, '*' )
Lösungen der Differentialgleichungen ohne Maßnahmen
Tatsächliche Daten Deutschland


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.

Szenario 1: r=0,1
Szenario 2: r=0,5
Szenario 3: r=0,8

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
Zeitabhängiger Erfassungsfaktor

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
Glättung der Daten täglich Geimpfter

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)
Linear vereinfachte Impffunktion


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[Bearbeiten]

1.Teil: Implementierung der Dirichlet-Randbedingung[Bearbeiten]

Zunä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:

Systemmatrix des Finiten Differenzen Verfahren
%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")
Quellfunktion mit den Dirichlet-Randbedingungen
Lösung des Systems



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.

Quellfunktion mit unterschiedlichen aber konstanten Dirichlet-Randbedingungen





Lösung des Systems mit unterschiedlichen Dirichlet Bedingungen


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

Quellfunktion mit gemischten Randbedingungen


Lösung des Finite Differenzen Verfahren mit gemischten Randbedingungen


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.

Quellfunktion mit Neumann-Bedingungen auf allen Rändern





Lösung des Finiten Differenzen Verfahren mit Neumann-Bedinung auf allen Rändern


Tutorium 5 - Räumliche Epidemieausbreitung mit Reaktionsdiffusionsprozess[Bearbeiten]

Ausgangsparameter[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[Bearbeiten]

Diese Ergänzungen wurden bei dem Octave-Tutorial vorgenommen, um die Implementation der Gruppe bzgl. der verwendeten Octave-Befehle nachvollziehbar zu machen.

Skripte in Octave[Bearbeiten]

Referenzieren Sie hier die von Ihnen verwendeten Skripte und Bibliotheken in Octave oder R/Studio

Literatur[Bearbeiten]

Notieren Sie hier die von Ihnen verwendete Literatur

  • Quelle 1