← Back to branch summary

~ubuntu-branches/debian/squeeze/maxima/squeeze

« back to all changes in this revision

Viewing changes to share/numeric/dblint_1.dem

  • Committer: Bazaar Package Importer
  • Author(s): Camm Maguire
  • Date: 2006-10-18 14:52:42 UTC
  • mto: (1.1.5 upstream)
  • mto: This revision was merged to the branch mainline in revision 4.
  • Revision ID: james.westby@ubuntu.com-20061018145242-vzyrm5hmxr8kiosf
Import upstream version 5.10.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
share/numeric/dblint_1.dem
5
5
FLOATDEFUNK and QUANC8 */
6
6
 
7
7
/* Get the fasl file for DBLINT */ 
8
 
LOAD('DBLINT); 
 
8
load('dblint); 
9
9
/* Get the fasl file for QUANC8 */
10
 
LOAD('QQ);
 
10
load('qq);
11
11
 
12
12
/* Use DBLINT to get the double integral of exp(x-y^2) over the region
13
13
bounded by y=1 and y=2+x^(3/2) and x=0 to x=1 */
14
14
 
15
15
/* Define the integrand as a function of two variables */
16
16
 
17
 
F(X,Y):=(MODE_DECLARE([X,Y],FLOAT),EXP(X-Y^2)); 
 
17
f(x,y):=(mode_declare([x,y],float),exp(x-y^2)); 
18
18
 
19
19
/* Define the lower and upper limits on the inner (y in this case)
20
20
integral as a function of the outer variable (x in this case) */
21
21
 
22
 
R(X):=(MODE_DECLARE(X,FLOAT),1.0); 
23
 
S1(X):=(MODE_DECLARE(X,FLOAT),2.0+X^(3/2));
 
22
r(x):=(mode_declare(x,float),1.0); 
 
23
s1(x):=(mode_declare(x,float),2.0+x^(3/2));
24
24
 
25
25
/* Now translate these functions for the sake of efficiency */
26
26
 
27
 
TRANSLATE(F,R,S1); 
 
27
translate(f,r,s1); 
28
28
 
29
29
/* Call the DBLINT function with quoted arguments for function names, and
30
30
floating point values for the endpoints of the outer (x) integration
31
31
*/ 
32
32
 
33
 
DBLINT_ANSWER:DBLINT('F,'R,'S1,0.0,1.0); 
 
33
dblint_answer:dblint('f,'r,'s1,0.0,1.0); 
34
34
 
35
35
/* Now generate the exact integral over y using RISCH */
36
36
 
37
 
INTY:RISCH(EXP(X-Y^2),Y); 
 
37
inty:risch(exp(x-y^2),y); 
38
38
 
39
39
/* Now get the integrand of the x integral */
40
40
 
41
 
XINT:EV(INTY,Y:2+X^(3/2))-EV(INTY,Y:1);
 
41
xint:ev(inty,y:2+x^(3/2))-ev(inty,y:1);
42
42
 
43
43
/* Try to do the x integral exactly */
44
44
 
45
 
RISCH(XINT,X);
46
 
RADCAN(%);
47
 
%,NOUNS;
 
45
risch(xint,x);
 
46
radcan(%);
 
47
%,nouns;
48
48
/* Still no luck with closed-form */
49
49
/* The integral over x can't be done in closed-form, so let's use a
50
50
one-dimensional integrator to do it numerically. First, we must
51
51
define the integrand as a floating point function. FLOATDEFUNK
52
52
does this very conveniently on the expression */
53
53
 
54
 
FLOATDEFUNK(fUN,[X],XINT);
 
54
floatdefunk(fun,[x],xint);
55
55
 
56
56
/* Now translate the function fun */
57
 
TRANSLATE(FUN);
 
57
translate(fun);
58
58
 
59
59
/* Now call QUANC8 on the translated function (see SHARE1;QQ USAGE for
60
60
details) */
61
61
 
62
 
QUANC8_ANSWER:QUANC8(FUN,0.0,1.0);
 
62
quanc8_answer:quanc8(fun,0.0,1.0);
63
63
 
64
64
/* Compare the answers from DBLINT and QUANC8 */
65
65
 
66
 
(DBLINT_ANSWER-QUANC8_ANSWER)/DBLINT_ANSWER;
 
66
(dblint_answer-quanc8_answer)/dblint_answer;
67
67
 
68
68
/* It appears to be small enough relative error for most uses */
69
69
"end of dblint.dm1"$