Overview of python functions from helmholtz.selfinteractions

Overview of python functions from biosspheres.helmholtz.selfinteractions.

import biosspheres.helmholtz.selfinteractions as selfinteractions

import numpy as np
import scipy.sparse as sparse
import matplotlib.pyplot as plt
from matplotlib import cm, colors

Preliminaries

Geometry

\(N\) disjoint spheres in \(\mathbb{R}^3\):

• 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:=\{ \mathbf{x} \in \mathbb{R}^3:||\mathbf{x}-\mathbf{p_j}||_2<\tilde{r}_j \}\).

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

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

Free space fundamental solution for the Helmholtz equation

The free space fundamental solution of the Helmholtz equation with \(k>0\), which satisfies the radiation condition is:

\[g_j\left(\mathbf{r},\mathbf{r'}\right):= \frac{e^{ik_j||\mathbf{r}-\mathbf{r}'||}}{4\pi ||\mathbf{r}-\mathbf{r}'||}, \quad \text{with } \mathbf{r}\not = \mathbf{r'}, \quad \text{and } -\left( \Delta + k_j^2 \right) g_j \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_0\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_0\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_j\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_j\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).\]

Boundary integral operators

For \(u\in C^\infty(\overline\Omega)\), Dirichlet and Neumann traces operators are defined as

\[\gamma_d u := u|_{\Gamma},\qquad \gamma_n u := \nabla u|_{\Gamma} \cdot \widehat{n},\]

where \(\widehat{n}\) is the exterior unit normal.

By density arguments, the definition of Dirichlet and Neumann traces operator can be extended to \(u_j \in H^1_{loc}(\Omega_j)\), with \(j \{0, ..., N \}\). We extend the notation as follows

\[\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}\).

Now, we recall the definition of the boundary integral operators:

\[\begin{split}\begin{aligned} V_{s,j}^0 &:= \frac{1}{2} \left( \gamma_d^{s} SL_{0j} + \gamma_d^{0s} SL_{0j} \right) , & V_{j}&:= \frac{1}{2} \left( \gamma_d^{0j} SL_{j} + \gamma_d^{j} SL_{j} \right) ,\\ K_{s,j}^0&:= \frac{1}{2} \left(\gamma_d^{s} DL_{0j} + \gamma_d^{0s} DL_{0j} \right) , &K_{j}&:= \frac{1}{2} \left(\gamma_d^{0j} DL_{j} + \gamma_d^{j} DL_{j} \right),\\ K^{*0}_{s,j}&:= \frac{1}{2} \left( - \gamma_n ^{s} SL_{0j} + \gamma_n ^{0s} SL_{0j} \right), & K^{*}_{j} &:= \frac{1}{2} \left( -\gamma_n ^{0j} SL_{j} + \gamma_n ^{j} SL_{j} \right), \\ W_{s,j}^0 &:= -\frac{1}{2} \left( - \gamma_n^{s} DL_{0j} + \gamma_n^{0s} 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_{s,j}^0 &: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_s), &V_{j}&: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_j),\\ W_{s,j}^0&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_s), &W_{j}&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_j) ,\\ K_{s,j}^0&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_s), &K_{j}&: H^{\frac{1}{2}}(\Gamma_j) \rightarrow H^{\frac{1}{2}}(\Gamma_j) ,\\ K^{*0}_{s,j}&: H^{-\frac{1}{2}}(\Gamma_j) \rightarrow H^{-\frac{1}{2}}(\Gamma_s), &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_{s,j}^0 &= \gamma_d^{0{s}} SL_{0j}, & V_{j}&= \gamma_d^{j} SL_{j},\\ W_{s,j}^0 &=- \gamma_n^{0s} DL_{0j}, & W_{j} &=- \gamma_n^{j} DL_{j}, \\ K_{s,j}^0&= \gamma_d^{0s} DL_{0j}, s \not=j, & K^{*0}_{s,j} &= \gamma_n ^{0s} SL_{0j}, s \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}\]

Calderón operator

\[\begin{split}A_{s,j}^0 = \begin{bmatrix} - K^0_{s,j} & V^0_{s,j} \\ W^0_{s,j} & K^{*0}_{s,j}\end{bmatrix}, \quad A_{j,j}= \begin{bmatrix} - K_{j,j} & V_{j,j} \\ W_{j,j} & K^{*}_{j,j}\end{bmatrix}.\end{split}\]

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}\]

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 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.

Complex spherical harmonics

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

