Simpson's integration on lograthmic radial grid

Simpson’s rule integration gives more accurate integration results fo data set’s that contains smooth curves using quadratic functions and some simple algebra. Here, I will show that on lograthmic grid, one can easily use the simpson’s rule with a modified weight factor which turns out to be the position of the grid points.

First, assume we have a bunch of data on a lograthmic radial grid and the grid can be expressed as:

\[R_i = R_0 e^{i \cdot dx},\]

where $i$ starts from 1 to N and $dx$ is the constant step size used to generate this radial grid. Converting to the continuous form:

\[R(x)= R_0 e^{x},\]

On this grid, we have a function $f(R)$ and we want to integrate it from $R=a$ to $R=b$ to get the area underneath $A$:

\[\begin{aligned} A &= \int_a^b f(R) dR \\ &= \int_a^b f(R) d(R_0 e^{x}) \\ &= \int_{ln(\frac{a}{R_0})}^{ln(\frac{b}{R_0})} [f(R) R] dx \\ \end{aligned}\]

Converting to discrete form:

\[\begin{aligned} A &=\int_a^b (f(R) R) dx \\ & = \sum_{i=0}^N f(R_i) R_i dx \end{aligned}\]

Since $dx$ is constant, this discrete form of integration can be performed using the simpson’s rule by treating $f(R_i)R_i$ as $F_i$. Let’s say we want to integrate $F_i$ to get the area betwee points $i=0,2$ (with an aditional point $i=1$ in the middle), simpson’s rule tells us:

\[A_{1} = \frac{dx}{3}\left[F_0+4 F_1+F_2\right]\]

or more generally:

\[A_{i} = \frac{dx}{3}\left[F_{i-1}+4 F_i+F_{i+1}\right]\]

And if we want to integrate the whole dataset, we can use:

\[\text{Sum} = \sum_{i=1,3,5,7,...} A_i\]

In code, assume we have a data set containing the position $R$ and the value of the function $f$ at these points, we can get the integration of this whole dataset using the following code:

def simpson(data):
    # note this code only works with odd number of grid points
    assert len(data)%2 == 1
    # note: on log grid, R_i/R_{i-1} = exp(dx)
    dx = np.log(data[1,0]/data[1-1,0])
    rab = data[:,0]*dx/3
    rho = data[:,1]
    f3 = rho[0]*rab[0]
    mesh = len(rho)
    asum = 0
    
    for i in range(1,mesh,2):
        f1 = f3
        f2 = rho[i] * rab[i]
        f3 = rho[i + 1] * rab[i + 1]
        asum = asum + f1 + 4.0 * f2 + f3
    
    return(asum)

data=np.loadtxt('AEWFC_Mg_3s.dat')
# we are using atomic wfc, its squared should be normalized
data[:,1] = data[:,1]**2
print(simpson(data))

This can be generalized to any grid that can be expressed using a uniform grid with transformation $R(x)$:

\[\begin{aligned} A &=\int_a^b (f(R) R') dx \\ & = \sum_{i=0}^N f(R_i) R'_i dx \end{aligned}\]

As long as we have the derivative of $R$, we can use the simpson’s rule to get the integral. Here, because $R$ is the exponential function, we get exactly the same $R$ back for $R’$.

A more general simpson’s rule of integration that allows uneven grid spacing (even irregular grid) can be found here:

\[\int_a^b f(x) d x=\sum_{i=0}^{N / 2-1} \frac{h_{2 i}+h_{2 i+1}}{6}\left[\left(2-\frac{h_{2 i+1}}{h_{2 i}}\right) f_{2 i}+\frac{\left(h_{2 i}+h_{2 i+1}\right)^2}{h_{2 i} h_{2 i+1}} f_{2 i+1}+\left(2-\frac{h_{2 i}}{h_{2 i+1}}\right) f_{2 i+2}\right]\]

However, this formulation is not exact equivalent with our previous result. Comparing the weight of the odd terms, our previous gives:

\[4*dx/3*[R_0*e^{dx}]\]

while this formula ($h_i = R_{i+1} - R_i$) gives:

\[R_0*(e^{2dx}-1)/6*(1+e^{dx})^2/e^{dx}\]

Analysis shows that they agree to the second order when $dx \to 0$.

The code to use this simpson’s rule:

import numpy as np

def simpson_uneven(data):
    # note this code only works with odd number of grid points
    assert len(data)%2 == 1
    # generate spacing
    h = [data[i,0]-data[i-1,0] for i in np.arange(1,len(data[:,0]))]
    mesh = len(data)
    asum = 0
    rho = data[:,1]
    
    for i in range(0,int((mesh-1)/2)):
        asum += (h[2*i]+h[2*i+1])/6 \
                * ((2-h[2*i+1]/h[2*i])*rho[2*i] \
                +(h[2*i]+h[2*i+1])**2/(h[2*i]*h[2*i+1])*rho[2*i+1] \
                +(2-h[2*i]/h[2*i+1])*rho[2*i+2])
    
    return(asum)

data=np.loadtxt('AEWFC_Mg_3s.dat')
# we are using atomic wfc, its squared should be normalized
data[:,1] = data[:,1]**2
print(simpson_uneven(data))

Finally, the normal integration method which assume constant $F$ for each intevral is calculated by:

import numpy as np

def step_intgration(data):
    drgrid=[data[i,0]-data[i-1,0] for i in np.arange(1,len(data[:,0]))]
    drgrid=np.insert(drgrid,0,data[0,0])
    integrate = drgrid*(data[:,1])
    result=integrate.sum()
		return(result)

data=np.loadtxt('AEWFC_Mg_3s.dat')
# we are using atomic wfc, its squared should be normalized
data[:,1] = data[:,1]**2
print(step_intgration(data))

These three methods give answers to a noramlized data set:

simpson:            0.9862851076279023
simpson_uneven:     0.9862850228267516
normal integration: 0.9721934311355657

Clearly, we see that the simpson based method is more accurate than the normal step like integration method while the two simpson’s methods differ from the seventh place after the decimal point for this dataset.

The simpson’s method is implemented in VASP’s radial.F and QE’s upflib/simpsn.f90.

The dataset (AEWFC_Mg_3s.dat) I’m using can be downloaded here.