Runge-Kutta (RK4) für System von Differentialgleichungen in Java
-
08-10-2019 - |
Frage
Dieses quation ist vor allem ein Ergebnis dieses Thread. Differentialgleichungen in Java
Im Grunde habe ich versucht, Jason S. beraten zu folgen und numerische Lösungen für Differentialgleichungen über Runge-Kutta-Verfahren (RK4).
Hallo an alle, Ich versuche, ein einfaches Simulationsprogramm von SIR-Epidemien Modell in Java zu erstellen. Grundsätzlich ist SIR definiert durch ein System von drei Differentialgleichungen:
S '(t) = - Lamda (t) * S (t)
I '(t) = lamda (t) * S (t) - gamma (t) * I (t)
R '(t) = gamma (t) * I (t)
S - anfällige Menschen, I - Infizierten, R - erholt Menschen. lamda (t) = [c * x * I (t)] / N (T) c - Anzahl der Kontakte, x - Infektiosität (Wahrscheinlichkeit krank nach dem Kontakt mit Kranken zu bekommen), N (t.) - Gesamtbevölkerung (die konstant)
gamma (t) = 1 / Dauer der Krankheit (konstant)
Nach dem ersten nicht sehr erfolgreichem Versuch, ich habe versucht, diese Gleichungen mit Runge-Kutta zu lösen, und dieser Versuch im folgenden Code resultierenden:
package test;
public class Main {
public static void main(String[] args) {
double[] S = new double[N+1];
double[] I = new double[N+1];
double[] R = new double[N+1];
S[0] = 99;
I[0] = 1;
R[0] = 0;
int steps = 100;
double h = 1.0 / steps;
double k1, k2, k3, k4;
double x, y;
double m, n;
double k, l;
for (int i = 0; i < 100; i++) {
y = 0;
for (int j = 0; j < steps; j++) {
x = j * h;
k1 = h * dSdt(x, y, S[j], I[j]);
k2 = h * dSdt(x+h/2, y +k1/2, S[j], I[j]);
k3 = h * dSdt(x+h/2, y+k2/2, S[j], I[j]);
k4 = h * dSdt(x+h, y + k3, S[j], I[j]);
y += k1/6+k2/3+k3/3+k4/6;
}
S[i+1] = S[i] + y;
n = 0;
for (int j = 0; j < steps; j++) {
m = j * h;
k1 = h * dIdt(m, n, S[j], I[j]);
k2 = h * dIdt(m+h/2, n +k1/2, S[j], I[j]);
k3 = h * dIdt(m+h/2, n+k2/2, S[j], I[j]);
k4 = h * dIdt(m+h, n + k3, S[j], I[j]);
n += k1/6+k2/3+k3/3+k4/6;
}
I[i+1] = I[0] + n;
l = 0;
for (int j = 0; j < steps; j++) {
k = j * h;
k1 = h * dRdt(k, l, I[j]);
k2 = h * dRdt(k+h/2, l +k1/2, I[j]);
k3 = h * dRdt(k+h/2, l+k2/2, I[j]);
k4 = h * dRdt(k+h, l + k3, I[j]);
l += k1/6+k2/3+k3/3+k4/6;
}
R[i+1] = R[i] + l;
}
for (int i = 0; i < 100; i ++) {
System.out.println(S[i] + " " + I[i] + " " + R[i]);
}
}
public static double dSdt(double x, double y, double s, double i) {
return (- c * x * i / N) * s;
}
public static double dIdt(double x, double y, double s, double i) {
return (c * x * i / N) * s - g * i;
}
public static double dRdt(double x, double y, double i) {
return g*i;
}
private static int N = 100;
private static int c = 5;
private static double x = 0.5;
private static double g = (double) 1 / x;
}
Dies scheint nicht zu arbeiten, weil Anzahl von Kranken (I) sollte zuerst increse und dann spezielle Massagen und Behandlungen auf etwa 0, und die Anzahl der wiedergewonnenen Menschen sollte unbedingt erhöhen. Gesamtzahl der Kranken + gesund + erholt 100 sein sollte, aber mein Code produziert einige seltsame Ergebnisse:
99.0 1.0 0.0
98.9997525 0.9802475 0.03960495
98.99877716805084 0.9613703819491656 0.09843730763898331
98.99661215494893 0.9433326554629141 0.1761363183872249
98.99281287394908 0.9261002702516101 0.2723573345404987
98.98695085435723 0.9096410034385773 0.3867711707625441
98.97861266355956 0.8939243545756241 0.5190634940761019
98.96739889250752 0.8789214477384787 0.6689342463444292
98.95292320009872 0.8646049401404658 0.8360970974155659
98.93481141227473 0.8509489367528628 1.0202789272217598
98.91270067200323 0.8379289104653137 1.22121933523726
98.8862386366277 0.8255216273600343 1.438670175799961
98.8550827193552 0.8137050767097959 1.672395117896858
Ich kann keinen Fehler finden, bitte angeben! Vielen Dank im Voraus!
Lösung
Keine wirkliche Programmierung Problem, das ich finden, aber ich werde auf jeden Fall beantworten.
Von einem kurzen Blick würde ich versuchen, zwei Dinge: Angenommen, Ihre Zeiteinheit Tage, im Moment scheinen Sie zu evaluieren, was die Situation für Tag 1 (korrigiert mich wenn ich falsch bin hier). Für den Fall, Sie präsentieren Ich glaube, Sie die Entwicklung über mehrere Tage wissen wollen. So haben Sie Ihre Anzahl der Schleifen zu erhöhen oder vielleicht Ihren Zeitschrittes (aber mit, dass vorsichtig sein)
Zweitens scheinen Sie hier einen Fehler zu haben: c * x * i / N ... sollte das nicht sein (c * x * i) / N? Überprüfen Sie, ob das einen Unterschied macht. Und ich glaube, Sie durch die Tatsache, zu überprüfen, dass S '+ I' + R‘sollte = 0 ...
Noch einmal, ich habe das nicht sehr tief überprüfen, sondern nur noch einen Blick und lassen Sie mich wissen, ob es irgendetwas ändert.