Over a million developers have joined DZone.
{{announcement.body}}
{{announcement.title}}

DZone's Guide to

# Simpson's Rule

·
Free Resource

Comment (0)

Save
{{ articles[0].views | formatCount}} Views
``````
from math import log, cos
def simpson(f, a, b, n):
"Approximate the definite integral of f from a to b by Simpson's rule."

if n % 2 != 0:
print "Ups: n must be even!"
return -1

h  = (float(b) - a)/n

si = 0.0
sp = 0.0

for i in range(1, n, 2):
xk = a + i*h
si += f(xk)

for i in range(2, n, 2):
xk = a + i*h
sp += f(xk)

s = 2*sp + 4*si + f(a) + f(b)

return (h/3)*s

# def f(x):
#     return x**4
def f(x):
return (log(x+1)**(1+cos(x)))

print simpson(f, 3, 4, 2)
exit()

ni = 50
nf = 1000000
n = ni
a = -20
b = 0
s = []
qc = []
ec = []
t = 0.0
div = 0.0
i = 0.0

# Calcular S para diferentes valores de n [n(i+1) = ni * 2]
# para poder ter S (h), S' (h'=h/2), S'' (h''=h/4)
while n < nf:
t = simpson(f, a, b, n)
s.append(t)
n = n * 2

# Calcular quociente (S'-S)/(S''-S') [que deve ser aprox igual a 16]
# e erro em que e = (S''-S')/15
for i in range(0, len(s)-2):
div = (s[i+2] - s[i+1]) # S''-S'
if div == 0:
break

t = (s[i+1] - s[i]) / div # (S'-S)/div com div = (S''-S')<=>(S'-S)/(S''-S')
qc.append(t)
t = div / 15
ec.append(t)

#qc ~= 16 implica validade => podemos usar os valores do erro: e => qualidade
for i in range(0, len(qc)):
print "S=%.12f, S'=%.12f, S''=%.12f\n\t=> qc=%.12f, e=%.12f\n" \
% (s[i], s[i+1], s[i+2], qc[i], ec[i])

``````
Topics:

Comment (0)

Save
{{ articles[0].views | formatCount}} Views

Opinions expressed by DZone contributors are their own.