share/numeric/dblint_1.dem

   1 /* a macro for floatdefunk, needed until optimu is used */
   2 load("fltdfnk.mc")$
   3
   4 /* Example of use of double integrator (DBLINT), part 2 */
   5
   6 /* Get the fasl file for DBLINT */
   7 load('dblint);
   8
   9 /* Use DBLINT to get the double integral of exp(x-y^2) over the region
  10 bounded by y=1 and y=2+x^(3/2) and x=0 to x=1 */
  11
  12 /* Define the integrand as a function of two variables */
  13
  14 f(x,y):=(mode_declare([x,y],float),exp(x-y^2));
  15
  16 /* Define the lower and upper limits on the inner (y in this case)
  17 integral as a function of the outer variable (x in this case) */
  18
  19 r(x):=(mode_declare(x,float),1.0);
  20 s1(x):=(mode_declare(x,float),2.0+x^(3/2));
  21
  22 /* Now translate these functions for the sake of efficiency */
  23
  24 translate(f,r,s1);
  25
  26 /* Call the DBLINT function with quoted arguments for function names, and
  27 floating point values for the endpoints of the outer (x) integration
  28 */
  29
  30 dblint_answer:dblint('f,'r,'s1,0.0,1.0);
  31
  32 /* Now generate the exact integral over y using RISCH */
  33
  34 inty:risch(exp(x-y^2),y);
  35
  36 /* Now get the integrand of the x integral */
  37
  38 xint:ev(inty,y:2+x^(3/2))-ev(inty,y:1);
  39
  40 /* Try to do the x integral exactly */
  41
  42 risch(xint,x);
  43 radcan(%);
  44 %,nouns;
  45 /* Still no luck with closed-form */
  46 "end of dblint.dm1"$