\[Y_{l,m}\left(\theta,\varphi\right) := \sqrt{ \frac{\left(2l+1\right)\left(l-m\right)!}{4 \pi \left(l+m\right)!}} P_l^{m} \left(\cos\left(\theta\right)\right) e^{i m \varphi},\]

\[Y_{l,-m}\left(\theta,\varphi\right) := (-1)^m\overline{Y}_{l,m}\left(\theta,\varphi\right),\]

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) \overline{\xi\left(\theta,\varphi\right) }\sin\left(\theta\right) d\theta d\varphi,\]

with \(\overline{\xi\left(\theta,\varphi\right)}\) the conjugate of \(\xi\left(\theta,\varphi\right)\).

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

Reference systems

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)\).

Complex spherical harmonics and Helmholtz free space fundamental solution

Let \(\mathbf{r}\), \(\mathbf{r}'\) be vectors, whose spherical coordinates in the reference system \(j\) are \(\left(r_j,\theta_j, \varphi_j \right)\) and \(\left(r_j',\theta_j', \varphi_j' \right)\), respectively. For \(r_j>r_j'\) we have

\[g_j\left(\mathbf{r}, \mathbf{r}'\right) = i k_j \sum_{l=0}^{\infty} j_l(k_j r_j') h_l^{(1)}(k_j r_j) \sum_{m=-l}^l Y_{l,m,j}\left(\mathbf{r}\right) \overline{Y_{l,m,j}}\left(\mathbf{r}'\right),\]

with \(j_l\) and \(h_l^{(1)}\) the spherical Bessel function and the spherical Hankel function of first kind. This series and its term by term first derivatives with respect to \(r_j\) or \(r_j'\) are absolutely and uniformly convergent on compact subsets with \(r_j>r_j'\).

Let be \(j\in \{1,\ldots,N\}\), and \(\tilde{r}_j\) the radius of \(\Omega_j\).

• Let \(\mathbf{r} \in \Omega_0\) and \(\mathbf{r}=(r,\theta,\varphi )\) be expressed in the \(j\)th spherical coordinate system, then

  \[\begin{split}\begin{aligned} DL_{0j} ( Y_{l,m,j}) (\mathbf{r}) &= -i \left( k_0 \tilde{r}_j \right)^2 j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 r) Y_{l,m}(\theta, \varphi ),\\ SL_{0j} ( Y_{l,m,j}) (\mathbf{r}) &= i k_0 \tilde{r}_j^2 j_l (k_0 \tilde{r}_j) h_l^{(1)}(k_0 r) Y_{l,m}(\theta, \varphi ). \end{aligned}\end{split}\]

• Let \(\mathbf{r} \in \Omega_j\) and \(\mathbf{r}=(r,\theta,\varphi )\) expressed in the \(j\)th spherical coordinate system, then

  \[\begin{split}\begin{aligned} DL_{j} ( Y_{l,m,j}) (\mathbf{r}) &= i \left( k_j \tilde{r}_j \right)^2 h_l^{(1)'}(k_j \tilde{r}_j) j_l (k_j r) Y_{l,m}(\theta, \varphi ),\\ SL_{j} ( Y_{l,m,j}) (\mathbf{r}) &= i k_j \tilde{r}_j^2 h_l^{(1)}(k_j \tilde{r}_j) j_l (k_j r) Y_{l,m}(\theta, \varphi ). \end{aligned}\end{split}\]

The diagonal forms of the BIOs are:

\[\begin{split}\begin{aligned} V_{j,j}^0\left( Y_{l,m,j} \right) &= i k_0 \tilde{r}_j^2 j_l (k_0 \tilde{r}_j) h_l^{(1)}(k_0 k_0 \tilde{r}_j) Y_{l,m,j}, & V_{j}\left( Y_{l,m,j} \right)&= i k_0 \tilde{r}_j^2 j_l (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) Y_{l,m,j}, \\ K_{j,j}^0\left( Y_{l,m,j} \right) &= -\frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j)) Y_{l,m,j} , & K_{j}\left( Y_{l,m,j} \right) &= \frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j)) Y_{l,m,j},\\ K^{*0}_{j,j}\left( Y_{l,m,j} \right) &= -\frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j))Y_{l,m,j}, & K^*_{j}\left( Y_{l,m,j} \right) &= \frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j))Y_{l,m,j},\\ W_{j,j}^0\left( Y_{l,m,j} \right) &= i k_0^3 \tilde{r}_j^2 j_l' (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j) Y_{l,m,j}, & W_{j}\left( Y_{l,m,j} \right)&= i k_0^3 \tilde{r}_j^2 j_l' (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j) Y_{l,m,j}. \end{aligned}\end{split}\]

