1
$\begingroup$

I would like to ask for a help from community in the numerical computation of an integral and in the preparation of a contour plot. I have already done some work in this direction, but everything is so slow and the plot is not smooth as I was expecting.

At first, we define the 4 compiled helper functions. I believe, there is not much that can be done on this step

f6c=Compile[{{x,_Complex},{a,_Real}},-((2 a x)/(π (a^2+x^2)^2)),CompilationTarget->"C"]
g6c=Compile[{{x,_Complex},{a,_Real}},(a^2-x^2)/(π (a^2+x^2)^2),CompilationTarget->"C"]
f7c=Compile[{{x,_Complex},{a,_Real}},(8 a^5)/(3 π (a^2+x^2)^3),CompilationTarget->"C"]
g7c=Compile[{{x,_Complex},{a,_Real}},(x (15 a^4+10 a^2 x^2+3 x^4))/(3 a^6 π (1+x^2/a^2)^3),CompilationTarget->"C"]

Now, comes a more difficult function definition, which is also compiled

dEYc=Compile[{{w,_Real},{c0,_Real},{c2,_Real}},
      Module[{α,rs,p,a,η,iDFx,iDFy},
       p=1.88104430514;
       a=0.1;
       η=0.00001;
       iDFx=1.+π/2 (1.-c2-c0) p(g7c[w+I η-p,a]-g7c[w+I η+p,a])+π/4 c2 p(g7c[w+I η-2p,a]-g7c[w+I η+2p,a])-c0 π p^2 g6c[w+I η,a];
       iDFy=(-(π/2)(1.-c2-c0)p(f7c[w+I η-p,a]-f7c[w+I η+p,a])-π/4 c2 p(f7c[w+I η-2p,a]-f7c[w+I η+2p,a])+c0 π p^2 f6c[w+I η,a]) Sign[Im[w]+η];
       Im[1./(iDFx+I iDFy)] 
      ],
      CompilationTarget->"C"]

Next, we define a wrapper function to make sure that integral is not evaluated symbolically

fInt[w_?NumericQ,x_?NumericQ,y_?NumericQ]:=Module[{p=1.88104430514},2/(π p^2) w dEYc[w,x,y]]

Finally, we define the numerical 1d integral of this function from 0 to infinity

fSum[x_?NumericQ,y_?NumericQ]:=NIntegrate[fInt[w,x,y],{w,0,∞},Method->{Automatic,"SymbolicProcessing"->0},MinRecursion->6,MaxRecursion->21]

I think, this part may profit from optimization. As a test, we have

fSum[0.00045,0.2]
Out[8]= 0.999391

which is returned almost instantly. Unfortunately the contour plot takes a lot of time (5 sec.) and looks ugly

ContourPlot[fSum[x,y]==1.0,{x,5 10^-5,5 10^-4},{y,10^-2,2 10^-1},FrameLabel->{"Subscript[c, 0]","Subscript[c, 2]"},ImageSize->400,PlotRange->All,PlotPoints->5]//Timing

Contour plot of a complicated integral

I am expecting a single line here, almost linear dependence. My question is, how can I optimize these calculations?

$\endgroup$
2
$\begingroup$

It looks like there is "numerically" only a region fSum[x, y]==1

RegionPlot[fSum[x, y] > 0.999, {x, 5 10^-5, 5 10^-4}, {y, 10^-2, 2 10^-1},PlotRange -> All ]

enter image description here

It seems to be an "edgecase" ( see comments, thanks @HenrikSchumacher )

ContourPlot[ fSum[x, y] , {x, 5 10^-5, 5 10^-4}, {y, 10^-2, 2 10^-1},PlotRange -> All , Contours -> {0.9, 0.99, 0.999, .9999},MaxRecursion -> 4]

enter image description here

$\endgroup$
  • $\begingroup$ Amazing, why such a difference in results? Have you change any other options? How long does it takes? $\endgroup$ – yarchik Aug 1 at 13:51
  • $\begingroup$ @UlrichNeumann Either that or there are no solutions at all. The plot Plot3D[ {fSum[x, y], 1}, {x, 5 10^-5, 5 10^-4}, {y, 10^-2, 2 10^-1}, PlotRange -> {1 - 0.1, 1 + 0.01}, PlotPoints -> 50 ] tells me that this is really an edge-case. ContourPlot has only a chance to create a good plot when the level-1-surface intersects the graph of the function transversally. Here, it does not this is clear from the plot. I guess the question whether there is a nonempty intersection or not can only be settled by further analysis. $\endgroup$ – Henrik Schumacher Aug 1 at 14:00
  • $\begingroup$ @HenrikSchumacher Yes it seems to be an edge case! $\endgroup$ – Ulrich Neumann Aug 1 at 14:01
  • $\begingroup$ @yarchik Evaluation takes around 80s! $\endgroup$ – Ulrich Neumann Aug 1 at 14:04
  • $\begingroup$ Yes, it was the edge-case. My ContourPlot failed because $f=1$ is not a 1D contour, but a 2D region. $\endgroup$ – yarchik Aug 1 at 16:28

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.