2
$\begingroup$

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.

$\endgroup$

1 Answer 1

0
$\begingroup$

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];
F4[\[Delta]F_?NumericQ]:=NIntegrate[f2[\[Gamma],\[Delta]F](Abs[F3[\[Theta]t,\[Phi]f,\[Delta]F]]) ,{\[Theta]t,0,\[Delta]F},{\[Phi]f,0,2\[Pi]},{\[Gamma],0,\[Delta]F},{\[Beta],0,2\[Pi]},MaxRecursion->1,AccuracyGoal->2,Method->"GlobalAdaptive"];
F4[0.2]//AbsoluteTiming 
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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