Runge-Kutta (RK4) für System von Differentialgleichungen in Java

Question

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

zu implementieren

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!

Answer 1

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.