v_jj_azimuthal_symmetry

This function returns a numpy array, which we will call v. Its length is big_l + 1. Each entry contains

\[v[l] = ( V_{j,j}^0 Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = ( V_{j,j} Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = i \tilde{r}_j^4 k_0 j_l(k_0\tilde{r}_j) h_l^{(1)}(k_0\tilde{r}_j).\]

Docstring

print(selfinteractions.v_jj_azimuthal_symmetry.__doc__)

    Returns a numpy array with the boundary integral operator V_{j,j}
    with Helmholtz kernel evaluated and tested with spherical harmonics
    of order 0 in a sphere of radius r.

    Notes
    -----
    v[l] = ( V_{j,j} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r).
    = i r**4 k j_l(kr) h_l(kr)
    for each l such that 0 <= l <= big_l, and with
    ( ., .)_L^2(surface sphere radius r): inner product where indicated.
    Y_l,0: spherical harmonic degree l, order 0.
    j_l: spherical Bessel function.
    h_l: spherical Hanckel function of first kind.

    Notice that the expression could also be computed as:
    0.5 i r pi J_{l+0.5}(k r) H_{l+0.5}(k r)
    with J_{nu} the Bessel function and
    H_{nu} the Hanckel function of first kind.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.

    Returns
    -------
    v : np.ndarray
        with complex values. Length (big_l+1). See the section notes for
        the ordering of the array.

Example

Definition of parameters:

r = 2.5
big_l = 50
k = 7.

Calling of the function:

v = selfinteractions.v_jj_azimuthal_symmetry(big_l, r, k)

Plot of the values:

plt.figure()
plt.plot(
    np.real(v),
    marker='x',
    label='Real part of V'
)
plt.plot(
    np.imag(v),
    marker='x',
    label='Imaginary part of V'
)
plt.title('V')
plt.xlabel('l')
plt.legend()
plt.show()

k_0_jj_azimuthal_symmetry

This function returns a numpy array, which we will call data_k. Its length is big_l + 1. Each entry contains

\[\begin{split}\begin{aligned} data_k[l] &= ( K_{j,j}^0 Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = ( K_{j,j}^{0*} Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})}\\ &= -\frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j)) Y_{l,0} \end{aligned}\end{split}\]

where \(\mathbb{S_j}\) is the surface of the sphere \(j\) and \(r_j\) its radius.

Docstring

print(selfinteractions.k_0_jj_azimuthal_symmetry.__doc__)

    Returns a numpy array with the boundary integral operator K_{j,j}^0,
    with normals from out to in of the sphere, with Helmholtz kernel
    evaluated and tested with spherical harmonics of order 0.

    Notes
    -----
    k_0_jj[l] = (K_{j,j}^0 Y_l,0 ; Y_l,0)_L^2(surface sphere radius r).
    = -1j * k**2 * r**4 * 0.5
      * (j_l(r k) * h_l_d(r k) + j_l_d(r k) * h_l(r k))
    for each l such that 0 <= l <= big_l, and with
    ( ., .)_L^2(surface sphere radius r): inner product where indicated.
    Y_l,0: spherical harmonic degree l, order 0.
    j_l: spherical Bessel function, and j_l_d its derivative.
    h_l: spherical Hanckel function of first kind, and j_l_d its
         derivative.

    Notice that in this specific case
    ( K_{j,j}^0 Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r)
    = ( K_{j,j}^{*0} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r)
    where K_{j,j}^{*0} has normals from out to in of the sphere, and
    with Helmholtz kernel.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.

    Returns
    -------
    k_0 : np.ndarray
        with complex values. Length (big_l+1). See the section notes for
        the ordering of the array.

Example

Definition of parameters:

r = 2.5
big_l = 50
k = 7.

Calling of the function:

k_0_jj = selfinteractions.k_0_jj_azimuthal_symmetry(big_l, r, k)

Plotting the values:

plt.figure()
plt.plot(
    np.real(k_0_jj),
    marker='x',
    label='Real part of K'
)
plt.plot(
    np.imag(k_0_jj),
    marker='x',
    label='Imaginary part of K'
)
plt.title('K')
plt.xlabel('l')
plt.legend()
plt.show()

