# how to reduce the time to calculate triple integral

i want to evaluate the following integral, but it takes 15 min to do it. i must reduce the time to calculate it, i have to do a table and create a graph, but takes too long. anyone have an idea ?

can i use mathlink to import the FORTRAN to do it ?

minroot[g_?NumericQ, b_?NumericQ] :=
Module[{rts, y},
rts = y /. Solve[g^2 - b^2*g^2  y - 4 y^6 + 4 y^3 == 0, y];
rts = Select[rts, With[{nval = N[#, 100]}, Im[nval] == 0 && nval > 0] &];
Min[rts]^(1/2)];

aA[g_?NumberQ, b_?NumberQ] :=
Pi - 2 b g NIntegrate[ 1/Sqrt[g^2 - b^2*g^2  y^2 - 4*y^12 + 4  y^6],
{y, 0, minroot[g, b] - 0.00001}];

qQ[g_?NumberQ] := NIntegrate[2*(1 - Cos[aA[g, b]]) b, {b, 0, 10}];

o[T_] = (1/T^3) NIntegrate[(g^5*qQ[g]) / E^(g^2/T), {g, 0, 50}];

o[4]


thanks for any help

• Try taking the Solve outside of minroot! May 3, 2013 at 15:00
• @PlatoManiac That might help slightly but it seems not to be a major bottleneck. Likewise replacing it with NSolve (keeping it in the loop) might help slightly. May 3, 2013 at 15:14
• I had been playing with various option settings and got reasonable speed that way. But then I sawthe response from @Philipp. Since my variations were no better I simply gave that one an upvote and left it at that. As for other ideas, I'm largely at a loss. About the only thing that comes to mind is to try to recast the underlying problem you wish to solve in a way that is less computation-hungry. barring that, maybe try to farm out the computations to several processors. May 6, 2013 at 14:54
• @DanielLichtblau my computer is octacore, an i7 processor ... how can i use the 8 nucleos ?? May 6, 2013 at 15:40
• I'd go with ParallelTable maybe. Check documntation on parallel operations in Mathematica. May 6, 2013 at 16:25

As always there are several ways to improve the speed

• Options Use MaxRecursions and MaxPoints
• Method Try using different methods to obtain quickest solution.
• Precalculate Use Block or Module to have some intermediate results only calculated once when required.
• Analyze Very general advice: use debug features as AbsoluteTiming at several places to see, where time is lost.

In your case I could improve speed by using

Method -> {Automatic, "SymbolicProcessing" -> 0}


and by controlling MaxRecursions. How few recursions you use depends on how exact you need to have the results:

aA[g_?NumberQ, b_?NumberQ, i_] :=   Pi - 2 b g NIntegrate[
1/Sqrt[g^2 - b^2*g^2 y^2 - 4*y^12 + 4 y^6], {y, 0,
minroot[g, b] - 0.00001}, MaxRecursion -> i, Method -> {Automatic, "SymbolicProcessing" -> 0}];
qQ[g_?NumberQ, i_] :=   NIntegrate[2*(1 - Cos[aA[g, b, i]]) b, {b, 0, 10},   MaxRecursion -> i, Method -> {Automatic, "SymbolicProcessing" -> 0}];


Lets see, what i we need:

Table[qQ[5, i] // Timing, {i, 1, 6}]
{{0.076005, 0.983541}, {0.116007, 0.854045}, {0.176011,   0.807208}, {0.176011, 0.782326}, {0.552034, 0.762594}, {0.632040,   0.762594}}


Looks like 6 is enough. Then:

o[T_, i_, j_] := (1/T^3) NIntegrate[(g^5*qQ[g, i])/E^(g^2/T), {g, 0, 50},
MaxRecursion -> j, Method -> {Automatic, "SymbolicProcessing" -> 0}];


And we see that

Table[o[5, 6, j] // Timing, {j, 1, 3}]
{{25.529596, 0.871541}, {40.394524, 0.868401}, {53.767361, 0.868459}}


j=1 is sufficient and calculation takes 26 seconds!

Kinda hacky, but it is very important to keep in mind how exact you really need the results.

• @Phillipp , thanks a lot for your help. But my graph have 1600 points ... and takes over 40 hours to do that .. Do you Know fortran ? can i minimize ( more ) the time ? any ideas ? May 5, 2013 at 15:52
• You can try to take that solve out of minroot so that you only have to solve this once. If that doesnt work, use Nsolve instead and use options that lead to speed and not more than the required accuracy. And then use ParrallelTable (after DistributeDefinitions["Global`"]) to create the Table with the values. And have a look especially the bold items 3 and 4 in my above post. May 7, 2013 at 8:02
• thanks for your help ... but can you rewrite the code using Parallel Table for me ? I try but I cant speed up , and I dont know nothing about using block or module. May 7, 2013 at 9:07
• I cannot, I am busy. ParallelTable is easy though: your code, then DistributeDefinitions and then instead using table to calculate all values you need use paralleltable. May 8, 2013 at 9:18