/***************************************************************************/ /* */ /* ACM Algorithm 60. ROMBERG INTEGRATION. */ /* */ /* F. L. Bauer, Gutenberg University, Mainz, Germany. */ /* */ /***************************************************************************/ /* Translated to PL/I by R. A. Vowels, 16 February 2006. */ /* ROMBERG is the value of the integral of the */ /* function FCT between the limits lgr and rgr, calculated by the */ /* algorithm of Romberg with an error term of the order */ /* 2*ord+2, ord >= 0. Computation time will be roughly be doubled */ /* when ord is increased by 1. */ ROMBERG: procedure (fct, lgr, rgr, ord) returns(float (18)); declare (lgr, rgr) float (18), ord fixed binary, fct entry (float (18)) returns (float (18)); /* declare (t(20), L, u, m) float, (f, h, j, n) fixed binary; /* DR PL/I */ declare (t(ord+1), L, u, m) float (18), (f, h, j, n) fixed binary; /* Full PL/I */ L = rgr - lgr; t(1) = (fct(lgr) + fct(rgr))/2; n = 1; do h = 1 to ord; u = 0; m = L/(2*n); do j = 1 to 2*n-1 by 2; u = u + fct(lgr+j*m); end; t(h+1) = (u/n + t(h))/2; f = 1; do j = h to 1 by -1; f = 4 * f; t(j) = t(j+1) + (t(j+1) - t(j)) / (f - 1); end; n = 2*n; end; return (t(1)*L); end ROMBERG;