share/contrib/Eulix/Eulix_Table_T5.mac

   1 load("Eulix.mac")$
   2
   3 /* An AIDS model with root finding
   4
   5   x' = p - sx*x -rx*x*v +g(v)*x            x(0)= 1000    x: healthy T-cells   per mm^3
   6   y' = rx*x*v -sy*y -g(v)*y                y(0)= 0       y: infected T-cells  per mm^3
   7   v' = n*g(v)*y +f(y) -sv*v -(rx+rv)*x*v   v(0)= 1       v: AIDS viruses      per mm^3
   8
   9   with  f(v):= 20*v/(v+1)  and  g(v):= v*g_inf/(vg+v)
  10
  11   n     : per time unit, g(v)*y infected T-cells die and produce n*g(v)*y
  12           new AIDS viruses [1000]
  13   p     : production of healthy T-cells by the thymus [10]
  14   rx    : there are rx*v new infected T-celss  [2.4E-5]
  15   rv    : rv*x*v viruses are kill by the immune system [7.4E-4]
  16   sx    : death rate of healthy T-cells [0.02]
  17   sy    : (natural) death rate of infected T-cells [0.265]
  18   sv    : (natural) death rate of the AIDS viruses [0]
  19   f(v)  : reproduction rate of viruses by non-T-celss
  20   g(v)  : cell multiplication rate of T-cells for fighting the AIDS virus
  21   g_inf : maximum cell multiplication rate of the AIDS viruses [0.01]
  22   vg    : number of viruses with a cell multiplication rate of half of
  23           the maximum cell multiplication rate [100]
  24   x == y[1] y == y[2]  v == y[3]
  25 */
  26
  27 g(v):= v/(100*(v+100))$
  28 f(v):= 20.0*v/(v+1.0)$
  29
  30 Rhs(t,y):= matrix([10.0-0.02*y[1,1]-2.4E-5*y[1,1]*y[3,1]+g(y[3,1])*y[1,1]],
  31                   [2.4E-5*y[1,1]*y[3,1]-0.265*y[2,1]-g(y[3,1])*y[2,1]],
  32                   [1000.0*g(y[3,1])*y[2,1]+f(y[3,1])-7.64E-4*y[1,1]*y[3,1]])$
  33
  34 mass_matrix: ident(3)$
  35
  36
  37 define(Rhs_t('t,'y),diff(Rhs('t,'y),'t))$
  38
  39 gen_jacobian(F,xx,Fdim)::= block([i,n:ev(Fdim),J,X:ev(xx)],
  40   local(_y,mynumer),
  41   mynumer: if fpprec <= 16 then 'float else 'bfloat, declare(mynumer,evfun),
  42   J: genmatrix(lambda([i,j],ev(diff(F(X,_y)[i,1],_y[j,1]),diff,mynumer)),n,n),
  43    buildq([J,t:X],lambda([t,_y],J))
  44 )$
  45
  46 Rhs_jac:gen_jacobian(Rhs,t,length(Rhs('t,'y)))$
  47
  48 root_1(x,y):= block(local(y),y[3,1]-3000)$
  49 define(root_1_y(x,y),transpose(matrix([0,0,1])))$
  50 root_1_x(x,y):= 0$
  51
  52 compile(Rhs,Rhs_t)$
  53 t0:0.0$
  54 y0:transpose(matrix([1000.0,0.0,1.0]))$
  55 t_end: 12*365+4$  /* 12 years */
  56 delta_t:7$
  57 Roots:[[root_1,root_1_y,root_1_x]]$
  58
  59 Start:elapsed_real_time()$
  60 debugmode:true$
  61
  62 Logging:false$
  63
  64 [Root_at,[tlist,ylist]]:Eulix_Table(t0,y0,delta_t,t_end,Rhs,Rhs_t,Rhs_jac,Roots,
  65                                      'mass_matrix=mass_matrix,check_parameters=true,
  66                                      combined_t_y_list=false,logging=Logging)$
  67 if first(Root_at) > 0 then (
  68   print("I've used ",elapsed_real_time()-Start,"seconds for computing the AIDS model"),
  69   print("Solution triggered by condition",first(Root_at)," at time ",second(Root_at),
  70         "  and state:",third(Root_at)),
  71   t1:second(Root_at),
  72   y1init:third(Root_at),
  73   [tlist2,ylist2]: Eulix_Table(t1,y1init,delta_t,t_end,Rhs,Rhs_t,Rhs_jac,[],
  74                                logging=true,combined_t_y_list=false),
  75
  76   tlist:append(tlist,tlist2),
  77   ylist:append(ylist,ylist2),
  78   tlist:tlist/365.25,  t1:t1/365.25,
  79   cut(z):=max(z,0.1),
  80   plot2d([[discrete,tlist,map(lambda([z],cut(z[1,1])),ylist)],
  81           [discrete,tlist,map(lambda([z],cut(z[2,1])),ylist)],
  82           [discrete,tlist,map(lambda([z],cut(z[3,1])),ylist)],
  83           [discrete,[t1,t1],[0.1,1E5]]], [logy,true],
  84          [legend,"healthy","infected","AIDS-viruses","root location"],[title,"An AIDS model"]),
  85   read("enter ; or $ to quit")
  86 ) else print("no root condition has been triggered")$
  87
  88
  89