How would I calculate the numerical value of this integral / function in python?

47 Views Asked by At

I am trying compute numerically the value of this function:

$ F(x) = x \, \int_x^{+\infty} K_{5/3} (z) \, dz $

Where $ K_{5/3} (z) $ is the modified Bessel function of the second kind and of order $ \frac{5}{3} $. This is the only thing I've been able to come up with:

from scipy.special import kv
from scipy.integrate import trapezoid
import numpy as np
import matplotlib.pyplot as plt

K = kv(5/3, np.linspace(0.1, 10, 1000))

F = []
for x in np.linspace(0, 10, 1000):
    F.append(x * trapezoid(K, np.linspace(x, 100, 1000)))
plt.plot(np.linspace(0, 10, 1000), F)

The result is certainly incorrect. The plot should look like the diagrams in page 23 of this pdf, shown both in linear and logarithmic axes.

I really can't come up with anything else. I think the problem might be that the upper limit in python is not $ +\infty $, but I have found no way to implement integration with infinite limits in python. Any help?

1

There are 1 best solutions below

0
On

I shall give you the analytical solution is order you be able to check your calculations.

$$F(x)=x \left(\frac{27 \sqrt{3} \pi x^{8/3} \, _1F_2\left(\frac{4}{3};\frac{7}{3},\frac{8}{3};\frac{x^2}{4}\right)}{160\ 2^{2/3} \Gamma \left(-\frac{1}{3}\right)}-\frac{2^{2/3} \sqrt{3} \pi \, _1F_2\left(-\frac{1}{3};-\frac{2}{3},\frac{2}{3};\frac{x^2}{4}\right)}{x^{2/3} \Gamma \left(-\frac{2}{3}\right)}-\frac{\pi }{\sqrt{3}}\right)$$

Being blind (and then unable to produce even decent plots), I just give you below a table to check your results $$\left( \begin{array}{cc} x & F(x) \\ 0.00 & 0.000000 \\ 0.25 & 0.915800\\ 0.50 & 0.870819 \\ 0.75 & 0.765267 \\ 1.00 & 0.651423 \\ 1.25 & 0.544867 \\ 1.50 & 0.450640 \\ 1.75 & 0.369775 \\ 2.00 & 0.301636 \\ 2.25 & 0.244926 \\ 2.50 & 0.198145 \\ 2.75 & 0.159812 \\ 3.00 & 0.128566 \\ 3.25 & 0.103203 \\ 3.50 & 0.082687 \\ 3.75 & 0.066140 \\ 4.00 & 0.052827 \\ 4.25 & 0.042139 \\ 4.50 & 0.033574 \\ 4.75 & 0.026722 \\ 5.00 & 0.021248 \end{array} \right)$$