Asymptotic Computation of Nodes and Weights

Code development
1

I have recently been working on implementing an asymptotic formula to help calculate the nodes and weights for the numerical evaluation of integrals.

After implementing Gauss-Legendre, the benefits have become clear, with runtime decreasing several fold when considering higher orders (n>10,000)

Code is present in the Higher Order Integration Repo within the Equadratures organization Github Page.

Currently working on graphically displaying error in calculating nodes and weights when considering different functions.

Below is the integration error for the following function, with the first graph being the Asymptotic computation, and the second using Equadratures.

F(x) = \int_{-1}^{1}f(x)dx=\int_{-1}^{1}xsin(x^2)dx \\ F(x) = \Bigg[-\frac{cos(x^2)}{2}\Bigg]_{-1}^{\;\,1} = 0

True value of F=0

image
image1216×1865 80.3 KB

As evident by the scale of the graphs, the asymptotic implementation reduces error approximately 10 fold, and this is already present with n <= 10,000. The asymptotic implementation allowed me to compute up to orders as high as 1,000,000 with a relatively short runtime. But, I was unable to compute the error using Equadratuers using higher order as the runtime was scaling exponentially, it is likely however that the error would increase even further as n increases towards 10e6.

Current Code for computing the nodes using the asymptotic implementation is shown here, with the weights and other functions available within the Higher Order Integration Repo linked below.

def legpts_nodes(a, u, a2, u2, n, vn):

F0 = a
F1 = 1/8 * (u * a - 1) / a

if n < 1e4:
    a3 = a ** 3
    F2 = 1/384 * (6 * a2 * (1 + u2) + 25 - u * (31 * u2 + 33) * a3) / a3
else:
    F2 = 0

if n < 1e3:
    u4 = u ** 4
    a5 = a ** 5
    R30 = u * (2595 + 6350 * u2 + 3779 * u4) / 15360
    R31 = -(31 * u2 + 11) / 1024
    R32 = u / 512
    R33 = -25 / 3072
    R35 = -1073 / 5120
    F3 = R30 + R35 / a5 + (1 + u2) * (R31 / a + R32 / a2 + R33 / a3)
else:
    F3 = 0

t = F0 + F1 * vn ** 2 + F2 * vn ** 4 + F3 * vn ** 6
x = np.cos(t)

return x, t

It uses terms up to order 3, which are calculated using the following formulas presented in “Iteration-Free Computation of Gauss-Legendre Quadrature Nodes and Weights” as referenced below.

\begin{aligned} F_0(x,u) &= x \\ F_1(x,u) &= \frac{1}{8}\frac{ux-1}{x} \\ F_2(x,u) & =\frac{1}{384}\frac{6x^2(1+u^2)+25-u(31u^2+33)x^3}{x^3} \end{aligned}

Similarly the weights are derived from the same paper, and up to order 2 are presented below:

\begin{aligned} W_0(x,u) &= 1 \\ W_1(x,u) &= \frac{1}{8}\frac{ux+x^2-1}{x^2} \\ W_2(x,u) & =\frac{1}{384}\frac{81-31ux-(3-6u^2)x^2+6ux^3-(27 + 84u^2 + 56u^4)x^4}{x^4} \end{aligned}

The current implementations for the asymptotic computation of nodes and weights only use terms up to the third order as previously stated (one more than shown in the above formulas as it is too big to include), however, data present within the referenced paper suggests there is merit in including higher order terms, potentially reducing error by many orders of magnitude.

Referenced Paper:
https://core.ac.uk/display/55887025?utm_source=pdf&utm_medium=banner&utm_campaign=pdf-decoration-v1
Higher Order Integration Repo:

2 Likes
2

Hi @RoshanAmir, what about when you are trying to integrate

\int f(x) \rho(x) dx

where \rho(x) is a Beta (Jacobi measure) distribution?