0

この引用は主にこのスレッドの結果です：Javaの微分方程式。
基本的に、私はJason S.のアドバイスに従い、ルンゲクッタ法（RK4）を介して微分方程式の数値解法を実装しようとしました。

みなさん、こんにちは。JavaでSIRエピデミックモデルの簡単なシミュレーションプログラムを作成しようとしています。基本的に、SIRは次の3つの微分方程式のシステムによって定義されます
。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 /病気の期間（一定）

最初の試みはあまり成功しませんでしたが、ルンゲクッタ法でこの方程式を解こうとしましたが、この試みの結果、次のコードが生成されました。 

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

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

簡単に見てみると、2つのことを試してみます。時間単位が日であると仮定すると、現時点では、1日目以降の状況を評価しているようです（ここで間違っている場合は訂正してください）。あなたが提示している場合、私はあなたが数日にわたる進化を知りたいと思うと思います。したがって、ループの数またはおそらくタイムステップを増やす必要があります（ただし、それに注意してください）

第二に、あなたはここで間違いをしているようです：c * x * i / N ...それは（c * x * i）/ Nであるべきではありませんか？それが違いを生むかどうかを確認してください。そして、S'+I'+R'が=0でなければならないという事実によって確認できると思います。

繰り返しになりますが、私はこれをあまり深くチェックしませんでしたが、それが何かを変えるかどうかを見て、私に知らせてください。

于 2010-12-05T18:53:35.133 に答える