# Speeding up the nested NIntegration consisted of Exp' and Bessel' functions

I got stuck with optimizing the nested numerical integration of function with a high number of Bessel functions (first and second kinds) summed and several exponents. I will appreciate any suggestions on how one can speed this integration up.

I am aiming at calculation of the nested integral F4 which can be written in the following general form: $$F4(x,\delta_F)= \int \int \int \int f_2(\gamma, \delta_F)\Big( \underbrace{\int f_1 (\theta,\phi,\gamma_F,\beta_F, x, \delta_F)\ \mathrm{d}\gamma_F\ \mathrm{d}\beta_F }_{F3(\theta,\phi, x,\delta_F)}\Big) \mathrm{d}\gamma \ \mathrm{d} \beta\ \mathrm{d} \theta \ \mathrm{d}\phi .$$ The complete routine, including integration and deduction the bulky function f1, is presented below:

(*function construction*)
Pn[j_]:=LegendreP[j,1,Cos[\[Theta]]]/Sin[\[Theta]];
Tn[j_]:=D[LegendreP[j,1,Cos[\[Theta]]],\[Theta]];
Fn[j_,ar_]:=ar (2SphericalBesselJ[j,ar]-I SphericalBesselY[j,ar]);
Xn[j_,ar_]:=ar (I SphericalBesselY[j,ar]);
DFn[j_,ar_]:=D[ar(2SphericalBesselJ[j,ar]-I SphericalBesselY[j,ar]),ar];
DXn[j_,ar_]:=D[ar (I SphericalBesselY[j,ar]),ar];
an[n_,m_,x_]=Sqrt[3](m Fn[n,mx]DFn[n,x]+Fn[n,x]DFn[n,mx])/(m Fn[n,mx]DXn[n,x]-Xn[n,x]DFn[n,mx])/.mx->m x;
bn[n_,m_,x_]=Sqrt[3](Fn[n,mx]DFn[n,x]-m Fn[n,x]DFn[n,mx])/(Fn[n,mx]DXn[n,x]-m Xn[n,x]DFn[n,mx])/.mx->m x;
ScIntDist=Exp[Sin[\[Gamma]]^2/Sin[\[Delta]]^2]Cos[\[Phi]] Sqrt[Cos[\[Gamma]]] Sin[\[Gamma]]Sin[\[Theta]t]Sum[(2 n+1)/(n(n+1)) (an[n,m,x]Pn[n]+bn[n,m,x]Tn[n]),{n,1,2}];
MT={Xx Cos[\[Beta]]^2 Cos[\[Gamma]]+Xx Sin[\[Beta]]^2+Cos[\[Beta]] (Yy (-1+Cos[\[Gamma]]) Sin[\[Beta]]-Zz Sin[\[Gamma]]),1/2 (Yy+Yy Cos[\[Gamma]]+2 (Yy Cos[2 \[Beta]]-Xx Sin[2 \[Beta]]) Sin[\[Gamma]/2]^2-2 Zz Sin[\[Beta]] Sin[\[Gamma]]),Zz Cos[\[Gamma]]+(Xx Cos[\[Beta]]+Yy Sin[\[Beta]]) Sin[\[Gamma]]};
Ar=TransformedField["Spherical"->"Cartesian",Re[ScIntDist],{r,\[Theta],\[Phi]}->{xx,yy,zz}]/.xx->MT[[1]]/.yy->MT[[2]]/.zz->MT[[3]];
ScIntDistTransformedr=TransformedField["Cartesian"->"Spherical",Ar,{Xx,Yy,Zz}->{r,\[Theta]t,\[Phi]f}];
Ai=TransformedField["Spherical"->"Cartesian",Im[ScIntDist],{r,\[Theta],\[Phi]}->{xx,yy,zz}]/.xx->MT[[1]]/.yy->MT[[2]]/.zz->MT[[3]];
ScIntDistTransformedi=TransformedField["Cartesian"->"Spherical",Ai,{Xx,Yy,Zz}->{r,\[Theta]t,\[Phi]f}];

