runge-kutta（RK4）用于Java中微分方程的系统

Question

此述说主要是该线程的结果： Java中的微分方程.
基本上，我尝试遵循Jason S.建议并通过Runge-Kutta方法（RK4）实现微分方程的数值解决方案。

大家好，我正在尝试在Java中创建一个简单的Sir-Epidemics模型模拟程序。基本上，SIR由三个微分方程的系统定义：
s'（t）= -lamda（t） * s（t）
i'（t）= lamda（t） * s（t）-gamma（t） * i（t）
r'（t）= gamma（t） * i（t）
S-易感人，我 - 感染的人，r-恢复的人。 lamda（t）= [c * x * i（t）] / n（t）c-联系人数量，x-感染性（与病人接触后生病的概率），n（t） - 总人口（其中是恒定的）。
伽马（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）应该首先增加，然后将其提高到大约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天后的情况（如果我在这里错了，请问我的情况）。对于您提出的情况，我认为您想知道几天的演变。因此，您必须增加循环的数量或时间段（但要小心）

其次，您似乎在这里有一个错误：c * x * i / n ...应该不是（c * x * i） / n吗？检查是否有所作为。而且我认为您可以检查S' + i' + r'应该= 0 ...

再一次，我没有非常深入地检查一下，而只是看一看，让我知道它是否改变了任何东西。