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
|
/* 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);
"end of dblint.dem"$
|
