Рунге-Кутта (RK4) для системы дифференциальных уравнений в Java

Question

Это место в основном является результатом этого потока: Дифференциальные уравнения в Java.
По сути, я пытался следить за Джейсоном С. Консультировать и реализовывать числовые решения для дифференциальных уравнений через метод Runge-Kutta (RK4).

Привет всем, я пытаюсь создать простую программу моделирования модели SIR-эпидемии в Java. По сути, сэр определяется системой трех дифференциальных уравнений:
S '(T) = - Lamda (T) * S (t)
I '(t) = ламда (т) * s (t) - гамма (т) * i (t)
R '(t) = гамма (т) * i (t)
S - восприимчивые люди, я - зараженные люди, R - восстановленные люди. Ламда (t) = [c * x * i (t)] / n (t) c - количество контактов, X - инфекционность (вероятность заболевать после контакта с больным человеком), n (t) - общая популяция (который является постоянным).
Гамма (т) = 1 / продолжительность болезни (постоянная)

После первой не очень успешной попытки я пытался решить эту уравнения с Runge-Kutta, и эта попытка привела к следующему коду:

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

Это, похоже, не работает, потому что количество больных людей (I) должно сначала скорее, а затем Dicrase до примерно 0, а количество восстановленных людей должно строго увеличить. Общее количество больных + здоровых + восстановленных должно быть 100, но мой код производит некоторые странные результаты:

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  

Я не могу найти ошибку, пожалуйста, посоветуйте! Спасибо заранее!

Answer 1

Не реальная проблема программирования я нахожу, но я все равно отвечу.

Из быстрого взгляда я бы попробую две вещи: предполагая, что ваш момент времени - это дни, в данный момент, когда вы, кажется, оцениваете, какова ситуация после 1 дня (поправьте меня, если я ошибаюсь здесь). Для этого вы представляете, я думаю, вы хотите знать эволюцию в течение нескольких дней. Таким образом, вы должны увеличить количество петель или, возможно, ваш Timestep (но будьте осторожны с этим)

Во-вторых, вам, кажется, ошибаюсь здесь: C * X * I / N ... Должен ли это не быть (c * x * i) / n? Проверьте, имеет ли это значение. И я думаю, что вы можете проверить то, что S '+ I' + R 'должен = 0 ...

Еще раз я не проверил это очень глубоко, но просто посмотрите и дайте мне знать, если это меняет что-нибудь.