Changing working precision to 10 in `NIntegrate` introducing error messages and giving bigger error than precision 5

Question

I have two numerical integrals as follows:

u[t_, p_] := Cos[t/2]*Exp[I*p/2];
v[t_, p_] := Sin[t/2]*Exp[-I*p/2];
r[t1_, p1_, t2_, p2_] := 2*Abs[u[t1, p2]*v[t2, p2] - u[t2, p2]*v[t1, p1]];

intk2[k_, twoS_, L_, wp_] := 
  NIntegrate[(r[t1, p1, t2, p2]/2)^(2*twoS - 2*L + k)*(2*
       Sqrt[twoS/2])^k*Abs[u[t1, p1]]^(2*L)*Abs[u[t2, p2]]^(2*L)*
    Exp[-Sqrt[twoS/2]*r[t1, p1, t2, p2]]*Sin[t1]*Sin[t2], {t1, 0, 
    Pi}, {t2, 0, Pi}, {p1, 0, 2*Pi}, {p2, 0, 2*Pi}, 
   WorkingPrecision -> wp];

norm[twoS_, L_, wp_] := 
  NIntegrate[(r[t1, p1, t2, p2]/2)^(2*twoS - 2*L)*
    Abs[u[t1, p1]]^(2*L)*Abs[u[t2, p2]]^(2*L)*Sin[t1]*Sin[t2], {t1, 0,
     Pi}, {t2, 0, Pi}, {p1, 0, 2*Pi}, {p2, 0, 2*Pi}, 
   WorkingPrecision -> wp];

I need this for a set of parameters, and when I run for example Table[intk2[i, 9, j, 5], {i, 0, 5}, {j, 0, 5}] and Table[norm[9, j, 5], {j, 0, 5}], with working precision 5, there is no problem. However I can check the accuracy against the exact answer (a long story involving several other functions, see below), and would like it to be better.

So I try to use better working precision and run Table[intk2[i, 9, j, 10], {i, 0, 5}, {j, 0, 5}] and Table[norm[9, j, 10], {j, 0, 5}]. This time however I get error messages like Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

Also I can check that the accuracy is lower than before.

This seems counter intuitive. How can better working precision give worse answers?

.

.

Here's how I can check the answers:

I define

i96 = Table[intk2[i, 9, j, 10], {i, 0, 5}, {j, 0, 5}];
n = Table[norm[9, j, 10], {j, 0, 5}];

and

VLLL[twoS_, J_] := 
  2*Binomial[2*twoS - 2*J, twoS - J]*
   Binomial[2*twoS + 2*J + 2, 
     twoS + J + 1]/(Sqrt[twoS/2]*Binomial[2*twoS + 2, twoS + 1]^2);
V[l_, L_, twoS_] := 
  Sum[(-1)^(twoS + L)*(2*l + 1)^2*SixJSymbol[{L, l, l}, {k, l, l}]*
     ThreeJSymbol[{l, -twoS/2}, {k, 0}, {l, twoS/2}]^2, {k, 0, 2*l}]/
   Sqrt[twoS/2];

and solve the simultaneous equations

Solve[{VLLL[9, 
     0] + (C0*i96[[1, 1]] + C1*i96[[2, 1]] + C2*i96[[3, 1]] + 
       C3*i96[[4, 1]] + C4*i96[[5, 1]] + C5*i96[[6, 1]])/n[[1]] == 
   V[5.5, 0, 9], 
  VLLL[9, 1] + (C0*i96[[1, 2]] + C1*i96[[2, 2]] + C2*i96[[3, 2]] + 
       C3*i96[[4, 2]] + C4*i96[[5, 2]] + C5*i96[[6, 2]])/n[[2]] == 
   V[5.5, 1, 9], 
  VLLL[9, 2] + (C0*i96[[1, 3]] + C1*i96[[2, 3]] + C2*i96[[3, 3]] + 
       C3*i96[[4, 3]] + C4*i96[[5, 3]] + C5*i96[[6, 3]])/n[[3]] == 
   V[5.5, 2, 9], 
  VLLL[9, 3] + (C0*i96[[1, 4]] + C1*i96[[2, 4]] + C2*i96[[3, 4]] + 
       C3*i96[[4, 4]] + C4*i96[[5, 4]] + C5*i96[[6, 4]])/n[[4]] == 
   V[5.5, 3, 9], 
  VLLL[9, 4] + (C0*i96[[1, 5]] + C1*i96[[2, 5]] + C2*i96[[3, 5]] + 
       C3*i96[[4, 5]] + C4*i96[[5, 5]] + C5*i96[[6, 5]])/n[[5]] == 
   V[5.5, 4, 9], 
  VLLL[9, 5] + (C0*i96[[1, 6]] + C1*i96[[2, 6]] + C2*i96[[3, 6]] + 
       C3*i96[[4, 6]] + C4*i96[[5, 6]] + C5*i96[[6, 6]])/n[[6]] == 
   V[5.5, 5, 9]}, {C0, C1, C2, C3, C4, C5}]

