archive/share/trash/qq.dem

   1
   2 /* Numerical Integration Demo for QUANC8*/
   3
   4 /* SETUP_AUTOLOAD(QQ,QUANC8); */
   5
   6 /* show the time and the garbage collection time, and use the
   7    GC-demon to make better use of flonum space */
   8
   9 showtime:'all;
  10 /* IF STATUS(FEATURE,ITS)=TRUE THEN LOAD("GCDEMN"); */
  11
  12 /* Sample highly oscillatory problem */
  13
  14 g(x):=(mode_declare(x,float),sin(1.0/x));
  15
  16 /* inefficient way to use QUANC8 */
  17
  18 quanc8(g(x),x,0.01,0.1);
  19
  20 /* give G as the thing to compile, then a semi-colon to read it in */
  21 compile('g);
  22
  23 /* see how much faster now, using fast version of QUANC8 */
  24
  25 quanc8('g,0.01,0.1);
  26
  27 /* note that romberg doesn't easily manage oscillatory behaviors */
  28
  29 errcatch(romberg('g,0.01,0.1));
  30
  31 /* a function with two large and narrow peaks,
  32    so hard to integrate with a fixed step-size */
  33
  34 h(x):=(mode_declare(x,float),1.0/(.0001+(x-.3)^2)+1.0/(.0025+(x-.9)^2));
  35 /* give H and a semi-colon in response */
  36 compile('h);
  37
  38 /* the exact answer for integrate(h(x),x) is
  39    100.0*atan(100.0*x-30.0)+20.0*atan(20.0*x-18.0) */
  40 /* the exact answer for integrate(h(x),x,0.0,1.0) is 361.847622.
  41    note that the romberg result is more accurate in this case */
  42 romberg('h,0.0,1.0);
  43 /* but at a large time cost compared to adaptive */
  44 quanc8('h,0.0,1.0);
  45 /* reduce QUANC8's relative error tolerance by a factor of 10 */
  46 quanc8_relerr:1.0E-5;
  47 /* does not take much longer, and gives much better result now */
  48 quanc8('h,0.0,1.0);
  49 /* put QUANC8_RELERR back to 1.0e-4 */
  50 quanc8_relerr:1.0E-4$
  51 /* the exact answer for integrate(h(x),x,0.0,10.0) is 372.33606.
  52    while too tough for fixed step size method */
  53 errcatch(romberg('h,0.0,10.0));
  54 /*QUANC8 is still super speedy even now */
  55 quanc8('h,0.0,10.0);
  56
  57 /* double integral of sample function exp(x)*sin(y)
  58    from x=0.0 to 2.0, y=0.0 to 1.0 */
  59 /* you must be sure that the quanc8 source is loaded
  60    (load("qq");) before you  translate a call to quanc8 in a function,
  61    else an error in macro expansion will occur when it can't figure
  62    out if the 3 or 4 arg version is used */
  63
  64 /* the exact answer for
  65    integrate(integrate(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0) is 2.93703434.
  66    integrate over y to get the final answer */
  67 quanc8(quanc8(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0);
  68
  69 p():=quanc8(lambda([y],mode_declare(y,float),
  70                    quanc8(lambda([x],mode_declare(x,float),
  71                                  exp(x)*sin(y)),
  72                           0.0,2.0)),0.0,1.0);
  73
  74 p();
  75 translate(p);
  76 p();
  77 /* give p and a semi-colon */
  78 compile('p);
  79 p();