2
$\begingroup$

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

Improve this question
$\endgroup$
5
  • $\begingroup$ Try taking the Solve outside of minroot! $\endgroup$
    – PlatoManiac
    Commented May 3, 2013 at 15:00
  • $\begingroup$ @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. $\endgroup$
    – Daniel Lichtblau
    Commented May 3, 2013 at 15:14
  • $\begingroup$ 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. $\endgroup$
    – Daniel Lichtblau
    Commented May 6, 2013 at 14:54
  • $\begingroup$ @DanielLichtblau my computer is octacore, an i7 processor ... how can i use the 8 nucleos ?? $\endgroup$
    – Lucas G Leite F Pollito
    Commented May 6, 2013 at 15:40
  • $\begingroup$ I'd go with ParallelTable maybe. Check documntation on parallel operations in Mathematica. $\endgroup$
    – Daniel Lichtblau
    Commented May 6, 2013 at 16:25

1 Answer 1

Reset to default
7
$\begingroup$

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.
  • Read Read stackexchange and find answers like: Methods to speed up numerical NDSolve, NIntegrate,

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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ @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 ? $\endgroup$
    – Lucas G Leite F Pollito
    Commented May 5, 2013 at 15:52
  • $\begingroup$ 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. $\endgroup$
    – Philipp
    Commented May 7, 2013 at 8:02
  • $\begingroup$ 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. $\endgroup$
    – Lucas G Leite F Pollito
    Commented May 7, 2013 at 9:07
  • $\begingroup$ 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. $\endgroup$
    – Philipp
    Commented May 8, 2013 at 9:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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