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

2020年9月27日 22点热度 0条评论

该引用主要是以下线程的结果:Differential Equations in Java
基本上,我尝试遵循Jason S.的建议,并通过Runge-Kutta方法(RK4)对微分方程实施数值解。

大家好,
我正在尝试用Java创建一个简单的SIR流行病模型仿真程序。
基本上,SIR由三个微分方程式的系统定义:
S'(t)=-lamda(t)* S(t)
I'(t)= lamda(t)* S(t)-γ(t)* I(t)
R'(t)=伽玛(t)* I(t)
S-易感人群,I-感染者,R-康复者。
λ(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  

我找不到错误,请指教!提前谢谢了!

解决方案如下:

我发现不是一个真正的编程问题,但是我还是会回答。

快速浏览一下,我将尝试两件事:
假设您的时间单位是天,此刻您似乎正在评估第一天后的情况(如果我错了,请纠正我)。对于您要介绍的情况,我想您想了解几天内的演变。因此,您必须增加循环数或时间步长(但要小心一点)

其次,您似乎在这里有一个错误:c * x * i / N ...不应该是(c * x * i)/ N吗?检查是否有区别。我认为您可以通过以下事实进行检验:S'+ I'+ R'= 0 ...

再一次,我没有对此进行非常深入的检查,只是看一看,让我知道它是否有任何改变。