How to remove error in NIntegrate

Question

I want to solev NIntegrate in two dimension. But when i run the code i got Throw::nocatch: Uncaught Throw[-HolonomicDifferentialRootReduceDumpy[NIntegrateLevinRuleDumpx]+(HolonomicDifferentialRootReduceDumpy^[Prime])[NIntegrateLevinRuleDumpx],NIntegrateLevinRuleDumpFastLookupHolonomicDifferentialEquation] this error. How can i overcome this error. I also attach my code there with Image code: In1:= Clear["Global`*"]

In[2]:= \[Gamma]21 = 0.1 ; \[Gamma]31 = 0.1;  \[Gamma]41 = 1;  w = 1; \
kp = 1; \[Beta] = 1;  \[CapitalDelta]2 = 0;  \[CapitalDelta]3 = 0 ;  \
\[CapitalDelta]1 = 1; 

In[3]:= \[Phi] = \[Pi]/6; \[CapitalDelta]33 = \[CapitalDelta]3 + kp*v; 

In[4]:=   \[CapitalDelta]11 =  \[CapitalDelta]1 + 
  kp*v ;  \[CapitalDelta]22 = \[CapitalDelta]2  + 
  kp*v;   d = 1; M1 = 4; N1 = 4 ; \[Lambda] = 0.25;

In[5]:= \[Phi] = \[Pi]/6;  L = 6 ;

In[6]:=  m1 = 0; n1 = 0;

In[7]:= \[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1 ; r = 1; l = 3;

In[8]:= \[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)] (r/w)^l;

In[9]:= \[CapitalOmega]3 = \[CapitalOmega]f1*
  Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22;  \[CapitalLambda] = 1; 

In[10]:= \[CapitalDelta]13 = \[CapitalDelta]11  - \[CapitalDelta]33 ; 

In[11]:= \[CapitalOmega]2 = \[CapitalOmega]20*(Sin[\[Pi] x/\
\[CapitalLambda]] + Sin[\[Pi] y/\[CapitalLambda]]);

In[12]:= \[Chi]p = (
  I*(\[Gamma]21  - I*\[CapitalDelta]12)*(\[Gamma]31 -  
     I*\[CapitalDelta]13))/((\[Gamma]41  - 
       I*   \[CapitalDelta]11)*(\[Gamma]21  - 
       I*\[CapitalDelta]12)*(\[Gamma]31  - 
       I*\[CapitalDelta]13) + (\[Gamma]31  - 
       I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21  - 
       I*\[CapitalDelta]12)* \[CapitalOmega]3^2);

In[13]:= a1 = ComplexExpand[Im[\[Chi]p]];

In[14]:= aa[v_] := 1/(d*Sqrt[2*\[Pi]])*Exp[-v^2/(2*d^2)];

In[15]:= aa[.1]

Out[15]= 0.396953

In[16]:= aaa[x_, y_] := NIntegrate[a1*aa[v], {v, -10, 10}];

In[17]:= b1 = ComplexExpand[Re[\[Chi]p]];

In[18]:= cc1[x_, y_] := NIntegrate[b1*aa[v], {v, -10, 10}];

In[19]:= tt[x_, y_] := Exp[-aaa[x, y]*L];

In[20]:= TT1[x_, y_] := Exp[-aaa[x, y]*L + I*cc1[x, y]];

In[21]:= q[\[Theta]1_, x_, \[Theta]2_, y_] := 
  Exp[-aaa[x, y]*L + I*cc1[x, y]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*x*Sin[\[Theta]1])/\[Lambda]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*y*Sin[\[Theta]2])/\[Lambda]];

In[22]:= dd[\[Theta]1_, \[Theta]2_] := 
 NIntegrate[
  q[\[Theta]1, x, \[Theta]2, y], {x, 0, 0.001}, {y, 0, 0.001}]


In[23]:= y[\[Theta]1_, \[Theta]2_] := (Abs[dd[\[Theta]1, \[Theta]2]])^2

