0

这个方程主要是这个线程的结果：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 - 易感人群，I - 感染者，R - 康复者。lamda(t) = [c * x * I(t)] / N(T) c - 接触人数，x - 传染性（与病人接触后生病的概率），N(t) - 总人口（其中是恒定的）。
gamma(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

找不到错误，请指教！提前谢谢了！

4

1 回答 1

1

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

快速浏览一下，我会尝试两件事：假设您的时间单位是天，此时您似乎正在评估第 1 天之后的情况（如果我在这里错了，请纠正我）。对于你提出的案例，我想你想知道几天的演变。所以你必须增加你的循环次数或者你的时间步长（但要小心）

其次，您这里似乎有一个错误：c * x * i / N ...不应该是（c * x * i）/ N吗？检查这是否有所作为。而且我认为您可以通过以下事实来检查 S' + I' + R' 应该 = 0 ...

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

于 2010-12-05T18:53:35.133 回答