java - Java 中微分方程组的 Runge-Kutta (RK4)

Question

这个公式主要是这个线程的结果:Differential Equations in Java .
基本上，我尝试遵循 Jason S. 的建议，并通过 Runge-Kutta 方法 (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)

在第一次不太成功的尝试之后，我尝试使用 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)应该首先增加，然后减少到大约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

这不是我发现的真正的编程问题，但无论如何我都会回答。

快速浏览一下，我会尝试两件事:假设你的时间单位是天，目前你似乎正在评估第一天之后的情况(如果我错了，请纠正我)。对于您所介绍的案例，我认为您想了解几天内的演变情况。所以你必须增加循环次数或者时间步长(但要小心)

其次，你似乎在这里犯了一个错误:c * x * i/N...不应该是(c*x*i)/N吗？检查这是否有影响。我认为你可以通过以下事实进行检查:S' + I' + R' 应该 = 0...

再一次，我没有深入检查这一点，只是看一下，如果它改变了什么，请告诉我。