In[24]:= s[\[Theta]1_, \[Theta]2_] := (Sin[(
    M1*Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]])^2/(
  M1^2*(Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]])^2)*(Sin[(
    N1*Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]])^2/(
  N1^2*(Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]])^2)

In[25]:= z[\[Theta]1_, \[Theta]2_] := 
 y[\[Theta]1, \[Theta]2]*s[\[Theta]1, \[Theta]2]


In[26]:= DensityPlot[
 z[\[Theta]1, \[Theta]2], {\[Theta]1, -0.6, 0.6}, {\[Theta]2, -0.6, 
  0.6}, PlotRange -> All, PlotPoints -> 150, 
 ColorFunction -> "Rainbow", 
 PlotLegends -> 
  Placed[BarLegend[Automatic, LegendMarkerSize -> 280], Right], 
 ImageSize -> 300, Background -> Transparent, FrameStyle -> Black, 
 LabelStyle -> {Black, FontSize -> 18}, RotateLabel -> True]

Welcome to the Mathematica Stack Exchange. Please fix the code to remove In[*] entries. Also provide the definitions for x and v. Thanks. — Syed, May 28 at 6:41

Alex Trounev · Accepted Answer · 2022-05-29 15:12:14Z

With using ?NumericQ we have

\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1; 
\[Phi] = Pi/6; \[CapitalDelta]33 = \[CapitalDelta]3 + kp*v; 
\[CapitalDelta]11 = \[CapitalDelta]1 + 
  kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 + 
  kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25; 
\[Phi] = Pi/6; L = 6; 
m1 = 0; n1 = 0; 
\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 3; 
\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l; 
\[CapitalOmega]3 = \[CapitalOmega]f1*
  Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1; 
\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33; 
\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] + 
     Sin[Pi*(y/\[CapitalLambda])]); 
\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 - 
       I*\[CapitalDelta]13))/((\[Gamma]41 - 
        I*\[CapitalDelta]11)*(\[Gamma]21 - 
        I*\[CapitalDelta]12)*(\[Gamma]31 - 
        I*\[CapitalDelta]13) + (\[Gamma]31 - 
        I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 - 
        I*\[CapitalDelta]12)*\[CapitalOmega]3^2); 
a1 = ComplexExpand[Im[\[Chi]p]]; 
aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)]; 
aaa[x1_?NumericQ, y1_?NumericQ] := 
  NIntegrate[a1*aa[v] /. {x -> x1, y -> y1}, {v, -10, 10}]; 
b1 = ComplexExpand[Re[\[Chi]p]]; 
cc1[x1_?NumericQ, y1_?NumericQ] := 
  NIntegrate[b1*aa[v] /. {x -> x1, y -> y1}, {v, -10, 10}]; 
tt[x_, y_] := Exp[(-aaa[x, y])*L]; 
TT1[x_, y_] := Exp[(-aaa[x, y])*L + I*cc1[x, y]]; 
q[\[Theta]1_, x_, \[Theta]2_, y_] := 
  Exp[(-aaa[x, y])*L + I*cc1[x, y]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*x*Sin[\[Theta]1])/\[Lambda]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*y*Sin[\[Theta]2])/\[Lambda]]; 
dd[\[Theta]1_?NumericQ, \[Theta]2_?NumericQ] := 
 NIntegrate[
  q[\[Theta]1, x, \[Theta]2, y], {x, 0, 0.001}, {y, 0, 0.001}]
