SHARE
TWEET

Untitled

a guest Oct 6th, 2015 17 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. p1:=-1.31*^-2;
  2. p2:=-1.35*^-2;
  3. p3:=2.9*^-6;
  4. h:=136;
  5. n:=.13;
  6. a:=3.11;
  7. beta:=1;
  8. gamma:=5.36*^-6;
  9. G0:=200;
  10. X0:=0;
  11. I10:=0;
  12. I20:=0;
  13. Gb:=100;
  14. Ib:=10;
  15. tmax:=1000;
  16. func:={(p1-X)*G-p1*Gb,p2*X+p3*(I1-Ib),a*Max[0,I2]-n*I1,beta*gamma*(G-h)-n*I2};
  17. yinit:={G0,X0,I10,I20};
  18. y:={G,X,I1,I2};
  19. step:=.1
  20. 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};
  21. RungeKutta[func_List,yinit_List,y_List,step_]:=Module[{k1,k2,k3,k4},k1=step N[func/.MapThread[Rule,{y,yinit}]];
  22. k2=step N[func/.MapThread[Rule,{y,k1/2+yinit}]];
  23. k3=step N[func/.MapThread[Rule,{y,k2/2+yinit}]];
  24. k4=step N[func/.MapThread[Rule,{y,k3+yinit}]];
  25. yinit+Total[{k1,2 k2,2 k3,k4}]/6]
  26. solutions:=NestList[RungeKutta[func,yinit,y,step]&,N[yinit],Round[tmax/step]]
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top