from numpy import * def simpson(x, f, F=[]): """ Use composite Simpsons Rule quadrature based on the subdivision xi Input: xi: location of nodes, defining the subdivision of the interval (including the endpoints) f: function to integrate F: antiderivative (for computing errors) Output: I = the numerical approximation to the integral on each subinterval E = error (absolute value) on each subinterval (only if F is known) """ n = len(x) # number of nodes assert(n % 2 == 1) # should be an odd number h = x[1]-x[0] # according to the book's definition, 1 panel in Simpson = 2h # number of panels m = int((n-1)/2) I = zeros(m) E = zeros(m) for i in range(m): I_simpson = h/3*(f(x[2*i])+4*f(x[2*i+1])+f(x[2*i+2])) if F: I_exact = F(x[2*i+2])-F(x[2*i]) E[i] = abs(I_simpson-I_exact) I[i] = I_simpson return I, E # Test Simpsons Rule if __name__ == "__main__": # Test N = 201 a = 0 b = 1 x = linspace(a,b,N) f = lambda x: sqrt(x) F = lambda x: x*sqrt(x)*2./3. I, E = simpson(x, f, F) Iq= sum(I) print("N = %6d, numerical = %e, error = %e" % (N,Iq,F(b)-F(a)-Iq)) m = 32 N = 2*32+1 # (a) a,b=0,1 x = linspace(a,b,N) f = lambda x: x**3 fprim = lambda x: 3*x*x g = lambda x: sqrt(1+fprim(x)**2) I,_=simpson(x, g) print('(a): arclength ~ %e' % sum(I)) # (b) a,b=0,pi/4 x = linspace(a,b,N) f = lambda x: tan(x) fprim = lambda x: 1/cos(x)**2 g = lambda x: sqrt(1+fprim(x)**2) I,_=simpson(x, g) print('(b): arclength ~ %e' % sum(I)) # (c) a,b=0,1 x = linspace(a,b,N) f = lambda x: atan(x) fprim = lambda x: 1/(1+x*x) g = lambda x: sqrt(1+fprim(x)**2) I,_=simpson(x, g) print('(c): arclength ~ %e' % sum(I))