Runge-Kutta (RK4) pour système d'équations différentielles en Java

Question

Cette quation est principalement le résultat de ce fil. équations différentielles en Java
En fait, je l'ai essayé de suivre Jason S. conseiller et de mettre en œuvre des solutions numériques aux équations différentielles par méthode de Runge-Kutta (RK4).

  

Bonjour tout le monde,   Je suis en train de créer un programme de simulation simple du modèle SIR-épidémies en java.   Fondamentalement, SIR est défini par un système de trois équations différentielles:
  S '(t) = - lamda (t) * S (t)
  I '(t) = lamda (t) * S (t) - gamma (t) * I (t)
  R '(t) = gamma (t) * I (t)
  S - personnes sensibles, I - personnes infectées, R - personnes récupérées.   lamda (t) = [c * x * I (t)] / N (T)   c - nombre de contacts, x - infectiosité (probabilité de tomber malade après un contact avec une personne malade), N (t) -. population totale (qui est constante)
  gamma (t) = 1 / durée de la maladie (constante)

Après la première tentative pas très réussie, j'ai essayé de résoudre cette équation avec Runge-Kutta, et cette tentative entraînant le code suivant:

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

Cela ne semble pas fonctionner, car le nombre de personnes malades (I) doit d'abord increse, puis dicrease à environ 0, et le nombre de personnes récupérées devrait strictement augmenter. Nombre total de malades + santé + récupéré devrait être 100, mais mon code produit des résultats étranges:

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  

Je ne peux pas trouver une erreur, s'il vous plaît conseiller! Un grand merci à l'avance!

Answer 1

Pas un vrai problème de programmation, je trouve, mais je vais y répondre de toute façon.

D'un rapide coup d'oeil, je voudrais essayer deux choses: En supposant que votre unité de temps est le jour au moment où vous semblez évaluer quelle est la situation après jour 1 (corrigez-moi si je me trompe). Pour le cas que vous présentez, je pense que vous voulez connaître l'évolution sur plusieurs jours. Donc, vous devez augmenter votre nombre de boucles ou peut-être votre timestep (mais attention à cela)

Deuxièmement, vous semblez avoir une erreur ici: c * x * i / N ... si cela ne pas être (c * x * i) / N? Vérifiez si cela fait une différence. Et je pense que vous pouvez vérifier par le fait que S '+ I' + R »= 0 devrait ...

Encore une fois, je ne suis pas le vérifier très profondément, mais juste un coup d'oeil et laissez-moi savoir si ça change quoi que ce soit.