Java中微分方程组的Runge-Kutta (RK4)
这个表达式主要是这个线程：Differential Equations in Java的结果。
基本上，我尝试遵循Jason .的建议，并通过Runge方法(RK4)实现微分方程的数值解。
大家好，我试着用java创建一个简单的SIR流行病模型模拟程序.基本上，SIR是由三个微分方程组定义的：
S'(t) =- lamda(t) * S(t)
I'(t) = lamda(t) * S(t) -γ(T)* I(t)
R'(t) =伽马(T)* I(t)
易感的人，感染的人，康复的人。lamda(t) =c*x* I(t) / N(T) c-接触次数，x-传染性(接触病人后患病的概率)，N(t) -总人数(不变)。
伽马(T)=1/病程(常数)
在进行了一次不太成功的尝试之后，我尝试用Runge来求解这个方程，这一尝试产生了以下代码：
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 我找不到错误，请指教！事先非常感谢！
回答 1
Stack Overflow用户
发布于 2010-12-05 18:53:35
我发现的不是一个真正的编程问题，但我还是会回答的。
简单地看一下，我会尝试两件事:假设你的时间单位是几天，现在你似乎在评估第一天之后的情况(如果我错了，请纠正我)。对于你提出的案件，我想你想知道几天来的演变情况。所以你必须增加你的循环数，或者你的时间步数(但是要小心)
其次，你好像搞错了:C*x*I/N.那不是(c*x*i)/N吗？看看这会不会有什么区别。我认为你可以证明S‘+ I’+ R‘应该=0.
再一次，我没有很深入地检查这个问题，只是看了一下，如果有什么变化，请告诉我。
https://stackoverflow.com/questions/4360160
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号