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.

  • $\begingroup$ Thanks for your prompt and helpful answer. I took your advice and got it to work. However, I think there may be 1-2 mistakes in your answer. When I made the change of variables $(r_1,r_0,t)\to(r_1,u,s)=(r_1, \frac{r_0}{r_1},\frac{t-r_1}{r_0})$ and took the Jacobian, I got $dt\,dr_0\,dr_1=r_0 r_1 ds\,du\,dr_1 = ur_1^2 ds\,du\,dr_1$, and the extra factors of $u$ and $r_1^2$ seems to be missing from your integral. I am not sure about the $\frac{3}{128}$ either. I did the integral that I got and I got answers that agreed nicely with numerical simulations. It is possible I made a mistake... $\endgroup$ – Stefan Smith Feb 25 '13 at 17:10
  • $\begingroup$ ... since I've never taken or taught a good course in advanced calculus. Please note that I accepted your answer. $\endgroup$ – Stefan Smith Feb 25 '13 at 17:12
  • $\begingroup$ @StefanSmith Thanks for pointing out the omissions. I actually had the Jacobian written in my post (in one-variable form, since I was really thinking in terms of making one substitution at a time), but forgot to include it in the integral (1). Fixed now. $\endgroup$ – user53153 Feb 25 '13 at 17:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.