1 /*
2 * $Id: romberg.i,v 1.3 2010-08-25 14:07:16 thiebaut Exp $
3 * Romberg integrator, after qromb in Numerical Recipes (Press, et.al.)
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func romberg(function, a, b, epsilon, notvector=)
12 /* DOCUMENT integral= romberg(function, a, b)
13 or integral= romberg(function, a, b, epsilon)
14 returns the integral of FUNCTION(x) from A to B. If EPSILON is
15 given, Simpson's rule is refined until that fractional accuracy
16 is obtained. EPSILON defaults to 1.e-6.
17
18 If the notvector= keyword is supplied and non-zero, then FUNCTION
19 may not be called with a list of x values to return a list of
20 results. By default, FUNCTION is assumed to be a vector function.
21
22 If the function is not very smooth, simpson may work better.
23
24 SEE ALSO: simpson, max_doublings
25 */
26 {
27 local _trapezoid_i, _trapezoid_s;
28 if (!epsilon || epsilon<0.) epsilon= 1.e-6;
29 a= double(a);
30 b= double(b);
31 s= array(0.0, 5);
32 h= [1.0, 0.25, 0.0625, 0.015625, 0.00390625];
33 for (i=1 ; i<=max_doublings ; ++i) {
34 ss= trapezoid(function, a, b, i, notvector);
35 if (i >= 5) {
36 s(5)= ss;
37 c= d= s; /* Neville algorithm tableau */
38 ns= 4; /* one less than smallest h, always last */
39 for (m=1 ; m<5 ; ++m) {
40 m5= 5-m;
41 ho= h(1:m5);
42 hp= h(m+1:5);
43 w= (c(2:m5+1)-d(1:m5))/(ho-hp);
44 d(1:m5)= hp*w;
45 c(1:m5)= ho*w;
46 dss= (2*ns<m5)? c(ns+1) : d(ns--);
47 ss+= dss;
48 }
49 if (abs(dss) <= epsilon*abs(ss)) return ss;
50 /* extrapolation to h=0 always uses last 5 points */
51 s(1:4)= s(2:5);
52 h*= 0.25;
53 } else {
54 s(i)= ss;
55 }
56 }
57 error, "no convergence after "+pr1(2^i)+" function evaluations";
58 }
59
60 func simpson(function, a, b, epsilon, notvector=)
61 /* DOCUMENT integral= simpson(function, a, b)
62 or integral= simpson(function, a, b, epsilon)
63 returns the integral of FUNCTION(x) from A to B. If EPSILON is
64 given, Simpson's rule is refined until that fractional accuracy
65 is obtained. EPSILON defaults to 1.e-6.
66
67 If the notvector= keyword is supplied and non-zero, then FUNCTION
68 may not be called with a list of x values to return a list of
69 results. By default, FUNCTION is assumed to be a vector function.
70
71 If the function is very smooth, romberg may work better.
72
73 SEE ALSO: romberg, max_doublings
74 */
75 {
76 local _trapezoid_i, _trapezoid_s;
77 if (!epsilon || epsilon<0.) epsilon= 1.e-6;
78 a= double(a);
79 b= double(b);
80 ost= os= -1.e-30;
81 for (i=1 ; i<=max_doublings ; ++i) {
82 st= trapezoid(function, a, b, i, notvector);
83 s= (4.*st - ost)/3.;
84 if (abs(s-os) <= epsilon*abs(os)) return s;
85 os= s;
86 ost= st;
87 }
88 error, "no convergence after "+pr1(2^i)+" function evaluations";
89 }
90
91 local max_doublings;
92 /* DOCUMENT max_doublings= 20
93 is the maximum number of times romberg or simpson will split the
94 integration interval in half. Default is 20.
95 */
96 max_doublings= 20;
97
trapezoid(function,a,b,n,notvector)98 func trapezoid(function, a, b, n, notvector)
99 {
100 extern _trapezoid_i, _trapezoid_s;
101 if (n==1) {
102 _trapezoid_i= 1;
103 return _trapezoid_s= 0.5*(b-a)*(function(a)+function(b));
104 }
105 dx= (b-a)/_trapezoid_i;
106 if (!notvector) {
107 /* conserve memory-- dont try more than 8192 points at a time */
108 if (_trapezoid_i <= 8192) {
109 s= sum(function(span(a,b,_trapezoid_i+1)(zcen)));
110 } else {
111 x= a+(indgen(8192)-0.5)*dx;
112 s= sum(function(x));
113 dx2= 8192*dx;
114 for (i=8193 ; i<=_trapezoid_i ; i+=8192) s+= sum(function((x+=dx2)));
115 }
116 } else {
117 x= a+0.5*dx;
118 s= function(x);
119 for (i=2 ; i<=_trapezoid_i ; ++i) s+= function((x+=dx));
120 }
121 _trapezoid_i*= 2;
122 return _trapezoid_s= 0.5*(_trapezoid_s + dx*s);
123 }
124