|
|
|
SIR Problem lösen mit 4- Stufen Runge- Kutta |
|
| marcella |
Forum-Newbie
|
 |
Beiträge: 2
|
|
|
|
Anmeldedatum: 15.11.14
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
 |
|
Verfasst am: 15.11.2014, 15:35
Titel: SIR Problem lösen mit 4- Stufen Runge- Kutta
|
 |
| |
 |
 |
Hallo ihr Lieben,
ich möchte gerne das SIR Epidemiemodell mit dem 4- stufigen Runge Kutta Verfahren lösen. Dabei geht es um folgende Differenzialgleichungen:
(I) d/dt u(t) = - Ru(t)v(t)
(II) d/dt v(t) = Ru(t)v(t)- v(t)
(III) d/dt w(t) = v(t)
R ist dabei eine Konstante. Zudem soll gelten u(t)+v(t)+w(t) = 1.
Ich hab das Verfahren nun zunächst mit ODE45 gelöst. Klappt wunderbar. Leider sollen wir nun das Runge Kutta Verfahren selber programmieren und das mag nicht so 100% klappen. Die Gleichungen u(t) und v(t) werden richtig geplottet, nur w(t) wird im Vergleich zu meiner ersten Lösung (mit der vordefinierten ODE45 Funktion) nicht richtig dagestellt (der Steigung ist etwas zu steil und somit der letzte Wert zu hoch, so dass u(t)+v(t)+w(t)=1 nicht mehr erfüllt ist).
Hier meine erste (richtige) Lösung:
%SIR Modell_ ODE45
figure();
r=1.6;
[t,Y]=ode45(@(t,y)[-1*r*y(1)*y(2); (r*y(1)-1)*y(2);y(2)],[0,25],[0.99;0.01;0]);
plot(t,Y)
title ('Krankheitsverlauf')
xlabel ('t');
ylabel ('Anteil der Population');
legend ('Gesund', 'Krank', 'Immun');
%Ende Code 1
Hier meine fehlerhafte zweite Lösung:
%4stufiges Runge Kutta Verfahren zum Lösen der Differenzialgleichung
h= 0.5; % Schrittweite
R= 1.6; % Reproduktionsrate, konstant
x= 0:h:25; % Simulation für 25 Tage
y = zeros(1,length(x));
z = zeros (1,length(x));
w = zeros (1,length(x));
y(1) = 0.99;
z(1) = 0.01;
w(1) = 0.0;
F_xyz = @(t,u,v) -R*u*v; %Differenzialgleichungssystem
G_xyz = @(t,u,v) R*u*v-v;
H_xz = @(t,v) v
for i=1:(length(x)-1) %Runge Kutta Verfahren
k_1 = F_xyz(x(i),y(i),z(i));
l_1 = G_xyz (x(i),y(i),z(i));
m_1 = H_xz (x(i),z(i));
k_2 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_1), (z(i)+0.5*h*l_1));
l_2 = G_xyz ((x(i)+0.5*h),(y(i)+0.5*h*k_1),(z(i)+0.5*h*l_1));
m_2 = H_xz ((x(i)+0.5*h), (z(i)+0.5*h*m_1));
k_3 = F_xyz((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*l_2));
l_3 = G_xyz ((x(i)+0.5*h),(y(i)+0.5*h*k_2),(z(i)+0.5*h*l_2));
m_3 = H_xz ((x(i)+0.5*h), (z(i)+0.5*h*m_2));
k_4 = F_xyz((x(i)+h),(y(i)+h*k_3),(z(i)+h*l_3));
l_4 = G_xyz ((x(i)+h),(y(i)+h*k_3),(z(i)+h*l_3));
m_4 = H_xz ((x(i)+h), (z(i)+h*m_3));
y(i+1) = y(i) + (1/6)*h*(k_1+2*k_2+2*k_3+k_4);
z(i+1) = z(i) + (1/6)*h*(l_1+2*l_2+2*l_3+l_4);
w(i+1) = w(i) + (1/6)*h*(m_1+2*m_2+2*m_3+m_4);
end
plot (x,y,x,z,x,w)
title('Krankheitsverlauf');
xlabel('Tage');
ylabel('Population');
legend('Gesund', 'Infiziert','Immun');
%Ende Code 2
Vielleicht sieht jemand den Fehler  Ich wäre euch super, super dankbar!!
Liebe Grüße,
Marcella
p.s. die beiden plots hänge ich an
| Beschreibung: |
|
Download |
| Dateiname: |
Bildrichtig.png |
| Dateigröße: |
5.25 KB |
| Heruntergeladen: |
311 mal |
| Beschreibung: |
| hier die falsche Trajektorie für w(t) |
|
Download |
| Dateiname: |
Bildfalsch.png |
| Dateigröße: |
5.23 KB |
| Heruntergeladen: |
290 mal |
|
|
|
|
|
|
| Harald |
Forum-Meister
|
|
Beiträge: 24.499
|
|
|
|
Anmeldedatum: 26.03.09
|
|
|
|
Wohnort: Nähe München
|
|
|
|
Version: ab 2017b
|
|
|
|
|
|
|
Verfasst am: 15.11.2014, 17:28
Titel:
|
|
Hallo,
wenn ich das richtig sehe, ist es nicht falsch, sondern nur ungenau.
Dein RK hat ja im Gegensatz zu ode45 keine Schrittweitensteuerung. Du musst die Schrittweite also klein genug wählen. Bei h = 0.05 sieht es gut aus.
Grüße,
Harald
|
|
|
|
| marcella |
Themenstarter
Forum-Newbie
|
|
Beiträge: 2
|
|
|
|
Anmeldedatum: 15.11.14
|
|
|
|
Wohnort: ---
|
|
|
|
Version: ---
|
|
|
|
|
|
|
Verfasst am: 15.11.2014, 17:37
Titel:
|
|
Lieber Harald!
Oh man  super, ja so sieht es gut aus.
Vielen, vielen Dank!
Marcella
|
|
|
|
|
|
|
Einstellungen und Berechtigungen
|
|
Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
|
|
Impressum
| Nutzungsbedingungen
| Datenschutz
| FAQ
|  RSS
Hosted by:
Copyright © 2007 - 2025
goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
|
|
