share/contrib/Eulix/Eulix_Step_T_Stiff.mac

   1 load("Eulix.mac")$
   2
   3 /*
   4   y1' = -0.04*y1 + 1E4*y2*y3             y1(0)= 1
   5   y2' =  0.04*y1 - 1E4*y2*y3 -3E7*y2^2   y2(0)= 0
   6   y3' =                       3E7*y2^2   y3(0)= 0
   7 */
   8
   9 Rhs(t,y):= matrix([ -0.04*y[1,1] + 1E4*y[2,1]*y[3,1]],
  10                   [  0.04*y[1,1] - 1E4*y[2,1]*y[3,1] -3E7*y[2,1]^2],
  11                   [                       3E7*y[2,1]^2])$
  12
  13 define(Rhs_Time('t,'y),diff(Rhs('t,'y),'t))$
  14
  15 gen_jacobian(F,xx,Fdim)::= block([i,n:ev(Fdim),J,X:ev(xx)],
  16   local(_y,mynumer),
  17   mynumer: if fpprec <= 16 then 'float else 'bfloat, declare(mynumer,evfun),
  18   J: genmatrix(lambda([i,j],ev(diff(F(X,_y)[i,1],_y[j,1]),diff,mynumer)),n,n),
  19    buildq([J,t:X],lambda([t,_y],J))
  20 )$
  21
  22 Rhs_Jac:gen_jacobian(Rhs,t,3)$
  23
  24 compile(Rhs,Rhs_Time)$
  25
  26
  27 t:0$
  28 y0:matrix([1],[0],[0])$
  29 atol:1e-6$  rtol:1e-4$
  30 h:0.1$
  31 mass_matrix: ident(3)$ /* must be either 'false or a square mass matrix
  32                        = ident(dim) for an ODE with trivial mass matrix */
  33 me:6$
  34
  35 y: copymatrix(y0)$
  36 t_end: 0.3;
  37 Check_Parm:true$
  38 Start:elapsed_real_time()$
  39 dense_out:true$
  40 debugmode:true$
  41
  42 while t < t_end do (
  43   h: min(h,t_end-t),
  44   [DO_DQ,tn,hn,me,failed]: Eulix_Step(y,t,Rhs,Rhs_Time,Rhs_Jac,h,me, dense_output=dense_out,
  45                                       absolute_tolerance=atol,relative_tolerance=rtol,
  46                                       'mass_matrix=mass_matrix,check_parameters=Check_Parm,
  47                                       logging=true),
  48   Check_Parm: false,
  49   h: tn-t,  /* <<<< in case of reject within Eulix */
  50   if failed then (
  51     printf(true,"Eulix failed at t: ~7,4f~%",t),
  52     throw('failed)
  53   ),
  54   if dense_out then (  /* dense_output=true */
  55     printf(true,"t: ~7,4f hn: ~10,8f  me:~2d y1=~13,6e  y2=~13,6e y3=~13,6e~%",
  56                  tn,      hn,     me,  y[1,1],    y[2,1],   y[3,1]),
  57     y:copymatrix(DO_DQ[1])
  58   ) else y:copymatrix(DO_DQ),
  59   printf(true,"t: ~7,4f hn: ~10,8f  me:~2d y1=~13,6e  y2=~13,6e y3=~13,6e~%",
  60                tn,      hn,     me,  y[1,1],    y[2,1],   y[3,1]),
  61   t: tn,  h: hn
  62 )$
  63
  64