pbeckman
commented
1 year ago If we run the following MATLAB code, which is an equivalent test using the implementation of the Hale and Townsend 2012 algorithm in Chebfun, we see much better results, even using many nodes with a singularity $\alpha = -0.99$
a = -0.99;
ns = 2.^ (1:16);
I_true = 2^(a+1)/(a+1);
Is_quad = zeros(length(ns), 1);
for i=1:length(ns)
[~, w] = jacpts(ns(i), a, 0);
Is_quad(i) = sum(w);
end
rel_errs = (I_true - Is_quad) / I_true;n relative error
----------------------
2 -4.516e-15
4 1.552e-15
8 4.145e-13
16 7.042e-13
32 1.437e-12
64 7.099e-12
128 -3.387e-15
256 3.810e-15
512 0.000e+00
1024 -3.528e-15
2048 -1.411e-15
4096 -2.258e-15
8192 4.234e-16
16384 -7.056e-16
32768 4.234e-16
65536 1.694e-15This seems to suggest an issue in the Julia implementation.
hyrodium
dlfivefifty
While using some large Gauss-Jacobi rules, I found that for highly singular integrands ($\alpha < -0.5$ or so), I'm losing digits in the resulting integrals. As a minimal example, consider $$\int_{-1}^1 (1-x)^\alpha dx = \frac{2^{\alpha+1}}{\alpha+1}.$$ For a range of $n$'s, we approximate the above integral by summing the weights of an $n$-point Gauss-Jacobi rule.
which yields the following errors
Is this expected behavior? Or perhaps related to this issue posted by @JamesCBremerJr?
Figure 4.5 in Hale and Townsend 2012 (on which I believe this routine is based?) shows 13 to 14 digits of accuracy in the weights for $\alpha = 2, \beta = -0.75$ up to $10^4$ nodes. But if we rerun the above code with $\alpha = -0.75$, we get only 10 digits in the integral for $n = 2^{16}$.
If we consider larger rules with $\alpha = -0.75$, we get only 8 digits for $n=2^{20}$. And if we take a harsher singularity, e.g. $\alpha = -0.99$, we get only 6 digits for $n=2^{16}$. So it looks like this loss of accuracy eventually (but long before any nodes coincide to machine precision) appears for any $\alpha < -0.5$, and becomes more severe as $\alpha \to -1$. Are such cases known limitations?