# function [integral_value,success] = composite_Newton_Cotes(f, a, b, m, eps) # Accompanying program for the text # # Classical and Modern Numerical Analysis: # Theory, Methods and Practice # by Azmy S. Ackleh, Edward J. Allen, # R. Baker Kearfott, and Padmanabhan Seshaiyer # # (Taylor and Francis / CRC Press, 2009) # # I = composite_Newton_Cotes(f, a, b, m, eps) computes an approximation to # the integral of f over the interval a to b, using composite Newton-Cotes # integration, according to the explanation in sections 6.3.1 and 6.3.2 of # the text. f is a string giving the name of an "m" file containing the # function to be integrated, a is the left end point of the integration # interval, b is the right end point of the integration interval, and m is # the number of points in the Newton-Cotes formula. (Closed Newton-Cotes # formulas are used.) The weights for various m are stored in this # routine, with valid values of m between 2 and 11. Computations stop when # the approximation with 2n subintervals differs by less than eps from # the approximation with n subintervals, provided that can be achieved with # less 15 doublings (2^{15} = 32768 subintervals. 'success' is set to # false if the tolerance cannot be met, and is set to 'true' if it can. # # For example, suppose the m-file f_sin_x_over_xsqp1.m contains # # function [y] = f_sin_x_over_xsqp1(x) # y = sin(x) / (x^2 + 1); # return # # Then, the following sequence will return an approximation to the integral # of sin(x) / (x^{2} + 1) over the interval [0,10], using Simpson's rule, # and with a stopping criterion for the composite integration process of # eps = 0.000001: # clear # a = 0; b = 10; m = 3; # [integral_value, success] = composite_Newton_Cotes('f_sin_x_over_xsqp1',... # a, b, m, 0.000001) # JR note: Did not do literal translation. The original .m file did not make any use of vectorization. # Here I vectorize over the panels. Shortens the code and makes it run faster. # I also changed the weights to a more straightforward form. import numpy as np def composite_Newton_Cotes(f, a, b, m, eps): if m == 2: # Trapezoidal rule w = (1/2) * np.array( [1, 1] ) elif m == 3: # Simpson's rule w = (1/3) * np.array( [1, 4, 1] ) elif m == 4: # The 3/8 rule w = (3/8) * np.array( [1, 3, 3, 1] ) elif m == 5: # Boole's rule w = (2/45) * np.array( [7, 32, 12, 32, 7] ) elif m == 6: w = (5/288) * np.array( [19, 75, 50, 50, 75, 19] ) elif m == 7: w = (1/140) * np.array( [41, 216, 27, 272, 27, 216, 41] ) elif m == 8: w = (7/17280) * np.array( [751, 3577, 1323, 2989, 2989, 1323, 3577, 751] ) elif m == 9: # Note weights of opposite signs, bad for roundoff error -- w = (4/14175) * np.array( [989, 5888, -928, 10496, -4540, 10496, -928, 5888, 989] ) elif m == 10: w = (9/89600) * np.array( [2857, 15741, 1080, 19344, 5778, 5778, 19344, 1080, 15741, 2857] ) elif m == 11: # Expect this formula to have very bad roundoff properties -- w = (5/299376) * np.array( [16067, 106300, -48525, 272400, -260550, 427368, -260550, 272400, -48525, 106300, 16067] ) else: print('Error in composite_Newton_Cotes:') print('The if m',m,'points is not implemented.') integral_value = (b-a)*feval(f,(a+b)/2) return False, integral_value n = 1 # number of panels error = 2*eps integral_estimate = (b-a)*(f(a)+f(b))/2 for n_doubling in range(1,16): n = 2**n_doubling x = np.linspace(a,b,n,endpoint=False) # left endpoints of panels H = (b-a)/n # width of panel xoff = np.linspace(0,H,m) h = xoff[1] - xoff[0] # width of subinterval in panel hw = h*w # scaled weights s = 0 for k in range(m): s += hw[k] * f(x + xoff[k]).sum() error_estimate = abs(s-integral_estimate) integral_estimate = s print('{:9.0f} {:12.6f} {:12.6f}'.format(n,s,error_estimate) ) if error_estimate <= eps: return True, s print('composite_Newton_Cotes could not achieve accuracy eps =', eps) return False, s if __name__ == '__main__': from f_sin_x_over_xsqp1 import f_sin_x_over_xsqp1 #def f_sin_x_over_xsqp1(x): return 0*x + 1 # test on very easy function (constant at 1) a = 0; b = 10; for m in range(2, 12): print() print('m =',m) success,integral_value = composite_Newton_Cotes(f_sin_x_over_xsqp1, a, b, m, 0.000001) if success: print('success:',integral_value)