I have written a sparse quadrature routine to integrate multidimensional Gaussian integrals. I'm trying to show convergence plots in my thesis and demonstrate that the code works correctly. I am in the process of benchmarking this using multidimensional integrals with closed forms. With normal numerical schemes one would assume the relative error to decreases continuously a with higher number of quadrature points. However the smolyak algorithm selects only those tensor products obeying the following criteria:
Where $M$ is the number of dimensions in the integral, $\ell$ is the order and $\|i\|_1$ the summation of the number of quadrature points used for each of the $M$ one-dimensional sets. This results in a sparse tensor product. As we can see when large $\ell$ is used the lower order tensor products are dismissed, particularly when $\ell>>M$. Hence would one still expect the error to decrease continuously with a higher order $\ell$? Within my code I am seeing that for higher $\ell$, the error plateaus, as $M$ is increased this plateaus moves such that the accuracy is decreased. For example when $M$ is 1, the error plateaus at 1E-16 i.e. double precision limit as expected. But when $M=5$ the plateau is at around 1E-12, this seem quite large and I've tested this for several different multidimensional integrals. Is this correct? For $\ell=35$ we have
So we are losing a lot of tensor products. Am I being to pedantic? I will only every require the first four significant figures due to input data.