Features implemented

b̃p(p; k=1)

Closure for the circular wave function in polar coordinates.

The closure captures the values of the mode number $p$ as well as the wavenumber $k$ and does not include any normalization. The function can then be evaluated at any $(r,θ)$ point.

The circular waves in polar coordinates are defined by

\[b̃_p := (r,θ) ↦ J_{p}(kr) e^{ı p θ}, \qquad ∀ p ∈ \mathbb{Z}.\]

βp(p; k=1)

Computation of the normalization constant of the circular waves.

The circular wave is defined by b̃p and the normalization is done using the $k$-weighted $H^1$ norm.

By definition

\[ \| b̃_p \|_{H^{1}}^2 = \| b̃_p \|_{L^{2}}^2 + k^{-2} \| \nabla b̃_p \|_{L^{2}}^2.\]

Using integration by parts,

\[ \| \nabla b̃_p \|_{L^{2}}^2 = k^2 \| b̃_p \|_{L^{2}}^2 + (\partial_n b̃_p , b̃_p).\]

On the one hand

\[ \| b̃_p \|_{L^{2}}^2 = \pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k)).\]

On the other hand

\[ (\partial_n b̃_p , b̃_p) = \pi k (J_{p-1}(k) - J_{p+1}(k))J_p(k).\]

Hence

\[ \| b̃_p \|_{H^{1}}^2 = 2\pi (J_p^2(k) - J_{p-1}(k)J_{p+1}(k)) + \pi k^{-1} (J_{p-1}(k) - J_{p+1}(k))J_p(k).\]

bp(p; k=1)

Closure for the normalized circular waves.

The closure captures the values of the mode number $p$ and the wavenumber $k$. The normalization is the default normalization defined in the function βp and is precomputed once for all. The function can then be evaluated at any $(r,θ)$ point.

The normalized circular waves are defined by

\[b_p := β_p b̃_p, \qquad ∀ p ∈ \mathbb{Z}.\]

solution_surrogate(U; k=1)

Computes a solution surrogate from a set of coefficients U.

The parameter U is a vector of coefficients $u_p$ for $p$ ranging from $-P$ to $P$ (hence of size $2P+1$). The solution surrogate is then

\[\mathbf{x} ↦ \sum_{|p| \leq P} u_p b_{p}(\mathbf{x}).\]

ãp(p)

Closure for the Herglotz polynomial.

The closure captures the values of the mode number $p$ and does not include any normalization. The function can then be evaluated at any $(ζ,φ)$ point in the complex strip. Here $ζ$ is the evanescence parameter and $φ$ is the angle indicating the direction of propagation.

The Herglotz polynomials are defined by

\[ã_p := (ζ,φ) ↦ e^{p (ζ + ıφ)}, \qquad ∀ p ∈ \mathbb{Z}.\]

wz(; k=1, z=1/4)

Closure to define the weight function.

The closure captures the values of the wavenumber $k$ and a parameter $z$. The function can then be evaluated at any $ζ$ point.

The weight function is defined as

\[w := ζ ↦ e^{z |ζ| - k \sinh|ζ|}, \qquad ∀ ζ ∈ \mathbb{R}.\]

αp(p, logweight)

Computation of the normalization constant of the Herglotz polynomial.

The Herglotz polynomial is defined by $ã_p$, see the function ãp. The normalization constant is defined by

\[α_p^{-2} := 2π \int_{\mathbb{R}} |ã_p|^{2} w^2 \mathrm{d}\zeta,\]

where $w$ is the weight function. The weight function $w$ can be taken to be the function wz.

The implementation relies on the convenience functions ϵsupport and adaptive_G_quad.

ap(p, logweight)

Closure for the normalized Herglotz polynomial.

The closure captures the values of the mode number $p$. The Herglotz polynomial is defined by $ã_p$, see the function ãp. The normalization is given by the second argument, as defined in the documentation of the function αp and is precomputed once for all. The function can then be evaluated at any $(ζ,φ)$ point in the complex strip. Here $ζ$ is the evanescence parameter and $φ$ is the angle indicating the direction of propagation.

The normalized Herglotz polynomials are defined by

\[a_p := α_p ã_p, \qquad ∀ p ∈ \mathbb{Z}.\]

gpw(ζ, φ; k=1)

Closure that defines a generalized plane wave in polar coordinates.

The closure captures the values of the evanescence parameter $ζ$, the angle $φ$ and the wavenumber $k$. It does not include any normalization. The function can then be evaluated at any $(r,θ)$ point.

