Runge-Kutta (RK4) for system of differential equations in Java
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.967398开发者_开发知识库89250752 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!
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.
精彩评论