De Runge-Kutta (RK4) para el sistema de ecuaciones diferenciales en Java
-
08-10-2019 - |
Pregunta
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!
Solución
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.