# Multiple integral over a disc

I would need some help on this integration problem:

$$I=\int_0^{2\pi}\int_0^{R}\int_0^{2\pi}\int_0^{R}\exp(-a\ r_{12}) \ r_1 \ r_2 \,\mathrm{d}r_1\,\mathrm{d}\phi_1\,\mathrm{d}r_2\,\mathrm{d}\phi_2$$

Here, $r_1, r_2, \phi_1$ and $\phi_2$ are the polar coordinates for the integral on a disc with radius R. And $r_{12}$ is the distance between the coordinates, which can be calculated by

$$r_{12}=\sqrt{r_1^2+r_2^2-2r_1r_2\cos(\phi_1-\phi_2)}.$$

Below a plot of the problem for $a=1$, $R=1$ and fixed coordinates of the first point at $(0,0)$:

• I believe the expression for $r_{12}$ is missing a factor of 2. shouldn't it be: $r_{12}=\sqrt{r_1^2+r_2^2-\color{red}2\,r_1r_2\cos(\phi_1-\phi_2)}$? Jun 11, 2014 at 22:04
• Due to the radial symmetry, you can perform the outermost angular integral to get a factor $2\pi$. I wouldn't be surprised if it turns out to be impossible to make any progress beyond that. Jun 12, 2014 at 7:14
• Thanks David, I corrected this error! Jun 12, 2014 at 7:48
• Thanks Joriki for the hint. I know, the integral of the exponential function looks really simple, but I already spend some time with it. Perhaps I am wrong, but can I really solve the outermost angular integral that easy? The angular cooridinate is still included in $\cos(\phi_1-\phi_2)$. Jun 12, 2014 at 7:53
• @AnkaiosArgo First change variables in the innermost angular integral by $\phi_1^\prime=\phi_1-\phi_2$. Then the integral $\int_{0}^{2\pi}d\phi_1\,f(\phi_1-\phi_2)$ becomes $\int_{-\phi_2}^{2\pi-\phi_2}d\phi_1^\prime\,f(\phi_1^\prime)=\int_{0}^{2\pi}d \phi _1^\prime\,f(\phi_1^\prime)$. Jun 12, 2014 at 8:13

This integral has a closed form in terms of special functions. First, start with the integral in the form $$I(\alpha, R) = \int_{D(0,R)}da\int_{D(0,R)}db e^{-\alpha|a-b|}$$ with parameters $\alpha$, $R$, and differentiate with respect to $R$: \begin{aligned} \frac{dI}{dR}(\alpha, R) &= \lim_{\epsilon\to0} \frac{2}{\epsilon}\int_{D(0,R)}da\int_{D(0,R+\epsilon)\setminus D(0,R)}db\, e^{-\alpha|a-b|} \\&= 4\pi R\int_{D(0,R)}da\, e^{-\alpha|a - R e_1|} \\&= 4\pi R^3\int_{D(0,1)}da\, e^{-\alpha R|a-e_1|}. \end{aligned} Here $e_1=(1,0)$ is a unit vector on the unit circle, and the integral simplified because $\int da\,e^{-\alpha|a-b|}$ is independent of where $b$ is on the circle $\partial D(0,R)$ of radius $R$. The factor of $2$ came from the two occurrences of $D(0,R)$ that are symmetric, and $2\pi R$ is the area of the strip $D(0,R+\epsilon)\setminus D(0,R)$ divided by $\epsilon$.

To evaluate this integral, write $a$ in polar coordinates, but with the origin at the point $(1,0)$, so that the unit disk has the form $$\pi/2 \leq \phi \leq 3\pi/2, \qquad 0 \leq r \leq -2\cos\phi = \rho(\phi).$$ Then the integral for $\frac{dI}{dR}$ becomes \begin{aligned} \frac{dI}{dR} &= 4\pi R^3\int_{\pi/2}^{3\pi/2}d\phi\int_0^{-2\cos\phi}r\,dr\, e^{-\alpha R r} \\&= 4\pi R^3\int_{\pi/2}^{3\pi/2}d\phi\frac{1-(1+\alpha R\rho)e^{-\alpha R\rho}}{\alpha^2 R^2}. \end{aligned} The integral over $\phi$ can be done with the help of computer algebra: $$\frac{dI}{dR} = \frac{4\pi^2 R}{\alpha^2}\big(1-I_0(2\alpha R)+2\alpha R I_1(2\alpha R)-2\alpha R \mathbf{L}_{-1}(2\alpha R)+\mathbf{L}_0(2\alpha R)\big),$$ where $I$ are the modified Bessel functions, and $\mathbf{L}$ are the modified Struve functions.

