Why can't Wolfram Alpha compute this integral?

Question

$$\int_0^\infty \int_0^\infty \lambda_1\lambda_2 \mid{x - y} \mid e^{-\lambda_1x - \lambda_2y} dy dx$$

$$= \int_0^{\infty} e^{-\lambda_1x}\bigg[ \int_0^x(x - y)e^{-\lambda_2y}dy + \int_x^\infty(y - x)e^{-\lambda_2y}dy\bigg] dx$$

It seems like a straight forward but tedious integral to compute. Is there a way I can input this so that Wolfram is less confused by $x$ and $y$ being treated as variables and constants in different situations? Also Is there a way to specify that $\lambda_1$ and $\lambda_2$ are positive constants?

I replaced $\lambda_1$ and $\lambda_2$ with $\pi$ and $e$ and it gave me an answer.

Answer 1

The problem is that you're not telling WolframAlpha that $\lambda_1$ and $\lambda_2$ are positive. If you don't tell it this, it assumes they could be any real number (or even complex), so it doesn't know if the integral is convergent. To fix this, use the Assumptions option in Integrate like this:

Integrate[a*b*Abs[x-y]Exp[-a*x-b*y], {x, 0, Infinity}, {y, 0, Infinity}, Assumptions -> a>0 && b>0]

And it will give the correct answer: $(a^2+b^2)/[ab(a+b)] = a^{-1} + b^{-1} - 2(a+b)^{-1}$.

Answer 2

Regardless of the capabilities of WA, $\mathbb{R}^+\times \mathbb{R}^+$ can be partitioned as $\{0\leq x\leq y\}\cup\{0\leq y\leq x\}$ and your integral (assuming $\lambda_j>0$) can be written as $\lambda_1 \lambda_2$ times $$ \int_{0}^{+\infty}\int_{x}^{+\infty}(y-x)e^{-\lambda_1 x-\lambda_2 y}\,dy\,dx+\int_{0}^{+\infty}\int_{y}^{+\infty}(x-y)e^{-\lambda_1 x-\lambda_2 y}\,dx\,dy $$ or $$ \int_{0}^{+\infty}\int_{0}^{+\infty} t e^{-\lambda_1 x-\lambda_2(x+t)}\,dt\,dx+\int_{0}^{+\infty}\int_{0}^{+\infty} t e^{-\lambda_1(x+t)-\lambda_2 x}\,dt\,dx $$ or $$ \int_{0}^{+\infty}\int_{0}^{+\infty}t \left(e^{-\lambda_1 t}+e^{-\lambda_2 t }\right)e^{-(\lambda_1+\lambda_2)x}\,dx\,dt $$ or $$ \frac{1}{(\lambda_1+\lambda_2)}\left(\frac{1}{\lambda_1^2}+\frac{1}{\lambda_2^2}\right).$$

Answer 3

I don't know what to do about Wolfram but by adapting the method in this answer you can find that the PDF of $Z=|X-Y|$ is given by $$\frac{\lambda_1 \lambda_2}{\lambda_1+\lambda_2} [e^{-\lambda_1 z}+e^{-\lambda_2 z}]$$ which is a weighted sum of PDFs of two exponential RVs with rates $\lambda_1$ and $\lambda_2$ respectively. Thus, the expectation can be quickly computed as, $$\mathbb{E}(Z)=\frac{1}{\lambda_1+\lambda_2}[\lambda_2 \mathbb{E}(X)+\lambda_1 \mathbb{E}(Y)],$$ where $X\sim Exp(\lambda_1)$ and $Y\sim Exp(\lambda_2)$, yielding $$\mathbb{E}(|X-Y|)=\frac{\lambda_1^2+\lambda_2^2}{\lambda_1 \lambda_2 (\lambda_1+\lambda_2)}.$$

Running this in R with,

n <- 50000
lambda <- 1/3
mu <- 1/5
x <- rexp(n, rate = lambda)
y <- rexp(n, rate = mu)
z <- abs(x-y)
exact  <- (mu^2+lambda^2)/((mu*lambda)*(mu+lambda))
sample <- mean(z)
print(data.frame(sample, exact))

gives $4.250225$ and $4.25$ for example.