Definite Integral Approximation using the Trapezoidal Method

Question

I wrote a program to calculate the value of Definite Integral of a function from a to b. It used the trapezoidal approximation method where the function area is divided into 'n' trapeziums and area of each is added.

def calculate_area(f, a, b, n):
    # Trapezoidal approximation: divide area of function into n trapeziums and add areas. 
    # Large n value means better approximation

    # Finds the value of a definite integral from lower_lim (a) to upper_lim (b)
    # integral f(x) from lower_lim to upper_lim is =
    # 0.5 * w * ( f(a) + f(b) ) + <----- part_1
    # w * SUMMATION(x varies from 1 to n-1): f(a + x * w) <----- part_2
    # where w = (b - a) / n

    w = (b - a) / float(n) 
    part_1 = 0.5 * w * ( f(a) + f(b) )
    part_2 = 0
    for x in range(1, n):
        part_2 += f(a + x * w)
    part_2 *= w
    total_area = part_1 + part_2
    return total_area

For example, this can be used to calculate the area of a circle (or \$\pi\$ if the radius of the circle is 1):

Consider the circle \$x^2 + y^2 = 1\$

Area of this circle = \$\pi\$

Expressed as a function of x is: f(x) = (1 - x^2) ^ 0.5

Integrating this from 0 to 1 gives us a quarter of a circle whose area is \$\frac\pi 4\$ (or the approximate value if we use approximate integration).

On multiplying 4, we get the value of \$\pi\$.

Pass in the following as arguments :

f = lambda x: (1 - x**2) ** 0.5; a = 0; b = 1

Now I face the following issues:

• Code is not very Pythonic.
• Code is not very efficient and fast.
• For better precision of PI value, the value of n must be increased. Python doesn't comply when n takes higher values (say 100 million). n = 10 million worked and I got a value of PI that was accurate upto 8 digits after the decimal.
• I feel that using floating point like datatypes which provide higher precision will increase precision of PI value. Correct me if I'm wrong.

I know that this is a very poor method of approximating the value of \$\pi\$ (this is the best I could come up with high school level knowledge). Let me know when I'm going to reach the limits. That is, when I will no longer be able to get noticeably more precise value of \$\pi\$.

Answer 1

If you are writing mathematical code in Python, then it is worth looking at NumPy, which is a library implementing (among other things) fast arithmetic operations on arrays of floating-point numbers.

In NumPy you'd write:

import numpy as np

def calculate_area(f, a, b, n):
    """Return an approximation to the definite integral of f from a to b
    using the trapezium rule with n intervals.

    """
    w = (b - a) / n                # width of the intervals
    x = np.linspace(a, b, n + 1)   # x coords of endpoints of intervals
    y = f(x)
    return w * (y.sum() - (y[0] + y[-1]) / 2)

Here np.linspace returns an array of numbers evenly spaced between two endpoints. For example, if we are asked to integrate between 0 and 1 with 4 intervals, then the endpoints of the intervals are:

>>> np.linspace(0, 1, 5)
array([0.  , 0.25, 0.5 , 0.75, 1.  ])

Under the assumption that f is a straightforward arithmetic function of its argument, then f(x) computes an array of the corresponding y values for each x. For example, if we have

def f(x):
    return x ** 2 + x + 1

Then:

>>> x = np.linspace(0, 1, 5)
>>> f(x)
array([1.    , 1.3125, 1.75  , 2.3125, 3.    ])

However, NumPy already has a built-in function np.trapz that computes an integral using the trapezium rule. Using this function, we could write:

def calculate_area(f, a, b, n):
    """Return an approximation to the definite integral of f from a to b
    using the trapezium rule with n intervals.

    """
    x = np.linspace(a, b, n + 1)   # x coords of endpoints of intervals
    return np.trapz(f(x), x)

Answer 2

1. You can use sum with a comprehension to create part_2.
2. You can move all the maths together if you think it looks nicer.

def calculate_area(f, a, b, n)
    w = (b - a) / n
    return (
        0.5 * (f(a) + f(b))
        + sum(f(a + x * w) for x in range(1, n))
    ) * w

Answer 3

If you want more accuracy than 8 or so digits, you probably will need a different method. The first obvious choice is Simpson's method, which approximates regions with parabolas instead of trapezoids. Better yet, check out this wikipedia article with a bunch of formulas for Pi. https://en.wikipedia.org/wiki/Approximations_of_%CF%80. Specifically of interest is https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_algorithm which shows a uses the following pseudocode gives the method

def calculate_pi(n):
    a = 1
    b = .5**.5
    t = .25
    p = 1
    for _ in range(n):
        a, b, t, p = a/2 + b/2, (a*b)**.5, t - p * (a/2 - b/2)**2, 2*p
    return (a + b)**2/(4*t)

calculate_pi(3)=3.14159265358979323 (with a precise enough data type this should probably use gmp for better accuracy).

Even with python's limited float data type, this yields 3.141592653589794 which is just the limit of the datatype.

Using gmpy2 to provide variable precision floats, we obtain the following

import gmpy2
from gmpy2 import mpfr

def calculate_pi(n):
    gmpy2.get_context().precision = 10*2**n
    a, b, t, p = mpfr(1), mpfr(.5)**mpfr(.5), mpfr(.25), mpfr(1)

    for _ in range(n):
        a, b, t, p = a/2 + b/2, (a*b)**.5, t - p * (a/2 - b/2)**2, 2*p
    return (a + b)**2/(4*t)

With n=6, this produces the first 171 digits in 0.00032782 seconds