Runge-Kutta（RK4）用於Java中的微分方程組

這個定義大多是這個線程的結果：Differential Equations in Java。
基本上，我試圖按照Jason S.的建議，並通過Runge-Kutta方法（RK4）實現微分方程的數值解。Runge-Kutta（RK4）用於Java中的微分方程組

大家好，我想在java中創建一個簡單的SIR流行病模型模擬程序。（t）= lamda（t）* S（t）
基本上，SIR由三個微分方程組定義： γ-（t）* I（t）
R'（t）= gamma（t）* I（t）
S-易感人羣，I-感染人羣，R-恢復人羣。（t）= [c * x * I（t）]/N（T） c - 接觸次數，x - 感染性（與病人接觸後感染的可能性），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; 
}

這似乎並沒有工作，因爲生病的人（我）的數量應該先適當增加，然後dicrease約0，恢復人的數量應嚴格增加。我的代碼產生了一些奇怪的結果：

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

我找不到錯誤，請指教！提前謝謝了！

來源

2010-12-05 Ypsilon IV

我還沒有深入瞭解你的代碼，但我是否正確地說現在你從時間t = 0到t = 1以0.01的增量移動？如果您嘗試在較大的時間段內運行系統，同時仍採用相同大小的步驟，會發生什麼情況？（說，採取500步） – Bart 2010-12-05 18:39:08

我找不到真正的編程問題，但無論如何我都會回答。

從一個快速的樣子我會嘗試兩件事：假設你的時間單位是幾天，此刻你似乎正在評估第一天後的情況（糾正我，如果我在這裏錯了）。對於你提出的案例，我想你想知道幾天的演變。所以你必須增加你的循環次數或者你的時間步（但要小心）

其次，你似乎有一個錯誤在這裏：c * x * i/N ...應該不是（c * X * I）/ N？檢查是否有所作爲。我認爲你可以通過S'+ I'+ R'應該= 0的事實來檢查...

再一次，我沒有檢查這個很深，但只是看看，讓我知道它是否改變一切。

來源

2010-12-05 18:53:35 Bart

Runge-Kutta（RK4）用於Java中的微分方程組

回答

相關問題