Combining Lebedev and Gauss chebyshev quadratures

Question

My problem is that I need to compute 3D integrals, and I need to combine a Gauss-Chebyshev (for the radial part + variable substitution) and a Lebedev quadrature (for the solid angle/unit sphere surface part).

The first part is to process the quadrature points :

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
pi = np.pi

##########################
# Quadrature weights/points
##########################

# for radial component
Nr = 30
ii  = np.arange(1,Nr+1).reshape(-1,1) 
sin_i = np.sin(ii*pi/(Nr+1))
cos_i = np.cos(ii*pi/(Nr+1))

# Chebyshev quadrature points on [-1,1]
Mu_quadrature = (Nr+1.0-2.0*ii)/(Nr+1.0) + 2.0/pi*(1.0 + 2.0/3.0*sin_i**2)*cos_i*sin_i

# Chebyshev quadrature points on [0, +infty [
# Bragg-slater Radius for Hydrogen
rm = 0.35 # Angstrom
# variable substitution for the radial component 
R_quadrature = rm*(1+Mu_quadrature)/(1-Mu_quadrature)

# Chebyshev quadrature weights
Wr_quadrature = 16.0/3.0/(Nr+1)*sin_i**4

##########################
# Quadrature for solid angle / angular coordinates
##########################

# file containing integration points (theta/phi in deg) and associated weights
fn = 'lebedenev_quadrature.txt'
Data = np.loadtxt(fn)
ThetaPhi_quadrature = Data[:,:2]

# convert to rad
ThetaPhi_quadrature = ThetaPhi_quadrature*pi/180.0
Wangular_quadrature = Data[:,2:3]

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
        np.tile(R_quadrature,(len(ThetaPhi_quadrature),1)),
        np.repeat(ThetaPhi_quadrature,len(R_quadrature),axis=0)   
    ],
    axis=1
)

R_int   = X_spherical_int[:,0:1]
Theta_int = X_spherical_int[:,1:2]
Phi_int = X_spherical_int[:,2:3]

Mu_int  = (R_int - rm)/(R_int + rm)

Wint = np.concatenate([
        np.tile(Wr_quadrature,(len(Wangular_quadrature),1)),
        np.repeat(Wangular_quadrature,len(Wr_quadrature),axis=0)   
    ],
    axis=1
)
Wint = np.prod(Wint,axis=1,keepdims=True)

X_cartesian_int = np.zeros_like(X_spherical_int)
R, Theta, Phi = R_int, Theta_int, Phi_int
Mu = (R - rm)/(R + rm)
X_cartesian_int[:,0:1] = R*np.cos(Theta)*np.sin(Phi)
X_cartesian_int[:,1:2] = R*np.sin(Theta)*np.sin(Phi)
X_cartesian_int[:,2:3] = R*np.cos(Phi)

The file containing the Lebedev quadrature points can be found here. At the bottom of the page as text files.

Then comes the tests, I tested the Chebyshev quadrature with the Gauss integral :

#################################
# Testing integration
#################################
def f(X):
    return np.exp(-np.sum(X,axis=1,keepdims=True)**2)

# solution
I1 = np.sqrt(pi)/2

# quadrature
I2 = np.einsum(
    'ij,ij',
    f(R_quadrature)*2.0*rm/(1-Mu_quadrature)**2,
    Wr_quadrature
)

