Compare numerial integral calculation#

quadrature
trapezoidal fomula
simpson’s rule

[1]:

import numpy as np

[2]:

def quad(func, lim=(0, 1), num=100, cof=1.0):
    """quadrature integral approximation"""
    x = np.linspace(*lim, num + 1) * cof
    k = np.ones(x.size) * (x[1] - x[0])
    return k.dot(func(x))


def trape(func, lim=(0, 1), num=100, cof=1.0):
    """trapezoidal integral approximation"""
    x = np.linspace(*lim, num + 1) * cof
    k = np.ones(x.size) * (x[1] - x[0])
    k[0] *= 0.5
    k[-1] *= 0.5
    return k.dot(func(x))


def simp(func, lim=(0, 1), num=100, cof=1.0):
    """simpson integral approximation"""
    x = np.linspace(*lim, num + 1) * cof
    k = np.ones(x.size) * (x[1] - x[0]) / 3.0
    k[1:-1:2] *= 4.0
    k[2:-2:2] *= 2.0
    return k.dot(func(x))


def sin_integ(lim=(0, 1), cof=1.0):
    """sin integral"""
    return np.cos(cof * lim[0]) - np.cos(cof * lim[1])


def error(ref, true):
    """relative error"""
    return 100 * abs(ref - true) / true

Numerical Integration of \(f(x) = x^2\)#

The number of descritazation is 100 points along to x axis (0, 1)

[3]:

true = 1.0 / 3.0  # analytical solution
def func(x):
    return x ** 2
qua = quad(func)
trap = trape(func)
simpson = simp(func)

[4]:

print("quadrature relative error  : {:.2e} %".format(error(qua, true)))
print("trapezoidal relative error : {:.2e} %".format(error(trap, true)))
print("simpsons relative error    : {:.2e} %".format(error(simpson, true)))

quadrature relative error  : 1.51e+00 %
trapezoidal relative error : 5.00e-03 %
simpsons relative error    : 3.33e-14 %

Numerical Integration of \(f(x) = \sin x\)#

The number of descritazation is 100 points along to x, the range of which is (0, \(\pi\))

[5]:

true = 2
def func(x):
    return np.sin(x)
qua = quad(func, lim=(0, np.pi))
trap = trape(func, lim=(0, np.pi))
simpson = simp(func, lim=(0, np.pi))

[6]:

print("quadrature relative error  : {:.2e} %".format(error(qua, true)))
print("trapezoidal relative error : {:.2e} %".format(error(trap, true)))
print("simpsons relative error    : {:.2e} %".format(error(simpson, true)))

quadrature relative error  : 8.22e-03 %
trapezoidal relative error : 8.22e-03 %
simpsons relative error    : 5.41e-07 %