y[\[Theta]1_, \[Theta]2_] := Abs[dd[\[Theta]1, \[Theta]2]]^2
s[\[Theta]1_, \[Theta]2_] := (Sin[(M1*Pi*\[CapitalLambda]*
         Sin[\[Theta]1])/\[Lambda]]^2/(M1^2*
      Sin[(Pi*\[CapitalLambda]*
           Sin[\[Theta]1])/\[Lambda]]^2))*(Sin[(N1*
         Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2/
        (N1^2*Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2))
z[\[Theta]1_, \[Theta]2_] := 
 y[\[Theta]1, \[Theta]2]*s[\[Theta]1, \[Theta]2]

DensityPlot[
 Evaluate[z[\[Theta]1, \[Theta]2]], {\[Theta]1, -0.6, 
  0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 50, 
 ColorFunction -> "Rainbow", 
   PlotLegends -> 
  Placed[BarLegend[Automatic, LegendMarkerSize -> 280], Right], 
 ImageSize -> 300, Background -> Transparent, FrameStyle -> Black, 
   LabelStyle -> {Black, FontSize -> 18}, RotateLabel -> True]

Update 1. In a case of integration limits {v, -Infinity, Infinity} and {x, 0, 1}, {y, 0, 1} we modified code as follows (we also update function s to exclude singularities)

\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1; 
\[Phi] = Pi/6; \[CapitalDelta]33 = \[CapitalDelta]3 + kp*v; 
\[CapitalDelta]11 = \[CapitalDelta]1 + 
  kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 + 
  kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25; 
\[Phi] = Pi/6; L = 6; 
m1 = 0; n1 = 0; 
\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 3; 
\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l; 
\[CapitalOmega]3 = \[CapitalOmega]f1*
  Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1; 
\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33; 
\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] + 
     Sin[Pi*(y/\[CapitalLambda])]); 
\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 - 
       I*\[CapitalDelta]13))/((\[Gamma]41 - 
        I*\[CapitalDelta]11)*(\[Gamma]21 - 
        I*\[CapitalDelta]12)*(\[Gamma]31 - 
        I*\[CapitalDelta]13) + (\[Gamma]31 - 
        I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 - 
        I*\[CapitalDelta]12)*\[CapitalOmega]3^2); 
a1 = ComplexExpand[Im[\[Chi]p]]; 
aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)]; 

aaa[x1_?NumericQ, y1_?NumericQ] := 
  NIntegrate[
   a1*aa[v] /. {x -> x1, y -> y1}, {v, -Infinity, Infinity}]; 
b1 = ComplexExpand[Re[\[Chi]p]]; 
cc1[x1_?NumericQ, y1_?NumericQ] := 
  NIntegrate[
   b1*aa[v] /. {x -> x1, y -> y1}, {v, -Infinity, Infinity}]; 

tt[x_, y_] := Exp[(-aaa[x, y])*L]; 
TT1[x_, y_] := Exp[(-aaa[x, y])*L + I*cc1[x, y]]; 
q[\[Theta]1_, x_, \[Theta]2_, y_] := 
  Exp[(-aaa[x, y])*L + I*cc1[x, y]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*x*Sin[\[Theta]1])/\[Lambda]]*
   Exp[(-2*Pi*I*\[CapitalLambda]*y*Sin[\[Theta]2])/\[Lambda]]; 

dd[\[Theta]1_?NumericQ, \[Theta]2_?NumericQ] := 
 NIntegrate[q[\[Theta]1, x, \[Theta]2, y], {x, 0, 1}, {y, 0, 1}, 
  PrecisionGoal -> 2, AccuracyGoal -> 2]