(*nested integration computation*)
f1[\[Theta]tF_?NumericQ,\[Phi]fF_?NumericQ,\[Gamma]F_?NumericQ,\[Beta]F_?NumericQ,xF_?NumericQ,\[Delta]F_?NumericQ]:=ScIntDistTransformedr+I ScIntDistTransformedi/.\[Delta]->\[Delta]F/.\[Theta]t->\[Theta]tF/.\[Phi]f->\[Phi]fF/.\[Gamma]->\[Gamma]F/.\[Beta]->\[Beta]F/.x->xF/.r->1/.m->3/2;
f2[\[Gamma]F_?NumericQ,\[Delta]F_?NumericQ]:= Exp[Sin[\[Gamma]F]^2/Sin[\[Delta]F]^2]Sqrt[Cos[\[Gamma]F]];
F3[\[Theta]t_?NumericQ,\[Phi]f_?NumericQ,xF_?NumericQ,\[Delta]F_?NumericQ]:=NIntegrate[f1[\[Theta]t,\[Phi]f,\[Gamma]F,\[Beta]F,xF,\[Delta]F],{\[Gamma]F,0,\[Delta]F},{\[Beta]F,0,0.2\[Pi]},MaxRecursion->1,AccuracyGoal->2];
F4[xF_?NumericQ,\[Delta]F_?NumericQ]:=NIntegrate[f2[\[Gamma],\[Delta]F](Abs[F3[\[Theta]t,\[Phi]f,xF,\[Delta]F]]) ,{\[Theta]t,0,\[Delta]F},{\[Phi]f,0,2\[Pi]},{\[Gamma],0,\[Delta]F},{\[Beta],0,2\[Pi]},MaxRecursion->1,AccuracyGoal->2];

(*computation run, a probe point (200, 0.2)*)
F4[200,0.2]//AbsoluteTiming
Out[97]= {29.3105,0.0194077}


I have tried all the Strategies attributed to NIntegration, and found the best one is "GlobalAdaptive". In the end computations take about 48 hours. Therefore, in this post, I attended to keep the main constructs persisting in the computations but save your time for test runs. In fact, for F3 and F4 I normally use AccuracyGoal-> 3 and MaxRecursion ->3 with large terms in the summ for ScIntDist.

I was desperate to find any further optimization steps, and will appreciate any of your suggestions.

Following the discussion given in the post some, I guess, evident steps for the experienced users were done. In fact It was possible to speed up the code in approx. 20 times.

One was necessary a) avoid substitutions ./x-> exactly in the functions definitions. b) Try to define all the Bessel functions numerically prior to sending them to functions. Two modified lines was enough:

FuncAfterSubst=N[ScIntDistTransformedr+I ScIntDistTransformedi/.\[Delta]->\[Delta]F/.\[Theta]t->\[Theta]tF/.\[Phi]f->\[Phi]fF/.\[Gamma]->\[Gamma]F/.\[Beta]->\[Beta]F/.x->200/.r->1/.m->3/2];
f1[\[Theta]tF_?NumericQ,\[Phi]fF_?NumericQ,\[Gamma]F_?NumericQ,\[Beta]F_?NumericQ,\[Delta]F_?NumericQ]:=Evaluate[FuncAfterSubst];


Here is the code:

