Gelöste Aufgaben/COVI

Aufgabenstellung

Als Ingenieure können wir die COVID-19-Epidemie genauso modellieren, wie technische Systeme.

Wir brauchen dazu "Koordinaten" - also Messgrößen - die die Anzahl der Individuen erfassen, die

• ansteckbar "a",
• infiziert "i" - mit und ohne Sympthome - und
• genesen "r"

sind. Für jede dieser Koordinaten müssen wir nun "Bewegungsgleichungen" - also z.B. Differentialbeziehungen in der Zeit - hinschreiben.

Das nennen wir "Modellbildung".

🧨 Dies ist kein zuverlässiges Prognosewerkzweug:

Hier geht es um die Modellierung der zeitlichen Zusammenhänge, die Interpretation der Zustandsgrößen des Systems und das Verständnis der Abläufe bei der Verbreitung des Virus in Abhängigkeit von zentralen Systemparametern. Das Modell ist nicht dafür gemacht, um daraus quantitative Schlussfolgerungen zu ziehen.

Numerische Lösung der COVID-19 Pandemmie als Anfangswertproblem.

Diese algebraischen und Differentialgleichungen sind dabei nicht das Ergebnis von Gleichgewichtsbeziehungen wie in der Technischen Mechanik. Wir begnügen uns statt dessen damit, Phänomene der Epidemie mit unserer Mathematik zu erfassen.

Und so gehen wir vor:

Wir modellieren eine Grundgesamtheit von n₀ = 80 Millionen Individuen.

Die Pandemie soll bei einer Durchseuchung von α=50% zum Stillstand kommen - die Individuen sind dann so weit voneinander entfernt, dass eine Ansteckung unwahrscheinlich ist.

Zum Zeitpunkt t₀=0 seien von insgesamt i₀ Individuen mit

$i_{0}=1000$

infiziert.

Lösung mit Matlab

Header

Wir lösen hier das Anfangswertproblem zu nichtlinearen Bewegungsgleichungen.

Diese Gleichungen haben wir nicht - wir entwickeln sie ad-hoc und achten nur darauf, zentrale Phänomene abzubilden. Eine Abbildung der realen Zahlen ist nicht beabsichtigt.

%*******************************************************
%* matlab script                                       *
%* version: R2020a                                     *
%* author: Andreas Baumgart                            *
%* last updated: 2020-03-25                            *
%* ref: Technische Mechanik mit Computer               *
%* description: COVID-19 Simulation                    *
%*              no spatial resoultion                  *
%*******************************************************

Declarations

Wir verwenden diese Parameter aus einer Tabelle

EXCEL-Datei mit Parametern

und erklären unten, wie wir sie einsetzten.

cd 'C:\Users\abs384\OneDrive\Confluence Sources\COVI'
addpath('.\functions')
 
% load parameters
parameter = readtable('parameter.xlsx');
for row=1:length(parameter.name)
   eval(strcat(char(parameter.name(row)),'=', num2str(parameter.zahl(row))))
end

Equations of Motion

Die Bewegungsgleichungen schreiben wir für die Anzahl der Individuen an,

• die infiziert sind (Koordinate "i") und
• die genesen sind (Koordinate "r").

Die Anzahl der Individuen, die sich noch anstecken können, ergibt sich dann aus

$a(t)=n_{0}-i(t)-r(t)$

Unsere "Minimal-Koordinaten" des Systems sind also i und r.

Jetzt lassen wir unserer Fantasie freien Lauf:

function dydt = covid19(t,y,dt,DS,n,a0,lambda)
%   COVID19 Evaluate the COVID19 ODEs
%
%   last updated: 2020-03-24
%   author: Andreas Baumgart
 
% present state of pandemy
q = y(end,:);
% individuals left to be infected
a = n - sum(q);
 
% assign i=0 for t<0
if length(y)>DS
    i = y(end-DS,1);
else
    i=0;
end
% recovery rate
r = dt*lambda(2)*i;
      % infected                   % recovered
