Java の微分方程式系の Runge-Kutta (RK4)

Question

この引用は主にこのスレッドの結果です。 Java の微分方程式.
基本的に、私は Jason S をフォローしようとしました。ルンゲ・クッタ法 (RK4) による微分方程式の数値解法をアドバイスし、実装します。

みなさん、こんにちは。私はJavaでSir-Epidemicsモデルのシンプルなシミュレーションプログラムを作成しようとしています。基本的に、SIR は 3 つの微分方程式系によって定義されます。
S'(t) = - ラムダ(t) * S(t)
I'(t) = ラムダ(t) * S(t) - ガンマ(t) * I(t)
R'(t) = ガンマ(t) * I(t)
S - 感染しやすい人、I - 感染した人、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

私が見つけた本当のプログラミングの問題ではありませんが、とにかく答えます。

ざっと見たところ、次の 2 つのことを試してみます。時間単位が日であると仮定すると、現時点では 1 日目以降の状況を評価しているようです (ここで間違っている場合は修正してください)。あなたが提示しているケースでは、数日間の推移を知りたいと考えています。したがって、ループの数を増やすか、タイムステップを増やす必要があります (ただし、それには注意が必要です)。

第二に、ここで間違いがあるようです。c * x * i / N...それは (c*x*i)/N ではないでしょうか?それが違いを生むかどうかを確認してください。そして、S' + I' + R' = 0であるはずであるという事実によって確認できると思います...

もう一度言いますが、私はこれをあまり深くチェックしませんでしたが、ただ見て、何か変化があればお知らせください。