k_1_jj_azimuthal_symmetry

This function returns a numpy array, which we will call k_j. Its length is big_l + 1. Each entry contains

\[\begin{split}\begin{aligned} k_j[l] &= ( K_{j,j} Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = ( K_{j,j}^* Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})}\\ &= \frac{i}{2} (k_0 \tilde{r}_j)^2 (j_l' (k_0 \tilde{r}_j) h_l^{(1)}(k_0 \tilde{r}_j) + j_l (k_0 \tilde{r}_j) h_l^{(1)'}(k_0 \tilde{r}_j)) Y_{l,0} \end{aligned}\end{split}\]

where \(\mathbb{S_j}\) is the surface of the sphere \(j\) and \(r_j\) its radius.

Docstring

print(selfinteractions.k_1_jj_azimuthal_symmetry.__doc__)

    Returns a numpy array with the boundary integral operator K_{j,j},
    with normals from in to out of the sphere, with Helmholtz kernel
    evaluated and tested with spherical harmonics of order 0.

    Notes
    -----
    k_jj[l] = ( K_{j,j} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r).
    = 1j * k**2 * r**4 * 0.5
      * (j_l(r k) * h_l_d(r k) + j_l_d(r k) * h_l(r k))
    for each l such that 0 <= l <= big_l, and with
    ( ., .)_L^2(surface sphere radius r): inner product where indicated.
    Y_l,0: spherical harmonic degree l, order 0.
    j_l: spherical Bessel function, and j_l_d its derivative.
    h_l: spherical Hanckel function of first kind, and j_l_d its
         derivative.

    Notice that in this specific case
    ( K_{j,j} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r)
    = ( K_{j,j}^{*} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r)
    where K_{j,j}^{*} has normals from in to out of the sphere, and with
    Helmholtz kernel.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.

    Returns
    -------
    k_jj : np.ndarray
        with complex values, length (big_l+1). See the section notes for
        the ordering of the array.

Example

Definition of parameters:

r = 2.5
big_l = 50
k = 7.

Calling of the function:

k_1_jj = selfinteractions.k_1_jj_azimuthal_symmetry(big_l, r, k)

Plotting the values:

plt.figure()
plt.plot(
    np.real(k_1_jj),
    marker='x',
    label='Real part of K'
)
plt.plot(
    np.imag(k_1_jj),
    marker='x',
    label='Imaginary part of K'
)
plt.title('K')
plt.xlabel('l')
plt.legend()
plt.show()

w_jj_azimuthal_symmetry

This function returns a numpy array, which we will call w. Its length is big_l + 1. Each entry contains

\[w[l] = ( W_{j,j}^0 Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = ( W_{j,j} Y_{l,0} , Y_{l,0} )_{L^2(\mathbb{S_j})} = i \tilde{r}_j^4 k_0 j_l(k_0\tilde{r}_j) h_l^{(1)}(k_0\tilde{r}_j).\]

Docstring

print(selfinteractions.w_jj_azimuthal_symmetry.__doc__)

    Returns a numpy array with the boundary integral operator W_{j,j}
    with Helmholtz kernel evaluated and tested with spherical harmonics
    of order 0.

    Notes
    -----
    w[l] = ( W_{j,j} Y_l,0 ; Y_l,0 )_L^2(surface sphere radius r).
    = -1j k**3 * r**4 * j_l_d(k r) h_l_d(k r)
    for each l such that 0 <= l <= big_l, and with
    ( ., .)_L^2(surface sphere radius r): inner product where indicated.
    Y_l,0: spherical harmonic degree l, order 0.
    j_l_d: derivative of the spherical Bessel function.
    h_l_d: derivative of the spherical Hanckel function of first kind.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.

    Returns
    -------
    w: np.ndarray
        with complex values, length (big_l+1). See the section notes for
        the ordering of the array.

Example

Definition of parameters:

r = 2.5
big_l = 50
k = 7.

Calling of the function:

w = selfinteractions.w_jj_azimuthal_symmetry(big_l, r, k)

Plotting of the values:

plt.figure()
plt.plot(
    np.real(w),
    marker='x',
    label='Real part of W'
)
plt.plot(
    np.imag(w),
    marker='x',
    label='Imaginary part of W'
)
plt.title('W')
plt.xlabel('l')
plt.legend()
plt.show()

bio_jj

Docstring

print(selfinteractions.bio_jj.__doc__)

    Returns a numpy array with the corresponding boundary integral
    operator from the function bio_azimuthal with Helmholtz kernel
    evaluated and tested with spherical harmonics of all orders.

    Notes
    -----
    D[l*(2l+1) + m] = (D_{j,j}Y_l,m;Y_l,m)_L^2(surface sphere radius r).
    which in this case is equal to
    = D[l*(2l+1)]
    for each l such that 0 <= l <= big_l, and with
    ( ., .)_L^2(surface sphere radius r): inner product where indicated.
    D: boundary integral operator corresponding to the function
       bio_azimuthal
    Y_l,m: spherical harmonic degree l, order m.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.
    bio_azimuthal : Callable[[int, float, float]
        python function for computing the operator for m = 0. Must be
        one of the functions written before this one.

    Returns
    -------
    np.ndarray
        with complex values, length (big_l+1)**2. See the section notes
        for the ordering of the array.

Example

Definition of the parameters:

r = 2.5
big_l = 3
k = 7.
operator = selfinteractions.w_jj_azimuthal_symmetry

Calling of the function:

full = selfinteractions.bio_jj(big_l, r, k, operator)

print(full)

[-11.59514193 -3.3131677j   13.61117226-39.00171561j
  13.61117226-39.00171561j  13.61117226-39.00171561j
 -17.04878169 -8.25398472j -17.04878169 -8.25398472j
 -17.04878169 -8.25398472j -17.04878169 -8.25398472j
 -17.04878169 -8.25398472j  20.40932274-28.34521185j
  20.40932274-28.34521185j  20.40932274-28.34521185j
  20.40932274-28.34521185j  20.40932274-28.34521185j
  20.40932274-28.34521185j  20.40932274-28.34521185j]

a_0j_matrix

Docstring

print(selfinteractions.a_0j_matrix.__doc__)

    Returns a numpy array that represents the following matrix boundary
    integral operator
    A_{j,j}^0 = [ -K_{j,j}^0 , V_{j,j}^0  ]
                [  W_{j,j}^0 , K*_{j,j}^0 ]
    with Helmholtz kernel evaluated and tested with spherical harmonics
    of order 0 if azimuthal = True, or all orders if azimuthal = False.

    Each block is a diagonal matrix.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.
    azimuthal : bool
        Default True.

    Returns
    -------
    matrix_0j : np.ndarray
        with complex values.
        If azimuthal = True
            Shape (2*(big_l+1), 2*(big_l+1))
        Else
            Shape (2*(big_l+1)**2, 2*(big_l+1)**2)

    See Also
    --------
    a_0j_linear_operator
    v_jj_azimuthal_symmetry
    k_0_jj_azimuthal_symmetry
    w_jj_azimuthal_symmetry
    bio_jj

Example

Definition of the parameters:

r = 2.5
big_l = 2
k = 7.

Calling of the function with azimuthal = True:

a_0j = selfinteractions.a_0j_matrix(big_l, r, k) # azimuthal = True

Plotting of the values:

fig = plt.figure(figsize=plt.figaspect(.5), dpi=68., layout='constrained')
ax_1 = fig.add_subplot(121)
ax_1.set_title("Real part of a_0j")
ax_1.imshow(np.real(a_0j), vmin=np.min(np.real(a_0j)), vmax=np.max(np.real(a_0j)))
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.real(a_0j)), vmax=np.max(np.real(a_0j)))), ax=ax_1, orientation="horizontal", shrink=0.75)
ax_1 = fig.add_subplot(122)
ax_1.imshow(np.imag(a_0j), vmin=np.min(np.imag(a_0j)), vmax=np.max(np.imag(a_0j)))
ax_1.set_title("Imaginary part of a_0j")
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.imag(a_0j)), vmax=np.max(np.imag(a_0j)))), ax=ax_1, orientation="horizontal", shrink=0.75)

<matplotlib.colorbar.Colorbar at 0x7dd95f9848f0>

Calling of the function with azimuthal = False:

a_0j = selfinteractions.a_0j_matrix(big_l, r, k, azimuthal=False)

Plotting of the values:

