/* Example of use of double integrator (DBLINT) */

/* F must be a function of two variables, R and S must be functions of
one variable, and A and B must be floating point numbers (all of the 
functions must be translated or compiled prior to using DBLINT) */

/* Get the fasl file containing DBLINT(F,R,S,A,B) */ 

batchload("dblint"); 

/* To get the double integral of exp(x-y^2) over the region bounded by
y=1 and y=2+x and x=0 to x=1 */

/* Define the integrand as a function of two variables */

F(X,Y):=(MODE_DECLARE([X,Y],FLOAT),EXP(X-Y^2)); 

/* Define the lower and upper limits on the inner (y in this case)
integral as a function of the outer variable (x in this case) */

R(X):=(MODE_DECLARE(X,FLOAT),1.0); 
S(X):=(MODE_DECLARE(X,FLOAT),2.0+X);

/* Now translate these functions for the sake of efficiency */

TRANSLATE(F,R,S); 

/* Call the DBLINT function with quoted arguments for function names, and
floating point values for the endpoints of the outer (x) integration
*/ 

DBLINT_ANS:DBLINT('F,'R,'S,0.0,1.0); 

/* Compare with the exact answer */

INTY:RISCH(EXP(X-Y^2),Y); 
XINT:EV(INTY,Y:2+X)-EV(INTY,Y:1);
DINT:RISCH(XINT,X); 
ANS:EV(DINT,X:1.0)-EV(DINT,X:0.0),NUMER; 

/* Relative error check here */
(ANS-DBLINT_ANS)/ANS;
/* Quite reasonable */
/* Of course, the DBLINT routine will still work even when we choose an
area over which the closed-form integral fails to be expressible in
terms of standard transcendental functions */ 

S1(X):=(MODE_DECLARE(X,FLOAT),2.0+X^(3/2));
TRANSLATE(S1);
DBLINT('F,'R,'S1,0.0,1.0); 

/* To see what happens when this integral is attempted in closed-form,
you may BATCH(DBLINT,DEMO1,DSK,SHARE1), and that will also show sample use
of FLOATDEFUNK and QUANC8 */
"end of dblint.dem"$