I have a function sigma[qr,q0] of two variables with a sharp delta-like peak that cannot be numerically integrated by NIntegrate, like this:

enter image description here

Since the upper and lower regions excluding the peak are flat enough, I want to integrate them separately up to an epsilon around the peak by NIntegrate. The peak itself can be approximated locally in that epsilon strip around the peak (e.g. peak position $\pm$ 15 q0).

Therefore I have defined an interpolating function PeakFit[qr] like:

Peaklist = Table[{qr, PeakPosition[qr]}, {qr, 0, 950}];
PeakInterpol = Interpolation[Peaklist, InterpolationOrder -> 1];
PeakFit[p_?NumericQ] := PeakInterpol[p];

that returns me the position of the peak along q0 for a given qr. This function works just perfect.

Now I'm getting to my actual problem: To integrate e.g. the lower part in the interval {qr, 0, 400} and {q0, 0, peak position -15} I want to create a region in {qr,q0} over that Mathematica should integrate sigma[qr,q0]. I create the region like this:

rgn1 = 
ImplicitRegion[0 <= qr <= 400 && q0 <= (PeakFit[qr] - 15) && 0 <= q0 <= 400
, {qr,q0}];

which does exactly what I intended:

enter image description here

But giving that region to the Plot3D command (I use the Plot as a cross check to see whether Mathematica chooses the right region before directly handing it to NIntegrate) Mathematica returns an error:

In[58]:=Plot3D[sigma[aa, bb], Element[{aa, bb}, rgn1], PlotPoints -> 150, 
 AxesLabel -> Automatic]

Out[58]= Plot3D[sigma[aa, bb], Element[{aa, bb}, rgn1], PlotPoints -> 150, 
 AxesLabel -> Automatic]



enter image description here

Surprisingly it works for simpler regions not including the interpolating function PeakFit[p] as upper limit. Furthermore Mathematica can't give me the dimension of rgn1

In[74]:= RegionDimension[rgn1]
Out[74]= RegionDimension[ImplicitRegion[0 <= qr <= 400 && q0 <= -15 + PeakFit[qr] && 0 <= q0 <= 400, {qr, q0}]]

To me it seems to be an issue with the interpolating function even though Mathematica is able to plot the region that I wanted.

Does anyone have a clue why Mathematica is able to create the right region but is not able to use it in the Plot3D command (or NIntegrate)?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.