ClearAll["Global*"]
Pn[j_]:=LegendreP[j,1,Cos[\[Theta]]]/Sin[\[Theta]];
Tn[j_]:=D[LegendreP[j,1,Cos[\[Theta]]],\[Theta]];
Fn[j_,ar_]:=ar (2SphericalBesselJ[j,ar]-I SphericalBesselY[j,ar]);
Xn[j_,ar_]:=ar (I SphericalBesselY[j,ar]);
DFn[j_,ar_]:=D[ar(2SphericalBesselJ[j,ar]-I SphericalBesselY[j,ar]),ar];
DXn[j_,ar_]:=D[ar (I SphericalBesselY[j,ar]),ar];
an[n_,m_,x_]=Sqrt[3](m Fn[n,mx]DFn[n,x]+Fn[n,x]DFn[n,mx])/(m Fn[n,mx]DXn[n,x]-Xn[n,x]DFn[n,mx])/.mx->m x;
bn[n_,m_,x_]=Sqrt[3](Fn[n,mx]DFn[n,x]-m Fn[n,x]DFn[n,mx])/(Fn[n,mx]DXn[n,x]-m Xn[n,x]DFn[n,mx])/.mx->m x;
ScIntDist=Exp[Sin[\[Gamma]]^2/Sin[\[Delta]]^2]Cos[\[Phi]] Sqrt[Cos[\[Gamma]]] Sin[\[Gamma]]Sin[\[Theta]t]Sum[(2 n+1)/(n(n+1)) (an[n,m,x]Pn[n]+bn[n,m,x]Tn[n]),{n,1,2}];
MT={Xx Cos[\[Beta]]^2 Cos[\[Gamma]]+Xx Sin[\[Beta]]^2+Cos[\[Beta]] (Yy (-1+Cos[\[Gamma]]) Sin[\[Beta]]-Zz Sin[\[Gamma]]),1/2 (Yy+Yy Cos[\[Gamma]]+2 (Yy Cos[2 \[Beta]]-Xx Sin[2 \[Beta]]) Sin[\[Gamma]/2]^2-2 Zz Sin[\[Beta]] Sin[\[Gamma]]),Zz Cos[\[Gamma]]+(Xx Cos[\[Beta]]+Yy Sin[\[Beta]]) Sin[\[Gamma]]};
Ar=TransformedField["Spherical"->"Cartesian",Re[ScIntDist],{r,\[Theta],\[Phi]}->{xx,yy,zz}]/.xx->MT[[1]]/.yy->MT[[2]]/.zz->MT[[3]];
ScIntDistTransformedr=TransformedField["Cartesian"->"Spherical",Ar,{Xx,Yy,Zz}->{r,\[Theta]t,\[Phi]f}];
Ai=TransformedField["Spherical"->"Cartesian",Im[ScIntDist],{r,\[Theta],\[Phi]}->{xx,yy,zz}]/.xx->MT[[1]]/.yy->MT[[2]]/.zz->MT[[3]];
ScIntDistTransformedi=TransformedField["Cartesian"->"Spherical",Ai,{Xx,Yy,Zz}->{r,\[Theta]t,\[Phi]f}];
FuncAfterSubst=N[ScIntDistTransformedr+I ScIntDistTransformedi/.\[Delta]->\[Delta]F/.\[Theta]t->\[Theta]tF/.\[Phi]f->\[Phi]fF/.\[Gamma]->\[Gamma]F/.\[Beta]->\[Beta]F/.x->200/.r->1/.m->3/2];

(*nested integration computation*)
f1[\[Theta]tF_?NumericQ,\[Phi]fF_?NumericQ,\[Gamma]F_?NumericQ,\[Beta]F_?NumericQ,\[Delta]F_?NumericQ]:=Evaluate[FuncAfterSubst];
f2[\[Gamma]F_?NumericQ,\[Delta]F_?NumericQ]:= Exp[Sin[\[Gamma]F]^2/Sin[\[Delta]F]^2]Sqrt[Cos[\[Gamma]F]];
F3[\[Theta]t_?NumericQ,\[Phi]f_?NumericQ,\[Delta]F_?NumericQ]:=NIntegrate[f1[\[Theta]t,\[Phi]f,\[Gamma]F,\[Beta]F,\[Delta]F],{\[Gamma]F,0,\[Delta]F},{\[Beta]F,0,0.2\[Pi]},MaxRecursion->1,AccuracyGoal->2];
`