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?
