Java の微分方程式系の Runge-Kutta (RK4)
-
08-10-2019 - |
質問
この引用は主にこのスレッドの結果です。 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
間違いが見つからないのでアドバイスお願いします!よろしくお願いします!
解決
私が見つけた本当のプログラミングの問題ではありませんが、とにかく答えます。
ざっと見たところ、次の 2 つのことを試してみます。時間単位が日であると仮定すると、現時点では 1 日目以降の状況を評価しているようです (ここで間違っている場合は修正してください)。あなたが提示しているケースでは、数日間の推移を知りたいと考えています。したがって、ループの数を増やすか、タイムステップを増やす必要があります (ただし、それには注意が必要です)。
第二に、ここで間違いがあるようです。c * x * i / N...それは (c*x*i)/N ではないでしょうか?それが違いを生むかどうかを確認してください。そして、S' + I' + R' = 0であるはずであるという事実によって確認できると思います...
もう一度言いますが、私はこれをあまり深くチェックしませんでしたが、ただ見て、何か変化があればお知らせください。