share/odepack/rtest_dlsode.mac

   1 (load("dlsode"),0);
   2 0$
   3 (closeto(value,compare,tol):=
   4   block(
   5     [abse],
   6     abse:abs(value-compare),if(abse<tol) then true else abse),
   7     done);
   8 done$
   9 y: [1.0, 0.0, 0.0];
  10 [1.0,0.0,0.0]$
  11 t: 0.0;
  12 0.0$
  13 rtol : 1e-4;
  14 1.0*10^-4$
  15 atol: [1e-6, 1e-10, 1e-6];
  16 [1.0*10^-6,1.0*10^-10,1.0*10^-6]$
  17 istate: 1;
  18 1$
  19 mf: 21;
  20 21$
  21 f1: -.04*y1 + 1e4*y2*y3;
  22 10000.0*y2*y3-0.04*y1$
  23 f3: 3e7*y2*y2;
  24 3.0*10^7*y2^2$
  25 fex: [f1, -f1-f3, f3];
  26 [10000.0*y2*y3-0.04*y1,-10000.0*y2*y3-3.0*10^7*y2^2+0.04*y1,3.0*10^7*y2^2]$
  27 t:0.0;
  28 0.0$
  29 tout:.4;
  30 0.4$
  31 result :[];
  32 []$
  33 (state : dlsode_init(fex, ['t,y1,y2,y3], mf),0);
  34 0$
  35 for k : 1 thru 12 do
  36   block([],
  37     result: dlsode_step(y, t, tout, rtol, atol, istate, state),
  38     printf(true, "At t = ~12,4,2e   y = ~{~14,6,2e~}~%", result[1], result[2]),
  39     istate : result[3],
  40     tout : tout * 10);
  41 done$
  42 closeto(result[1],4.0*10^10,1e-10);
  43 true$
  44 closeto(result[2][1],-7.170560923964532*10^-8,1e-10);
  45 true$
  46 closeto(result[2][2],-2.868223874302304*10^-13,1e-10);
  47 true$
  48 closeto(result[2][3],1.000000071705913,1e-10);
  49 true$
  50 closeto(result[3],2,1e-10);
  51 true$
  52 (kill(closeto),0)$
  53 0$