0
$\begingroup$

I'm using Mathematica to integrate this function:

w[x_, z_] := 
 E^x/(E^x + 1)^2 Log[(E^(z^2/(4 x)) + E^-x)/(E^(z^2/(4 x)) - 1)]
W[z_] := NIntegrate[w[x, z], {x, 0, \[Infinity]}]

where z > 0 is a positive parameter. Analytically, I know that the integrated is non-negative, so the integral itself should be non-negative. However, I'm getting some incorrect results due to a lack of accuracy:

W[100]

enter image description here enter image description here

-6.05023*10^-20

My question: How can I get a reliable result? Thanks!

$\endgroup$

2 Answers 2

3
$\begingroup$
w[x_, z_] := E^x/(E^x + 1)^2 Log[(E^(z^2/(4 x)) + E^-x)/(E^(z^2/(4 x)) - 1)]

For larger z the integrand is very small and high precision is required in the integration.

Plot3D[w[x, z], {x, 0, 40}, {z, 0, 30}, WorkingPrecision -> 15,
 PlotPoints -> 50,
 MaxRecursion -> 5,
 ClippingStyle -> None]

enter image description here

W2[z_?NumericQ] := NIntegrate[w[x, z], {x, 0, ∞},
  WorkingPrecision -> 40,
  Method -> "LocalAdaptive"]

W2[100]

(* 4.679853458969239635780655689865016458810*10^-43 *)
$\endgroup$
1
  • $\begingroup$ Thank you so much! $\endgroup$ Nov 11, 2021 at 22:33
2
$\begingroup$

Simplifying and rearranging the integrand, along with high precision, solves the problem in this case:

With[{z = 100}, 
     NIntegrate[1/4 Sech[x/2]^2 Log[(1 + E^-x E^(-(z^2/(4 x))))/(1 - E^(-(z^2/(4 x))))],
                {x, 0, ∞}, WorkingPrecision -> 25]]
   4.679853735636909286562544*10^-43

In particular, the important part here was to try to express the integrand in terms of functions that decay at infinity, like $\exp(-x)$ or $\operatorname{sech} x$.

$\endgroup$
1
  • $\begingroup$ Point taken, thanks a lot. $\endgroup$ May 21, 2022 at 22:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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