Let $f=\frac{1}{|x|},x\in\mathbb{R^3}$ and $\Omega=[-b,b]^3$. How to construct a quadrature scheme to solve 
$$
\int\_\Omega f\phi\psi dx\quad ?
$$ 
where $\phi\psi$ is smooth function.


I know there exists a transformation called Duffy transformation which can eliminate the singularity at $x=0$ by using the jacobian. Indeed, let $$\begin{aligned}T&=\{(x\_1,x\_2,x\_3);0\leq x\_1\leq 1,0\leq x\_2\leq x\_1,0\leq x\_3\leq x\_1\}\\ \hat{T}&=\{(\hat{x}\_1,\hat{x}\_2,\hat{x}\_3);0\leq \hat{x}\_1\leq 1,0\leq\hat{x\_2}\leq1,0\leq\hat{x\_3}\leq 1\}\end{aligned}$$


Define $F:\hat{T}\rightarrow T$ as 
$$
F(\hat{x}\_1,\hat{x}\_2,\hat{x}\_3)=(\hat{x}\_1,\hat{x}\_1\hat{x}\_2,\hat{x}\_1\hat{x}\_3)
$$
Then 
$$
\begin{aligned}
\int\_{T}\frac{1}{|x|}\phi\psi dx&=\int\_{\hat{T}}\frac{1}{\sqrt{x\_1^2+x\_2^2+x\_3^2}}\phi\cdot\psi(x)dx
=\int\_{\hat{T}}\frac{\phi\cdot\psi(\hat{x\_1},\hat{x}\_1\hat{x}\_2,\hat{x}\_1\hat{x}\_3)}{\sqrt{\hat{x}\_1^2+\hat{x}\_1^2\hat{x}\_2^2+\hat{x}\_1^2\hat{x}\_3^2}}\hat{x}\_1^2d\hat{x}\\
&=\int\_\hat{T}\frac{\hat{x}\_1\phi\cdot\psi(\hat{x\_1},\hat{x}\_1\hat{x}\_2,\hat{x}\_1\hat{x}\_3)}{\sqrt{1+\hat{x}\_2^2+\hat{x}\_3^2}}d\hat{x}
\end{aligned}
$$ 
which can be solved by usual quadrature scheme.


On the other hand, someone can also avoid the singularity by omitting a $\varepsilon$ cube which contains the origin.


However, these are not enough for me when I want to study the convergence rate of numerical integration on $\Omega$. I hope the error can be bounded by $O(h)\Vert \phi\Vert\_{1,2,\Omega}\Vert \psi\Vert\_{1,2,\Omega}$


Is there any other way to integrate a $L^2$ function or a function with singularity at vertex directly with the error bounded by $L^2$ norm?


Thank you very much.



How to integrate the $L^2$ function $1/|x|$ numerically
Following Michael Renardy's suggestion, write $g=\phi\psi$, and write
$$
g(x)=g(0)+x\cdot G(x) ,
$$
where $G$ is a smooth function. Then the integral becomes
$$
\int\_\Omega fg = g(0)\int\_\Omega f + \sum\_i\int\_\Omega \frac{x\_i}{|x|}G\_i(x)d^3x .
$$
The integrands in the latter integrals are bounded, and the first integral can be written as
$$
\int\_\Omega f = \int\_{B\_r}f + \int\_{\Omega\setminus B\_r}f ,
$$
where the first can be computed analytically, and the second has a smooth integrand. Here $B\_r$ is some ball centred at the origin.


