I am now trying to compute numerically the following integral. $$ \begin{split} L_1^s(\hat{\phi}_s)(r,\zeta,\theta_\zeta) &=\frac{1}{\sqrt{2}\pi} \int\limits_{0}^{2\pi} d\varphi \int\limits_0^{\pi} d\bar{\theta}_\zeta \int\limits_0^{\infty} d\bar{\zeta} \bar{\zeta}^2\sin\bar{\theta}_\zeta\cos\varphi F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})^{-\frac{1}{2}} \hat{\phi}_s(r,\bar{\zeta},\bar{\theta}_\zeta) \exp \left(-\bar{\zeta}^2 +\frac{F_2(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})}{F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})}\right)\\ \\ F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta}) &=\bar{\zeta}^2+\zeta^2-2\bar{\zeta}\zeta\cos\bar{\theta}_\zeta\cos\theta_\zeta-2\bar{\zeta}\zeta\sin\bar{\theta}_\zeta\sin\theta_\zeta\cos \varphi\\ F_2(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta}) &=\bar{\zeta}^2\zeta^2-\bar{\zeta}^2\zeta^2 (\cos\bar{\theta}_\zeta\cos\theta_\zeta+\sin\bar{\theta}_\zeta\sin\theta_\zeta\cos\varphi)^2 \end{split} $$
However, this integral has singularity because of $F_1^{-1/2}$.
When calculating this integral directly using the Gauss quadrature method, the results only approximately match those obtained through computation in Mathematica by approximately two decimal places. (I do not expect Mathematica to provide an exact solution, but I am currently using it as a benchmark.I am using GlobalAdaptive for mathematica's Nintegral option. However, when I use DoubleExponential, the value changes from the third significant digit, so I don't know which one to trust as a benchmark.)
I want to remove that singularity with a variable transformation, but I can't come up with a good way to do it. PS: When considering the problem of the singularity of this integral, it can be simplified by removing the unnecessary parts as follows. $$ \begin{split} L_1&=\int\limits_{0}^{2\pi} d\varphi \int\limits_0^{\pi} d\bar{\theta}_\zeta \int\limits_0^{\infty} d\bar{\zeta } \bar{\zeta}^2\sin\bar{\theta}_\zeta\cos\varphi F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})^{-\frac{ 1}{2}} \\ F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta}) &=\bar{\zeta}^2+\zeta^2-2\bar{\zeta }\zeta\cos\bar{\theta}_\zeta\cos\theta_\zeta-2\bar{\zeta}\zeta\sin\bar{\theta}_\zeta\sin\theta_\zeta\cos \varphi\\ \end{split} $$ I am not sure if this is necessary, but I will provide some background on the $L_1$.
This integral is part of the collision term in Boltzmann's equation.
$$ \zeta_i\frac{\partial\phi}{\partial x_i} =\frac{1}{k}[L_1(\phi)-L_2(\phi)-\nu(\zeta_i)\phi] $$ $L_1$ is expressed as follows ($\times$ is the cross product): $$ L_1(\phi)(\boldsymbol{x},\boldsymbol{\zeta})=\frac{1}{\sqrt{2}\pi}\int\limits_{\text{all }\bar{\boldsymbol{\zeta}}}\frac{1}{|\bar{\boldsymbol{\zeta}}-\boldsymbol{\zeta}|}\exp\left(-\bar{\boldsymbol{\zeta}}^2+\frac{|\bar{\boldsymbol{\zeta}}\times\boldsymbol{\zeta}|^2}{|\bar{\boldsymbol{\zeta}}-\boldsymbol{\zeta}|^2}\right)\phi(\boldsymbol{x},\bar{\boldsymbol{\zeta}})d\boldsymbol{\bar{\zeta}} $$ For $\phi(\boldsymbol{x},\boldsymbol{\zeta})$, transforming the space $\boldsymbol{x},\boldsymbol{\zeta}$to polar coordinates.
$$ L_1(\phi)(r,\theta,\psi,\zeta,\theta_\zeta,\psi_\zeta) = \frac{1}{\sqrt{2}\pi}\int\limits_0^{2\pi} d\bar{\psi}_\zeta \int\limits_0^{\pi} d\bar{\theta}_\zeta \int\limits_0^{\infty} d\bar{\zeta} \bar{\zeta}^2\sin\bar{\theta}_\zeta F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})^{-\frac{1}{2}} \phi(r,\psi,\theta,\bar{\zeta},\bar{\theta}_\zeta,\bar{\psi}_\zeta) \exp \left(-\bar{\zeta}^2+\frac{F_2(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})}{F_1(\bar{\boldsymbol{\zeta}},\boldsymbol{\zeta})}\right),\\ $$ and introducing a similarity solution $\phi=\sin\theta\zeta\sin\theta_\zeta\sin\psi_\zeta\phi(r,\zeta,\theta_\zeta)$ yields the following equation. $$ \begin{split} L_1(\sin{\theta}{\zeta}_\psi\phi_s)({r},{\zeta},{\theta}_\zeta) & = \frac{\sin{\theta}}{\sqrt{2}\pi} \int\limits_0^{2\pi} d\bar{\psi}_\zeta \int\limits_0^{\pi} d\bar{\theta}_\zeta \int\limits_0^{\infty} d\bar{\zeta} \bar{\zeta}^3\sin^2\bar{\theta}_\zeta\sin\bar{\psi}_\zeta F_1(\bar{\boldsymbol{\zeta}},{\boldsymbol{\zeta}})^{-\frac{1}{2}} \phi_s({r},\bar{\zeta},\bar{\theta}_\zeta) \exp\left (-\bar{\zeta}^2+\frac{F_2(\bar{\boldsymbol{\zeta}},{\boldsymbol{\zeta}})}{F_1(\bar{\boldsymbol{\zeta}},{\boldsymbol{\zeta}})}\right),\\ \\ F_1(\bar{\boldsymbol{\zeta}},{\boldsymbol{\zeta}}) & =\bar{\zeta}^2+{\zeta}^2-2\bar{\zeta}{\zeta}\cos\bar{\theta}_\zeta\cos{\theta}_\zeta-2\bar{\zeta}{\zeta}\sin\bar{\theta}_\zeta\sin{\theta}_\zeta\cos (\bar{\psi}_\zeta-{\psi}_\zeta)\\ F_2(\bar{\boldsymbol{\zeta}},{\boldsymbol{\zeta}}) & =\bar{\zeta}^2{\zeta}^2-\bar{\zeta}^2{\zeta}^2 (\cos\bar{\theta}_\zeta\cos{\theta}_\zeta+\sin\bar{\theta}_\zeta\sin{\theta}_\zeta\cos(\bar{\psi}_\zeta-{\psi}_\zeta))^2 \end{split} $$ where $\varphi=\bar{\psi}_\zeta-\psi_\zeta$ and using the $2\pi$ periodicity with respect to $\psi_\zeta$ of the integrand.Let $\hat{\phi}_s=\zeta\sin\theta_\zeta\phi_s$, to obtain the first equation.I want to emphasize here that $F_1$ is originally $|\bar{\boldsymbol{\zeta}}-\boldsymbol{\zeta}|^2$, and $F_2$ is originally $|\bar{\boldsymbol{\zeta}}\times\boldsymbol{\zeta}|$.
I am not used to asking questions yet, so I apologize if I am rude. Please point out any missing information. I will add it as soon as possible.
