program romb1 c c -- numerical integration by the romberg method c integer out real sum, upper, lower, tol external f common /inout/ out data lower/1.0/, upper/9.0/, tol/1.0e-5/ c out = 6 write(out,*) ' romberg integration' call romber(lower, upper, tol, sum, f) write(out, 104) sum stop 104 format(/' area =', f10.5/) end function f(x) c f = 1.0 / x return end subroutine romber(lower, upper, tol, ans, f) c c -- numerical integration by the romberg method c integer out, pieces, i, nx(16) integer nt, ii, n, nn, l, ntra, k, m, j, l2 real x, delta, lower, upper, sum, tol real c, fotom, t(136), ans common /inout/ out c pieces = 1 nx(1) = 1 delta = (upper - lower) / pieces c = (f(lower) + f(upper)) / 2.0 sum = c t(1) = delta * c n = 1 nn = 2 5 n = n + 1 fotom = 4.0 nx(n) = nn pieces = pieces * 2 l = pieces - 1 l2 = (l+1) / 2 delta = (upper - lower) / pieces c -- compute trapezoidal sum for 2**(n-1)+1 points do 10 ii = 1, l2 i = ii * 2 - 1 x = lower + delta * i sum = sum + f(x) 10 continue t(nn) = sum * delta ntra = nx(n-1) k = n-1 c -- compute n-th row ot t array do 20 m = 1, k j = nn + m nt = nx(n-1) + m - 1 t(j) = (fotom * t(j-1) - t(nt)) / (fotom - 1) fotom = fotom * 4 20 continue write(out, 101) pieces, t(nn), j, t(j) if (n .lt. 5) goto 40 if (t(nn+1) .eq. 0.0) goto 40 if (abs(t(ntra+1)-t(nn+1)) .le. abs(t(nn+1)*tol)) goto 99 if (abs(t(nn-1) - t(j)) .le. abs(t(j)*tol)) goto 99 if (n .gt. 15) goto 100 40 nn = j + 1 goto 5 99 ans = t(j) return 100 write(out,*) ' romberg error' return 101 format(1x, i7, f9.5, i5, f9.5) end