share/contrib/Eulix/Eulix_T8.mac

   1 load("Eulix.mac")$
   2
   3
   4 A:-(U^2+V^2)*sinh(u)/(cosh(u)^3-cosh(u));
   5 B:2*U*V*sinh(u)/cosh(u);
   6
   7 /*
   8 dA_du:(-V^2-U^2)*cosh(u)-(-V^2-U^2)*sinh(u)*(3*cosh(u)^2*sinh(u)-sinh(u))/(cosh(u)^3-cosh(u))^2;
   9 dA_du:-(U^2+V^2)*(cosh(u)-sinh(u)*(3*cosh(u)^2*sinh(u)-sinh(u))/(cosh(u)^3-cosh(u))^2);
  10 dB_du:2*U*V*(1-sinh(u)^2/cosh(u)^2);
  11 */
  12
  13
  14
  15 odes:[U,V,A,B]$
  16 funcs:[u,v,U,V]$
  17
  18 /* stop if  u-3e-2 = 0  OR  u+3e-2 = 0 */
  19
  20 interval:[s,0,0.4,0.002,[u-3e-2,u+3e-2]]$
  21 /* this doesn't find the root : interval:[s,0,0.4,0.01,[u-3e-2,u+3e-2]] */
  22
  23 initial:[1,0,20,5];
  24
  25 /* debug_mode(true); */
  26
  27
  28 [Root_at,eusol]:Eulix(odes,funcs,initial,interval,
  29                 absolute_tolerance=1E-8,combined_t_y_list=false,
  30                 check_parameters=true)$
  31 if first(Root_at) > 0 then
  32   print("root at ",second(Root_at),"   triggered by condition ",first(Root_at),
  33         "root solution=
  34 ",third(Root_at))
  35 else (
  36   print("no root condition was triggered"),
  37   [tlist,ylist]:eusol,
  38   print("Final Solution=
  39   ", ylist[length(tlist)])
  40 );