De Runge-Kutta (RK4) para el sistema de ecuaciones diferenciales en Java

Question

Esta quation es principalmente un resultado de este hilo:. Ecuaciones diferenciales en Java
Básicamente, he intentado seguir a Jason S. asesorar e implementar soluciones numéricas de las ecuaciones diferenciales a través del método de Runge-Kutta (RK4).

Hola a todos, Estoy tratando de crear un programa de simulación sencilla del modelo SIR-epidemias en Java. Básicamente, SIR se define por un sistema de tres ecuaciones diferenciales:
S '(t) = - lambda (t) * S (t)
I '(t) = lambda (t) * S (t) - gamma (t) * I (t)
R '(t) = gamma (t) * I (t)
S - personas susceptibles, I - las personas infectadas, R - la gente recuperado. lambda (t) = [c * x * I (t)] / N (T) c - número de contactos, x - infecciosidad (probabilidad a enfermar después del contacto con la persona enferma), N (t) -. población total (que es constante)
gamma (t) = 1 / duración de la enfermedad (constante)

Después del primer intento no tuvo mucho éxito, he tratado de resolver este ecuaciones con Runge-Kutta, y este intento resulta en el siguiente código:

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;
}

Esto no parece funcionar, porque el número de personas enfermas (I) debe primero increse, y luego dicrease a aproximadamente 0, y el número de personas recuperadas debería aumentar de manera estricta. Número total de enfermos + + sana recuperado debe ser 100, pero mi código produce algunos resultados extraños:

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  

No puedo encontrar un error, por favor avise! Muchas gracias de antemano!

Answer 1

No es un problema de programación real que encuentro, pero voy a responder de todos modos.

De un vistazo rápido me gustaría probar dos cosas: Asumiendo que su unidad de tiempo se día, en el momento en que parece estar evaluando cuál es la situación después del día 1 (corrígeme si estoy equivocado aquí). Para el caso de que esté presentando Creo que quiere conocer la evolución a lo largo de varios días. Así que hay que aumentar el número de bucles o tal vez su paso de tiempo (pero tenga cuidado con eso)

En segundo lugar, parece que tienes un error aquí: c * x * i / N ... que no debería ser (c * x * i) / N? Compruebe si eso hace la diferencia. Y creo que se puede comprobar por el hecho de que S '+ I' + R' debería = 0 ...

Una vez más, no mira esto muy profundamente, pero sólo echar un vistazo y quiero saber si cambia nada.