1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79

/* Numerical Integration Demo for QUANC8*/
/* SETUP_AUTOLOAD(QQ,QUANC8); */
/* show the time and the garbage collection time, and use the
GCdemon to make better use of flonum space */
SHOWTIME:'ALL;
IF STATUS(FEATURE,ITS)=TRUE THEN LOAD("GCDEMN");
/* Sample highly oscillatory problem */
G(X):=(MODE_DECLARE(X,FLOAT),SIN(1.0/X));
/* Inefficient way to use QUANC8 */
QUANC8(G(x),x,0.01,0.1);
/* give G as the thing to compile, then a semicolon to read it in */
compile();
/* see how much faster now, using fast version of QUANC8 */
QUANC8('G,0.01,0.1);
/* note that romberg doesn't easily manage oscillatory behaviors */
ERRCATCH(ROMBERG('G,0.01,0.1));
/* a function with two large and narrow peaks,
so hard to integrate with a fixed stepsize */
H(X):=(MODE_DECLARE(X,FLOAT),1.0/(.0001+(X.3)^2)+1.0/(.0025+(X.9)^2));
/* give H and a semicolon in response */
compile();
/* the exact answer for integrate(h(x),x) is
100.0*atan(100.0*x30.0)+20.0*atan(20.0*x18.0) */
/* the exact answer for integrate(h(x),x,0.0,1.0) is 361.847622.
note that the romberg result is more accurate in this case */
ROMBERG('H,0.0,1.0);
/* but at a large time cost compared to adaptive */
QUANC8('H,0.0,1.0);
/* reduce QUANC8's relative error tolerance by a factor of 10 */
QUANC8_RELERR:1.0E5;
/* does not take much longer, and gives much better result now */
QUANC8('H,0.0,1.0);
/* put QUANC8_RELERR back to 1.0e4 */
QUANC8_RELERR:1.0E4$
/* the exact answer for integrate(h(x),x,0.0,10.0) is 372.33606.
while too tough for fixed step size method */
ERRCATCH(ROMBERG('H,0.0,10.0));
/*QUANC8 is still super speedy even now */
QUANC8('H,0.0,10.0);
/* double integral of sample function exp(x)*sin(y)
from x=0.0 to 2.0, y=0.0 to 1.0 */
/* you must be sure that the quanc8 source is loaded
(load("qq");) before you translate a call to quanc8 in a function,
else an error in macro expansion will occur when it can't figure
out if the 3 or 4 arg version is used */
/* the exact answer for
integrate(integrate(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0) is 2.93703434.
integrate over y to get the final answer */
QUANC8(quanc8(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0);
p():=QUANC8(lambda([y],mode_declare(y,float),
quanc8(lambda([x],mode_declare(x,float),
exp(x)*sin(y)),
0.0,2.0)),0.0,1.0);
p();
translate(p);
p();
/* give p and a semicolon */
compile();
p();