print(f'Testing Gauss-Chebyshev quadrature with Gauss Integral :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

In case, you are woundering with there is the

2.0*rm/(1-Mu_quadrature)**2

This is because the Chebyshev quadrature is on $\mu \in [-1,1]$, the actual integration on $r\in R⁺$ and the variable substitution is $r = r_m \frac{1+\mu}{1-\mu}$

The test of Lebedev with the surface of the unit sphere :

def f(Theta,Phi):
    return np.ones_like(Theta)

# solution
I1 = 4*pi #np.exp(-1)*4*pi

# quadrature
I2 = 4.0*pi*np.einsum(
    'ij,ij',
    f(ThetaPhi_quadrature[:,1:2],None),
    Wangular_quadrature
)
 
print(f'Testing Ledenev quadratude with Integral over unit sphere surface :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

Now I am trying to compute the volume of a sphere, but with no good results :

#""" Unit Sphere
def f(X):
    R = norm(X,axis=1,keepdims=True)
    return np.where(R<=1,1,0)
I1 = 4/3*pi
#"""

""" Testing with Gauss integral
def f(X):
    return np.exp(-norm(X,axis=1,keepdims=True)**2)
I1 = 4*pi*np.sqrt(pi)/2
#"""

I2 = 4*pi*np.einsum(
    'ij,ij',
    f(X_cartesian_int)*2.0*rm/(1-Mu_int)**2,
    Wint
)

print(f'Testing Ledenv + Chebyshev qudrature with Integral over unit sphere volume :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

The output is :

Testing Gauss-Chebyshev quadrature with Gauss Integral :
 I1 = 0.8862269254527579 | I2 = 0.8862570079725882 | I1-I2 = -3.0082519830276766e-05

Testing Ledenev quadratude with Integral over unit sphere surface :
 I1 = 12.566370614359172 | I2 = 12.566370614362953 | I1-I2 = -3.780087354243733e-12

Testing Ledenv + Chebyshev qudrature with Integral over unit sphere volume :
 I1 = 4.1887902047863905 | I2 = 12.623540989301834 | I1-I2 = -8.434750784515444

I also tried to only use spherical coordinates, but it did not work.

The full code :

import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import norm
pi = np.pi

##########################
# Quadrature weights/points
##########################

# for radial component
Nr = 30
ii  = np.arange(1,Nr+1).reshape(-1,1) 
sin_i = np.sin(ii*pi/(Nr+1))
cos_i = np.cos(ii*pi/(Nr+1))

# Chebyshev quadrature points on [-1,1]
Mu_quadrature = (Nr+1.0-2.0*ii)/(Nr+1.0) + 2.0/pi*(1.0 + 2.0/3.0*sin_i**2)*cos_i*sin_i

# Chebyshev quadrature points on [0, +infty [
# Bragg-slater Radius for Hydrogen
rm = 0.35 # Angstrom
# variable substitution for the radial component 
R_quadrature = rm*(1+Mu_quadrature)/(1-Mu_quadrature)

# Chebyshev quadrature weights
Wr_quadrature = 16.0/3.0/(Nr+1)*sin_i**4

##########################
# Quadrature for solid angle / angular coordinates
##########################

# file containing integration points (theta/phi in deg) and associated weights
fn = 'lebedenev_quadrature.txt'
Data = np.loadtxt(fn)
ThetaPhi_quadrature = Data[:,:2]

# convert to rad
ThetaPhi_quadrature = ThetaPhi_quadrature*pi/180.0
Wangular_quadrature = Data[:,2:3]

# make a tensor product of the Quadrature points
X_spherical_int = np.concatenate([
        np.tile(R_quadrature,(len(ThetaPhi_quadrature),1)),
        np.repeat(ThetaPhi_quadrature,len(R_quadrature),axis=0)   
    ],
    axis=1
)

R_int   = X_spherical_int[:,0:1]
Theta_int = X_spherical_int[:,1:2]
Phi_int = X_spherical_int[:,2:3]

Mu_int  = (R_int - rm)/(R_int + rm)

Wint = np.concatenate([
        np.tile(Wr_quadrature,(len(Wangular_quadrature),1)),
        np.repeat(Wangular_quadrature,len(Wr_quadrature),axis=0)   
    ],
    axis=1
)
Wint = np.prod(Wint,axis=1,keepdims=True)

X_cartesian_int = np.zeros_like(X_spherical_int)
R, Theta, Phi = R_int, Theta_int, Phi_int
Mu = (R - rm)/(R + rm)
X_cartesian_int[:,0:1] = R*np.cos(Theta)*np.sin(Phi)
X_cartesian_int[:,1:2] = R*np.sin(Theta)*np.sin(Phi)
X_cartesian_int[:,2:3] = R*np.cos(Phi)

#################################
# Testing integration
#################################
def f(X):
    return np.exp(-np.sum(X,axis=1,keepdims=True)**2)

# solution
I1 = np.sqrt(pi)/2

# quadrature
I2 = np.einsum(
    'ij,ij',
    f(R_quadrature)*2.0*rm/(1-Mu_quadrature)**2,
    Wr_quadrature
)

print(f'Testing Gauss-Chebyshev quadrature with Gauss Integral :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

def f(Theta,Phi):
    return np.ones_like(Theta)

# solution
I1 = 4*pi #np.exp(-1)*4*pi

# quadrature
I2 = 4.0*pi*np.einsum(
    'ij,ij',
    f(ThetaPhi_quadrature[:,1:2],None),
    Wangular_quadrature
)
 
print(f'Testing Ledenev quadratude with Integral over unit sphere surface :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')

#""" Unit Sphere
def f(X):
    R = norm(X,axis=1,keepdims=True)
    return np.where(R<=1,1,0)
I1 = 4/3*pi
#"""

""" Testing with Gauss integral
def f(X):
    return np.exp(-norm(X,axis=1,keepdims=True)**2)
I1 = 4*pi*np.sqrt(pi)/2
#"""

I2 = 4*pi*np.einsum(
    'ij,ij',
    f(X_cartesian_int)*2.0*rm/(1-Mu_int)**2,
    Wint
)

print(f'Testing Ledenv + Chebyshev qudrature with Integral over unit sphere volume :\n I1 = {I1} | I2 = {I2} | I1-I2 = {I1-I2}\n')
```