# ====================================================================== # ch07_quanc8.py # # Adaptive numerical integration using QUANC8. # # The program evaluates # # pi # / # | sin(x) dx = 2 # / # 0 # # with the adaptive 8-panel Newton-Cotes QUANC8 algorithm. # # QUANC8 repeatedly subdivides intervals and compares successive # Newton-Cotes estimates to control the integration error. # # Based on the QUANC8 routine presented by # G. E. Forsythe, M. A. Malcolm, and C. B. Moler, # Computer Methods for Mathematical Computations, # Prentice-Hall, 1977. # # Python version prepared for the companion website, 2026. # Matplotlib is used only for visualization. # ====================================================================== import math import matplotlib.pyplot as plt def f(x): # ------------------------------------------------------------------ # Function to be integrated. # # To use this program for another problem, replace only the line # defining y below and change the integration limits in main() # if needed. # # Input: # x - integration variable # # Output: # y - value of the integrand # ------------------------------------------------------------------ y = math.sin(x) return y def quanc8(fun, a, b, abserr, relerr): # ------------------------------------------------------------------ # Adaptive integration using the 8-panel Newton-Cotes QUANC8 # algorithm. # # Input: # fun - function to integrate # a - lower integration limit # b - upper integration limit; b may be less than a # abserr - requested absolute error tolerance, abserr >= 0 # relerr - requested relative error tolerance, relerr >= 0 # # Output: # result - numerical approximation to the integral # errest - estimated magnitude of the numerical error # nfun - number of function evaluations # flag - reliability indicator # # Reliability flag: # flag = 0 indicates that the requested tolerance was probably met. # A nonzero value indicates that one or more intervals did not # converge normally or that the function-evaluation limit was # approached. # # Method: # QUANC8 is an automatic adaptive quadrature routine based on an # 8-panel Newton-Cotes formula. The interval is subdivided where # the difference between successive estimates is too large. # # The local error estimate is based on # # |Q_new - Q_old| / 1023. # # Previously computed function values are reused whenever an # interval is subdivided. # # Based on the QUANC8 routine presented by # G. E. Forsythe, M. A. Malcolm, and C. B. Moler, # Computer Methods for Mathematical Computations, # Prentice-Hall, 1977. # ------------------------------------------------------------------ if abserr < 0.0 or relerr < 0.0: raise ValueError("quanc8: error tolerances must be non-negative.") # Initialize outputs so the zero-length interval case is defined. result = 0.0 errest = 0.0 flag = 0.0 nfun = 0 if a == b: return result, errest, nfun, flag # ---- Storage ------------------------------------------------------ qright = [0.0] * 33 fval = [0.0] * 17 x = [0.0] * 17 fsave = [[0.0] * 31 for _ in range(9)] xsave = [[0.0] * 31 for _ in range(9)] # ---- General initialization -------------------------------------- levmin = 1 levmax = 30 levout = 6 nomax = 5000 # Trouble section is entered when nfun approaches this limit. nofin = nomax - 8 * (levmax - levout + 128) # Newton-Cotes coefficients. w0 = 3956.0 / 14175.0 w1 = 23552.0 / 14175.0 w2 = -3712.0 / 14175.0 w3 = 41984.0 / 14175.0 w4 = -18160.0 / 14175.0 cor11 = 0.0 area = 0.0 # ---- Initialize first interval ----------------------------------- lev = 0 nim = 1 x0 = a x[16] = b qprev = 0.0 f0 = fun(x0) stone = (b - a) / 16.0 x[8] = (x0 + x[16]) / 2.0 x[4] = (x0 + x[8]) / 2.0 x[12] = (x[8] + x[16]) / 2.0 x[2] = (x0 + x[4]) / 2.0 x[6] = (x[4] + x[8]) / 2.0 x[10] = (x[8] + x[12]) / 2.0 x[14] = (x[12] + x[16]) / 2.0 for j in range(2, 17, 2): fval[j] = fun(x[j]) nfun = 9 # ---- Main adaptive loop ------------------------------------------ while nfun <= nomax: # Complete the 17-point grid on the current interval. x[1] = (x0 + x[2]) / 2.0 fval[1] = fun(x[1]) for j in range(3, 16, 2): x[j] = (x[j - 1] + x[j + 1]) / 2.0 fval[j] = fun(x[j]) nfun += 8 step = (x[16] - x0) / 16.0 # Newton-Cotes estimates on the left and right halves. qleft = ( w0 * (f0 + fval[8]) + w1 * (fval[1] + fval[7]) + w2 * (fval[2] + fval[6]) + w3 * (fval[3] + fval[5]) + w4 * fval[4] ) * step qright[lev + 1] = ( w0 * (fval[8] + fval[16]) + w1 * (fval[9] + fval[15]) + w2 * (fval[10] + fval[14]) + w3 * (fval[11] + fval[13]) + w4 * fval[12] ) * step qnow = qleft + qright[lev + 1] qdiff = qnow - qprev area += qdiff # ---- Local convergence test ---------------------------------- esterr = abs(qdiff) / 1023.0 tolerr = max(abserr, relerr * abs(area)) tolerr *= step / stone if lev < levmin: key = 1 elif lev >= levmax: key = 2 elif nfun > nofin: key = 3 elif esterr <= tolerr: key = 4 else: key = 1 if key == 1: # No convergence: subdivide the current interval. nim = 2 * nim lev += 1 # Save the right half for later. for i in range(1, 9): fsave[i][lev] = fval[i + 8] xsave[i][lev] = x[i + 8] # Continue immediately with the left half. qprev = qleft for i in range(1, 9): j = -i fval[2 * j + 18] = fval[j + 9] x[2 * j + 18] = x[j + 9] continue elif key == 2: # Maximum subdivision level reached. flag += 1.0 elif key == 3: # Function-evaluation limit is being approached. nofin = 2 * nofin levmax = levout flag += (b - x0) / (b - a) # key == 4 means the current interval satisfies the test. # ---- Accept current interval --------------------------------- result += qnow errest += esterr cor11 += qdiff / 1023.0 # Locate the next interval that still has to be processed. while nim % 2 != 0: nim //= 2 lev -= 1 nim += 1 if lev <= 0: break # Restore saved data for the next interval. qprev = qright[lev] x0 = x[16] f0 = fval[16] for i in range(1, 9): fval[2 * i] = fsave[i][lev] x[2 * i] = xsave[i][lev] # ---- Final correction and error estimate ------------------------- result += cor11 if errest == 0.0: return result, errest, nfun, flag # Make sure errest is not below the representable roundoff level. temp = abs(result) + errest while temp == abs(result): errest *= 2.0 temp = abs(result) + errest return result, errest, nfun, flag def plot_integrand(a, b): """Plot the test integrand before calling QUANC8.""" nplot = 1000 xplot = [a + i * (b - a) / nplot for i in range(nplot + 1)] yplot = [f(x) for x in xplot] plt.figure() plt.plot(xplot, yplot, linewidth=1.5) plt.xlabel("x") plt.ylabel("f(x)") plt.title("Integrand for QUANC8 adaptive integration") plt.grid(True) def main(): # ---- Problem setup ------------------------------------------------ a = 0.0 b = math.pi abserr = 0.0 relerr = 1.0e-8 exact = 2.0 plot_integrand(a, b) result, errest, nfun, flag = quanc8(f, a, b, abserr, relerr) print("QUANC8 adaptive integration") print("Integral of sin(x) from 0 to pi") print(f"Exact value = {exact:15.7e}") print(f"Calculated value = {result:15.7e}") print(f"Estimated error = {errest:12.4e}") print(f"Function calls = {nfun:12d}") print(f"Flag = {flag:12.5f}") plt.show() if __name__ == "__main__": main()