Definite Integral Approximation using the Trapezoidal Method 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) ) + &lt;----- part_1
    # w * SUMMATION(x varies from 1 to n-1): f(a + x * w) &lt;----- 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:

&gt;&gt;&gt; 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:

&gt;&gt;&gt; x = np.linspace(0, 1, 5)
&gt;&gt;&gt; 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)