Numerical Integration

Trapezoidal rule

This rule is used to find the approximate value of a definite integral. In trapezoidal rule, the region under the curve of a given function is considered to be a trapezoid and the area under the curve is calculated. It follows that \begin{equation} \int_{a}^{b}f(x)dx \approx (b-a)\left[\frac{f(a)+f(b)}{2}\right] \end{equation} For more accurate result, the domain of the graph is divided into \(n\) segments of equal size so that grid spacing or segment size is \(h=\frac{b-a}{n}\). Therefore the approximate value of the integral can be written as \begin{equation} \int_{a}^{b}f(x)dx = \frac{(b-a)}{2n}\left[f(a)+2\sum_{i=1}^{n-1}f(a+ih)+f(b)\right] \end{equation}

Simpson’s 1/3rd rule

\begin{align} &I(a,b)==\frac{1}{3}h\left[f(a) + f(b) + 4\sum_{k = 1\dots N-1}f(a+kh) + 2\sum_{k =2 \dots N-2}f(a+kh) \right]\\ &h=\frac{(b-a)}{N} \end{align}

Algorithm: Numerical Integration using Trapezoidal and Simpson’s Rule

Algorithm Overview

We will approximate definite integrals of a given function f(x) over the interval [a, b] using:

  • Trapezoidal Rule

  • Simpson’s 1/3 Rule

We test the integration on two types of functions:

  1. A rational function: $ f(x) = \frac{1}{1 + x^2} $

  2. A Gaussian-like function multiplied by a linear term

Trapezoidal Rule Algorithm

Input:

  • Function f

  • Interval [a, b]

  • Number of subintervals N

Steps:

  1. Compute step size: ( h = \frac{b - a}{N} )

  2. Generate N+1 equally spaced x values in [a, b]

  3. Evaluate f(x) at all points

  4. Compute:

    \[I = h \left[ \frac{f(x_0) + f(x_N)}{2} + \sum_{i=1}^{N-1} f(x_i) \right]\]

  5. Return the integral I

Simpson’s 1/3 Rule Algorithm

Input:

  • Function f

  • Interval [a, b]

  • Number of subintervals N (must be even)

Steps:

  1. If N is odd, increment it by 1 to make it even

  2. Compute step size: ( h = \frac{b - a}{N} )

  3. Generate N+1 equally spaced x values in [a, b]

  4. Evaluate f(x) at all points

  5. Compute:

    \[I = \frac{h}{3} \left[ f(x_0) + f(x_N) + 4\sum_{\text{odd } i} f(x_i) + 2\sum_{\text{even } i} f(x_i) \right]\]

  6. Return the integral I

[1]:

# Defining function to evaluate  by Trapezoidal Rule
def trapzInt(f, a, b, N):
    x = np.linspace(a, b, N+1)
    h = (b-a)/N
    y = f(x)
    I = (sum(y[1:-1]) + (y[0] + y[-1])/2)*h
    return I

# Defining function to evaluate  by Simpson's 1/3rd Rule
def simp(f, a, b, N):
    if N%2 == 1:
        N += 1
    x = np.linspace(a, b, N+1)
    h = (b-a)/N
    y = f(x)
    I = (y[0] + y[-1] + 4*sum(y[1:-1:2]) + 2*sum(y[2:-2:2]))*(h/3)
    return I

Evaluate:

\[\int_{-\infty}^{\infty}\frac{1}{1+x^2}dx\]
[2]:

import numpy as np

f = lambda x: 1/(1 + x**2)
a, b = -10, 10
N    = 100
value = trapzInt(f, a, b, N)
print("The value of the integral by Trapezoidal rule = "+str(round(value, 3)))

value = simp(f, a, b, N)
print("The value of the integral by Simpson's 1/3rd rule = "+str(round(value, 3)))

The value of the integral by Trapezoidal rule = 2.942
The value of the integral by Simpson's 1/3rd rule = 2.942

Evaluate

\[\frac{1}{\sqrt{2\pi \sigma^2}}\int_{-\infty}^{\infty}\exp\left(-\frac{(x - 2)^2}{2\sigma^2}\right)(x + 3) dx\]
[3]:

def f(x):
    int = (1/np.sqrt(2*np.pi*sigma**2))*np.exp(-(x-2)**2/(2*sigma**2))*(x + 3)
    return int

a, b = -10, 10
N    = 1000
h    = (b - a)/N
sigma = 1



value_trapz = trapzInt(f, a, b, N)
value_simp = simp(f, a, b, N)

print("The value of the integral by Trapezoidal rule = "+str(round(value_trapz, 3)))
print("The value of the integral by Simpson's 1/3rd rule = "+str(round(value_simp, 3)))

The value of the integral by Trapezoidal rule = 5.0
The value of the integral by Simpson's 1/3rd rule = 5.0