share/numeric/rtest_romberg.mac

   1 /* (kill (all), load (romberg), float_approx_equal_tolerance : rombergtol, 0); */
   2 (load (romberg), float_approx_equal_tolerance : rombergtol, 0);
   3 0;
   4
   5 (I : [2*x - -(log((4 + %e)/(2*%pi)))*(((4 + %e)/(2*%pi))^x), x, -1, 0],
   6  apply (romberg, I));
   7 ''(first (apply (quad_qags, I)));
   8
   9 (I : [2*x - cos((%e + %pi)*x), x, 0, 1],
  10  apply (romberg, I));
  11 ''(apply (integrate, I),
  12  ev (%%, numer, %enumer));
  13
  14 (expr : x^2 - 5,
  15  romberg (expr, x, 0, 10));
  16 ''(integrate (x^2 - 5, x, 0, 10),
  17  ev (%%, numer));
  18
  19 /* verify that romberg expression is returned for non-numeric expression or bounds
  20  */
  21 (kill (a, b), romberg (x^a - 5, x, 0, 10));
  22 romberg (x^a - 5, x, 0.0, 10.0);
  23
  24 romberg (x^3 - 5, x, a, b);
  25 romberg (x^3 - 5, x, a, b);
  26
  27 /* verify that romberg nested inside another function call is OK */
  28 (quad_qags (romberg (sin(a)*cos(x) - a*x, x, 0, 10), a, 1, 3),
  29  first (%%));
  30 ''(integrate (integrate (sin(a)*cos(x) - a*x, x, 0, 10), a, 1, 3),
  31  ev (%%, numer));
  32
  33 romberg (romberg (a^2 - x, a, 0, x) - 7, x, 0, 100);
  34 7999300.0;
  35
  36 /* verify that symbolic function name is OK */
  37 (foo (a) := 3^a - 5,
  38  bar : foo,
  39  romberg (bar, 0, 10));
  40 ''(integrate (bar (a), a, 0, 10),
  41  ev (%%, numer));
  42
  43 /* misc examples */
  44
  45 (expr : t - (297 * exp ((1000000 * t) / 33) - 330) / 10000000,
  46  romberg (expr, t, 1e-9, 0.003));
  47 ''(float ((-9801*%e^(1000/11)+9801*%e^(1/33000)+9197999933999/200000) /10000000000000));
  48
  49 (expr : 6096 * tan((2 * atan(c/(2 * fl))) / r) / (tan((1/60) * (%pi/180))),
  50  ev (romberg (expr-6096, fl, 1, 10), c=7.176, r=3264));
  51 ''(first (ev (quad_qags (expr-6096, fl, 1, 10), c=7.176, r=3264)));
  52
  53 /* more fun with evaluation */
  54
  55 (g (a) := romberg (f (x, a), x, 0, 200),
  56  f (x, a) := x^a - 5,
  57  0);
  58 0;
  59
  60 g (0.5);
  61 ''(integrate (f (x, 1/2), x, 0, 200), ev (%%, numer));
  62
  63 expr : g (z + z);
  64 romberg (x^(2 * z) - 5, x, 0.0, 200.0);
  65
  66 ''(at (expr, z=0.25));
  67 ''(integrate (f (x, 1/2), x, 0, 200), ev (%%, numer));
  68
  69 (quad_qags (g (z), z, 1, 3),
  70  first (%%));
  71 ''(assume (z > 0),
  72  integrate (f (x, z), x, 0, 200),
  73  quad_qags (%%, z, 1, 3),
  74  first (%%));
  75
  76 /* verify that symbolic constants (in integrand or limits of integration) evaluate to numbers */
  77
  78 (f(x) := cos (sin (%e) + %e^%pi + %pi^%e) * sin(x) - x/2, 0);
  79 0;
  80
  81 [romberg (cos (sin (%e) + %e^%pi + %pi^%e) * sin(x) - x/2, x, 0.1, %pi),
  82  romberg (f(x), x, 0.1, %pi),
  83  romberg (f, 0.1, %pi)];
  84 ''(integrate (f (x), x, a, b),
  85  ev (%%, a=0.1, b=%pi, numer, %enumer),
  86  [%%, %%, %%]);
  87
  88 [romberg (f, 1/(%pi*%e), 2*%pi*sin(%e)),
  89  romberg (f, log(%pi), %e^%pi),
  90  romberg (f, sin(%e), %pi^%e),
  91  romberg (f, exp(1/5), exp(cos(%e + %pi))),
  92  romberg (f, cos(exp(2))/10, 10*cos(exp(2)))];
  93 ''([integrate (f (x), x, 1/(%pi*%e), 2*%pi*sin(%e)),
  94  integrate (f (x), x, log(%pi), %e^%pi),
  95  integrate (f (x), x, sin(%e), %pi^%e),
  96  integrate (f (x), x, exp(1/5), exp(cos(%e + %pi))),
  97  integrate (f (x), x, cos(exp(2))/10, 10*cos(exp(2)))],
  98  ev (%%, numer, %enumer));
  99
 100 /* adapted from the mailing list 2007/06/10
 101  * charfun2 copied from the interpol share package
 102  */
 103
 104 /* KNOWN FAILURE: romberg computes a result but relative error is greater than rombergtol */
 105
 106 (charfun2 (z, l1, l2) := charfun (l1 <= z and z < l2),
 107  e1 : (-.329*x^3+.494*x^2 +.559*x+.117),
 108  e2 : (.215*x^3-1.94*x^2 +4.85*x-2.77),
 109  e3 : (.0933*x^3-1.02*x^2+2.56*x-.866),
 110  e4 : (.0195*x^3-.581*x^2+1.67*x-.275),
 111  e5 : (.00117*x^3-.498*x^2 +1.55*x -.213),
 112  expr : e1 *charfun2(x,minf,1.0)
 113     +e2 *charfun2(x,2.5,inf)
 114     +e3 *charfun2(x,2.0,2.5)
 115     +e4 *charfun2(x,1.5,2.0)
 116     +e5 *charfun2(x,1.0,1.5),
 117  romberg (expr, x, 0, 4));
 118 ''(ev (integrate (e1, x,0,1)
 119  + integrate (e2, x,5/2,4)
 120  + integrate (e3, x,2,5/2)
 121  + integrate (e4, x,3/2,2)
 122  + integrate (e5, x,1,3/2), keepfloat));
 123
 124 (reset(float_approx_equal_tolerance),0);
 125 0$