How to intergrate a funtion with high accuracy

Question

I am trying to integrate the product of functions of sin and cos. The results I am getting match when the result has a relatively large magnitude (e.g. 1e-2~1e-4). But the results with small magnitude (e.g. 1e-11~1e-32) have large relative error and are sometimes of opposite sign.

import numpy as np
import scipy.integrate as spi

# Define the integrand function
def f0(x,n):
    if n == 1:
        return (1-x)
    elif n == 2:
        return x
    else:
        return np.sin((n - 2) * np.pi * x)


def f00(x, m, n):
    return f0(x,m)*f0(x,n)

# Define the dimensions of the tensor
NMHT_0 = 12

# Create the tensor
I00 = np.zeros((NMHT_0, NMHT_0))

# Perform the integration
for m in range(NMHT_0):
    for n in range(NMHT_0):
        result1, error1 = spi.quad(f00, 0, 1, args=(m, n))
        
        
        I00[m, n] = result1      

I tried to integrate first with the scipy.integrate.quad_vec function, but all the small numbers are not matching. I replaced the quad_vec function with the scipy.integrate.quad function, and it managed to get some of the small numbers to match, but the majority of the small numbers still do not match. Any way to get more accurate results, or should I try a different library?

Answer 1

The results are already accurate. You could try using the epsrel and epsarg arguments of scipy.integrate.quad to try to get those small magnitude numbers closer to zero, but they are theoretically zero, so what you're seeing is a little numerical noise.

As suggested in the comments, mpmath.quad can be used for arbitrary precision numerical integration in Python (when you need it).

import numpy as np
from mpmath import mp
mp.dps = 100  # set the precision

# Define the integrand function
def f0(x, n):
    if n == 1:
        return (mp.one - x)
    elif n == 2:
        return x
    else:
        return mp.sin((n - 2) * mp.pi * x)

def f00(x, m, n):
    return f0(x, m) * f0(x, n)

# Define the dimensions of the tensor
NMHT_0 = 12

# Create the tensor
I00 = np.zeros((NMHT_0, NMHT_0))

# Perform the integration
for m in range(NMHT_0):
    for n in range(NMHT_0):        
        I00[m, n] = mp.quad(lambda x: f00(x, m, n), (0, 1), method='tanh-sinh')
        # I00[m, n] = mp.quad(lambda x: f00(x, m, n), (0, 1), method='gauss-legendre')

To be more confident that the precision is high enough, you can take a look for agreement between the two methods of mpmath.quad and check for convergence as mp.dps is increased.

But as @jared notes, this does not require numerical integration. We can use sympy.integrate (or look at a calculus textbook) to see that those tiny values are really zero.

from sympy import integrate, symbols, sin, pi
m, n = symbols('m, n', positive=True, integer=True)
x = symbols('x')
integrate(sin(m*pi*x) * sin(n*pi*x), (x, 0, 1))