Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- p1:=-1.31*^-2;
- p2:=-1.35*^-2;
- p3:=2.9*^-6;
- h:=136;
- n:=.13;
- a:=3.11;
- beta:=1;
- gamma:=5.36*^-6;
- G0:=200;
- X0:=0;
- I10:=0;
- I20:=0;
- Gb:=100;
- Ib:=10;
- tmax:=1000;
- func:={(p1-X)*G-p1*Gb,p2*X+p3*(I1-Ib),a*Max[0,I2]-n*I1,beta*gamma*(G-h)-n*I2};
- yinit:={G0,X0,I10,I20};
- y:={G,X,I1,I2};
- step:=.1
- NDvec:={G'[t]==(p1-X[t])*G[t]-p1*Gb,X'[t]==p2*X[t]+p3*(I1[t]-Ib),I1'[t]==a*Max[0,I2[t]]-n*I1[t],I2'[t]==beta*gamma*(G[t]-h)-n*I2[t],G[0]==G0,X[0]==X0,I1[0]==I10,I2[0]==I20};
- RungeKutta[func_List,yinit_List,y_List,step_]:=Module[{k1,k2,k3,k4},k1=step N[func/.MapThread[Rule,{y,yinit}]];
- k2=step N[func/.MapThread[Rule,{y,k1/2+yinit}]];
- k3=step N[func/.MapThread[Rule,{y,k2/2+yinit}]];
- k4=step N[func/.MapThread[Rule,{y,k3+yinit}]];
- yinit+Total[{k1,2 k2,2 k3,k4}]/6]
- solutions:=NestList[RungeKutta[func,yinit,y,step]&,N[yinit],Round[tmax/step]]
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement