Solving a Laplace transmission problem with the Multiple Traces Formulation

Geometry

We consider \(N\) disjoint spheres in \(\mathbb{R}^3\), with the following notation:

• Position vectors of the sphere centers: \(\mathbf{p_j} \in \mathbb{R}^3\), with \(j \in \{1, ..., N\}\).

• Radii: \(\tilde{r}_j \in \mathbb{R}^+\), with \(j \in \{1, ..., N\}\).

• Interior of sphere \(j\): \(\Omega_j\), with \(j \in \{1, ..., N\}\).

• Boundary of sphere \(j\): \(\Gamma_j := \partial \Omega_j\).

• Exterior medium: \(\Omega_0\), defined as \(\mathbb{R}^3\) without the spheres and their boundaries.

Volume problem

We want to find \(u_j \in H^1_{loc}(\Omega_j)\), \(j \{0, ..., N \}\), such that they solve:

\[\text{div} (\sigma_j \nabla u_j) = 0 \quad \text{ in } \Omega_j, \text{ with } j \in \{0, ..., N \},\]

\[-\gamma_d^{0j} u_0 + \gamma_d^{j} u_j = \gamma_d^{0j} \phi_e \quad \text{ on } \Gamma_j, \text{ with } j \in \{1, ..., N \},\]

\[\sigma_0 \gamma_n^{0j} u_0 + \sigma_j \gamma_n^{j} u_j = - \sigma_0\gamma_n^{0j} \phi_e \quad \text{ on } \Gamma_j, \text{ with } j \in \{1, ..., N \},\]

\[\lim_{||\mathbf{x}||_2 \rightarrow \infty} u_0(\mathbf{x}) = \mathcal{O}\left(\frac{1}{||\mathbf{x}||_2}\right),\]

with

\[\gamma_d^{0j} u_0 := u_0|_{\Gamma_j},\]

\[\gamma_d^{j} u_j := u_j|_{\Gamma_j},\]

\[\gamma_n^{0j} u_0 := \nabla u_0|_{\Gamma_j} \cdot \widehat{n}_{0j},\]

\[\gamma_n^{j} u_j := \nabla u_j|_{\Gamma_j}\cdot \widehat{n}_{j},\]

where \(\widehat{n}_{j}\) is the exterior normal of \(\Omega_j\), with \(j\in \{ 1, ..., N\}\) and \(\widehat{n}_{0j}=-\widehat{n}_{j}\).

\(\sigma_j\), with \(j\in \{ 1, ..., N\}\), are positive parameters and

\[\phi_e(\mathbf{x}) = \frac{1}{4 \pi \sigma_0 ||\mathbf{x} - \mathbf{p}||_2},\]

with \(\mathbf{p} \in \Omega_0\).

Boundary formulation

Free space fundamental solution for the Laplace equation

The free space fundamental solution of the Laplace equation, which satisfies the radiation condition is:

\[g\left(\mathbf{r},\mathbf{r'}\right):= \frac{1}{4\pi ||\mathbf{r}-\mathbf{r}'||_2}, \quad \text{with } \mathbf{r}\not = \mathbf{r'}, \quad \text{and } -\Delta g \left(\mathbf{r},\mathbf{r'}\right)=\delta \left(\mathbf{r}-\mathbf{r'}\right).\]

Layer operators

Single and double layer operators defined for smooth densities:

\[DL_{0j} \left(\psi\right)\left(\mathbf{r}\right):=\int_{\Gamma_j} \psi\left(\mathbf{r}'\right) \nabla g\left(\mathbf{r},\mathbf{r'}\right) \cdot \widehat{\mathbf{n}}_{0j} \ dS',\]

\[SL_{0j} \left(\psi\right)\left(\mathbf{r}\right):=\int_{\Gamma_j} {\psi\left(\mathbf{r}'\right) g\left(\mathbf{r},\mathbf{r'}\right) dS'},\]

\[DL_j \left(\psi\right)\left(\mathbf{r}\right):=\int_{\Gamma_j} {\psi\left(\mathbf{r}'\right) \nabla g\left(\mathbf{r},\mathbf{r'}\right) \cdot \widehat{\mathbf{n}}_{j} \ dS'},\]

\[SL_j \left(\psi\right)\left(\mathbf{r}\right):=\int_{\Gamma_j} {\psi\left(\mathbf{r}' \right) g\left(\mathbf{r},\mathbf{r'}\right) dS'}\]

with the gradient being taken with respect to \(\mathbf{r}'\), \(\widehat{\mathbf{n}}_j\) being the exterior normal vector of \(\Omega_j\), and \(\widehat{\mathbf{n}}_j=-\widehat{\mathbf{n}}_{0j}\).

These operators are linear and continuous in the following Sobolev spaces:

\[DL_{0j}: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^1_{loc} \left(\mathbb{R}^3 \setminus \cup_{j=1}^{N}\Gamma_j\right),\]

\[SL_{0j}: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^1_{loc} \left(\mathbb{R}^3 \setminus \cup_{j=1}^{N}\Gamma_j\right) ,\]

\[DL_{j}: H^{\frac{1}{2}}(\Gamma_j)\rightarrow H^1_{loc} \left(\mathbb{R}^3 \setminus \cup_{j=1}^{N}\Gamma_j\right),\]

\[SL_{j}: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^1_{loc}\left(\mathbb{R}^3 \setminus \cup_{j=1}^{N}\Gamma_j\right).\]

Representation formula

The integral representation formulas for \(u_j \in H^1(\Omega_j)\), \(u_0 \in H^1_{loc}(\Omega_0)\) read

\[\begin{split}\begin{aligned} u_0 &= - \sum_{j=1}^N DL_{0j} \left(\gamma_d^{0j} u_0 \right) + \sum_{j=1}^N SL_{0j} \left(\gamma_n^{0j} u_0 \right) ,\\ u_j&= -DL_j \left(\gamma_d^{j} u_j \right) +SL_j \left(\gamma_n{j} u_j \right), \quad \forall j\in \{1, ..., N\}. \end{aligned}\end{split}\]

where \(u_j\) are zero-valued on the complement of \(\Omega_j\).

Boundary integral operators

We recall the definition of the boundary integral operators:

\[\begin{split}\begin{aligned} V_{i,j}^0 &:= \frac{1}{2} \left( \gamma_d^{i} SL_{0j} + \gamma_d^{0i} SL_{0j} \right) , & V_{j}&:= \frac{1}{2} \left( \gamma_d^{0j} SL_{j} + \gamma_d^{j} SL_{j} \right) ,\\ K_{i,j}^0&:= \frac{1}{2} \left(\gamma_d^{i} DL_{0j} + \gamma_d^{0i} DL_{0j} \right) , &K_{j}&:= \frac{1}{2} \left(\gamma_d^{0j} DL_{j} + \gamma_d^{j} DL_{j} \right),\\ K^{*0}_{i,j}&:= \frac{1}{2} \left( - \gamma_n ^{i} SL_{0j} + \gamma_n ^{0i} SL_{0j} \right), & K^{*}_{j} &:= \frac{1}{2} \left( -\gamma_n ^{0j} SL_{j} + \gamma_n ^{j} SL_{j} \right), \\ W_{i,j}^0 &:= -\frac{1}{2} \left( - \gamma_n^{i} DL_{0j} + \gamma_n^{0i} DL_{0j} \right) , & W_{j} &:=- \frac{1}{2} \left( -\gamma_n^{0j} DL_{j} + \gamma_n^{j} DL_{j} \right). \end{aligned}\end{split}\]

These operators are linear and continuous in the following Sobolev spaces:

\[\begin{split}\begin{aligned} V_{{i},j}^0 &: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_i), &V_{j}&: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_j),\\ W_{{i},j}^0&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_i), &W_{j}&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_j) ,\\ K_{{i},j}^0&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_i), &K_{j}&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_j) ,\\ K^{*0}_{{i},j}&: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_i), &K^*_{j}&: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_j). \end{aligned}\end{split}\]

Since the domains are smooth, the jump relations for the potentials across a closed boundary yield

\[\begin{split}\begin{aligned} V_{{i},j}^0 &= \gamma_d^{0{i}} SL_{0j}, & V_{j}&= \gamma_d^{j} SL_{j},\\ W_{{i},j}^0 &=- \gamma_n^{0{i}} DL_{0j}, & W_{j} &=- \gamma_n^{j} DL_{j}, \\ K_{{i},j}^0&= \gamma_d^{0{i}} DL_{0j}\text{ with } {i} \not=j, & K^{*0}_{{i},j} &= \gamma_n ^{0{i}} SL_{0j}\text{ with } {i} \not=j, \end{aligned}\end{split}\]

and

\[\begin{split}\begin{aligned} K_{j,j}^0(\psi) &= \frac{1}{2}\psi +\gamma_d^{0j} {DL_{0j}(\psi)} , &K_{j}(\psi) &= \frac{1}{2} \psi +\gamma_d^{j} {DL_{j}(\psi)} ,\\ K^{*0}_{j,j}(\psi) &= -\frac{1}{2} \psi + \gamma_n^{0j} {SL_{0j}(\psi)}, &K^*_{j}(\psi) &= -\frac{1}{2} \psi + \gamma_n^j {SL_{j}(\psi)}. \end{aligned}\end{split}\]

Multiple traces formulation

For \(j\in \{1,...,N\}\), we introduce the Cartesian product of Hilbert spaces \({\boldsymbol{H}_j} := H^{\frac{1}{2}}(\Gamma_j) \times H^{-\frac{1}{2}}(\Gamma_j)\), with norm \(\|\cdot \|_{\boldsymbol{H}_j}= \|\cdot \|_{H^{\frac{1}{2}}(\Gamma_j)} + \| \cdot \|_{H^{-\frac{1}{2}}(\Gamma_j)}\).

Let be \(\boldsymbol{\phi}\), \(\boldsymbol{\xi}\) \(\in {\boldsymbol{H}_j}\), such that \(\boldsymbol{\phi} = (\phi_D, \phi_N )\) and \(\boldsymbol{\xi} = (\xi_D, \xi_N )\). We introduce the cross-product over \(\Gamma_j\) (Section 2.2.1, HiptmairJerez2012) by \(\langle \boldsymbol{\phi}, \boldsymbol{\xi}\rangle_{\times,j} := \langle \phi_d, \xi_n \rangle_{j} + \langle \xi_d , \phi_n\rangle_{j},\) where for brevity we denote \(\langle \xi_d , \phi_n\rangle_{j}\) := \(\langle \xi_d, \phi_n\rangle_{H^{\frac{1}{2}}(\Gamma_j) \times H^{-\frac{1}{2}}(\Gamma_j)}\).

We define also the flip-sign operator \(\mathsf{X}_j:{\boldsymbol{H}_j}\rightarrow {\boldsymbol{H}_j}\), \(\gamma^{0j}:H^1_{loc}(\Delta,\Omega_0)\rightarrow {\boldsymbol{H}_j}\) and \(\gamma^{j}:H^1(\Delta,\Omega_j)\rightarrow {\boldsymbol{H}_j}\) as:

\[\begin{split}\mathsf{X}_j:= \begin{bmatrix} I & 0 \\[2mm] \displaystyle 0 & -\frac{\sigma_0}{\sigma_j} I \end{bmatrix}, \quad \gamma^{0j} := \begin{pmatrix} \gamma_d^{0j} \\[2mm] \gamma_n^{0j} \end{pmatrix} \quad \text{ and } \quad \gamma^{j} := \begin{pmatrix} \gamma_d^{j} \\[2mm] \gamma_n^{j} \end{pmatrix},\quad j \in \{1,..., N\},\end{split}\]

with \(I\) being the identity operator in the corresponding functional space, and for simplicity, we adopt the same notation for \(I\) in different spaces.

Then, we write Dirichlet and Neumamn boundary conditions succinctly as

\[\begin{split}\begin{aligned} -\mathsf{X}_j \gamma^{0j} u_0+ \gamma^j u_j &= \mathsf{X}_j \gamma^{0j} \phi_e , \\ \gamma^{0j} u_0 -\mathsf{X}_j^{-1}\gamma^j u_j &=- \gamma^{0j} \phi_e, \end{aligned}\end{split}\]

where both equations are equivalent.

Taking Dirichlet and Neumann traces of both representation formulas for \(u_0\) and \(u_j\), and rewriting in terms of BIOs, we obtain

\[\begin{split}\begin{aligned} \gamma_d^{0j} u_0 &= -\left(-\frac{1}{2}I\left(\gamma_d^{0j} u_0 \right)+\sum_{i=1}^n K_{j,i}^0 \left(\gamma_d^{0i} u_0 \right)\right) +\sum_{i=1}^n V_{j,i}^0 \left(\gamma_n^{0i} u_0 \right),\\ \gamma_n^{0j} u_0 &= \sum_{i=1}^n W_{j,i}^0 \left(\gamma_d^{0i} u_0 \right) +\left( \frac{1}{2}I\left(\gamma_n^{0j} u_0 \right)+\sum_{i=1}^n K_{j,i}^{*0} \left(\gamma_n^{0i} u_0 \right)\right),\\ \gamma_d^{j} u_j &= -\left(-\frac{1}{2}I\left(\gamma_d^{j} u_j \right) +K_j\left(\gamma_d^{j} u_j \right) \right) +V_j \left(\gamma_n^{j} u_j \right),\\ \gamma_n^{j} u_j &= W_j \left(\gamma_d^{j} u_j \right) +\left( \frac{1}{2}I\left(\gamma_n^{j} u_j \right) + K^*_j\left(\gamma_n^{j} u_j \right) \right) . \end{aligned}\end{split}\]

After some algebra, one can write

\[\gamma^{0j} u_0= 2\sum_{i=1}^N \mathsf{A}_{j,i}^0 \gamma^{0i}u_0, \quad \gamma^j u_j= 2\mathsf{A}_j \gamma^j u_j, \quad j \in \{1,...,N\},\]

with :math:`displaystyle mathsf{A}_{j,i}^0:= begin{bmatrix}

-K_{j,i}^0 & V_{j,i}^0 \ W_{j,i}^0 & K_{j,i}^{*0} \

end{bmatrix} text{ and }mathsf{A}_j:= begin{bmatrix}

-K_j & V_j \ W_j & K_j^* \

end{bmatrix}.`

By replacing \(\gamma^{0j} u_0\), \(\gamma^j u_j\) into the boundary conditions we obtain

\[\begin{split}\begin{aligned} 2\sum_{i=1}^n \mathsf{A}_{j,i}^0 \gamma^{0{i}} u_0 -\mathsf{X}_j^{-1}\gamma^j u_j &= - \gamma^{0j} \phi_e ,\\ -\mathsf{X}_j\gamma^{0j} u_0 +2\mathsf{A}_j \gamma^j u_j &= \mathsf{X}_j \gamma^{0j} \phi_e \quad \text{ on } \Gamma_j. \end{aligned}\end{split}\]

We define Cartesian product space of multiple traces \(\mathbb{H}:=\Pi_{j=1}^N \boldsymbol{H}_j\) and \(\mathbb{H}^{(2)} := \mathbb{H} \times \mathbb{H} = \Pi_{j=1}^N \boldsymbol{H}_j \times \Pi_{j=1}^N \boldsymbol{H}_j\), and the cross-product:

\[\displaystyle \langle \boldsymbol{\phi}, \boldsymbol{\xi } \rangle_{\times} = \sum_{j=1}^N \langle \boldsymbol{\phi}^{0j}, \boldsymbol{\xi}^{0j} \rangle_{\times,j} + \sum_{j=1}^N \langle \boldsymbol{\phi}^{j}, \boldsymbol{\xi}^{j} \rangle_{\times,j},\]

with \(\boldsymbol{\phi} = (\boldsymbol{\phi}^{01},..., \boldsymbol{\phi}^{0N}, \boldsymbol{\phi}^{1}, ..., \boldsymbol{\phi}^{N})\) and \(\boldsymbol{\xi} = (\boldsymbol{\xi}^{01},..., \boldsymbol{\xi}^{0N}, \boldsymbol{\xi}^{1}, ..., \boldsymbol{\xi}^{N})\).

Now, we can introduce the local Multiple Trace formulation (MTF) operator (Section 3.2.3, HiptmairJerez2012), \(\boldsymbol{\mathsf{M}}_N: \mathbb{H}^{(2)} \rightarrow \mathbb{H}^{(2)}\):

\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{M}}_N:= \begin{bmatrix} 2\boldsymbol{\mathsf{A}}_{0,N} & -\boldsymbol{\mathsf{X}}_N^{-1} \\ -\boldsymbol{\mathsf{X}}_{N} & 2\boldsymbol{\mathsf{A}}_{1,N} \end{bmatrix}, \text{ with } \boldsymbol{\mathsf{A}}_{0,N}:=\begin{bmatrix} \mathsf{A}_{1,1}^0 & \mathsf{A}_{1,2}^0 & ... & \mathsf{A}_{1,N}^0\\ \mathsf{A}_{2,1}^0 & \mathsf{A}_{2,2}^0 & ... & \mathsf{A}_{2,N}^0\\ \vdots & &\ddots & \vdots\\ \mathsf{A}_{N,1}^0 & \mathsf{A}_{N,2}^0 & ... & \mathsf{A}_{N,N}^0 \end{bmatrix}, \end{aligned}\end{split}\]

\[\begin{split}\begin{aligned} &\boldsymbol{\mathsf{A}}_{1,N}:=\begin{bmatrix} \mathsf{A}_{1} & 0 & ... & 0\\ 0 & \mathsf{A}_{2} & ... & 0\\ \vdots & &\ddots & \vdots\\ 0 & ... & 0& \mathsf{A}_N\\ \end{bmatrix} \text{ and } \boldsymbol{\mathsf{X}}_N:=\begin{bmatrix} \mathsf{X}_{1} & 0 & ... & 0\\ 0 & \mathsf{X}_{2} & ... & 0\\ \vdots & &\ddots & \vdots\\ 0 & ... & 0& \mathsf{X}_N \end{bmatrix}. \end{aligned}\end{split}\]

With the MTF operator, the transmission problem can be written as:

\[\begin{split}\begin{aligned} \boldsymbol{\mathsf{M}}_N\begin{pmatrix} \gamma^0_u \\ \gamma_u \end{pmatrix} = \begin{pmatrix} -\gamma^0_{\phi_e} \\ \boldsymbol{\mathsf{X}}_N \gamma^0_{\phi_e} \end{pmatrix}, \end{aligned}\end{split}\]

where we use the notation:

\[\begin{split}\begin{aligned} \gamma^0_u &:=\begin{pmatrix} \gamma^{01} u_0, \gamma^{02} u_0, ..., \gamma^{0N} u_0 \end{pmatrix}^t , &\gamma_u &:=\begin{pmatrix} \gamma^1 u_1 , \gamma^2 u_2, ..., \gamma^N u_N \end{pmatrix}^t,\\ \gamma^0_{\phi_e} &: =\begin{pmatrix} \gamma^{01} \phi_e, \gamma^{02} \phi_e, ..., \gamma^{0N} \phi_e \end{pmatrix}^t, \end{aligned}\end{split}\]

with superscript \(t\) denoting the transposition.

This is the system that it is solved.

Preliminaries of the discretization

Spherical coordinate system

A vector is written as \(\mathbf{r}=\left(r,\varphi,\theta\right)^t\), with \(r \in [0,\infty)\), \(\varphi \in [0,2\pi)\) and \(\theta \in [0,\pi]\), which in Cartesian coordinates is equivalent to \(\mathbf{r}=r\left(\sin \theta \cos \varphi,\sin \theta \sin \varphi,\cos \theta\right)^t\). The unitary vectors of the spherical coordinate system can be written in Cartesian coordinates as

\[\begin{split}\begin{aligned} &\widehat{\mathbf{e}}_r= (\sin \theta \cos \varphi , \sin \theta \sin \varphi , \cos \theta )^t, \\ &\widehat{\mathbf{e}}_\theta=(\cos \theta \cos \varphi , \cos \theta \sin \varphi , -\sin \theta )^t, \\ &\widehat{\mathbf{e}}_{\varphi}=(-\sin \varphi , \cos \varphi , 0 )^t. \end{aligned}\end{split}\]

Also, the gradient operator in spherical coordinates:

\[\begin{aligned} \nabla f&=\frac{\partial f}{\partial r} \widehat{\mathbf{e}}_r+\frac{1}{r} \frac{\partial f}{\partial \theta} \widehat{\mathbf{e}}_\theta+\frac{1}{r \sin \theta} \frac{\partial f}{\partial \varphi} \widehat{\mathbf{e}}_{\varphi} \end{aligned}\]

Associated Legendre functions

\(P_l^m\) are the associated Legendre functions of degree \(l\) and order \(m\) defined as:

\[P_{l}^m\left(x\right) := (-1)^m \left( 1- x^2\right)^{\frac{m}{2}} \frac{d^m}{dx^m}P_l(x), \quad \text{with} \quad P_{l}\left(x\right) := \frac{1}{2^ll!}\frac{d^l}{dx^l}(x^2-1)^l.\]

Here, the term \((-1)^m\) is the Condon-Shortley phase factor.

Real spherical harmonics

Real spherical harmonics of degree \(l\) and order \(m\) are defined using spherical coordinates:

\[\begin{split}\begin{aligned} Y_{l,m}\left(\theta,\varphi\right) &:= \sqrt{ (2-\delta_{m,0}) \frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\theta\right) \cos m \varphi , \text{ and}\\ Y_{l,-m}\left(\theta,\varphi\right) &:= \sqrt{ (2-\delta_{m,0})\frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\theta\right) \sin m \varphi , \end{aligned}\end{split}\]

with \(l\in \mathbb{N}_0\), \(m\in \mathbb{Z}\) such that \(0\leq m\leq l\). If \(m=0\), then \(\delta_{m,0}=1\), and it is zero otherwise.

Orthonormality of spherical harmonics

Spherical harmonics are dense in \(C(\mathbb{S}^2)\), with \(\mathbb{S}^2\) the surface of the unit sphere, and form a complete orthonormal system in \(L^2(\mathbb{S}^2)\) with respect to the internal product defined by:

\[\left( \psi , \xi \right)_{L^2(\mathbb{S}^2)} = \int_{0}^{2\pi}\int_{0}^{\pi} \psi\left(\theta,\varphi\right) \xi\left(\theta,\varphi\right) \sin\left(\theta\right) d\theta d\varphi.\]

They also are orthogonal in \(H^1(\mathbb{S}^2)\).

Discretization

Let be \(j\in \{1,...,N\}\). We define the reference system \(j\) as the one centered at \(\mathbf{p_j}\) with the same orientation that the reference system centered in the origin. Furthermore, we denote by \(Y_{l,m,j}\) the spherical harmonic \(Y_{l,m}\) centered in the origin of the reference system \(j\). Thus, if \(\left( r_j, \varphi_j, \theta_j \right)\) are the vector spherical coordinates of \(\mathbf{r_j}\) in the reference system \(j\), we have that \(Y_{l,m,j}\left(\mathbf{r}_j\right)=Y_{l,m}\left(\theta_j, \varphi_j\right)\).

For \(L \in \mathbb{N}_0\) and \(j\in \{1,...,N\}\), we define subspaces

\[\mathcal{Y}_L\left(\Gamma_j \right):= \text{span}\left\lbrace Y_{l,m,j}: l \in \mathbb{N}_0, m \in \mathbb{Z}, l \leq L, |m|\leq l \right\rbrace,\]

equipped with the \(L^2(\Gamma_j)\)-norm. Notice that the dimension of each subspace is \((L+1)^2\), and that the sequence of subspaces \(\lbrace \mathcal{Y}_L \left(\Gamma_j \right) \rbrace_{L \in \mathbb{N}_0}\) is dense in \(H^{\frac{1}{2}}(\Gamma_j)\) and in \(H^{-\frac{1}{2}}(\Gamma_j)\). This result justifies the discretization of all boundary Dirichlet and Neumann unknowns with spherical harmonics.

For \(j \in \lbrace 1, ..., N \rbrace\), we write \(u^L_{d,0j}\), \(u^L_{n,0j}\), \(u^L_{d,j}\) and \(u^L_{n,j}\) in \(\mathcal{Y}_L(\Gamma_j)\) for the approximations of \(\gamma_d^{0j} u_0\), \(\gamma_n^{0j} u_0\), \(\gamma_d^{j} u_j\) and \(\gamma_n^{j} u_j\), respectively. They can be written as the following series expansions:

\[\begin{split}\begin{aligned} u^L_{d,0j}&=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{d,0j} Y_{l,m,j}, & u^L_{n,0j} &=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{n,0j} Y_{l,m,j},\\ u^L_{d,j}&=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{d,j} Y_{l,m,j},& u^L_{d,j} &=\sum_{l=0}^{L} \sum_{m=-l}^l u^{l,m}_{n,j} Y_{l,m,j}, \end{aligned}\end{split}\]

with \(u^{l,m}_{d,0j}\), \(u^{l,m}_{n,0j}\), \(u^{l,m}_{d,j}\), and \(u^{l,m}_{n,j}\) being constants.

Solving

In the module biosspheres.formulations.mtf.solvertemplates there are two python functions to solve the presented problem. One uses a direct method to solve, and the other one an indirect one.

Import of the needed libraries

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm, colors
import biosspheres.formulations.mtf.solvertemplates as solver
import biosspheres.quadratures.sphere as quadratures
import biosspheres.miscella.spherearrangements as pos
import biosspheres.laplace.drawing as draw

/opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

Definition of needed variables

We will need:

• The number of spheres,

• The maximum degree of spherical harmonics used to discretize each trace.

• The degree parameter for the numerical quadrature.

• An array with the radii of the spheres.

• A list with the positions of the spheres.

• An array with the values of \(\sigma_j\).

• The vector \(\mathbf{p}\) of \(\phi_e\).

Definition of the number of spheres, it must be an integer:

n = 8

Definition of the maximum degree of spherical harmonics to be used, it must be an integer:

big_l = 15

Definition of the degree parameter for the numerical quadrature, it must be an integer:

big_l_c = 55

Definition of the array with the radii of the spheres, it must be a numpy array of one dimension. In this case we are going to define all radius equal to a parameter named r.

r = 0.875
radii = np.ones(n) * r

Definition of the list with the positions of the spheres, for this we use a existing routine:

d = 1.15
center_positions = pos.cube_vertex_positions(int(n**(1/3)), r, d)

We will plot the spheres to help visualizing them

quantity_theta_points, quantity_phi_points, weights, pre_vector = quadratures.gauss_legendre_trapezoidal_2d(50)
vector = pre_vector * r
fig = plt.figure()
ax_1 = fig.add_subplot(111, projection='3d')
for s in np.arange(0,n):
    ax_1.plot_surface(
        vector[0, :, :]+center_positions[s][0],
        vector[1, :, :]+center_positions[s][1],
        vector[2, :, :]+center_positions[s][2],
        rstride=1, cstride=1)
ax_1.set_xlabel('$x$')
ax_1.set_ylabel('$y$')
ax_1.set_zlabel('$z$')
ax_1.set_aspect('equal')
plt.show()

Definition of the numpy array containing \(\sigma_j\). In this example, we will make only \(\sigma_0\) different:

sigma_e = 1.75
sigma_i = 0.75
sigmas = np.ones(n+1) * sigma_i
sigmas[0] = sigma_e

Finally, the definition of the point source vector:

p0 = np.ones(3)*-5.
p0

array([-5., -5., -5.])

mtf_laplace_n_spheres_point_source_direct_solver

Now, we solve using a direct solver:

traces = solver.mtf_laplace_n_spheres_point_source_direct_solver(n, big_l, big_l_c, radii, center_positions, sigmas, p0)

The solver used is numpy.linalg.solve from the numpy library.

The following is a plot of the exterior dirichlet traces:

quantity_theta_points, quantity_phi_points, \
    weights, pre_vector, spherical_harmonics = \
    quadratures.gauss_legendre_trapezoidal_real_sh_mapping_2d(
        big_l, big_l_c)
vector = pre_vector * r
num = (big_l + 1)**2
surface_d0 = np.zeros((n, quantity_theta_points, quantity_phi_points))
for s in np.arange(0, n):
    surface_d0[s, :, :] = np.sum(
        traces[num*s:num*(s+1), np.newaxis, np.newaxis]
        * spherical_harmonics, axis=0)
