0
$\begingroup$

I am working in a triple integral, but the second integral (qQ) and the first too ( aA ) is a oscillatory funcion. The third integral ( os11 ) is like a gaussian, and is easy to integrate.

i am trying

\[Sigma] = 6/10;

minroot[gg_?NumericQ, bb_?NumericQ] := 
  Module[{b, g, rts, r}, b = Rationalize[bb, 0];
   g = Rationalize[gg, 0];
   rts = r /. 
     Solve[1 - (b/
           r)^2 - (g^-2)*(2/
             15*\[Sigma]^9 (1/(r - 1)^9 - 1/(r + 1)^9 - 
              9/(8 r)*(1/(r - 1)^8 - 
                 1/(r + 1)^8)) - \[Sigma]^3 (1/(r - 1)^3 - 
              1/(r + 1)^3 - (3/(2 r))*(1/(r - 1)^2 - 1/(r + 1)^2))) ==
        0, r];
   Max[Select[rts, Im[#] == 0 && # > 0 &]]];



aA[g_?NumberQ, b_?NumberQ] := 
  Pi - 2 b NIntegrate[
     1/(r^2*Sqrt[
         1 - (b/r)^2 - 
          g^(-2)*(2/
               15*\[Sigma]^9 (1/(r - 1)^9 - 1/(r + 1)^9 - 
                9/(8 r)*(1/(r - 1)^8 - 
                   1/(r + 1)^8)) - \[Sigma]^3 (1/(r - 1)^3 - 
                1/(r + 1)^3 - (3/(2 r))*(1/(r - 1)^2 - 
                   1/(r + 1)^2)))]), {r, minroot[g, b], \[Infinity]}, 
     Method -> {"LevinRule", 
       Method -> {"GaussKronrodRule", "Points" -> 11}}]; 

qQ[g_?NumberQ] := 
 NIntegrate[2*(1 - Cos[Re[aA[g, b]]]) b, {b, 0.001, 15}, Method -> {"LevinRule", 
 Method -> {"GaussKronrodRule", "Points" -> 11}}]

os11[T_, j_] := (1/T^3) NIntegrate[(g^5*Re[qQ[g]])/E^(g^2/T), {g,0.001, 5}, Method -> {"LevinRule",Method -> {"GaussKronrodRule", "Points" -> 11}}];

AbsoluteTiming[ParallelTable[os11[T, 5], {T, 0.1, 0.2, 0.1}]]

but appears to me the message that integral options should be ina built list ( something like that ), and takes HOURS and HOURS and DONT give me the final value os os11 = 5.44 and 4.68

any ideas ?

$\endgroup$
  • $\begingroup$ If you encounter a error message, it would be great to include it to your question. I suspect it is rather a warning about your highly oscillatory function...? $\endgroup$ – Lukas Nov 28 '15 at 17:56
1
$\begingroup$

This is rather a detailled comment than an answer. To me, it seems that you need to spend more time on characterizing your function(s). Putting a lot of AbsoluteTiming inside the evaluations reveals that the evaluation of qQ clearly is a bottleneck w.r.t. timing (so far, no comment about oscillatory character). What gave a speedup of a factor between 2-6 (depending on parameters) was using the Method Method -> {Automatic, "SymbolicProcessing" -> 0} instead of your LevinRule. However, evaluation is still slow and throws a lot of warnings regarding possible singularity, highly oscillatory behavior etc. So it is worth a shot to identify the problem.

data = ParallelTable[{b, g, 2*(1 - Cos[Re[aA[g, b]]]) b}, {b, 0.001, 
 15, 0.05}, {g, 0.001, 0.15, 0.005}]; // AbsoluteTiming
ListPlot3D[Flatten[data, 1], PlotRange -> All,  AxesLabel -> {"b", "g", "value"}]

yieldsenter image description here. You can clearly see that the function you are integrating within a call of qQ is well-behaved and smooth for almost all parameters. But unfortunately it is a) strongly peaked and b) additionally oscillatory along a single curve specified by certain combinations of g,b. Therefore, to get rid of the warnings you need to use a very high WorkingPrecision for the integration and additionally consider increasing MaxRecursion since there are also some convergence issues sometimes.

Also, I would recommend defining your functions like this f[x_]:=f[x]=... which avoids reevaluation. It does not help to speed up the code, but in case you want to debug after a first evaluation you get the results almost instantly since they are stored in memory from previous evaluation(s).

My advice is to separately investigate all functions you are integrating throughout your evaluation and see which one requires high WorkingPrecision because of a behavior similar to what is shown above. A first try (I aborted after couple of minutes) shows that for aA WorkingPrecision -> 20, MaxRecursion -> 20 seems to be sufficient, whereas you need WorkingPrecision -> 300, MaxRecursion -> 20 for qQ. Having said that, this will still need a lot of time to evaluate. Unfortunately I cannot help to speed up by orders of magnitude while achieving reliable results. My first attempt would be to see how it actually performs with the method named in the beginning together with a very high WorkingPrecision.

$\endgroup$
  • $\begingroup$ thanks for your advices. But how can i resolve the oscilations in [qQ] ? should i use a different integration method ?? Like Method -> {"GlobalAdaptive", Method -> "MultidimensionalRule"} or Levin Rule at some part of the problem? $\endgroup$ – Lucas G Leite F Pollito Nov 28 '15 at 22:37
  • $\begingroup$ and, how can i can define my functions different than this ? try to give me some simple example ... $\endgroup$ – Lucas G Leite F Pollito Nov 28 '15 at 22:38
  • $\begingroup$ @LukasGLeiteFPollito What do you mean with “avoid oscillations“? Do you mean avoid the warning? Because the oscillations seem to be part of your functions. So it's not possible to get rid of them... To get rid of the warning, do what the help suggests: increase working precision (as I said). I showed how you can define the functions. For instance, just write qQ[g_?NumberQ]=qQ[g]:=.... Then evaluate qQ[5] twice and observe the timings compared to those in your version. As said, this will not help to speed up the code but it will help if you want to analyze results/behavior later $\endgroup$ – Lukas Nov 29 '15 at 12:28
  • $\begingroup$ thanks again. But i spend about 10 hours and the integral still working, so i quite up. Do you now matlab ? mathematica takes too long for doing it. $\endgroup$ – Lucas G Leite F Pollito Nov 29 '15 at 13:52
  • $\begingroup$ @LucasGLeiteFPollito I know matlab, but not goodgood enough to advise you wwhat to do there. Also, it your function appears to be hard to integrate because of the strong oscillatory behavior at this one line. No idea how matlab can perform $\endgroup$ – Lukas Nov 29 '15 at 13:56

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.