y[\[Theta]1_, \[Theta]2_] := Abs[dd[\[Theta]1, \[Theta]2]]^2
s[\[Theta]1_, \[Theta]2_] := 
  If[Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]]^2 <= 10^-10, 
    1, (Sin[(M1*Pi*\[CapitalLambda]*
           Sin[\[Theta]1])/\[Lambda]]^2/(M1^2*
        Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]1])/\[Lambda]]^2))]*
   If[Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2 <= 10^-10,
     1, (Sin[(N1*Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2/
          (N1^2*
        Sin[(Pi*\[CapitalLambda]*Sin[\[Theta]2])/\[Lambda]]^2))];
z[\[Theta]1_, \[Theta]2_] := 
 y[\[Theta]1, \[Theta]2]*s[\[Theta]1, \[Theta]2]

To plot z we first prepare interpolation function

lst = Table[{\[Theta]1, \[Theta]2, {y[\[Theta]1, \[Theta]2]}}, {\
\[Theta]1, -0.6, 0.6, .1}, {\[Theta]2, -0.6, 0.6, .1}]

Y = Interpolation[Flatten[lst, 1]];

Finally we plot functions s, Y, z

{DensityPlot[
  Evaluate[s[\[Theta]1, \[Theta]2]], {\[Theta]1, -0.6, 
   0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 100, 
  ColorFunction -> "Rainbow", 
    PlotLegends -> Automatic], 
 DensityPlot[Y[t1, t2], {t1, -.6, .6}, {t2, -.6, .6}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  PlotLegends -> Automatic, PlotPoints -> 100], 
 DensityPlot[
  s[\[Theta]1, \[Theta]2] Y[\[Theta]1, \[Theta]2], {\[Theta]1, -0.6, 
   0.6}, {\[Theta]2, -0.6, 0.6}, PlotRange -> All, PlotPoints -> 100, 
  ColorFunction -> "TemperatureMap", 
    PlotLegends -> 
   Placed[BarLegend[Automatic, LegendMarkerSize -> 280], Right], 
  ImageSize -> 300, Background -> Transparent, FrameStyle -> Black, 
    LabelStyle -> {Black, FontSize -> 18}, RotateLabel -> True]}

Update 2. This code is answer to the question on this page

Clear["Global`*"]

plotset = {FrameStyle -> Directive[Thickness[0.004]], 
   TicksStyle -> Directive[Black, 18]};
plotset2 = 
  FrameTicksStyle -> {{Directive[Black, 18], 
     Directive[FontOpacity -> 0, FontSize -> 0]}, {Directive[Black, 
      18], Directive[FontOpacity -> 0, FontSize -> 0]}};

\[Gamma]21 = 0.1; \[Gamma]31 = 0.1; \[Gamma]41 = 1; w = 1; kp = 1; \
\[Beta] = 1; \[CapitalDelta]2 = 0; \[CapitalDelta]3 = 0; \
\[CapitalDelta]1 = 1;

\[CapitalDelta]33 = \[CapitalDelta]3 + kp*v;

\[CapitalDelta]11 = \[CapitalDelta]1 + 
  kp*v; \[CapitalDelta]22 = \[CapitalDelta]2 + 
  kp*v; d = 1; M1 = 4; N1 = 4; \[Lambda] = 0.25;

L = 6;

m1 = 0; n1 = 0;

\[CapitalOmega]1 = 1; \[CapitalOmega]20 = 1; r = 1; l = 1;

\[CapitalOmega]f1 = 2*\[CapitalOmega]1*Exp[-(r^2/w^2)]*(r/w)^l;

\[CapitalOmega]3 = \[CapitalOmega]f1*
  Cos[l*\[Phi]]; \[CapitalDelta]12 = \[CapitalDelta]11 - \
\[CapitalDelta]22; \[CapitalLambda] = 1;

\[CapitalDelta]13 = \[CapitalDelta]11 - \[CapitalDelta]33;

\[CapitalOmega]2 = \[CapitalOmega]20*(Sin[Pi*(x/\[CapitalLambda])] + 
     Sin[Pi*(y/\[CapitalLambda])]);

\[Chi]p = (I*(\[Gamma]21 - I*\[CapitalDelta]12)*(\[Gamma]31 - 
       I*\[CapitalDelta]13))/((\[Gamma]41 - 
        I*\[CapitalDelta]1)*(\[Gamma]21 - 
        I*\[CapitalDelta]12)*(\[Gamma]31 - 
        I*\[CapitalDelta]13) + (\[Gamma]31 - 
        I*\[CapitalDelta]13)*\[CapitalOmega]2^2 + (\[Gamma]21 - 
        I*\[CapitalDelta]12)*\[CapitalOmega]3^2);

a1 = ComplexExpand[Im[\[Chi]p]];

aa[v_] := (1/(d*Sqrt[2*Pi]))*Exp[-v^2/(2*d^2)];

aaa[x1_?NumericQ, y1_?NumericQ, phi_?NumericQ] := 
  NIntegrate[
   a1*aa[v] /. {x -> x1, y -> y1, \[Phi] -> phi}, {v, -Infinity, 
    Infinity}];

b1 = ComplexExpand[Re[\[Chi]p]];

cc1[x1_?NumericQ, y1_?NumericQ, phi_?NumericQ] := 
  NIntegrate[
   b1*aa[v] /. {x -> x1, y -> y1, \[Phi] -> phi}, {v, -Infinity, 
    Infinity}];

tt[x_, y_, \[Phi]_] := Exp[(-aaa[x, y, \[Phi]])*L];

TT1[x_, y_, \[Phi]_] := 
  Exp[(-aaa[x, y, \[Phi]])*L + I*cc1[x, y, \[Phi]]];

q[x_, y_, \[Phi]_] := 
  Exp[-aaa[x, y, \[Phi]]*L + I*cc1[x, y, \[Phi]]]*Exp[-2*Pi*I*m1*x]*
   Exp[-2*Pi*I*n1*y];


pp1[\[Phi]_?NumericQ] := 
  NIntegrate[q[x, y, \[Phi]], {x, 0, 1}, {y, 0, 1}, 
   PrecisionGoal -> 2, AccuracyGoal -> 2];

yy[\[Phi]_] := (Abs[pp1[\[Phi]]])^2;


Y = Interpolation[Table[{x, yy[x]}, {x, 0, 4, .1}]]

Visualization

p1 = Plot[Y[\[Phi]], {\[Phi], 0, 4}, PlotRange -> {All, All}, 
   PlotStyle -> {Blue, Thickness[0.007]}, GridLines -> Automatic, 
   Frame -> True, 
   FrameLabel -> {Style["", 18, Bold], Style["", 18, Bold]}, 
   PlotLegends -> 
    Placed[LineLegend[{"\!\(\*SubscriptBox[\(I\), \
\(p\)]\)(\!\(\*SubsuperscriptBox[\(\[Theta]\), \(x\), \
\(0\)]\),\!\(\*SubsuperscriptBox[\(\[Theta]\), \(y\), \(0\)]\))"}, 
      LegendMarkers -> Automatic, 
      LegendMarkerSize -> {{30, 25}}], {After, Top}], 
   Evaluate@plotset, Evaluate@plotset2, Axes -> True] // 
  AbsoluteTiming

Alex Trounev, Thank you so much for your answer. I have one more question: when I want to DensityPlot, it takes too long. Could you please suggest a way for me to speed up calculation? — abdul wahab, May 29 at 1:41
because I want the Integrate limit from, dd[[Theta]1_?NumericQ, [Theta]2_?NumericQ] := NIntegrate[ q[[Theta]1, x, [Theta]2, y], {x, 0, 1}, {y, 0,1}] — abdul wahab, May 29 at 1:45
I see, that my answer is only part of answer. Is these limits final or you suppose to extend it more? — Alex Trounev, May 29 at 2:18
Thanks for your reply, the limits for: aaa[x1_?NumericQ, y1_?NumericQ] := NIntegrate[a1*aa[v] /. {x -> x1, y -> y1}, {v, -Infinity, Infinity}]; cc1[x1_?NumericQ, y1_?NumericQ] := NIntegrate[b1*aa[v] /. {x -> x1, y -> y1}, {v, - Infinity, Infinity}]; and for dd[[Theta]1_?NumericQ, [Theta]2_?NumericQ] := NIntegrate[ q[[Theta]1, x, [Theta]2, y], {x, 0, 1}, {y, 0, 1}]. These are the final limits. For first case limits maybe less but in last case it must be ( {x, 0, 1}, {y, 0, 1}). — abdul wahab, May 29 at 2:25
What commands do you recommend I use to speed up my calculations? Because density plotting takes a long time. — abdul wahab, May 29 at 2:30

Stack Exchange Network

current community

your communities

more stack exchange communities

How to remove error in NIntegrate

1 Answer 1

Your Answer

Not the answer you're looking for? Browse other questions tagged plotting numerical-integration or ask your own question.

Hot Network Questions

How to remove error in NIntegrate

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged plotting numerical-integration or ask your own question.

Related

Hot Network Questions