share/contrib/Eulix/Eulix_T2.mac

   1 load("Eulix.mac")$
   2 load("rkf45")$
   3
   4 /* This is a 1-d example,
   5    the right-hand side Rhs must return a 1x1 matrix
   6    and the initial value must be 1x1 matrix, as well
   7
   8    given a differentiable function  stat(t)  and lamba > 0
   9    this is the Robin test equation
  10
  11    y'(t)= -Lambda*(y-stat(t)) + diff(stat(t),t)   y(0)= y0
  12 */
  13
  14 stat(t):= sin(t)$
  15 Lambda: 1000$
  16 debugmode:true$
  17 y0:1$
  18 solution(t,y0):= stat(t)+exp(-Lambda*t)*(y0-stat(0))$
  19 timer(rkf45,Eulix)$
  20 print("rkf45 errors:")$
  21 RR:rkf45(-Lambda*(y-stat(t)) + diff(stat(t),t),y,1,[t,0,1])$
  22 for i:1 thru length(RR) do block([ty:RR[i],tt,yt],
  23   tt:first(ty),yt:second(ty),
  24   printf(true,"at t= ~7,4f  err= ~10,3e~%",tt,(yt-solution(tt,y0)))
  25 );
  26
  27 print("Eulix errors:")$
  28 ER:Eulix(-Lambda*(y-stat(t)) + diff(stat(t),t),y,1,[t,0,1])$
  29 for i:1 thru length(ER) do block([ty:ER[i],tt,yt],
  30   tt:first(ty),yt:second(ty),
  31   printf(true,"at t= ~7,4f  err= ~10,3e~%",tt,(yt-solution(tt,y0)))
  32 );
  33 timer_info();  /* rkf45: 0.487 sec   Eulix: 0.906 sec */