Рунге-Кутта (RK4) для системы дифференциальных уравнений в Java - PullRequest
0 голосов
/ 05 декабря 2010

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

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

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

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) должно сначала увеличиваться, а затем уменьшаться примерно до 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  

Не могу найти ошибку, пожалуйста, сообщите! Большое спасибо заранее!

1 Ответ

1 голос
/ 05 декабря 2010

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

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

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

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

...