File: qq.dem

package info (click to toggle)
maxima 5.6-17
  • links: PTS
  • area: main
  • in suites: woody
  • size: 30,572 kB
  • ctags: 47,715
  • sloc: ansic: 154,079; lisp: 147,553; asm: 45,843; tcl: 16,744; sh: 11,057; makefile: 7,198; perl: 1,842; sed: 334; fortran: 24; awk: 5
file content (79 lines) | stat: -rw-r--r-- 2,464 bytes parent folder | download
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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

/* Numerical Integration Demo for QUANC8*/

/* SETUP_AUTOLOAD(QQ,QUANC8); */

/* show the time and the garbage collection time, and use the 
   GC-demon to make better use of flonum space */

SHOWTIME:'ALL;
IF STATUS(FEATURE,ITS)=TRUE THEN LOAD("GCDEMN");

/* Sample highly oscillatory problem */

G(X):=(MODE_DECLARE(X,FLOAT),SIN(1.0/X));

/* Inefficient way to use QUANC8 */

QUANC8(G(x),x,0.01,0.1);

/* give G as the thing to compile, then a semi-colon to read it in */
compile();

/* see how much faster now, using fast version of QUANC8 */

QUANC8('G,0.01,0.1);

/* note that romberg doesn't easily manage oscillatory behaviors */

ERRCATCH(ROMBERG('G,0.01,0.1));

/* a function with two large and narrow peaks,
   so hard to integrate with a fixed step-size */

H(X):=(MODE_DECLARE(X,FLOAT),1.0/(.0001+(X-.3)^2)+1.0/(.0025+(X-.9)^2));
/* give H and a semi-colon in response */
compile();

/* the exact answer for integrate(h(x),x) is
   100.0*atan(100.0*x-30.0)+20.0*atan(20.0*x-18.0) */
/* the exact answer for integrate(h(x),x,0.0,1.0) is 361.847622.
   note that the romberg result is more accurate in this case */
ROMBERG('H,0.0,1.0);
/* but at a large time cost compared to adaptive */
QUANC8('H,0.0,1.0);
/* reduce QUANC8's relative error tolerance by a factor of 10 */
QUANC8_RELERR:1.0E-5;
/* does not take much longer, and gives much better result now */
QUANC8('H,0.0,1.0);
/* put QUANC8_RELERR back to 1.0e-4 */
QUANC8_RELERR:1.0E-4$
/* the exact answer for integrate(h(x),x,0.0,10.0) is 372.33606.
   while too tough for fixed step size method */
ERRCATCH(ROMBERG('H,0.0,10.0));
/*QUANC8 is still super speedy even now */
QUANC8('H,0.0,10.0);

/* double integral of sample function exp(x)*sin(y) 
   from x=0.0 to 2.0, y=0.0 to 1.0 */
/* you must be sure that the quanc8 source is loaded
   (load("qq");) before you  translate a call to quanc8 in a function, 
   else an error in macro expansion will occur when it can't figure 
   out if the 3 or 4 arg version is used */

/* the exact answer for 
   integrate(integrate(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0) is 2.93703434.
   integrate over y to get the final answer */
QUANC8(quanc8(exp(x)*sin(y),x,0.0,2.0),y,0.0,1.0);

p():=QUANC8(lambda([y],mode_declare(y,float),
                   quanc8(lambda([x],mode_declare(x,float),
                                 exp(x)*sin(y)),
                          0.0,2.0)),0.0,1.0);

p();
translate(p);
p();
/* give p and a semi-colon */
compile();
p();