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 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
    leading_x=a+i

    while (a<=leading_x<=b) or (a>=leading_x>=b):
        area=f((trailing_x+leading_x)/2)*i
        cumulative_area+=area

        leading_x+=i
        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
    leading_x=a+i

    while (a<=leading_x<=b) or (a>=leading_x>=b):
        area=(f(trailing_x)+f(leading_x)))*i/2
        cumulative_area+=area

        leading_x+=i
        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.