15
$\begingroup$

Computational Science people:

I originally posted this question at Math Stack Exchange and someone commented that I might get "much better" answers here:

I am a novice at numerical methods and Matlab. I am attempting to evaluate the following sum of two triple integrals (it can obviously be written more simply, but you still cannot evaluate it symbolically (?)). I am having trouble getting the $\LaTeX$ to work here, so I reluctantly broke it up into pieces here: I want to find the sum of

$$\frac{2}{((1/0.3) - 1)^2}\left(\int_1^{1/0.3}\int_1^{r_1}\int_0^{r_1-r_0}F_1(r_0,r_1,t)\exp(-\frac{(0.3)^2 t^2}{4})\,dt\,dr_0\,dr_1 \right),$$

and

$$\frac{2}{((1/0.3) - 1)^2}\left(\int_1^{1/0.3}\int_1^{r_1}\int_{r_1-r_0}^{r_1+r_0} F_2(r_0,r_1,t)\exp(-\frac{(0.3)^2 t^2}{4})\,dt\,dr_0\,dr_1 \right),$$

where

$$F_1(r_0,r_1,t)=\frac{t^2 r_0^3*(0.3)^3}{2r_1^3\sqrt{\pi}}$$

and

$$F_2(r_0,r_1,t)=\frac{(0.3)^3\pi^{3/2}(r_0+r_1-t)^4 (t^2+2t(r_0+r_1)-3(r_1-r_0)^2)^2}{288(\frac{4}{3}\pi r_0^3)(\frac{4}{3}\pi r_1^3)}.$$

EDIT (March 2 2013): Someone responded that they got Mathematica to do the integrals symbolically. I just attempted to do this (with simplified versions of the integrals) and Mathematica could only do the outer two of the first one, and stalled on the second one. I would appreciate some help. Here is what I did.:

I attempted to evaluate

$$\int_1^2 \int_1^{r_2} \int_0^{r_2-r_1} \frac{r_1^3 t^2 \exp(-t^2)}{r_2^3}\,dt\,dr_1\,dr_2$$ via

Integrate[r1^3/r2^3*t^2*Exp (-t^2), {t, 0, r2 - r1}, {r1, 1, r2}, {r2, 1, 2}]

and Mathematica returns (I had trouble with the $\LaTeX$ here because the result is long. I broke it into two equations. if anyone knows a good way to display this please tell me):

$$\int_1^2 \frac{1}{64r2^2} e^{-1-r2^2}(2e^{2r2}(25+r2(19+2r2(1+r2)))-$$

$$e^{1+r2^2}(32r2(2+r2^2)) +\sqrt{\pi}(11+4r2^2(9+r2^2))\operatorname{Erf}[1-r2])\,dr2.$$

Then I tried to evaluate

$$\int_1^2\int_1^{r_2}\int_{r_2-r_1}^{r_2+r_1} \ldots \qquad \qquad \qquad $$

$$\ldots\frac{\exp(-t^2)(r_1+r_2-t)^4(t^2+2t(r_1+r_2)-3(r_2-r_1)^2)^2}{r_1^3 r_2^3}\,dt\,dr_\,dr_2$$

using

Integrate[(r1 + r2 - t)^4*(t^2 + 2*t*(r1 + r2) - 3*(r2 - r1)^2)^2* Exp[-t^2]/r1^3/r2^3, {r2, 1, 2}, {r1, 1, r2}, {t, r2-r1, r2 + r1}]

just now, and Mathematica has not returned an answer after about half an hour (but I am having computer network problems right now, and they may be to blame).

[END OF MARCH 2 EDIT]

I used Matlab's "triplequad" command, with no extra options. I handled the variable limits of integration by means of heaviside functions, because I didn't know any other way to do it. Matlab gave me $0.007164820144202$.

I know Matlab is good software, but I have heard that numerical triple integrals are hard to do accurately, and mathematicians are supposed to be skeptical, so I want some way to verify the accuracy of this answer. The integrals give the expected value of a certain experiment (if anyone wants, I can edit this question to describe the experiment): I implemented the experiment in Matlab using appropriately randomly generated numbers, a million times, and averaged the results. I repeated this process four times. Here are the results (I apologize if I have used the word "trial" improperly):

Trial 1: $0.007133292603256$

Trial 2: $0.007120455071989$

Trial 3: $0.007062595022049$

Trial 4: $0.007154940168452$

Trial 5: $0.007215000289130$

Although each trial used a million samples, the simulation values only agree in the first significant digit. They are not close enough to each to each other for me to determine whether the numerical triple integral is accurate.

So can anyone tell me whether I can trust the result of "triplequad" here, and under what circumstances one can trust it in general?

