java 微分方程求解_怎样用java解下面的微分方程组

dI(t)/dt=G(t)*I(t-1/30)*S(t)-0.05*I(t)………………(1)dX(t)/dt=(0.06/360000)*S(t)*(I(t)+Y(t))…………(2)dY(t)/dt=0.05*I(t)…………………………………(3)G(t)=G

dI(t)/dt=G(t)*I(t-1/30)*S(t)-0.05*I(t)………………(1)dX(t)/dt=(0.06/360000)*S(t)*(I(t)+Y(t))…………(2)dY(t)/dt=0.05*I(t)…………………………………(3)G(t)=G(0)*[(1-I(t))/360000]^3……………(4)V(t)=360000-X(t)-Y(t)=I(t)+S(t)…………………(5)其中未知量有：X(t),Y(t),I(t),G(t),S(t),V(t)I(0)=I0=1V(0)=V0=360000由公式(5)知道：S(0)=V(0)-I(0)=359999X(0)=Y(0)=0G(0)=k/Ω=358/2^32并且：当t<0时，I(t)=0要求：步长=1/300 t=0--800------------------------------------------μ=0.006/ V0，λ=0.005，G(0)=258/(2^32)我的程序如下，但是结果不对：public class func{ public static void main(String[] args) { double X[]=new double[801]; double Y[]=new double[801]; double I[]=new double[801]; double G[]=new double[801]; double S[]=new double[801]; double V[]=new double[801]; double tmpI=1; double tmpV=360000; double tmpS=359999; double tmpX=0; double tmpY=0; double tmpG=358/Math.pow(2,32); double h=1.0/300; I[0]=tmpI; V[0]=tmpV; S[0]=tmpS; X[0]=tmpX; Y[0]=tmpY; G[0]=tmpG; double tmpIt[]=new double[10]; for(int t=0;t<800;t++) { for(int i=0;i<300;i++) { tmpG=G[0]*Math.pow(1-tmpI/V[0],3); tmpV=V[0]-tmpX-tmpY; tmpS=V[0]-tmpX-tmpY-tmpI; if(t==0&i<10){ tmpIt[i]=tmpI; tmpI=tmpI-h*0.005*tmpI; } else { tmpI=tmpI+h*(tmpG*tmpIt[i%10]*tmpS-0.005*tmpI); tmpIt[i%10]=tmpI; } tmpX=tmpX+h*(0.006/V[0])*tmpS*(tmpI+tmpY); tmpY=tmpY+h*0.005*tmpI; } I[t+1]=tmpI; X[t+1]=tmpX; Y[t+1]=tmpY; G[t+1]=tmpG; V[t+1]=tmpV; S[t+1]=tmpS; } for(int i=0;i<=800;i++) { intln("t="+i+" "+"I="+I[i]); intln(" "+"X="+X[i]); intln(" "+"Y="+Y[i]); intln(" "+"S="+S[i]); intln(" "+"G="+G[i]); intln(" "+"V="+V[i]); } /*for(int i=0;i<=800;i++) { if(I[i+1]<=I[i]) intln(i); }*/ }}

展开

全部