ルンゲクッタ（RK4）は、Java

Question

における微分方程式のシステムのためのこのquationは、主にこのスレッドの結果である：Differential Equations in Java。
基本的には、Jason S. adviseに従って、Runge-Kuttaメソッド（RK4）を使って微分方程式の数値解を実装しようとしました。

こんにちは皆さん、 私はjavaでSIR-流行モデルのシンプルなシミュレーションプログラムを作成しようとしています。 基本的に、SIRは、3つの微分方程式のシステムによって定義される：
S '（t）= - lamda（t）* S（t）
I'（t）= lamda（t）* Sガンマ（トン）* I（t）は
R '（T）=γ（t）は* I（t）は
S - 影響を受けやすい人々、I - 感染した人々、Rは - 人々を回復しました。 ラムダ（T）= C *第X * I（T）]/N（T） C - コンタクトの数、X - infectiveness（病人との接触後に病気する確率）、N（T） - 総人口（これは一定である）。
γ（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）の数は最初increse、次いで約0にdicrease、回収した人の数は、厳密に増加するはずである必要がありますので、これは、動作するようには思えません。回復病気+健康+の総数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 ... * x * i）/ N？違いがあるかどうかを確認してください。そして、私はあなたがS '+ I' + R '= 0であることを確認することができると思います。

もう一度、私はこれを非常に深くチェックしませんでしたが、何かを変える。