Numerical Simpsons rule
Posted March 08, 2013 at 06:18 PM | categories: integration, math | tags:
A more accurate numerical integration than the trapezoid method is Simpson's rule. The syntax is similar to trapz, but the method is in scipy.integrate.
import numpy as np from scipy.integrate import simps, romb a = 0.0; b = np.pi / 4.0; N = 10 # this is the number of intervals x = np.linspace(a, b, N) y = np.cos(x) t = np.trapz(y, x) s = simps(y, x) a = np.sin(b) - np.sin(a) print print 'trapz = {0} ({1:%} error)'.format(t, (t - a)/a) print 'simps = {0} ({1:%} error)'.format(s, (s - a)/a) print 'analy = {0}'.format(a)
>>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> >>> trapz = 0.70665798038 (-0.063470% error) simps = 0.707058914216 (-0.006769% error) analy = 0.707106781187
You can see the Simpson's method is more accurate than the trapezoid method.
Copyright (C) 2013 by John Kitchin. See the License for information about copying.