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\$.


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


>>> 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)
  • \$\begingroup\$ Using an external library would definitely make it much easier and efficient. Is there any library which you know that deals with higher precision floating point numbers. Plus, I can't compute the integral for larger values of n (say 100 million). Any solutions? \$\endgroup\$ – PotCoder May 11 '18 at 15:47
  • \$\begingroup\$ @PotCoder: (1) NumPy supports the "long double" data type, which varies according to the platform: on typical x86-64 machines, this is 80 bits. (2) Why would you want to make n so large? \$\endgroup\$ – Gareth Rees May 11 '18 at 16:22
  • \$\begingroup\$ Higher the n value. More accurate the integral value is. The actual purpose of this was to compute the value of pi. I'll add further details in the question. \$\endgroup\$ – PotCoder May 11 '18 at 18:00
  • \$\begingroup\$ Made some major revision. Please Check \$\endgroup\$ – PotCoder May 11 '18 at 18:47
  • 2
    \$\begingroup\$ @PotCoder - I have restored the original intent of the question to be about the more generic Trapezoidal Approximation, rather than specifically approximating \$\pi\$. Your edits went too far and changed the sense of the question too much, making other answers look incomplete. See What you may and may not do after receiving answers \$\endgroup\$ – rolfl May 12 '18 at 12:27
  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
  • 2
    \$\begingroup\$ It appears too clustered. By the way, thanks for fixing all the typos. Can't believe I made so many mistakes while typing. \$\endgroup\$ – PotCoder May 11 '18 at 9:38

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

  • \$\begingroup\$ Pretty much solves half of the updated problem. I wanted this and the solution given by Gareth Rees to be the accepted answer. I think I should've just made another question instead of completely changing the question. \$\endgroup\$ – PotCoder May 12 '18 at 18:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.