-
Notifications
You must be signed in to change notification settings - Fork 190
/
Simpson.py
61 lines (51 loc) · 1.56 KB
/
Simpson.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def Simpson(f, a, b, n=500):
"""
Return the approximation of the integral of f
from a to b using Simpson's rule with n intervals.
"""
if a > b:
print 'Error: a=%g > b=%g' % (a, b)
return None
# Check that n is even
if n % 2 != 0:
print 'Error: n=%d is not an even integer!' % n
n = n+1 # make n even
h = (b - a)/float(n)
sum1 = 0
for i in range(1, n/2 + 1):
sum1 += f(a + (2*i-1)*h)
sum2 = 0
for i in range(1, n/2):
sum2 += f(a + 2*i*h)
integral = (b-a)/(3*n)*(f(a) + f(b) + 4*sum1 + 2*sum2)
return integral
def h(x):
return (3./2)*sin(x)**3
from math import sin, pi
def application():
print 'Integral of 1.5*sin^3 from 0 to pi:'
for n in 2, 6, 12, 100, 500:
approx = Simpson(h, 0, pi, n)
print 'n=%3d, approx=%18.15f, error=%9.2E' % \
(n, approx, 2-approx)
def test_Simpson():
"""Test exact integration of quadratic polynomials."""
a = 1.5
b = 2.0
n = 8
g = lambda x: 3*x**2 - 7*x + 2.5 # test integrand
G = lambda x: x**3 - 3.5*x**2 + 2.5*x # integral of g
exact = G(b) - G(a)
approx = Simpson(g, a, b, n)
success = abs(exact - approx) < 1E-14
msg = 'Simpson: %g, exact: %g' % (approx, exact)
assert success, msg
def experiments():
"""Some other tests with known answers."""
from math import sin, pi, exp, log
n = 200
print 'sin from 0 to pi:', Simpson(sin, 0, pi, n)
print 'exp from 0 to log(3):', Simpson(exp, 0, log(3), n)
test_Simpson()
application()
experiments()