# finding definite integrals in Python using Trapezium rule

I wrote this program which finds the definite integral of a function. Where in the program could I optimize this code:

def trapezium(f,n,a,b):
h=(b-a)/float(n)

area = (0.5)*h
sum_y = (f(0)+f(b))

i=a
while i<b:
print i
sum_y += 2*f(i)
i += h

area *= sum_y

return area

def f(x):
return x ** 2

print trapezium(f, 10000, 0, 5)


First, there is an error where you have:

sum_y = (f(0)+f(b))


f(0) should be f(a). It doesn't matter on your example, because you start with 0, but would otherwise.

Another error is that you add f(a) 3x instead of just once. i=a+h should be the line before while.

Your while loop makes O(n) multiplications and twice as many additions. Instead you should have something like:

i = a+h
part_sum = 0
while i<b:
part_sum += f(i)
i += h
sum_y += 2*part_sum


Still same number of additions, but only one multiplication.

Using list comprehensions might be a bit faster, but you'd spend too much memory at scale where it matters.

print significantly slows your function. If you actually need to print this, then store intermediate results in array and then print them before final result with print "\n".join(intermediates).

def trapezium(f,n,a,b):


I recommend a docstring explaining the arguments

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


I recommend adding from __future__ import division so that division always results in a float rather then doing this.

    area = (0.5)*h


Drop the (): area = 0.5*h

    sum_y = (f(0)+f(b))


Drop the (): sum_y = f(0) + f(b)

    i=a
while i<b:
print i
sum_y += 2*f(i)
i += h


Should you be printing here? I assume not. You should use a for loop

for i in xrange(a,b,h):
sum_y += 2 * f(i)


Or better yet a generator expression:

sum_y += sum(2*f(i) for i in xrange(a,b,h))

...

    area *= sum_y

return area


Why so many blank lines?

def f(x):
return x ** 2

print trapezium(f, 10000, 0, 5)


If you are interested in getting more speed, you should look at numpy.

• I had to use a while loop because range doesn't accept a float value(here, h needs to be float). Know of any way around that? Thanks for all the other advice too, I'll keep them in mind. Commented Aug 28, 2011 at 8:42
• @infoquad, I should have known that. I'd solve it by using numpy which has arange which allows floats. Commented Aug 28, 2011 at 13:07