One suggestion I got at Math Stack Exchange was to try other software like Mathematica, Octave, Maple, and SciPy. Is this good advice? Do people actually do numerical work in Mathematica and Maple? Octave is kind of a Matlab clone, so can I assume it uses the same integration algorithms? I haven't even heard of SciPy before and would appreciate any opinions about it.


UPDATE: Someone from Math Stack Exchange did it in Maple and got $0.007163085468$. That is agreement to three significant figures. That is a good sign.

Also, I would appreciate suggestions on how to enter long, multi-line expression in $\LaTeX$ in Stack Exchange. Can you use the "aligned" environment here? I tried, and I couldn't get it to work.

$\endgroup$
  • 2
    $\begingroup$ Your simulation results are perfectly consistent with the numerical value returned by Matlab: their mean of $0.00713726$ is just $-1.11$ standard errors less than what Matlab returned. FWIW, Mathematica returns $0.00716308537\ldots$. It can also evaluate these integrals symbolically in terms of polynomials and error functions. $\endgroup$ – whuber Mar 1 '13 at 17:29
  • $\begingroup$ @whuber Thanks. I could swear I tried it symbolically in Maple and Maple couldn't do it. I'll try again in Maple, and if it doesn't work, I'll try it in Mathematica. BTW, I did a similar integral in Maple, and I got a huge symbolic answer. It appeared to be a sum and difference of very large numbers whose grand total was quite small. I suspect roundoff error was likely in the final answer. In a problem like this, should you use the symbolic answer, or just do the integral numerically? $\endgroup$ – Stefan Smith Mar 2 '13 at 0:46
  • $\begingroup$ Symbolic answers have the advantage of being combinations of functions that (often) can efficiently be computed to arbitrary precision. Usually, too, the symbolic solution also lends itself to rapid recomputation when parameters are varied. For these reasons it's often worthwhile to seek a symbolic solution. $\endgroup$ – whuber Mar 2 '13 at 4:34
  • $\begingroup$ @whuber :I tried doing some essentially equivalent integrals (changing some of the constants and removing some multiplicative constants) in Mathematica, and Mathematica could only do the outer two integrations of the first integral, and seems to have stalled on the second one. I posted my code and results above. $\endgroup$ – Stefan Smith Mar 2 '13 at 15:24
  • 1
    $\begingroup$ Re The March 2 edit: By reducing the triple integral symbolically to a single integral (in the first half of your integrals) you have accomplished a lot. The integrand is very nicely behaved and can be numerically integrated to extremely high precision within a fraction of a second. $\endgroup$ – whuber Mar 2 '13 at 20:31
9
$\begingroup$

First of all, it is not the software (or at least it shouldn't be) that determines the quality of the solution to a problem, it's the quality and appropriateness of the algorithm that is applied. You should check what algorithm is being used by triplequad in Matlab (I would guess it uses a nested adaptive Gaussian quadrature). And you should check what the requested tolerances are (required absolute and relative tolerance). Chances are that, by default, it only asks for $10^{-8}$ relative precision.

The answer coming from Maple is probably done by Computer Algebra and maybe it could find a closed solution which was then evaluated using double precision floating point. This has the advantage that you are not approximating the integral by a finite summation (and hence are introducing approximation errors) but the Computer Algebra System will find an expression for the integral which can then be evaluated. Of course, care must be taken when evaluating this expression (for round-off).

If you would want to do this with SciPy, you would also need to resort to nested adaptive Gaussian quadrature using the underlying Quadpack (Piessens et al.) routines. In Octave, you'll have the same approach. And I wouldn't be too surprised if Matlab also uses Quadpack as quadrature engine (since it is the reference).

$\endgroup$
  • $\begingroup$ @GretVdE : Thanks for the info. I tried evaluating the integral symbolically first, and Maple couldn't do it (so it was probably impossible, using standard functions), so I asked Maple to do it numerically. I don't know what algorithm it used. $\endgroup$ – Stefan Smith Mar 1 '13 at 12:21
  • $\begingroup$ @StefanSmith : you can find out by setting the infolevel in Maple: infolevel[`evalf/int`] := 4. Are you sure Mape can't find a closed solution? The integral doesn't seem to be too complicated. Could you make your Maple sheet public somewhere? $\endgroup$ – GertVdE Mar 1 '13 at 13:04
  • $\begingroup$ @StefanSmith : I would post the Maple code in the question above. $\endgroup$ – GertVdE Mar 1 '13 at 13:36
  • $\begingroup$ I can't get Maple to work on my system right now, but I attempted equivalent integrals in Mathematica, and Mathematica did only the inner two of the first triple integral, and stalled on the second triple integral. Please see the edited question. $\endgroup$ – Stefan Smith Mar 2 '13 at 15:14

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.