Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question

2 Answers 2

up vote 6 down vote accepted

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.

share|improve this answer
    
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 3: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 in ver. 8. I think you could try some methods in Sum` to work around this issue. –  Artes Mar 16 '13 at 3: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 3:50
    
@Artes: I have acknowledged your assistance in my answer. –  robjohn Mar 16 '13 at 4:02
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

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.