Is there a simple way to compute the following integral
$$I(y,\beta,\gamma)=\int_1^y \frac{x^2}{\sqrt{\left(x^2-1\right)\left(x^2-\gamma ^2\right) \left(x^2+\beta ^2\right) }}\mathrm{d}x$$
with $y>1$, $0<\gamma<1$ and $\beta>0$?
Integral involving ratios of polynomial and square roots
definite integralsimproper-integralsintegration
Related Solutions
3.471.9
I will start with 3.471.9, since it is much simpler. The key ingredient
Slater theorem, saying that $$ \int_0^\infty t^{\alpha-1} f(t) g\left(\frac{x}{t}\right) \mathrm{d} t = \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma + i \infty} \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) x^{-s} \mathrm{d} s = \mathcal{M}^{-1} \left\{ \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) \right\}(x) $$ where $\mathcal{M}_f(s) = \int_0^\infty t^{s-1} f(t) \mathrm{d} t$ is the Mellin transform of $f$, valid within some strip $\min < \Re(s) < \max $, and $\gamma$ is a real constant, such that $\gamma$ is within the strip of validity of the Mellin transform of $g$, and $\gamma + \Re(\alpha)$ in within the strip of the Mellin transform of $f$.
The proof is easy: $$ \begin{eqnarray} \int_0^\infty t^{\alpha-1} f(t) g\left(\frac{x}{t}\right) \mathrm{d} t &=& \int_0^\infty t^{\alpha-1} f(t) \left( \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \mathcal{M}_g(s) x^{-s} t^s \mathrm{d} s \right) \mathrm{d} t \\ &=& \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \left( \int_0^\infty t^{s+\alpha-1} f(t) \mathrm{d} t \right) \mathcal{M}_g(s) x^{-s} \mathrm{d} s \\ &=& \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+\infty} \mathcal{M}_f(s+\alpha) \mathcal{M}_g(s) x^{-s} \mathrm{d} s \end{eqnarray} $$
We now apply this theorem to 3.471.9, using Cahen-Mellin integral $\mathcal{M}_{\exp(-\bullet)}(s) = \Gamma(s) $: $$\begin{eqnarray} \int_0^\infty t^{\alpha-1} \exp\left(-y t\right) \exp\left( - \frac{x}{t}\right) \mathrm{d} t &=& \frac{1}{2 \pi i} \int_{\gamma - i \gamma}^{\gamma + i \infty} \Gamma(s+\alpha) \Gamma(s) y^{-\alpha -s} x^{-s} \mathrm{d}s \\ &=& y^{-\alpha} \left( \frac{1}{2 \pi i} \int_{\gamma - i \gamma}^{\gamma + i \infty} \Gamma(s+\alpha) \Gamma(s) \left( x y \right)^{-s} \mathrm{d}s \right) \\ &=& y^{-\alpha} \left( 2 (x y)^{\alpha/2} K_{\alpha}(2 \sqrt{x y} ) \right) \\ &=& 2 \left( \frac{x}{y} \right)^{\alpha/2} K_\alpha\left(2 \sqrt{x y}\right) \end{eqnarray} $$ where DMLF 10.43.19 was used.
6.635.3
A formula of kin is worked out in Bateman, Erdelyi et al, "Higher Transcendental Functions", chapter 7, section 7.7.6 on Macdonald's and Nochilson's formulas, formula (37). It says:
$$
\int_0^\infty \exp\left(-\frac{t}{2} - \frac{x^2+X^2}{2 t} \right) I_\nu \left( \frac{x X}{t} \right) \frac{\mathrm{d} t}{t} = \cases{2 I_\nu\left( x \right) K_\nu\left( X \right) & x < X \\ 2 I_\nu\left( X \right) K_\nu\left( x \right) & x > X }
$$
By doing an analytic continuation $x \to \mathrm{e}^{i \pi/2} x$ of the (the only possible) first branch, we get
$$
\int_0^\infty \exp\left(-\frac{t}{2} - \frac{X^2-x^2}{2 t} \right) J_\nu \left( \frac{x X}{t} \right) \frac{\mathrm{d} t}{t} = 2 J_\nu(x) K_\nu(X)
$$
Now, performing reparameterization $X = \sqrt{2 \alpha \left( \sqrt{\beta^2 + \gamma^2}+ \beta \right)}$, $x = \sqrt{2 \alpha \left( \sqrt{\beta^2 + \gamma^2} - \beta \right)}$, accompanied by a change of variables $t \to \left(2 \alpha t\right)^{-1}$ we arrive at 6.635.3.
Derivation of the point of departure formula relies on two auxiliary results: $$ \int_0^\infty J_{\nu}(x u) J_{\nu}(X u) \exp\left(-\frac{t u^2}{2}\right) u \mathrm{d} u = \frac{1}{t} \exp\left(-\frac{x^2+X^2}{2t}\right) I_\nu\left( \frac{x X}{t}\right) $$ which is obtained by expanding the exponential in series and integrating term-wise. Then, using the above integral representation in the formula under consideration, and carrying out simple integration with respect to $t$, we get $$ \int_0^\infty \frac{\mathrm{e}^{-t/2}}{t} \exp\left(-\frac{x^2+X^2}{2t}\right) I_\nu\left( \frac{x X}{t}\right) \mathrm{d} t = \int_0^\infty J_\nu(x u) J_\nu(X u) \frac{2 u \mathrm{d} u}{1+u^2} $$ The resulting integral is of Sonine-Gegenabuer-type, and equals to $2 I_\nu( \min(x,X)) K_\nu(\max(x,X))$.
I do not know in which context you faced this function (for mutual conveniency, I reuse the formula as @Mycroft wrote it). $$A(a) = \frac{2 a(a-1) K\left(\frac{a+1}{2 a+1}\right)+4(2 a+1) E\left(\frac{a+1}{2 a+1}\right)}{3 \sqrt{a (2 a+1)}}$$ We faced a very similar one years ago in thermodynamics for $a>1$ and, for obvious computing reasons, we develop it as a series $$A(a)=\sum_{n=0}^p \frac{\alpha_n}{\beta_n} a^{1-n}$$ The table below reproduces the values for the first $n$'s $$\left( \begin{array}{ccc} n & \alpha_n & \beta_n \\ 0 & \sqrt{2} \left(24 \pi ^2+\Gamma \left(\frac{1}{4}\right)^4\right) & 16 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right)^2 \\ 1 & \sqrt{2} \Gamma \left(\frac{1}{4}\right) \Gamma \left(\frac{5}{4}\right)^2+\pi \Gamma \left(\frac{7}{4}\right) & 2 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right) \\ 2 & -\sqrt{2} \left(3 \pi ^2+40 \Gamma \left(\frac{1}{4}\right) \Gamma \left(\frac{5}{4}\right)^3\right) & 96 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right)^2 \\ 3 & 5 \sqrt{2} \Gamma \left(\frac{1}{4}\right) \Gamma \left(\frac{5}{4}\right) & 512 \sqrt{\pi } \\ 4 & 3 \sqrt{2} \left(14 \pi ^2-15 \Gamma \left(\frac{1}{4}\right)^3 \Gamma \left(\frac{5}{4}\right)\right) & 10240 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right)^2 \\ 5 & \sqrt{2} \left(-77 \pi ^2+45 \Gamma \left(\frac{1}{4}\right)^3 \Gamma \left(\frac{5}{4}\right)\right) & 20480 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right)^2 \\ 6 & 3 \sqrt{2} \left(154 \pi ^2-65 \Gamma \left(\frac{1}{4}\right)^3 \Gamma \left(\frac{5}{4}\right)\right) & 163840 \sqrt{\pi } \Gamma \left(\frac{1}{4}\right)^2 \end{array} \right)$$
A few results for $a=10^k$
$$\left( \begin{array}{ccc} k & \text{approximation} & \text{exact} \\ 0 & 2.91216151929910 & 2.91258419032827 \\ 1 & 10.3487069573255 & 10.3487069573835 \\ 2 & 88.9615534510584 & 88.9615534510584 \\ 3 & 875.573856650175 & 875.573856650175 \\ 4 & 8741.74602244500 & 8741.74602244500 \\ 5 & 87403.4726014927 & 87403.4726014927 \\ 6 & 874020.738884157 & 874020.738884157 \\ 7 & 8740193.40176002 & 8740193.40176002 \end{array} \right)$$
Best Answer
Define the function $\mathcal{I}:\left(0,1\right)\times\mathbb{R}_{>0}\times\left(1,\infty\right)\rightarrow\mathbb{R}$ via the improper integral
$$\mathcal{I}{\left(p,q,y\right)}:=\int_{1}^{y}\mathrm{d}x\,\frac{2x^{2}}{\sqrt{\left(x^{2}-1\right)\left(x^{2}-p^{2}\right)\left(x^{2}+q^{2}\right)}}.$$
Suppose $\left(p,q,y\right)\in\left(0,1\right)\times\mathbb{R}_{>0}\times\left(1,\infty\right)$, and set $z:=y^{2}\land a:=1\land b:=p^{2}\land c:=0\land d:=-q^{2}$.
Note that $z>a>b>c>d$, and then also
$$0<\frac{\left(a-d\right)\left(b-c\right)}{\left(a-c\right)\left(b-d\right)}<1\land0<\frac{\left(b-d\right)\left(z-a\right)}{\left(a-d\right)\left(z-b\right)}<1.$$
The integral $\mathcal{I}$ can be transformed into an elliptic integral by means of a simple power substitution:
$$\begin{align} \mathcal{I}{\left(p,q,y\right)} &=\int_{1}^{y}\mathrm{d}x\,\frac{2x^{2}}{\sqrt{\left(x^{2}-1\right)\left(x^{2}-p^{2}\right)\left(x^{2}+q^{2}\right)}}\\ &=\int_{1}^{y^{2}}\mathrm{d}t\,\frac{t}{\sqrt{t}\sqrt{\left(t-1\right)\left(t-p^{2}\right)\left(t+q^{2}\right)}};~~~\small{\left[x=\sqrt{t}\right]}\\ &=\int_{a}^{z}\mathrm{d}t\,\frac{t}{\sqrt{\left(t-a\right)\left(t-b\right)\left(t-c\right)\left(t-d\right)}}.\\ \end{align}$$
Next, consider the linear fractional transformation given implicitly by the relation
$$\left(t-b\right)\left(1-u\right)=\left(a-b\right).$$
Then,
$$t=\frac{a-bu}{1-u}\implies dt=du\,\frac{\left(a-b\right)}{\left(1-u\right)^{2}},$$
and
$$\frac{t-a}{t-b}=u\land\frac{t-c}{t-b}=\frac{a-c}{a-b}-\frac{b-c}{a-b}u\land\frac{t-d}{t-b}=\frac{a-d}{a-b}-\frac{b-d}{a-b}u.$$
Applying this LFT to $\mathcal{I}$, we obtain
$$\begin{align} \mathcal{I}{\left(p,q,y\right)} &=\int_{a}^{z}\mathrm{d}t\,\frac{t}{\sqrt{\left(t-a\right)\left(t-b\right)\left(t-c\right)\left(t-d\right)}}\\ &=\int_{a}^{z}\mathrm{d}t\,\frac{t}{\left(t-b\right)^{2}\sqrt{\left(\frac{t-a}{t-b}\right)\left(\frac{t-c}{t-b}\right)\left(\frac{t-d}{t-b}\right)}}\\ &=\int_{0}^{\frac{z-a}{z-b}}\mathrm{d}u\,\frac{\left(a-b\right)}{\left(1-u\right)^{2}}\cdot\frac{a-bu}{1-u}\cdot\frac{\left(1-u\right)^{2}}{\left(a-b\right)^{2}}\\ &~~~~~\times\frac{1}{\sqrt{u\left(\frac{a-c}{a-b}-\frac{b-c}{a-b}u\right)\left(\frac{a-d}{a-b}-\frac{b-d}{a-b}u\right)}};~~~\small{\left[t=\frac{a-bu}{1-u}\right]}\\ &=\int_{0}^{\frac{z-a}{z-b}}\mathrm{d}u\,\frac{a-bu}{1-u}\cdot\frac{1}{\sqrt{u\left[\left(a-d\right)-\left(b-d\right)u\right]\left[\left(a-c\right)-\left(b-c\right)u\right]}}\\ &=\int_{0}^{\frac{\left(b-d\right)\left(z-a\right)}{\left(a-d\right)\left(z-b\right)}}\mathrm{d}v\,\frac{a-b\left(\frac{a-d}{b-d}\right)v}{1-\left(\frac{a-d}{b-d}\right)v}\\ &~~~~~\times\frac{1}{\sqrt{v\left(1-v\right)\left[\left(a-c\right)\left(b-d\right)-\left(b-c\right)\left(a-d\right)v\right]}};~~~\small{\left[u=\left(\frac{a-d}{b-d}\right)v\right]}\\ &=\frac{1}{\sqrt{\left(a-c\right)\left(b-d\right)}}\int_{0}^{\frac{\left(b-d\right)\left(z-a\right)}{\left(a-d\right)\left(z-b\right)}}\mathrm{d}v\,\left[\frac{a-b}{1-\left(\frac{a-d}{b-d}\right)v}+b\right]\frac{1}{\sqrt{v\left(1-v\right)\left[1-\frac{\left(a-d\right)\left(b-c\right)}{\left(a-c\right)\left(b-d\right)}v\right]}}\\ &=\frac{1}{\sqrt{\left(a-c\right)\left(b-d\right)}}\int_{0}^{\sin^{2}{\left(\varphi\right)}}\mathrm{d}v\,\left[\frac{a-b}{1-\eta\,v}+b\right]\frac{1}{\sqrt{v\left(1-v\right)\left(1-\kappa^{2}v\right)}}\\ &=\frac{2}{\sqrt{\left(a-c\right)\left(b-d\right)}}\int_{0}^{\sin{\left(\varphi\right)}}\mathrm{d}x\,\left[\frac{a-b}{1-\eta\,x^{2}}+b\right]\frac{1}{\sqrt{\left(1-x^{2}\right)\left(1-\kappa^{2}x^{2}\right)}};~~~\small{\left[v=x^{2}\right]}\\ &=\frac{2}{\sqrt{\left(a-c\right)\left(b-d\right)}}\left[\left(a-b\right)\Pi{\left(\varphi,\eta,\kappa\right)}+bF{\left(\varphi,\kappa\right)}\right],\\ \end{align}$$
where we've introduced the auxiliary parameters
$$\kappa:=\sqrt{\frac{\left(a-d\right)\left(b-c\right)}{\left(a-c\right)\left(b-d\right)}}\in\left(0,1\right),$$
$$\varphi:=\arcsin{\left(\sqrt{\frac{\left(b-d\right)\left(z-a\right)}{\left(a-d\right)\left(z-b\right)}}\right)}\in\left(0,\frac{\pi}{2}\right),$$
$$\eta:=\frac{a-d}{b-d}\in\left(1,\infty\right).$$
Note: even though $\eta>1$, we can still define $\Pi$ in the usual way without worrying about Cauchy principal values since $0<1-\eta\sin^{2}{\left(\varphi\right)}$.