surface_d0_max = np.max(np.abs(surface_d0))
surface_d0 = (surface_d0 + surface_d0_max) / (2. * surface_d0_max)

fig = plt.figure()
ax_1 = fig.add_subplot(111, projection='3d')
for s in np.arange(0,n):
    ax_1.plot_surface(
        vector[0, :, :]+center_positions[s][0],
        vector[1, :, :]+center_positions[s][1],
        vector[2, :, :]+center_positions[s][2],
        rstride=1, cstride=1,
        facecolors=cm.coolwarm(surface_d0[s, :, :]))
ax_1.set_xlabel('$x$')
ax_1.set_ylabel('$y$')
ax_1.set_zlabel('$z$')
ax_1.set_aspect('equal')
fig.colorbar(cm.ScalarMappable(
    norm=colors.CenteredNorm(halfrange=surface_d0_max), cmap=cm.coolwarm),
    ax=ax_1)
plt.show()

mtf_laplace_n_spheres_point_source_indirect_solver

Now, we solve using the solver provided by scipy.sparse.linalg.gmres. This time, we have to specify a tolerance:

tolerance = 10**(-10)

traces_2 = solver.mtf_laplace_n_spheres_point_source_indirect_solver(n, big_l, big_l_c, radii, center_positions, sigmas, p0, tolerance)

Let us see the difference between the solutions:

np.linalg.norm(traces-traces_2)

np.float64(6.452772922078449e-13)

In the notebook mtf_times.ipynb we measure the difference of time between routines (the iterative one for this example is faster).

Plots using the representation formula

Cut in a plane

Center of the plot:

center = np.asarray([1., 1., 1.])*(r+d/2)

Size of the window to be plotted:

horizontal = 10.
vertical = 10.

For the resolution of the image:

resolution = 50 # points per side
inter_horizontal = resolution
inter_vertical = resolution

If we do not want to plot the external function, we can use the following:

def zero(x: np.ndarray) -> float:
    return 0.

cut = 1
x1, y1, data = draw.draw_cut_representation_formula_n_sphere(cut, center,
                                                             horizontal,
                                                             vertical,
                                                             inter_horizontal,
                                                             inter_vertical,
                                                             traces,
                                                             radii,
                                                             center_positions,
                                                             big_l,
                                                             zero)
plt.figure()
plt.imshow(data, origin='lower',
           extent=[-horizontal / 2 + center[1], horizontal / 2 + center[1],
                   -vertical / 2 + center[2], vertical / 2 + center[2]],
           norm=colors.CenteredNorm()
           )
plt.xlabel('$y$')
plt.ylabel('$z$')
plt.title('$x = $'+str(center[0]))
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x77dc324539b0>

cut = 2
x1, y1, data = draw.draw_cut_representation_formula_n_sphere(cut, center,
                                                             horizontal,
                                                             vertical,
                                                             inter_horizontal,
                                                             inter_vertical,
                                                             traces,
                                                             radii,
                                                             center_positions,
                                                             big_l,
                                                             zero)
plt.figure()
plt.imshow(data, origin='lower',
           extent=[-horizontal / 2 + center[1], horizontal / 2 + center[1],
                   -vertical / 2 + center[2], vertical / 2 + center[2]],
           norm=colors.CenteredNorm()
           )
plt.xlabel('$x$')
plt.ylabel('$z$')
plt.title('$y = $'+str(center[1]))
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x77dc324e9a00>

cut = 3
x1, y1, data = draw.draw_cut_representation_formula_n_sphere(cut, center,
                                                             horizontal,
                                                             vertical,
                                                             inter_horizontal,
                                                             inter_vertical,
                                                             traces,
                                                             radii,
                                                             center_positions,
                                                             big_l,
                                                             zero)
plt.figure()
plt.imshow(data, origin='lower',
           extent=[-horizontal / 2 + center[1], horizontal / 2 + center[1],
                   -vertical / 2 + center[2], vertical / 2 + center[2]],
           norm=colors.CenteredNorm()
           )
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.title('$z = $'+str(center[2]))
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x77dc32450dd0>