dydt=[dt*lambda(1)*q(1)*(a-a0)-r,  r          ];
 
end

Solving

1: Die Bewegungsgleichung für r

Die Änderung der Anzahl der Genesenen ist proportional zu Anzahl der Infizierten. Man könnte man also

$\displaystyle {\frac {d}{dt}}r(t)=\lambda _{2}\cdot i(t)$

anschreiben. Allerdings genesen zum Zeitpunkt t die Personen, die sich zum Zeitpunkt t-ΔT infiziert hatten. Dazwischen liegt gerade der gesamte Krankheitsverlauf mit der Dauer ΔT. Also schreiben wir

$\displaystyle {\frac {d}{dt}}r(t)=\lambda _{2}\cdot i(t-\Delta T)$

und nennen ΔT die Genesungsdauer - oder in Begriffen der Schwingungslehre: Totzeit (engl.: "dead-time").

2: Die Bewegungsleichungen für i

Die Änderung der Anzahl der Infizierten steigt mit der Anzahl der Infizierten i selbst und mit der Anzahl der ansteckbaren Individuen a. Sie reduziert sich um die Änderung der Anzahl der genesenen Individuen r. Allerdings weiß man aus Erfahrung, dass eine Epedemie endet, wenn ein bestimmter Durchseuchungsgrad α erreicht ist. Also wählen wir

$\displaystyle {\frac {d}{dt}}i(t)=\lambda _{1}\;\;i(t)\cdot (a(t)-a_{0})-{\frac {d}{dt}}r(t)\;\;{\text{ mit }}a_{0}=\alpha \cdot n_{0}$.

Beim Auftragen des Produkts i∙(a-a0) unten erkennt man gut, wie die Funktion für i=0 und (a-a0)=0 Null liefert - also die Zuwachsrate der Infizierten Null ist, wenn niemand infiziert ist oder die Anzahl der ansteckbaren Individuen Null ist. Dazwischen steigt die Funktion steil an:

Die "Ansteckungsfunktion" des Modells

Damit haben wir zwei Bewegungsdifferentialgleichungen, die wir in die Matlab-Funktion "covid19" schrieben.

Die Bewegungsgleichungen machen wir dimensionslos mit

$t=\Delta T\cdot \tau$

und erhalten als Bewegungsgleichungen

$\displaystyle {\frac {d}{d\tau }}\left({\begin{array}{c}i(\tau )\\r(\tau )\end{array}}\right)=\left({\begin{array}{c}\Delta T\cdot \lambda _{1}\;\cdot \;i(\tau )\cdot (a(\tau )-a_{0})-\displaystyle {\frac {dr(\tau )}{d\tau }}\\\Delta T\cdot \lambda _{2}\;\cdot \;i(\tau -1)\end{array}}\right)$.

% limit of infection
a0= n*alpha

% initial values in Mio. individuals
y = [ 1E-3 0];

% time steps and "dates"
dt=0.1;
t=0:dt:50;

% dead-time given in integration steps
DS=int16(round(DT/dt));

% make lambda dim'less
lambda = DT*lambda;

% Eulers Method
for j=2:length(t)   
    y(end+1,:)=y(end,:) + covid19(t,y,dt,DS,n,a0,lambda);
end

Postprocessing

Die Ergebnisse plotten wir für verschiedene Werte von λ1:

... für λ₁ = 0.016

... für λ₁ = 0.012

Sie sehen hier, wie unterschiedliche Ansteckungsraten zu unterschiedlich hohen Maximalwerten in den Infizierten i zu einem bestimmten Zeitpunkt führen. Und warum es sinnvoll ist, dieses λ₁ niedrig zu halten.

% plot results
hold on
plot(t,y(:,1))
plot(t,y(:,2))
plot(t,n-y(:,1)-y(:,2))
hold off
legend('i(t)','r(t)','a(t)')
xlabel('t / tBez →') 
ylabel('Individuen / Mio →')