share/numeric/dblint.dem

   1 /* Example of use of double integrator (dblint) */
   2
   3 /* F must be a function of two variables, R and S must be functions of
   4 one variable, and A and B must be floating point numbers (all of the
   5 functions must be translated or compiled prior to using dblint) */
   6
   7 /* Get the fasl file containing dblint(F,R,S,A,B) */
   8
   9 batchload("dblint");
  10
  11 /* To get the double integral of exp(x-y^2) over the region bounded by
  12 y=1 and y=2+x and x=0 to x=1 */
  13
  14 /* Define the integrand as a function of two variables */
  15
  16 F(X,Y):=(mode_declare([X,Y],float),exp(X-Y^2));
  17
  18 /* Define the lower and upper limits on the inner (y in this case)
  19 integral as a function of the outer variable (x in this case) */
  20
  21 R(X):=(mode_declare(X,float),1.0);
  22 S(X):=(mode_declare(X,float),2.0+X);
  23
  24 /* Now translate these functions for the sake of efficiency */
  25
  26 translate(F,R,S);
  27
  28 /* Call the dblint function with quoted arguments for function names, and
  29 floating point values for the endpoints of the outer (x) integration
  30 */
  31
  32 dblint_ans:dblint('F,'R,'S,0.0,1.0);
  33
  34 /* Compare with the exact answer */
  35
  36 inty:risch(exp(X-Y^2),Y);
  37 xint:ev(inty,Y:2+X)-ev(inty,Y:1);
  38 dint:risch(xint,X);
  39 ans:ev(dint,X:1.0)-ev(dint,X:0.0),numer;
  40
  41 /* Relative error check here */
  42 (ans-dblint_ans)/ans;
  43 /* Quite reasonable */
  44 /* Of course, the dblint routine will still work even when we choose an
  45 area over which the closed-form integral fails to be expressible in
  46 terms of standard transcendental functions */
  47
  48 S1(X):=(mode_declare(X,float),2.0+X^(3/2));
  49 translate(S1);
  50 dblint('F,'R,'S1,0.0,1.0);
  51
  52 "end of dblint.dem"$