Implementation of Point Source

Junxiong Jia, 2020.12.02

Introduction

For solving inverse problems, we often meet point source when tacking adjoint equations, e.g., the following equation
$-\Delta u - e^{2q}u = -\frac{1}{\sigma^2}\sum_{j=1}^{K}(w(x)-d_j)\delta(x-x_j) \quad\text{ in }\Omega,$
where $\{x_j\}_{j=1}^{K}$ are some measurement points, $w(x)$ is usually some function obtained from forward equation. In the following, we explain that how to implement the point source term on the right-hand side of the above equation.

Measurement Matrix

Usually, for a function in $L^2(\Omega)$ ( $\Omega$ is some open bounded domain), it can be decomposed by some fixed orthonormal basis $\{ \varphi_{i} \}_{i=1}^{\infty}$ . In the finite element method, we project functions on a finite number of basis as follow
$u(x) \approx \sum_{i=1}^{N}u_i \varphi_i(x).$
And once we specify the coefficients $\{u_i\}_{i=1}^{N}$ , the function $u(x)$ will be specified. We introduce the following measurement matrix
$S = \left[\begin{matrix} \varphi_1(x_1) & \cdots & \varphi_N(x_1) \\ \vdots & & \vdots \\ \varphi_1(x_K) & \cdots & \varphi_N(x_K) \end{matrix}\right].$
The matrix has the following property
$\left[\begin{matrix} u(x_1) \\ \vdots \\ u(x_K) \end{matrix}\right] = S \left[\begin{matrix} u_1 \\ \vdots \\ u_N \end{matrix}\right].$

Implement of the Point Source

Now we can give specific suggestions for implementing the point source. For discretization, we usually introduce the following quantities
$f_k := \int_{\Omega}\sum_{j=1}^{K}(w(x)-d_j)\delta(x-x_j)\varphi_{k}(x)dx$
with $k = 1,2,\cdots, N$ . Actually, we can reduce $f_k$ by some simple calculations as follow:
$f_k = \sum_{j=1}^{K}(w(x_j)-d_j)\varphi_k(x_j).$
With these calculations, we find that
$\left[\begin{matrix} f_1 \\ \vdots \\ f_N \end{matrix}\right] = \left[\begin{matrix} \varphi_1(x_1) & \cdots & \varphi_1(x_K) \\ \vdots & & \vdots \\ \varphi_{N}(x_1) & \cdots & \varphi_{N}(x_K) \end{matrix}\right] \left[\begin{matrix} w(x_1)-d_1 \\ \vdots \\ w(x_K)-d_K \end{matrix}\right] = S^T \left[\begin{matrix} w(x_1)-d_1 \\ \vdots \\ w(x_K)-d_K \end{matrix}\right].$
Denote $d = (d_1,\cdots,d_K)^T$ , $f=(f_1,\cdots,f_N)^T$ and $w=(w_1,\cdots,w_N)^T$ (coefficients of the function $w$ ), we can reduce the above formula as follows:
$f = S^T(Sw - d),$
which is convenient for implementations by some numerical algorithm, e.g., finite element method.