share/contrib/Eulix/Eulix_Step_TP.mac

   1 load("Eulix.mac")$
   2
   3 /*
   4   y' =  2 * t * y * log(z)    y(0)=1
   5   z' = -2 * t * z * log(y)    z(0)=%e
   6
   7   exact solution  y(t)=exp(sin(t^2))  z(t)=exp(cos(t^2))
   8 */
   9
  10 Rhs(t,y):= matrix([ 2*t*y[1,1]*log(y[2,1])],
  11                   [-2*t*y[2,1]*log(y[1,1])] )$
  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,2)$
  23
  24 compile(Rhs,Rhs_Time)$
  25
  26 solution(t):= matrix([exp(sin(t^2))],[exp(cos(t^2))])$
  27 high_precision: true$
  28
  29 if high_precision then (
  30   fpprec:40,
  31   t:bfloat(0.1),
  32   y0:solution(t),
  33   atol:1e-30,  rtol:1e-28,
  34   h:1b-3,
  35   me:20,
  36   if true then
  37     mass_matrix:bfloat(ident(2)) /* must be either 'false or a square mass matrix
  38                                  = ident(dim) for an ODE with trivial mass matrix */
  39   else mass_matrix:false
  40
  41 ) else (
  42
  43   t:0.1,
  44   y0:solution(t),
  45   atol:1e-10,  rtol:1e-8,
  46   h:0.1,
  47   mass_matrix: ident(2), /* must be either 'false or a square mass matrix
  48                          = ident(dim) for an ODE with trivial mass matrix */
  49   me:6
  50 )$
  51
  52 y: copymatrix(y0)$
  53 t_end: sqrt(float(%pi));
  54 Check_Parm:true$
  55 Start:elapsed_real_time()$
  56 dense_out:true$
  57 debugmode:true$
  58
  59 while t < t_end do (
  60   h: min(h,t_end-t),
  61   [DO_DQ,tn,hn,me,failed]: Eulix_Step(y,t,Rhs,Rhs_Time,Rhs_Jac,h,me, dense_output=dense_out,
  62                                       absolute_tolerance=atol,relative_tolerance=rtol,
  63                                       'mass_matrix=mass_matrix),
  64   Check_Parm: false,
  65   h: tn-t,  /* <<<< in case of reject within Eulix */
  66   if failed then (
  67     printf(true,"Eulix failed at t: ~7,4f~%",t),
  68     throw('failed)
  69   ),
  70   if dense_out then (  /* dense_output=true */
  71     printf(true,"t: ~7,4f  mid interval         err: ~13,6e ~%",t+0.5*h,
  72                      mat_norm(Eulix_Dense_Out(-0.5*h,DO_DQ)-solution(t+0.5*h),'inf)),
  73
  74     y:copymatrix(DO_DQ[1])
  75   ) else y:copymatrix(DO_DQ),
  76   printf(true,"t: ~7,4f hn: ~10,8f  me:~2d err: ~13,6e ~%",tn,hn,me,
  77                      mat_norm(y-solution(tn),'inf)),
  78   t: tn,  h: hn
  79 )$
  80
  81