# Optimization of numerical integration

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


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

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 ]


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]


• Amazing, why such a difference in results? Have you change any other options? How long does it takes? – yarchik Aug 1 '19 at 13:51
• @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. – Henrik Schumacher Aug 1 '19 at 14:00
• @HenrikSchumacher Yes it seems to be an edge case! – Ulrich Neumann Aug 1 '19 at 14:01
• @yarchik Evaluation takes around 80s! – Ulrich Neumann Aug 1 '19 at 14:04
• Yes, it was the edge-case. My ContourPlot failed because $f=1$ is not a 1D contour, but a 2D region. – yarchik Aug 1 '19 at 16:28