Simpson rule

The trapezoidal rule uses straight-line segments for approximations, but we can use curves instead. Simpson』s rule fits intervals with quadratic curves.

We require three points to specify a quadratic (not just two, as with a line). So the fit is of a quadratic to pairs of adjacent slices. Again, suppose the integrand is f(x), and the spacing is h. Suppose the three points are at ?h, 0,+h. Then fitting a quadratic Ax2 + Bx + C through these points gives

and the area under the curve from ?h to h is approximated by the area under the quadratic, so

which is Simpson』s rule. It approximates the area under two adjacent slices. You need to have an even number of slices for this to work consistently, and the approximate value of the integral comes out to

The two sums are calculatedseparately, using range(1,N, 2) and range(2,N, 2) for their summation conditions.

code:

from __future__ import division # Python 2 compatibilitynndef simpson(f, a, b, n):n """Approximates the definite integral of f from a to b by then composite Simpsons rule, using n subintervals (with n even)"""nn if n % 2:n raise ValueError("n must be even (received n=%d)" % n)nn h = (b - a) / nn s = f(a) + f(b)nn for i in range(1, n, 2):n s += 4 * f(a + i * h)n for i in range(2, n-1, 2):n s += 2 * f(a + i * h)nn return s * h / 3nn# Demonstrate that the method is exact for polynomials up to 3rd ordernprint(simpson(lambda x:x**3, 0.0, 10.0, 2)) # 2500.0nprint(simpson(lambda x:x**3, 0.0, 10.0, 100000)) # 2500.0nnprint(simpson(lambda x:x**4, 0.0, 10.0, 2)) # 20833.3333333nprint(simpson(lambda x:x**4, 0.0, 10.0, 100000)) # 20000.0n

推薦閱讀：