7

I try this numerical summation (in two parts)

a = NSum[(HarmonicNumber[2 m])/m^3, {m, 1, Infinity}, 
      WorkingPrecision -> 100, PrecisionGoal -> 100];
N[Pi^3/24 Log[2]^2 + Log[2] Pi/16 Zeta[3] - Pi^5/960 - Pi/16 a, 100]

and this numerical integration

NIntegrate[x^2 Log[Sin[x]] Log[Cos[x]], {x, 0, Pi/2}, 
  WorkingPrecision -> 100, PrecisionGoal -> 100]

which are supposed to give the same result, and they do, but only to 25 places.

Obviously, at least one of the results is off. How can I increase the precision so that these agree past 25 places? If that can't be done, which of these is the more accurate?

If I try to evaluate the first quantity symbolically, I get

b = Sum[(HarmonicNumber[2 m])/m^3, {m, 1, Infinity}]
1/144 (π^4 + 72 (EulerGamma + Log[4]) Zeta[3] - 36 Sqrt[π]  
    (HypergeometricPFQRegularized^({0, 0, 0, 0, 0}, {0, 0, 0, 1}, 0))[
      {1, 1, 1,1, 3/2}, {2, 2, 2, 3/2}, 1])

However, N[b,20] never returns. The problem seems to be the evaluation of the derivative of HypergeometricPFQRegularized.

m_goldberg
  • 107,779
  • 16
  • 103
  • 257
robjohn
  • 1,071
  • 8
  • 18

2 Answers2

7

To make the result of NSum more precise you can use also the NSumTerms option (15 by default, see e.g. Numerical Evaluation of Sums and Products) appropriately increased. Let's try e.g. :

a1 = NSum[ HarmonicNumber[2 m]/m^3, {m, 1, Infinity}, WorkingPrecision -> 140, 
             PrecisionGoal -> 70, NSumTerms -> 2000]
 1.9746275368413284954203787248027995910222173561519748313727983181550691548

now compare

NIntegrate[ x^2 Log[Sin[x]] Log[Cos[x]], {x, 0, Pi/2}, WorkingPrecision -> 75, 
                                                       PrecisionGoal -> 70 ]
N[ Pi^3/24 Log[2]^2 + Log[2] Pi/16 Zeta[3] - Pi^5/960 - Pi/16 a1, 75]

% - %%
 0.0778219793722938643380943991911599389199168078241333818284167516820632583615
 0.07782197937229386433809439919115993891991680782413338182841675168206325836

0.*10^-75

That's pretty close.

Regarding your definition of b I had no problems with evaluating it, e.g.

N[b, 70]
1.974627536841328495420378724802799591022217356151974831372798318155069

You can see that the NSum result is really close to this value.

For some closely related problems with NSum see e.g. this question Precision differences.

Artes
  • 57,212
  • 12
  • 157
  • 245
  • Thank you very much! I knew there was some option that I hadn't found that I needed to set. I've had this problem before and I am glad to know how to fix it. As for the problem evaluating the derivative of HypergeometricPFQRegularized, it might be my version, which is 8.0.1.0. – robjohn Mar 16 '13 at 03:30
  • @robjohn I'm glad I could help. I had no problem in ver. 9.0.1 with N[b, 70], but trying to evaluate it in 8.0.4 its just taken a few minutes and I still have got no result. I'll look later if one can do it successfully inver. 8. I think you could try some methods inSum` to work around this issue. – Artes Mar 16 '13 at 03:46
  • @Artes: I will look at options in Sum to see what I can do. I tried simply evaluating the derivative of HypergeometricPFQRegularized by itself at the given point, and have left it running for over 30 minutes with no response. Thanks for the help. – robjohn Mar 16 '13 at 03:50
  • @Artes: I have acknowledged your assistance in my answer. – robjohn Mar 16 '13 at 04:02
  • This question seems like it is a duplicate of the question you linked to. Do you disagree? – Mr.Wizard Mar 16 '13 at 06:27
  • 2
    @Mr.Wizard Clearly we encounter the same issue here and to provide a higher precision we use the same option. On the other hand I haven't marked it as an exact duplicate since there is some kind of bug in ver.8 not present in ver.9, see comments above. I'll have to take a closer look if it can be explained to a larger extent and then I'll update my answer. – Artes Mar 16 '13 at 12:30
2

Of course, one should also remember that the Method options of NSum[] accept sub-options as well. For instance,

NSum[HarmonicNumber[2 m]/m^3, {m, 1, Infinity}, 
     Method -> {"EulerMaclaurin", "ExtraTerms" -> 50, 
                Method -> {NIntegrate, Method -> "DoubleExponential"}}, 
     NSumTerms -> 50, PrecisionGoal -> 90, VerifyConvergence -> False, 
     WorkingPrecision -> 120]
   1.9746275368413284954203787248027995910222173561519748313727983181550691

In particular, when using the Euler-Maclaurin method, I have found it helpful to use the double exponential quadrature scheme of Takahasi and Mori for evaluating the needed (often improper) integrals; I have thus done so here.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574