Runge-Kutta (RK4) per il sistema di equazioni differenziali in Java

Question

Questa quation è in gran parte il risultato di questa discussione:. Differential Equations in Java
Fondamentalmente, ho cercato di seguire Jason S. consulenza e di implementare soluzioni numeriche di equazioni differenziali tramite metodo di Runge-Kutta (RK4).

  

Ciao a tutti,   Sto cercando di creare un semplice programma di simulazione del SIR-epidemie modello in java.   Fondamentalmente, SIR è definito da un sistema di tre equazioni differenziali:
  S '(t) = - Lamda (t) * S (t)
  I '(t) = Lamda (t) * S (t) - gamma (t) * I (t)
  R '(t) = gamma (t) * I (t)
  S - persone sensibili, I - persone infette, R - recuperato persone.   Lamda (t) = [c * x * I (t)] / N (T)   c - numero di contatti, x - infettività (probabilità di ammalarsi dopo il contatto con la persona malata), N (t) -. popolazione totale (che è costante)
  gamma (t) = 1 / durata della malattia (costante)

Dopo il tentativo prima non molto successo, ho cercato di risolvere questo equazioni con Runge-Kutta, e questo tentativo con conseguente il seguente codice:

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

Questo non sembra funzionare, perché il numero di malati (I) deve prima increse, e quindi ad allontanare a circa 0, e il numero di persone recuperati dovrebbe rigorosamente aumentare. Numero totale di malati + sana + recuperato dovrebbe essere 100, ma il mio codice produce alcuni risultati strani:

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  

Non riesco a trovare un errore, si prega di avvisare! Molte grazie in anticipo!

Answer 1

Non è un vero problema di programmazione lo trovo, ma io a rispondere in ogni caso.

Da una rapida occhiata vorrei provare due cose: Assumendo che l'unità di tempo è giorni, in questo momento ti sembra di essere valutando qual è la situazione dopo il giorno 1 (mi corregga se sbaglio qui). Per il caso in cui si sta presentando Penso che si desidera conoscere l'evoluzione di diversi giorni. Quindi bisogna aumentare il numero di cicli o magari il passo temporale (ma attenzione a quello)

In secondo luogo, ti sembra di avere un errore qui: c * x * i / N ... Se questo non si (c * x * i) / N? Verificare se questo fa la differenza. E penso che è possibile controllare con il fatto che S '+ I' + R' dovrebbe = 0 ...

Ancora una volta, non ho controllare questo molto profondamente, ma basta dare un'occhiata e fatemi sapere se cambia qualcosa.