fig = plt.figure(figsize=plt.figaspect(.5), dpi=68., layout='constrained')
ax_1 = fig.add_subplot(121)
ax_1.set_title("Real part of a_0j")
ax_1.imshow(np.real(a_0j), vmin=np.min(np.real(a_0j)), vmax=np.max(np.real(a_0j)))
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.real(a_0j)), vmax=np.max(np.real(a_0j)))), ax=ax_1, orientation="horizontal", shrink=0.75)
ax_1 = fig.add_subplot(122)
ax_1.imshow(np.imag(a_0j), vmin=np.min(np.imag(a_0j)), vmax=np.max(np.imag(a_0j)))
ax_1.set_title("Imaginary part of a_0j")
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.imag(a_0j)), vmax=np.max(np.imag(a_0j)))), ax=ax_1, orientation="horizontal", shrink=0.75)

<matplotlib.colorbar.Colorbar at 0x7dd95f84c950>

a_j_matrix

Docstring

print(selfinteractions.a_j_matrix.__doc__)

    Returns a numpy array that represents the following matrix boundary
    integral operator
    A_{j,j}^0 = [ -K_{j,j} , V_{j,j}  ]
                [  W_{j,j} , K*_{j,j} ]
    with Helmholtz kernel evaluated and tested with spherical harmonics
    of order 0 if azimuthal = True, or all orders if azimuthal = False.

    Each block is a diagonal matrix.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.
    azimuthal : bool
        Default True.

    Returns
    -------
    matrix_j : np.ndarray
        with complex values.
        If azimuthal = True
            Shape (2*(big_l+1), 2*(big_l+1))
        Else
            Shape (2*(big_l+1)**2, 2*(big_l+1)**2)

    See Also
    --------
    a_j_linear_operator
    v_jj_azimuthal_symmetry
    k_1_jj_azimuthal_symmetry
    w_jj_azimuthal_symmetry
    bio_jj

Example

Definition of the parameters:

r = 2.5
big_l = 2
k = 7.

Calling of the function with azimuthal = True:

a_j = selfinteractions.a_j_matrix(big_l, r, k) # azimuthal = True

Plotting of the values:

fig = plt.figure(figsize=plt.figaspect(.5), dpi=68., layout='constrained')
ax_1 = fig.add_subplot(121)
ax_1.set_title("Real part of a_j")
ax_1.imshow(np.real(a_j), vmin=np.min(np.real(a_j)), vmax=np.max(np.real(a_j)))
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.real(a_j)), vmax=np.max(np.real(a_j)))), ax=ax_1, orientation="horizontal", shrink=0.75)
ax_1 = fig.add_subplot(122)
ax_1.imshow(np.imag(a_j), vmin=np.min(np.imag(a_j)), vmax=np.max(np.imag(a_j)))
ax_1.set_title("Imaginary part of a_j")
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.imag(a_j)), vmax=np.max(np.imag(a_j)))), ax=ax_1, orientation="horizontal", shrink=0.75)

<matplotlib.colorbar.Colorbar at 0x7dd95f9afef0>

Calling of the function with azimuthal = False:

a_j = selfinteractions.a_j_matrix(big_l, r, k, azimuthal=False)

Plotting of the values:

fig = plt.figure(figsize=plt.figaspect(.5), dpi=68., layout='constrained')
ax_1 = fig.add_subplot(121)
ax_1.set_title("Real part of a_j")
ax_1.imshow(np.real(a_j), vmin=np.min(np.real(a_j)), vmax=np.max(np.real(a_j)))
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.real(a_j)), vmax=np.max(np.real(a_j)))), ax=ax_1, orientation="horizontal", shrink=0.75)
ax_1 = fig.add_subplot(122)
ax_1.imshow(np.imag(a_j), vmin=np.min(np.imag(a_j)), vmax=np.max(np.imag(a_j)))
ax_1.set_title("Imaginary part of a_j")
fig.colorbar(cm.ScalarMappable(norm=colors.Normalize(vmin=np.min(np.imag(a_j)), vmax=np.max(np.imag(a_j)))), ax=ax_1, orientation="horizontal", shrink=0.75)

<matplotlib.colorbar.Colorbar at 0x7dd95f7d1490>

a_0j_linear_operator

Docstring

print(selfinteractions.a_0j_linear_operator.__doc__)

    Returns a scipy linear operator equivalent to the given by
    a_0j_matrix.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.
    azimuthal : bool
        Default True

    Returns
    -------
    linear_operator : sparse.linalg.LinearOperator
        a scipy linear operator.

    See Also
    --------
    a_0j_matrix

Examples

Simpler example

Definition of parameters:

r = 2.5
big_l = 2
k = 7.

Calling of the function with azimuthal = True:

a_0j = selfinteractions.a_0j_linear_operator(big_l, r, k) # azimuthal = True

Calling of the function with azimuthal = False:

a_0j = selfinteractions.a_0j_matrix(big_l, r, k, azimuthal = False)

Comparing with the matrix version

Definition of parameters:

big_l = 10
r = 0.3
k = 1.3
tole = 10.**-10

Initialization of a random vector as right hand side:

num = big_l + 1
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))

Calling of the functions to be compared

linear_operator = selfinteractions.a_0j_linear_operator(big_l, r, k) # azimuthal = True
matrix = selfinteractions.a_0j_matrix(big_l, r, k) # azimuthal = True

Solving of the iterative system:

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**3,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Solving of the direct system:

solution2 = np.linalg.solve(matrix, b)

Difference between solutions:

print(np.linalg.norm(solution2 - solution))

4.5723754933069785e-12

Which means that solving the numpy array from a_0j_matrix and the linear operator from a_0j_linear_operator can be used to solve the same linear system.

Let us see what happens with azimuthal=False.

Definition of parameters:

big_l = 10
r = 0.3
k = 1.3
tole = 10.**-5

Initialization of a random vector as right hand side:

num = (big_l + 1)**2
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))

Calling of the functions to be compared

linear_operator = selfinteractions.a_0j_linear_operator(big_l, r, k, False)
matrix = selfinteractions.a_0j_matrix(big_l, r, k, False)

Solving of the iterative system:

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**3,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Solving of the direct system:

solution2 = np.linalg.solve(matrix, b)

Difference between solutions:

print(np.linalg.norm(solution2 - solution))

2.5550504989755213e-10

Comparing the azimuthal version against the no azimuthal

For a right hand with azimuthal symmetry when extended by 0 to explicitly compute the result without considering the symmetry, the result should be the same with the azimuthal and the not azimuth versions.

Definition of parameters:

big_l = 10
r = 1.3
k = 2
tole = 10.**-7

Initialization of a random vector as right hand side, and its extension by zero:

import biosspheres.miscella.extensions as righthandsextensions
num = big_l + 1
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))
b_2_1 = righthandsextensions.azimuthal_trace_to_general_with_zeros(
    big_l, b[0:num])
b_2_2 = righthandsextensions.azimuthal_trace_to_general_with_zeros(
    big_l, b[num:2 * num])
b2 = np.concatenate((b_2_1, b_2_2))

Calling of the function with azimuthal=False:

linear_operator = selfinteractions.a_0j_linear_operator(big_l, r, k, False)

Solving of the iterative system:

num = (big_l + 1)**2

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b2,
    rtol=tole,
    restart=(2 * num)**2,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Calling of the function with azimuthal=True:

linear_operator = selfinteractions.a_0j_linear_operator(big_l, r, k) # azimuthal = True

Solving of the iterative system:

num = big_l + 1

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution2, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**2,
    callback=callback_function,
    callback_type='pr_norm',
)
print(info)

0

Comparison of the solutions:

solution2_1 = \
    righthandsextensions.azimuthal_trace_to_general_with_zeros(
        big_l, solution2[0:num])
solution2_2 = \
    righthandsextensions.azimuthal_trace_to_general_with_zeros(
        big_l, solution2[num:2 * num])
solution2 = np.concatenate((solution2_1, solution2_2))

print(np.linalg.norm(solution2 - solution))

7.940632618156166e-15

a_j_linear_operator

Docstring

print(selfinteractions.a_j_linear_operator.__doc__)

    Returns a scipy linear operator equivalent to the given by
    a_j_matrix.

    Parameters
    ----------
    big_l : int
        >= 0, max degree.
    r : float
        > 0, radius.
    k : float
        wave number.
    azimuthal : bool
        Default True.

    Returns
    -------
    linear_operator : sparse.linalg.LinearOperator
        a scipy linear operator.

    See Also
    --------
    a_j_matrix

Examples

Simple example

Definition of parameters:

r = 2.5
big_l = 2
k = 7.

Calling of the function with azimuthal=True:

a_0j = selfinteractions.a_0j_linear_operator(big_l, r, k) # azimuthal = True

Calling of the function with azimuthal=False:

a_0j = selfinteractions.a_0j_matrix(big_l, r, k, azimuthal = False)

Comparing with the matrix version

Definition of parameters:

big_l = 10
r = 0.3
k = 1.3
tole = 10.**-10

Initialization of a random vector as right hand side:

num = big_l + 1
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))

Calling of the functions to be compared

