The integral is probably divergent, or slowly convergent

Question

from scipy import integrate
import numpy as np
from mpmath import coth

DvDc = 6.5
dens = 5.65
Vs = 4.6e3
a = 3.3e-9
w0 = np.sqrt(2)*Vs/a
T = 50 
kb = 1.38064852e-5  # in eV`
j0 = (DvDc)**2 / ((2*np.pi)**2 *dens*Vs**5)

def func(x):
    return x*np.exp(-(x/w0)**2)*coth(x/(2*kb*T))
S = j0*integrate.quad(func, 0, np.inf)[0]
print(S)

Hi, so I have no experience using numerical integration, but I was wondering if anybody can help. Here I have some variables defined, and have a function that I want to integrate. However I get back "The integral is probably divergent, or slowly convergent." And pretty much get the same result despite what I put in the denominator of coth.

Does anybody know if I'm making a mistake or perhaps I just can't do this integral.

Thanks!

Answer 1

For positive <math>w</math> and <math>t</math>
<math>
\int_0^\infty xe^{-x^2/w^2}\coth(x/t) dx
=
t^2
\int_0^\infty xe^{-x^2(t/w)^2}\coth(x) dx</math>

<math>
=
t^2 
\int_0^\infty xe^{-x^2(t/w)^2}\frac{e^{2x}+1}{e^{2x}-1} dx


=
t^2
\int_0^\infty xe^{-x^2(t/w)^2}\frac{1+e^{-2x}}{1-e^{-2x}} dx

=
t^2 
\int_0^\infty xe^{-x^2(t/w)^2}[1+2e^{-2x}+2e^{-4x}+2e^{-6x}+\cdots] dx
.
</math>
As a shortcut $t\approx 0.00138064$, $w\approx 0.19713\times 10^{13}$, $\alpha\equiv w/t\approx 0.14278\times 10^{16}$.
All these integrals can be written in terms of the complementary error function:
<math>

=t^2
\big[
\frac{\alpha^2}{2}
+\sqrt\pi \alpha e^{\alpha^2} \mathrm{erfc} (\alpha)
+\sqrt\pi \alpha e^{4\alpha^2} \mathrm{erfc}(2 \alpha)
+\sqrt\pi \alpha e^{9\alpha^2} \mathrm{erfc}(3 \alpha)+\cdots
\big]
</math>
Because <math>\alpha</math> is large, one can use the factorial series representations
<math>
j \sqrt\pi \alpha e^{j^2\alpha^2} \mathrm{erfc}(j\alpha)\approx   1-\frac{1}{2}\frac{1}{(j^2\alpha^2+1)}
+\frac14\frac{1}{(j^2\alpha^2+1)(j^2\alpha^2+2)}



\approx
t^2
\big[
\frac{\alpha^2}{2}
+1-\frac12\frac{1}{\alpha^2+1}
+\frac12-\frac14\frac{1}{4\alpha^2+1}
\big]

So the result is basically <math>(t\alpha)^2/2=w^2/2 \approx 1.94307\times 10^{24}</math>.

Answer 2

There a some difficulties with your integrand. First, although the function is not divergent, it may look like it for the integration scheme. Just plot it:

That's because your w0 is 1971327996035.2234, so the exp(-x**2) suppression only happens very far to the right. Use variable substitution y = x / w0 to get it right. Second, the function is not defined at x=0 (because you divide by tanh) but can be continuously extended with 2 * kb * T.

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

DvDc = 6.5
dens = 5.65
Vs = 4.6e3
a = 3.3e-9
w0 = np.sqrt(2) * Vs / a
T = 50
kb = 1.38064852e-5  # in eV`
j0 = (DvDc) ** 2 / ((2 * np.pi) ** 2 * dens * Vs ** 5)


def func(x):
    return x * np.exp(-((x / w0) ** 2)) / np.tanh(x / (2 * kb * T))


def func2(y):
    return np.where(
        y > 0,
        y * w0 * np.exp(-(y ** 2)) / np.tanh(y * w0 / (2 * kb * T)),
        2 * kb * T
    )

val, err = integrate.quad(func2, 0, np.inf)
val *= w0  # y = x / w0  =>  dx = dy * w0
print(val)

1.943067033976125e+24

A plot of func2: