3
$\begingroup$

Can someone point me in the right direction? I 've been puzzling with this for weeks...

I have a huge function (of the form f(y,ξ)sin(ξx) ) to integrate in ξ in region (0, Infinity) and I use NIntegrate with the strategy "DoubleExponentialOscillatory" which is very fast and efficient:

U1X[x_, y_] := NIntegrate[f(y,ξ)*Sin[ξ*x], {ξ, 0, Infinity}, Method->"DoubleExponentialOscillatory", WorkingPrecision -> 40, MaxRecursion -> 40]

After this calculation I get a new function, say F(x,y) and I calculate a point to check the time. It is below 1/2 second:

F[1/100, 1/1000] // Timing
{0.4368028, 1.12454848442106131564112930702}

So I proceed to make a List plot like this:

ff = ParallelTable[{x, F[x, 1/1000]/P}, {x, 1/10000, 10, 1/30}] // Chop;
fff = MapThread[{-#1, -#2} &, {ff[[1 ;; All, 1]], ff[[1 ;; All, 2]]}];
f1 = Join[ff, fff];
ListPlot[{f1}]

The above procedure gives me the following errors:

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[\[Ellipsis], _SystemException].

Parallel`Developer`QueueRun::hmm: Received unexpected result SystemException[MemoryAllocationFailure] from KernelObject[2,local].

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[\[Ellipsis], _SystemException].

Parallel`Developer`QueueRun::hmm: Received unexpected result SystemException[MemoryAllocationFailure] from KernelObject[1,local].

In order to overcome this issue, I changed integration strategy to "ExtrapolatingOscillatory", but the time to calculate one point of F(x,y) rises to 85 seconds!

F[1/100, 1/1000] // Timing
{85.5353, 1.12454848442106131564112930702}

The strange thing is that now the above procedure to plot F(x,y) completes with no errors BUT it needs a huge amount of time (almost a day)!!! If I significantly reduce the estimation points to {x, 1/10000, 2, 1/5} it drops to an hour or so:

ff = ParallelTable[{x/h, μ2*U1X[x, 1/10000]/P}, {x, 1/10000, 2, 
 1/5}] // Chop;
fff = MapThread[{-#1, -#2} &, {ff[[1 ;; All, 1]], ff[[1 ;; All, 2]]}];
f1 = Join[ff, fff];
ListPlot[{f1}]

plot

Since I have to make 36 such plots and 36 (2D) contour plots with a similar to the above procedure, which need even more time to complete, it is obvious that "ExtrapolatingOscillatory" method is completely inefficient...

$\endgroup$
2
$\begingroup$

I think I have pinpointed the source of ListPlot's inability to complete the requested above plotting of the F(x,y). It was the highly oscillating behaviour of this function around zero. Also I needed to define that ξ > 0.

Nevertheless, I have managed to plot it for {x, 1/10000, 10, 1/10} but first I had to use FullSimplify to make it smaller and perhaps remove any complicated parts in it.

The result of the plot of the numerical integration of the function with the method: "DoubleExponentialOscillatory" is shown in the picture below. It came out rather quicly after all the simplifications in 4-5 minutes.

enter image description here

Then I did an integration of the same function with the method "ExtrapolatingOscillatory" and plotted the result for {x, 1/10000, 10, 1/10}. The plot is shown below.

enter image description here

As you can observe it is smoother at some points hence in my opinion the numerical integration with "ExtrapolatingOscillatory" method is more accurate. This plot took more than an hour to complete.

Mathematica also took about 20 minutes to complete FullSimplify of the function so I guess if accuracy is your objective you will accept the extra time expence.

Otherwise, using "DoubleExponentialOscillatory" in NIntegrate will produce very fast and acceptable results.

What I noticed today is that I used as integration limit the infinity for ξ e.g. {ξ,0,Infinity}. That was a worse slowing down factor than the simplification of the integrand. When I set a specific large value as the integration limit i.e. {ξ,0,1000000} the calculation finishes in minutes!!

And since the integration is numerical using a large number instead of infinity I guess is completely acceptable and justified accuracy wise.

$\endgroup$

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.