Home Manual Packages Global Index Keywords Quick Reference ``` /* ROMBERG.I Romberg integrator, after qromb in Numerical Recipes (Press, et.al.) \$Id\$ */ /* Copyright (c) 1994. The Regents of the University of California. All rights reserved. */ func romberg (function, a, b, epsilon, notvector=) /* DOCUMENT integral= romberg(function, a, b) or integral= romberg(function, a, b, epsilon) returns the integral of FUNCTION(x) from A to B. If EPSILON is given, Simpson's rule is refined until that fractional accuracy is obtained. EPSILON defaults to 1.e-6. If the notvector= keyword is supplied and non-zero, then FUNCTION may not be called with a list of x values to return a list of results. By default, FUNCTION is assumed to be a vector function. If the function is not very smooth, simpson may work better. SEE ALSO: simpson, max_doublings */ { if (!epsilon || epsilon<0.) epsilon= 1.e-6; a= double(a); b= double(b); s= array(0.0, 5); h= [1.0, 0.25, 0.0625, 0.015625, 0.00390625]; for (i=1 ; i<=max_doublings ; ++i) { ss= trapezoid(function, a, b, i, notvector); if (i >= 5) { s(5)= ss; c= d= s; /* Neville algorithm tableau */ ns= 4; /* one less than smallest h, always last */ for (m=1 ; m<5 ; ++m) { m5= 5-m; ho= h(1:m5); hp= h(m+1:5); w= (c(2:m5+1)-d(1:m5))/(ho-hp); d(1:m5)= hp*w; c(1:m5)= ho*w; dss= (2*ns