The generalized plane wave in polar coordinates is defined by

\[ϕ := (r,θ) ↦ e^{ı k \mathbf{d} ⋅ \mathbf{x}},\]

where

\[\mathbf{d} = (\cos[φ+ıζ], \sin[φ+ıζ]), \qquad\text{and}\qquad \mathbf{x} = (r \cos θ, r\sin θ).\]

approximation_set(Y, W; k=1)

Generate a set of generalized plane waves.

The parameters Y and W are respectively the sets of complex angles $(ζ_j,φ_j)$ and associated weights $ω_j$ corresponding to the generalized plane waves, see gpw

\[ϕ_j := (r,θ) ↦ e^{ı k \mathbf{d}_j ⋅ \mathbf{x}},\]

where

\[\mathbf{d}_j = (\cos[φ_j+ıζ_j], \sin[φ_j+ıζ_j]), \qquad\text{and}\qquad \mathbf{x} = (r \cos θ, r\sin θ).\]

Approximations are constructed in the form of

\[(r, θ) ↦ \sum_j ξ_j ω_j ϕ_j(r, θ),\]

where $(ξ_j)_j$ is the set of unknown coefficients.

The value of the wavenumber k defaults to $1$ and can be provided as an optional parameter.

approximation_set(N; k=1)

Generate a set of $N$ propagative plane waves with equispaced angles.

The equispaced angles are defined by

\[θ_n := 2π s / R, \qquad n = 0,…,N-1.\]

The value of the wavenumber k defaults to $1$ and can be provided as an optional parameter.

approximation_set(N, qs, smpl_type; k=1)
approximation_set(N, Q::Integer, smpl_type; k=1)

Generate a set of $N$ evanescent plane waves according to some sampling type.

The second parameter can be an integer $Q$ which controls the truncation parameter of the TruncKernel, or a set of integers qs to specify only some specific terms in the expansion. The parameter smpl_type should be a function name: see uniform_sampling, sobol_sampling and random_sampling.

The value of the wavenumber k defaults to $1$ and can be provided as an optional parameter.

RegularizedSVDPseudoInverse(A; ϵ=1e-8)

Structure able to solve linear system A using a regularized SVD.

The threshold ϵ controls the regularization. Let $σ_{\max}$ be the larger singular value. A singular value $σ$ such that $σ \leq ϵ σ_{\max}$ is approximated by zero when solving the least-squares problem.

condition_number(B::RegularizedSVDPseudoInverse)

Computes condition number of the matrix for which B was constructed.

The condition number is defined as the ratio of the largest over the smallest singular values of the matrix $A$.

\[κ := \frac{σ_{\max}}{σ_{\min}}.\]

solve_via_regularizedSVD(B::RegularizedSVDPseudoInverse, b)

Solve the linear system $Ax=b$ using the regularized SVD of $A$.

TruncKernel(P::Integer, logweight)
TruncKernel(ps::NTuple{N,Int64}, logweight) where N

Structure to hold precomputed normalization constants of Herglotz polynomials.

This is a convenience structure to hold precomputed normalization constants αp of Herglotz polynomials which are expensive to compute.

All the normalized Herglotz polynomials $a_p$ (see ap) for $p$ in the argument ps (noted $\mathcal{P}$) are precomputed and stored in the attribute as which is a Tuple. The logweight function is used to compute the normalization constants. If only an Integer $P$ is given as first argument, the terms in the kernel are the $a_p$ functions for $p$ ranging from $-P$ to $P$.

A kernel can be evaluated at (x,y) where x=(ζx,φx) and y=(ζy,φy) are points in the complex strip.

\[(\mathbf{x},\mathbf{y}) ↦ K(\mathbf{x}, \mathbf{y}) = \sum_{p \in \mathcal{P}} \overline{a_{p}(\mathbf{x})} a_{p}(\mathbf{y}).\]

kernel_diagonal(K::TruncKernel)

Closure that defines the evaluation function of the diagonal of a TruncKernel.

\[\mathbf{y} ↦ K(\mathbf{y}, \mathbf{y}) = \sum_{p \in \mathcal{P}} |a_{p}(\mathbf{y})|^2, \qquad\mathbf{y}=(ζ,φ).\]

probabilitydensityfunction(K::TruncKernel)

Probability density function defined by a TruncKernel.

