# Numerical integration of a function with several variables

I need to calculate the next integral

$$F_l(R)=\int_0^{\infty} {dk k^2j_l(k R)*2*\exp[-0.36k^2](k)\rho_l(k)(18/k^2)P_l(\cos{\beta})}$$

where ρl(k) is given by

$$\rho_l(k)=\int_{0}^{\infty}{dr r^2 j_l(k r)\left(\int_{0}^{\pi}{2\pi\rho(r,\theta)Y_{l0}(\theta)\sin{\theta}d\theta}\right)}$$

At the same time ρ(r,θ) is given by

$$\rho(r,\theta)=\rho_0/(1+\exp[(r-1.07Ad^{1/3}[1+b_2Y_{20}(\theta)+b_4Y_{40}(\theta)])/0.54]$$

In order to calculate this integral I put everything in a single function and then integrate all variables θ, r and k for given value of l. In all next equations Ad=178, Zd=80, b2=-0.113 and b4=-0.026; we have

Rd[th_, Ad_, b2_, b4_] := 1.07 Ad^(1/3)(1 + b2*SphericalHarmonicY[2,0,
th,0]+b4*SphericalHarmonicY[4, 0, th, 0])


and then the function

fd[r_, th_, Ad_, b2_, b4_] := 1/(1 + Exp[(r - Rd[th, Ad, b2, b4])/0.54])


The function Ddcd[r, Ad, Zd, th, b2, b4] (which is ρ(r,θ)) is given by

Ddcd[r_, th_, Ad_, Zd_, b2_, b4_] :=   Zd/(NIntegrate[      2 Pi*rr^2*fd[rr,
tth, Ad, b2, b4]*Sin[tth], {rr, 0,        Infinity}, {tth, 0, Pi},
Method -> {Automatic, "SymbolicProcessing" -> 0}])*   fd[r, th, Ad, b2, b4]


The integrand of Fl(R) is

fcd[R_, k_, r_, th_, Ad_, Zd_, b2_, b4_, l_] :=
k^2*SphericalBesselJ[l, k R]*(2*Exp[-0.36^2])*
(18/k^2)*(2 Pi*r^2*SphericalBesselJ[l, k r]*
Ddcd[r, Ad, Zd, th, b2, b4])*SphericalHarmonicY[l, 0, th, 0]*
Sin[th]


The integral I am trying to calculate is given by the code

FR0[R_, Ad_, Zd_, b2_, b4_, l_] =
Interpolation[
Table[{R,
NIntegrate[
fcd[R, k, r, th ,Ad ,Zd ,b2 ,b4 ,l], {k, 0, Infinity}, {r, 0, Infinity}, {th,
0, Pi}, Method ->
{Automatic, "SymbolicProcessing" -> 0}, PrecisionGoal -> 3, AccuracyGoal ->
2]}, {R, 0.01, 40, 0.1}]]
[R]


My problem is the integrals for l=0,2,4 are taking a lot of time (hours). How can I make the calculations faster ? Is there anyway to modify the integral so the calculations are faster ?Thank you in advance

• Have you tried Monte-Carlo methods? Method -> {"MonteCarlo", "MaxPoints" -> 10^6}. Also it would be helpful if you could write your constants and functions in a top-down order, and I'm getting an error Raw object ... cannot be used as an iterator – flinty Jul 5 '20 at 20:14
• In Ddcd you're integrating with respect to r, but you've passed in r as an argument. Which means r is getting replaced by a value and I get the message Raw object 178 cannot be used as an iterator. The same happens with th. If you replace these variables with different names from the arguments it fixes it: Ddcd[r_, th_, Ad_, Zd_, b2_, b4_] := Zd/(NIntegrate[ 2 Pi*rr^2*fd[rr, tth, Ad, b2, b4]*Sin[tth], {rr, 0, Infinity}, {tth, 0, Pi}, Method -> {Automatic, "SymbolicProcessing" -> 0}])* fd[r, th, Ad, b2, b4] – flinty Jul 6 '20 at 14:26
• Hi. Thank you for your answer. I do not have that mistake because my code is slightly different; the integral part is another function. I did not write it here because I thought it would be too long. Thank you again – Jhoan Perez Jul 6 '20 at 14:37