Then the integral $I(\alpha,R) = \int \frac{dI}{dR}$ can also be done by a computer in terms of Struve and hypergeometric functions, which yields $$I(\alpha,R) = \frac{2\pi R}{\alpha^3}\big( -2+\alpha \pi R-\alpha \pi R \,{}_0\mathbf{F}_1(2,\alpha^2R^2)+2 \alpha^3 \pi R^3 \,{}_0\mathbf{F}_1(3,\alpha^2 R^2)-2 \alpha \pi R \mathbf{L}_{-2}(2 a R)+\pi \mathbf{L}_1(2\alpha R)\big)$$ where ${}_0\mathbf{F}_1(a;z) = {}_0F_1(a;z)/\Gamma(a)$ is the regularized hypergeometric function.

This is the Mathematica expression for the above formula:

1/a^3 2 \[Pi] R (-2 + a \[Pi] R -
a \[Pi] R Hypergeometric0F1Regularized[2, a^2 R^2] +
2 a^3 \[Pi] R^3 Hypergeometric0F1Regularized[3, a^2 R^2] -
2 a \[Pi] R StruveL[-2, 2 a R] + \[Pi] StruveL[1, 2 a R])

EDIT Differentiating the integral. Let $D=D(0,R)$, $C=D(0,R+\epsilon)\setminus D(0,R)$. Then $$\left(\int_{D}+\int_C\,da\right)\left(\int_D+\int_C\,db\right) - \int_Dda\int_C db = \int_Cda\int_D db + \int_Dda\int_Cdb + \int_C\int_C\,da\,db.$$ The term $\int_C\int_C$ is on the order of $O(\epsilon^2)$ because the area of $C$ is $2\pi R\epsilon$, so it may be ignored.

Because the function $e^{-\alpha|a-b|}$ is symmetric in $a,b$, $$\int_Cda\int_D db + \int_Dda\int_Cdb = 2\int_Dda\int_Cdb.$$ Also, the integral $\int_D e^{-\alpha|a-b|}\,da$ depends only on the distance between $b$ and the origin, so picking a special value of $b=Re_1$ gives $$2\int_D da\int_C db\,e^{-\alpha|a-b|} = 2\mathop{\mathrm{area}}(C)\int_Dda\,e^{-\alpha|a-Re_1|}.$$

• Thanks Kirill, this looks very smart. However, I tried to verify the solution and the value I found was $31.55$ for $a=1$ and $R=1$. But perhaps I did a mistake, because I had to do it with matlab, because I didn't had mathematica at hand. Please find the code together with the other code above. Could you please verify it? The solution should be close to $17.43$. Jun 12, 2014 at 12:41
• @AnkaiosArgo In your code why do you multiply the sum by $(2\pi R^2)^2$? I think this should be $(\pi R^2)^2$ as in the Riemann sum you gave. My formula, and the original integral (with Mathematica's NIntegrate) at $a=1,R=1$ are both equal to $4.35588845190304646851\ldots$ ($\frac14$ times what you said). I have no access to Matlab right now, but when you evaluate the formula in Matlab, are you sure "struvem" is the modified Struve function $\mathbf{L}$? Jun 12, 2014 at 13:28
• @AnkaiosArgo By the way, Matlab has good support for numerical integration, you don't have to implement it yourself. Have a look at this page for example: matlab's default method will be usually faster and more accurate. Jun 12, 2014 at 13:32
• @AnkaiosArgo Also, it's the regularized hypergeometric function, so it should be hypergeom([],a,z)/gamma(a). Jun 12, 2014 at 13:35
• Thanks a lot. I corrected the factor 4 above and in the plot! Thus, your solution should be correct! Might be a problem of the matlab functions. I am not really sure, if struvem is what we need... and also thanks for the hint with hypergeom... but in this special case ($a=1$) it does not change the value for the code. I will try it in Mathematica. Would be great, if you code post your code for your formula. Thanks again for your help! Jun 12, 2014 at 13:47