Runge-Kutta (RK4) для системы дифференциальных уравнений в Java

Question

Этот результат в основном является результатом этой темы: Differential Equations in Java.
В принципе, я попытался следовать совету Джейсона С. и реализовать численные решения дифференциальных уравнений методом Рунге-Кутты (RK4).

Привет всем, Я пытаюсь создать простую симуляционную программу модели SIR-эпидемий в java. В принципе, SIR определяется системой трех дифференциальных уравнений:
S '(t) = - lamda (t) * S (t)
I' (t) = lamda (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) - общая популяция (который является постоянным).
гамма (т) = 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, а затем dicrease до 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 

Я не могу найти ошибку, пожалуйста, сообщите! Спасибо заранее!

Answer 1

Не знаю, какую проблему программирования я нахожу, но я все равно отвечу.

С быстрым взглядом я бы попробовал две вещи: Предполагая, что ваша единица времени - это дни, на данный момент вы, кажется, оцениваете, что происходит после 1-го дня (исправьте меня, если я ошибаюсь здесь). Для случая, который вы представляете, я думаю, что вы хотите знать эволюцию в течение нескольких дней. Поэтому вам нужно увеличить количество циклов или, возможно, ваш timestep (но будьте осторожны)

Во-вторых, вы, кажется, ошибаетесь здесь: c * x * i/N ... не должно быть (c * х * я)/N? Убедитесь, что это имеет значение. И я думаю, вы можете проверить, что S '+ I' + R 'должен = 0 ...

Еще раз, я не проверял это очень глубоко, но просто посмотрю и дайте мне знать, если это ничего не меняет.