Numerical Integration in Python

last updated Jan 5, 2017

Numerical integration aims to find the area under a curve without using analytical methods.

Rectangular

The rectangle rule states that:

In other words, the area under a curve $$f(x)$$ between a point $$a$$ and a point $$b$$ is roughly equal to the area of a rectangle with width $\lvert a-b\rvert$ and height $$\frac{f(a)+f(b)}{2}$$.

This approach is easily implemented in Python:

def rectint(f,a,b,rectangles):
cumulative_area=0

a=float(a)
b=float(b)
rectangles=float(rectangles)

i=(b-a)/rectangles

trailing_x=a

cumulative_area+=area

trailing_x+=i

return cumulative_area


Trapezoidal

The trapezoid rule is a small departure from the rectangle rule and tends to produce more accurate approximations:

The Python code is itself similar to that of the rectangle rule:

def trapint(f,a,b,trapezoids):
cumulative_area=0

a=float(a)
b=float(b)
trapezoids=float(trapezoids)

i=(b-a)/trapezoids

trailing_x=a

cumulative_area+=area

trailing_x+=i

return cumulative_area


In each code example, rules for negativity are preserved in line with the FTC. Predictably, the more rectangles or trapezoids requested, the lower the error. Indeed, the limit definition of an integral is based around infinite subintervals.

Comparison

Let’s test the functions with an integral which is easily evaluated analytically:

With rectint:

>>> def func(x):
...     return x**2-2*x+3
>>> rectint(func,-5,10,15):
343.75
>>> rectint(func,-5,10,100):
332.72353125
>>> rectint(func,-5,10,1000):
344.99971875
>>> rectint(func,-5,10,10000):
344.875517437
>>> rectint(func,-5,10,100000):
344.999999971


With trapint:

>>> def func(x):
...     return x**2-2*x+3
>>> trapint(func,-5,10,15):
347.5
>>> trapint(func,-5,10,100):
332.8070625
>>> trapint(func,-5,10,1000):
345.0005625
>>> trapint(func,-5,10,10000):
344.875525873
>>> trapint(func,-5,10,100000):
345.000000056


QED.