linear_operator = selfinteractions.a_j_linear_operator(big_l, r, k) # azimuthal = True
matrix = selfinteractions.a_j_matrix(big_l, r, k) # azimuthal = True

Solving of the iterative system:

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**3,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Solving of the direct system:

solution2 = np.linalg.solve(matrix, b)

Difference between solutions:

print(np.linalg.norm(solution2 - solution))

1.3355528626430198e-11

Which means that solving the numpy array from a_0j_matrix and the linear operator from a_0j_linear_operator can be used to solve the same linear system.

Let us see what happens with azimuthal=False.

Definition of parameters:

big_l = 10
r = 0.3
k = 1.3
tole = 10.**-5

Initialization of a random vector as right hand side:

num = (big_l + 1)**2
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))

Calling of the functions to be compared

linear_operator = selfinteractions.a_j_linear_operator(big_l, r, k, False)
matrix = selfinteractions.a_j_matrix(big_l, r, k, False)

Solving of the iterative system:

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**3,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Solving of the direct system:

solution2 = np.linalg.solve(matrix, b)

Difference between solutions:

print(np.linalg.norm(solution2 - solution))

1.7724454839951905e-10

Comparing the azimuthal version against the no azimuthal

For a right hand with azimuthal symmetry when extended by 0 to explicitly compute the result without considering the symmetry, the result should be the same with the azimuthal and the not azimuth versions.

Definition of parameters:

big_l = 10
r = 1.3
k = 2
tole = 10.**-7

Initialization of a random vector as right hand side, and its extension by zero:

import biosspheres.miscella.extensions as righthandsextensions
num = big_l + 1
b = np.random.random((2 * num)) + 1j * np.random.random((2 * num))
b_2_1 = righthandsextensions.azimuthal_trace_to_general_with_zeros(
    big_l, b[0:num])
b_2_2 = righthandsextensions.azimuthal_trace_to_general_with_zeros(
    big_l, b[num:2 * num])
b2 = np.concatenate((b_2_1, b_2_2))

Calling of the function with azimuthal=False:

linear_operator = selfinteractions.a_j_linear_operator(big_l, r, k, False)

Solving of the iterative system:

num = (big_l + 1)**2

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution, info = sparse.linalg.gmres(
    linear_operator, b2,
    rtol=tole,
    restart=(2 * num)**2,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Calling of the function with azimuthal=True:

linear_operator = selfinteractions.a_j_linear_operator(big_l, r, k) # azimuthal = True

Solving of the iterative system:

num = big_l + 1

norms = []

def callback_function(pr_norm):
    norms.append(pr_norm)
    return 0.

solution2, info = sparse.linalg.gmres(
    linear_operator, b,
    rtol=tole,
    restart=(2 * num)**2,
    callback=callback_function,
    callback_type='pr_norm')
print(info)

0

Comparison of the solutions:

solution2_1 = \
    righthandsextensions.azimuthal_trace_to_general_with_zeros(
        big_l, solution2[0:num])
solution2_2 = \
    righthandsextensions.azimuthal_trace_to_general_with_zeros(
        big_l, solution2[num:2 * num])
solution2 = np.concatenate((solution2_1, solution2_2))

print(np.linalg.norm(solution2 - solution))

9.130324111516091e-15

a_0_a_n_sparse_matrices

Example

Definition of parameters:

n = 2
big_l = 5
r = 1.3
radii = np.ones(n) * r
kii = np.ones(n+1)
k = 7.
kii[0] = k

Calling of the function with azimuthal = True:

sparse_big_a_0_self, sparse_big_a_n = selfinteractions.a_0_a_n_sparse_matrices(n, big_l, radii, kii) # azimuthal = False

Calling of the function with azimuthal = False:

sparse_big_a_0_self, sparse_big_a_n = selfinteractions.a_0_a_n_sparse_matrices(n, big_l, radii, kii, azimuthal = True)

reduced_a_sparse_matrix

Example

Definition of parameters:

n = 2
big_l = 5
r = 1.3
radii = np.ones(n) * r
kii = np.ones(n+1)
k = 7.
kii[0] = k
pi0 = 0.5
pii = np.ones(n) * pi0

Calling of the function with azimuthal = True:

reduced_a = selfinteractions.reduced_a_sparse_matrix(n, big_l, radii, kii, pii) # azimuthal = False

Calling of the function with azimuthal = False:

reduced_a = selfinteractions.reduced_a_sparse_matrix(n, big_l, radii, kii, pii, azimuthal=True)