I then define the answer {C1,...,C5} as c2 and make the function

ppVeff2[twoS_, L_, c_, K_, wp_] := 
  NIntegrate[((r[t1, p1, t2, p2]/2)^(2*twoS - 2*L - 1)*(2*
           Sqrt[twoS/2])^(-1) + 
       Sum[c[[k + 
           1]]*(r[t1, p1, t2, p2]/2)^(2*twoS - 2*L + k)*(2*
            Sqrt[twoS/2])^k*Exp[-Sqrt[twoS/2]*r[t1, p1, t2, p2]], {k, 
         0, K - 1}])*Abs[u[t1, p1]]^(2*L)*Abs[u[t2, p2]]^(2*L)*
     Sin[t1]*Sin[t2], {t1, 0, Pi}, {t2, 0, Pi}, {p1, 0, 2*Pi}, {p2, 0,
      2*Pi}, WorkingPrecision -> wp]/n[[L + 1]];

Finally I run ppVeff2[9, 0, c2, 6, 5], where I know the answer should be 0.270607. Using the answers with Working precision 10 I find 0.29043, while precision 5 gives 0.27696.

Increasing the number of coefficients in c2={C1,...} would also make this answer better, but I'd like to get as good accuracy as possible from the numerical integration as well.

Answer 1

Setting the option Method -> {"Multidimensional", "Generators" -> 9} takes care of most of the cases, at least as far as the warning messages go. There were a few exceptions, which suggests that intervention on a case by case basis may be necessary if further investigations are needed.

gen[2, 3] = 5;
gen[4, 0] = 6;
gen[4, 4] = 5;
gen[5, 0] = 5;
gen[5, 1] = 5;
gen[5, 3] = 5;
gen[5, 4] = 5;
gen[__] = 9;

intk2[k_, twoS_, L_, wp_] := 
  NIntegrate[(r[t1, p1, t2, p2]/2)^(2*twoS - 2*L + k)*(2*
       Sqrt[twoS/2])^k*Abs[u[t1, p1]]^(2*L)*Abs[u[t2, p2]]^(2*L)*
    Exp[-Sqrt[twoS/2]*r[t1, p1, t2, p2]]*Sin[t1]*Sin[t2], 
   {t1, 0, Pi}, {t2, 0, Pi}, {p1, 0, 2*Pi}, {p2, 0, 2*Pi}, 
   WorkingPrecision -> wp,
   Method -> {"Multidimensional", "Generators" -> gen[k, L]}];

norm[twoS_, L_, wp_] := 
  NIntegrate[(r[t1, p1, t2, p2]/2)^(2*twoS - 2*L)*
      Abs[u[t1, p1]]^(2*L)*Abs[u[t2, p2]]^(2*L)*Sin[t1]*Sin[t2],
   {t1, 0, Pi}, {t2, 0, Pi}, {p1, 0, 2*Pi}, {p2, 0, 2*Pi}, 
   WorkingPrecision -> wp,
   Method -> {"Multidimensional", "Generators" -> If[L == 0, 6, 9]}];

The value returned by ppVeff2[9, 0, c2, 6, 10] is 0.2769071357, a very modest improvement over a working precision of 5.

As to what's going on numerically, the computation is so slow, it discourages tracking down precisely what's going on. High dimensional integrals are particularly hard to evaluate, and four-dimensional ones are hard to visualize. But generally, errors sometimes cancel each other out, sometimes not. A change in the working precision causes a change in the precision goal, which leads to changes in the criteria for convergence as well as in sampling. That it happens once, that a precision of 5 gives a better result than a precision of 10, is not that surprising to me.