Finite difference methods for approximating Heaviside functions (Q1017609)

The paper presents two algorithms to evaluate the integral \(\mathcal{I}:=\int_\Omega f(\vec{x})d\vec{x}\), where \(\vec{x}\in \mathbb{R}^n\), \(\Omega=\{\vec{x}:\,u(\vec{x})>0\}\) and \(\partial\Omega=\{\vec{x}:\,u(\vec{x})=0\}\) is a compact manifold of codimension one. The method to be employed is motivated by the expression \(\mathcal{I}=\int_{\mathbb{R}^n} H(u(\vec{x}))f(\vec{x})d\vec{x}\), where \(H\) is the Heaviside function. The approach consists of approximating \(H\) by finite differencing its first few primitives, a technique already used by the author to approximate delta functions. A brief presentation of the algorithms is given below. Let \[ I(z)=\int_0^z H(\zeta)d\zeta\quad\text{and}\quad J(z)=\int_0^z I(\zeta)d\zeta. \] The following relationships are derived: \[ I(u)=\langle\nabla J(u),\nabla u \rangle /|\nabla u|^2,\quad H(u)=\langle\nabla I(u),\nabla u \rangle /|\nabla u|^2, \] where \(\langle \cdot,\cdot \rangle\) stands for the inner product. By discretizing \(H(u)\) the one-step algorithm \(FDMH_1\) is obtained which converges at a rate of \(\mathcal{O}(h^2)\) when \(u\) is smooth enough. By discretizing both relationships the two-step algorithm \(FDMH_2\) is derived which can converge at a rate of \(\mathcal{O}(h^3)\). These results are validated by means of some numerical examples.