\[(ζ, φ) ↦ \frac{1}{|\mathcal{P}|} K(\mathbf{y}, \mathbf{y}) = \frac{1}{|\mathcal{P}|} \sum_{p \in \mathcal{P}} |a_{p}(\mathbf{y})|^2, \qquad\mathbf{y}=(ζ,φ).\]

Notice that it is indeed a PDF since all the $a_p$ functions are orthonormal.

cumulativedensityfunction(K::TruncKernel)

Cumulative density function defined by a TruncKernel.

The probability density function is provided by the function probabilitydensityfunction.

\[ζ ↦ \frac{2π}{|\mathcal{P}|} \int_{-\infty}^{ζ} K(ζ,ζ) \mathrm{d}ζ = \frac{2π}{|\mathcal{P}|} \int_{-\infty}^{ζ} \sum_{p \in \mathcal{P}} |a_{p}|^2(ζ) \mathrm{d}ζ,\]

where, abusing notations, we used $K(ζ,ζ) = K(\mathbf{y},\mathbf{y})$ and $|a_{p}|^2(ζ) = |a_{p}|^2(\mathbf{y})$ for any $\mathbf{y} = (ζ,φ)$ since they both are quantities independent of $φ$.

The implementation relies on the convenience functions ϵsupport and adaptive_G_quad.

Sampler(K::TruncKernel)

Structure to compute sampling nodes and weights.

The nodes and weights are sampled according to the probability density function defined by a TruncKernel.

inversion(smpl::Sampler, u)

Computes the pre-image $ζ$ of u under the cumulative density function.

This function allows to perform Inversion Transform Sampling.

This is a univariate inversion function because the cdf is constant with respect to the other variable $φ$.

weight(smpl::Sampler, node)

Computes the weight associated to a node $(ζ,φ)$.

The weight at $\mathbf{y}=(ζ,φ)$ is defined by

\[ω(\mathbf{y}) := \sqrt{\frac{|\mathcal{P}|}{\sum_{p \in \mathcal{P}} |a_{p}(\mathbf{y})|^2}}, \qquad\mathbf{y}=(ζ,φ).\]

random_sampling(smpl::Sampler)
random_sampling(smpl::Sampler, M)

Draw M random samples according to the distribution stored in smpl.

uniform_sampling(smpl::Sampler, M)

Draw M deterministic samples according to the distribution stored in smpl.

This has a tensor-like structure, with the same number of samples in both directions.

sobol_sampling(smpl::Sampler, M)

Draw M quasi-random samples according to the distribution stored in smpl.

This function uses Sobol sequences computed according to the package QuasiMonteCarlo.jl.

samples_from_nodes(f, X)
samples_from_nodes(f::Vector, X)

Evaluate f at the sampling nodes X.

If f is a vector, a matrix is computed.

boundary_sampling_nodes(S)

The sampling nodes on the boundary of the unit disk.

The samplings nodes are then $(1, θ_s)$ with

\[θ_s := 2π s / S, \qquad s = 0,…,S-1.\]

number_of_boundary_sampling_nodes(M; η=2, P=0)

Number of sampling nodes necessary for an approximation set of dimension $M$.

The number of sampling points on the boundary of the unit disk is given by $S := \max(ηM, 2P+1)$. The (overdetermined) linear system that will be solved is then of size $S$ by $M$.

The oversampling parameter $η$ defaults to $2$ but can be provided by the user as an optional parameter η. The mode number $P$ defaults to $0$ but can be provided by the user as an optional parameter P.

Dirichlet_sampling(k, N, U; smpl_type=nothing, η=2, ϵ=1e-14, Q=(length(U)-1)//2)

Example of reconstruction of a solution surrogate via Dirichlet sampling.

The solution surrogate is defined by its vector of coefficients U. The approximation set of dimension N can be composed of propagative plane waves (default) or evanescent plane waves (depending on the value of smpl_type and Q the maximum mode number assumed to be in the target). The number of sampling points on the boundary can be controlled via the oversampling ration η. The amount of regularization in the SVD can be controlled by the regularization parameter ϵ.

function ϵsupport(f, ϵ; δ=1e-3)

Computes the ϵ-support of a function.

This is useful to compute quadratures of compactly ϵ-supported functions on unbounded domains.

adaptive_G_quad(f; quadrule=gausslegendre, a=-1, b=1)

Computes the integral of $f$ on $[a,b]$ using adaptive quadrature rule.

Here adaptivity is only understood with respect to the number of nodes. There is no-subdivision of the interval.