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中的微分方程组

回答

相关问题