Runge-Kutta (RK4) for system of differential equations in Java

Question

This quation is mostly a result of this thread: Differential Equations in Java.
Basically, I've tried to follow Jason S. advise and to implement numerical solutions to differential equations via Runge-Kutta method (RK4).

Hello everyone, I am trying to create a simple simulation program of SIR-epidemics model in java. Basically, SIR is defined by a system of three differential equations:
S'(t) = - lamda(t) * S(t)
I'(t) = lamda(t) * S(t) - gamma(t) * I(t)
R'(t) = gamma(t) * I(t)
S - susceptible people, I - infected people, R - recovered people. lamda(t) = [c * x * I(t)] / N(T) c - number of contacts, x - infectiveness (probability to get sick after contact with sick person), N(t) - total population (which is constant).
gamma(t) = 1 / duration of illness (constant)

After first not very successful attempt, I've tried to solve this equations with Runge-KUtta, and this attempt resulting in the following code:

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

This does not seem to work, because number of sick people (I) should first increse, and then dicrease to about 0, and number of recovered people should strictly increase. Total number of sick + healthy + recovered should be 100, but my code produces some strange results:

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  

I cannot find a mistake, please advise! Many thanks in advance!

Answer 1

Not a real programming problem I find, but I'll answer anyway.

From a quick look I would try two things: Assuming your time unit is days, at the moment you seem to be evaluating what the situation is after day 1 (correct me if I'm wrong here). For the case you're presenting I think you want to know the evolution over several days. So you have to increase your number of loops or perhaps your timestep (but be careful with that)

Second, you seem to have a mistake here: c * x * i / N... should that not be (c*x*i)/N? Check if that makes a difference. And I think you can check by the fact that S' + I' + R' should = 0...

Once again, I did not check this very deeply, but just have a look and let me know if it changes anything.