Bounty ending tomorrow:
Consider the following function defined by two infinite integrals: $$ F(\epsilon) = \int_0^\infty \frac{q^3 (1+q)e^{-2q}}{q^4+\epsilon^4} \, \mathrm{d} q + \int_0^\infty \frac{\, G_q \, e^{-q}}{q^2(q^4+\epsilon^4)} \, \mathrm{d} q $$ where $$ G_q := G \left( \left[ \left[\frac{1}{2} \right] , [\,]\right] , \left[ \left[ \frac{9}{2},\frac{5}{2} \right], \left[\frac{3}{2}\right] \right] ,\frac{q^2}{4} \right) \, , $$ with $G$ being the Meijer G-function. The goal is to study the behavior of $F$ around the singularity point, i.e. as $\epsilon\to 0$. Numerically, it can clearly be observed that $F$ scales logarithmically with $\epsilon$. I am wondering whether this behavior can be shown analytically via a rigorous analysis, e.g. using perturbation techniques. Your help or hints are highly appreciated and are most welcome.
Thanks H
Your Integral may be expressed as a Meijer G-Function of two Variables, which is discussed by Agarwal in 1964 for the first time. A more general form is the H-Fox-Function of two Variables which can be found in the standard publication of Mathai: https://www.researchgate.net/publication/266566090_The_H-function_Theory_and_Applications.
$$F_{1}\left( \epsilon \right) =\int_{0}^{\infty }\frac{q^{3}\left( 1+q\right) \exp \left( -2q\right) }{q^{4}+\epsilon ^{4}}dq$$
can be solved analytically in terms of two Meijer G-functions:
$$F_{1}\left( \sigma \right) \frac{1}{16\sqrt{2}\pi ^{\frac{3}{2}}\sigma ^{3}}% \left[ G_{1,5}^{5,1}\left( 16\,\sigma ^{4}\left\vert \begin{array}{c} -,-,\frac{3}{4},-,- \\ \frac{3}{4},1,\frac{5}{4},\frac{3}{2},\frac{7}{4}% \end{array}% \right. \right) +\frac{1}{4\,\sigma }G_{1,5}^{5,1}\left( 16\,\sigma ^{4}\left\vert \begin{array}{c} -,-,1,-,- \\ 1,1,\frac{5}{4},\frac{3}{2},\frac{7}{4}% \end{array}% \right. \right) \right] $$
where we use $\sigma =\frac{\epsilon }{4}$ to have the same notation later on for the second indefinite integral:
$$F_{2}\left( \epsilon \right) =\int_{0}^{\infty }\frac{1}{q^{2}\left( q^{4}+\epsilon ^{4}\right) }G_{1,3}^{2,1}\left( \frac{q^{2}}{4}\left\vert \begin{array}{c} -,\frac{1}{2},- \\ \frac{9}{2},\frac{5}{2},\frac{3}{2}% \end{array}% \right. \right) \exp \left( -q\right) dq$$
The first and the last function can be found in The wolfram functions side. They can be expressed in terms of their Meijer G-identities and then $F_{2}$ can be evaluated using the integration theorem for Meijer G-functions (http://functions.wolfram.com/HypergeometricFunctions/MeijerG/21/02/04/)
$$\int_{0}^{\infty }t^{\alpha -1}G_{p,q}^{m,n}\left( \left\vert \begin{array}{c} a_{1},...,a_{p} \\ b_{1},...,b_{q}% \end{array}% \right. \right) \times $$
$$G_{p_{1},q_{1}}^{m_{1},n_{1}}\left( x\,t\left\vert \begin{array}{c} a_{11},...,a_{1p_{1}} \\ b_{11},...,b_{1q_{1}}% \end{array}% \right. \right) \times G_{p_{2},q_{2}}^{m_{2},n_{2}}\left( y\,t\left\vert \begin{array}{c} a_{21},...,a_{2p_{2}} \\ b_{21},...,b_{2q_{2}}% \end{array}% \right. \right) dt$$ $$z^{-\alpha }G_{p,q,p_{1},q_{1},p_{2},q_{2},}^{m,n,m_{1},n_{1},m_{2},n_{2}}\left( \begin{array}{c} 1-\alpha -b_{1},...,1-\alpha -b_{q} \\ 1-\alpha -a_{1},...,1-\alpha -a_{p}% \end{array}% \left\vert \begin{array}{c} a_{11},...,a_{1p_{1}} \\ b_{11},...,b_{1q_{1}}% \end{array}% \right. \left\vert \begin{array}{c} a_{21},...,a_{2p_{2}} \\ b_{21},...,b_{2q_{2}}% \end{array}% \right. \left\vert \frac{x}{z},\frac{y}{z}% \begin{array}{c} \\ \end{array}% \right. \,\right) $$
The Meijer G-identities corresponding to the functions are: $$\exp \left( -q\right) =G_{0,1}^{1,0}\left( q\left\vert \begin{array}{c} -,- \\ 0,-% \end{array}% \right. \right) $$ $$\frac{1}{q^{2}\left( q^{4}+\epsilon ^{4}\right) }=\frac{1}{\epsilon ^{6}}% G_{1,1}^{1,1}\left( \left( \frac{q}{\epsilon }\right) ^{4}\left\vert \begin{array}{c} -\frac{1}{2},- \\ -\frac{1}{2},-% \end{array}% \right. \right) $$
After some algebraic manipulations and variable transformations, using (http://functions.wolfram.com/HypergeometricFunctions/MeijerG/17/02/03/) we conclude:
$$F_{2}\left( \sigma \right) =\frac{1}{32\sqrt{2}\pi ^{\frac{5}{2}}\sigma ^{6}}\int_{0}^{\infty }t^{\frac{1% }{4}}G_{1,1}^{1,1}\left( \frac{t}{\sigma ^{4}}\left\vert \begin{array}{c} -\frac{1}{2},- \\ -\frac{1}{2},-% \end{array}% \right. \right) \times G_{2,6}^{4,2}\left( t\left\vert \begin{array}{c} -,-,\frac{1}{4},\frac{3}{4},-,- \\ \frac{9}{4},\frac{11}{4},\frac{5}{4},\frac{7}{4},\frac{3}{4},\frac{5}{4}% \end{array}% \right. \right) $$
$$\times G_{0,4}^{4,0}\left( t\left\vert \begin{array}{c} -,-,-,- \\ 0,\frac{1}{2},\frac{1}{4},\frac{3}{4}% \end{array}% \right. \right) dt =$$
$$= \frac{1}{32\sqrt{2}\pi ^{\frac{5}{2}}\sigma ^{5}}% ~G_{1,1,2,6,0,4}^{1,1,4,2,4,0}\left( \begin{array}{c} \frac{5}{4},- \\ \frac{5}{4},-% \end{array}% \left\vert \begin{array}{c} -,-,\frac{1}{4},\frac{3}{4},-,- \\ \frac{9}{4},\frac{11}{4},\frac{5}{4},\frac{7}{4},\frac{3}{4},\frac{5}{4}% \end{array}% \right. \left\vert \begin{array}{c} -,-,-,- \\ 0,\frac{1}{2},\frac{1}{4},\frac{3}{4}% \end{array}% \right. \left\vert \sigma ^{4},\sigma ^{4}% \begin{array}{c} \\ \end{array}% \right. \,\right) $$
A similar example can be found here: https://www.researchgate.net/publication/269504875_Capacity_of_k_-_m_Shadowed_Fading_Channels. You may find also the source code for implementing the function in Mathematica in the publication above. Due to the lack of time, I can not check, if the solution of your Integral reduces to functions of only one Variable, since $\frac{x}{z}=\frac{y}{z}=\sigma ^{4}$. You may do this numerical, also I can provide the H-Fox-Function of one variable.