Math people:
I am trying to use Matlab's "triplequad" to numerically evaluate the triple integral
$$\int_0^1 \int_0^{r_1} \int_{r_1-r_0}^{r_1+r_0} \left[\frac{Ov(r_0,r_1,t)^2}{(\frac{4}{3}\pi r_0^3)(\frac{4}{3}\pi r_1^3)}\right] *\frac{t^2}{2\sqrt{\pi}}\exp(-t^2/4)\,dt\,dr_0\,dr_1,$$
where $$Ov(r_0,r_1,t)=\left[(1/12)\pi(r_0+r_1-t)^2(t^2+2t(r_0+r_1)-3(r_1-r_0)^2)\right]/t.$$ I am handling the variable limits of integration by means of heaviside functions. It is not obvious, but the expression in square brackets in the triple integral is between $0$ and $1$. Before I integrate, I combine $Ov(r_0,r_1,t)^2$ and $t^2$ so Matlab doesn't have to deal with small values of $t$ in a denominator.
Letting $r_0$ and $r_1$ start at $0$ doesn't work because Matlab thinks you are trying to divide by zero. I tried letting $r_0$ and $r_1$ start at $0.001$ and I obtained results that I am pretty sure were poor. Based on where this problem comes from (I can add to the question if anyone is curious), I ran a simulation that involved generating 1000000 pairs of points in $\mathbf{R}^3$ from a certain distribution and averaging a quantity computed from these points, and the simulation gave me $0.0036776...$. I rescaled the triple integral above so $r_0$ and $r_1$ went up to $100$ instead of $1$, and Matlab gave me $0.00343179...$, which is reasonably close to the simulation but not spectacular.
My question is, how do you deal with integrals like this in Matlab, where the integrand involves quotients of small numbers? I don't need anyone to solve the whole problem for me, but I'd like to know if there are any general principles that people use or options that you can use in Matlab's triplequad that will yield better results. I am a novice at numerical methods and Matlab. All the examples of Matlab's multiple integral commands I have seen online work out nice and don't help me.
STEFAN (STack Excange FAN)
Try to integrate over rectangular boxes whenever possible. More vaguely stated, look for substitutions that cancel "large integrand" with "thin region of integration".
In this case, this means making the substitutions $r_0=u r_1$ and $t=r_1+sr_0=r_1+su\,r_1$. Then $dr_0=u\,dr_1$ and $dt=ur_1\,ds$, and the integral becomes $$\int_0^1 \int_0^{1} \int_{-1}^{1} \left[\frac{F(u,r_1,s)^2}{(\frac{4}{3}\pi u^3 r_1^3)(\frac{4}{3}\pi r_1^3)}\right] \frac{1}{2\sqrt{\pi}}\exp(-r_1^2(1+su)^2/4)\,u\,r_1^2\,ds\,du\,dr_1 \tag1$$ where $$F(u,r_1,s)=\frac{\pi}{12}u^3 r_1^4(1-s)^2 (4s+s^2u+2su-3u+8) \tag2$$ The expression (1) simplifies to $$\frac{3}{128\sqrt{\pi}}\int_0^1 \int_0^{1} \int_{-1}^{1} G(u,r_1,s) \exp(-r_1^2(1+su)^2/4)\,ds\,du\,dr_1 \tag3$$ where $G(u,r_1,s)= u^4 r_1^4 (1-s)^4 (4s+s^2u+2su-3u+8)^2 $.
That's right: you get a polynomial times a Gaussian. No quotients of small numbers.