Double integral $\int_0^\infty \int_0^\infty \frac{e^{-(x+y)/2}x^j y^k}{ax+by+c} dx dy$

Question

In an ideal scenario I would like to have an analytic expression for the integral $$ \operatorname{I}\left(\,{a,b,c,j,k}\,\right) = \int_{0}^{\infty} \int_{0}^{\infty} \dfrac{\operatorname{e}^{-\left(\,{x + y}\,\right)/2}\,\,x^{j}\, y^{k}}{ax + by + c}\,{\rm d}x\,{\rm d}y $$ for positive real $a,b,c$ and non-negative integers $j,k$, which can then be evaluated for very large values of $c$. The reason I mention the latter circumstance is because while Mathematica ($13.1$ for sure) is able to take this integral for small values of $j$ and $k$, the resulting expression is numerically unstable for large values of $c$ (meaning that it's a product of terms exponentially large and small in $c$). For example, in the case $j = k = 0$ one gets: $$ \frac{2 e^{\frac{c}{2 b}} \operatorname{Ei}\left(-\frac{c}{2 b}\right)-2 e^{\frac{c}{2 a}} \operatorname{Ei}\left(-\frac{c}{2 a}\right)}{a-b} $$

Would appreciate any advice on how to:

1. Take the original integral analytically,
2. Efficiently evaluate $\operatorname{e}^\alpha \operatorname{Ei}(-\alpha)$ for large values of $\alpha$.

Given that Mathematica can take this integral for small $j$ and $k$, I would guess that there exists a recursive relation.

Answer 1

$\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}[1]{{\displaystyle #1}} \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\on}[1]{\operatorname{#1}} \newcommand{\pars}[1]{\left(\,{#1}\,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\sr}[2]{\,\,\,\stackrel{{#1}}{{#2}}\,\,\,} \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$ \begin{align} \on{I}\pars{a,b,c,j,k} & \equiv\color{#44f} {\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c} \,\dd x\,\dd y} \\[5mm] & \sr{\substack{\,\,\,\,x/c\ \mapsto\ x\\[0.5mm] \atop y/c\ \mapsto\ y\\[0.5mm] \mbox{}}}{=} c^{j + k + 1}\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-c\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + 1} \,\dd x\,\dd y \\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim} c^{j + k + 1}\pars{\int_{0}^{\infty}\expo{-cx/2}\,x^{j}\,\dd x} \pars{\int_{0}^{\infty}\expo{-cy/2}\,y^{k}\,\dd y}\tag{1}\label{1} \\[5mm] & = c^{j + k + 1}\ \bracks{\pars{2 \over c}^{j + 1}\Gamma\pars{j + 1}} \bracks{\pars{2 \over c}^{k + 1}\Gamma\pars{k + 1}} \end{align} (\ref{1}) involves the $\underline{Laplace's\ Method}$ because the main contribution to the integral, as $\ds{c \to \infty}$, occurs in a neighborhood of "$\ds{x = 0\ \mbox{and}\ y = 0}$". Therefore, \begin{align} \on{I}\pars{a,b,c,j,k} & \equiv\color{#44f} {\int_{0}^{\infty}\int_{0}^{\infty} {\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c} \,\dd x\,\dd y} \\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim} \bbx{\color{#44f}{2^{j + k + 2}\,\,\,\Gamma\pars{j + 1}\Gamma\pars{k + 1}} \over c} \end{align}

Answer 2

For the second question $$e^{\alpha } \,\text{Ei}(-\alpha )=\sum _{n=1}^{\infty } (-1)^n \,\frac{(n-1)!}{\alpha^n}$$ You will not need many terms for a very good approximation.

You could also use the $[k,k+1]$ Padé approximant $P_k$. For example $$P_2=-\frac{\alpha (3 \alpha +13)}{3 \alpha ^3+16 \alpha ^2+10 \alpha -4}$$ whose error is $\frac {44}{3 \alpha^6 }$.

For $\alpha=100$, this would give $$e^{100 } \,\text{Ei}(-100)=-\frac{7825}{790249}=-0.009901942299$$ while the exact value is $-0.009901942287$.

Using $P_2$ the relative error is smaller than $0.01$% as soon as $\alpha \geq 8$ and smaller than $0.001$% as soon as $\alpha \geq 14$.