share/contrib/Eulix/Eulix_T5.mac

   1 load("Eulix.mac")$
   2
   3 Start:elapsed_real_time()$
   4
   5 n:1000$ p:10$
   6 rx:2.4E-5$ rv:7.4E-4$
   7 sx:0.02$ sy:0.265$
   8 ginf:0.01$
   9 vg:100$ sv:0$
  10 g(v):= ginf*v/(v+100)$
  11 f(v):= 20.0*v/(v+1.0)$
  12
  13 debugmode(true)$
  14 /* :break Eulix$ */
  15
  16 Logging:false$
  17
  18 [Solution,[tlist,ylist]]: Eulix([p - sx*x -rx*x*v +g(v)*x,
  19                                  rx*x*v -sy*y -g(v)*y,
  20                                  n*g(v)*y +f(v) -sv*v -(rx+rv)*x*v],[x,y,v],
  21                                  [1000.0,0,1.0],[t,0.0,12*365+4,1,v-t-245.0],
  22                                  initial_step_size=1,absolute_tolerance=1e-8,logging=Logging,
  23                                  combined_t_y_list=false)$
  24 print("I've used ",elapsed_real_time()-Start,"seconds for computing the AIDS model")$
  25 print("Solution triggered by condition",first(Solution)," at time ",second(Solution),
  26       "  and state:",third(Solution))$
  27 t1:second(Solution)$
  28 y1init:third(Solution)$
  29 [tlist2,ylist2]: Eulix([p - sx*x -rx*x*v +g(v)*x,
  30                         rx*x*v -sy*y -g(v)*y,
  31                         n*g(v)*y +f(v) -sv*v -(rx+rv)*x*v],
  32                         [x,y,v],y1init,[t,t1,12*365+4,1],
  33                         absolute_tolerance=1e-8,logging=Logging,
  34                         combined_t_y_list=false)$
  35
  36 tlist:append(tlist,tlist2)$
  37 ylist:append(ylist,ylist2)$
  38 tlen:length(tlist)$
  39 print("tlist:",tlist[tlen],"   ylist:",ylist[tlen])$
  40 log10(z):= if z > 0 then log(z)/log(10) else 0$
  41 tlist:tlist/365.4$  t1:t1/365.4$
  42 plot2d([[discrete,tlist,map(lambda([z],log10(z[1])),ylist)],[discrete,tlist,map(lambda([z],log10(z[2])),ylist)],
  43         [discrete,tlist,map(lambda([z],log10(z[3])),ylist)],
  44         [discrete,[t1,t1],[0,5]]],
  45        [legend,"healthy","infected","AIDS-viruses","root location"],[title,"An AIDS model"])$
  46 read("enter ; or $ to quit");
  47