Singularity extraction for elliptic equations with coefficients with jumps and Dirac sources

Abstract

We explore the mathematical structure of the solution to an elliptic diffusion problem with point-wise Dirac sources. The conductivity parameter is space-varying, may have jumps and the Dirac sources may be located along the discontinuity curves of that parameter. The variational problem, issued by duality, is proven to be well posed using a sharp elliptic regularity result by Di-Giorgi [Mem. Accad. Sci. Torino, 3, 1957]. The paper is aimed at a key expansion into a split singular/regular contributions. The singular part is calculated by an explicit formula, while the regular correction can be computed as the solution to a standard variational Poisson problem. The latter can be successfully approximated by most of the numerical methods practiced nowadays. Some analytical examples are discussed at last to assess the minimality of the assumptions we use to establish our theoretical results.

Keywords

Elliptic Problems point-wise Dirac sources non-smooth conductivity singular/regular splitting

1. Introduction: Setting of the problem

Let Ω be a bounded connected domain in $R^{d}$ , $d = 2$ or 3, with a Lipschitz-continuous boundary $Γ = \partial Ω$ . The generic point in Ω is denoted by $x$ and the generic point on Γ by $τ$ , while $n$ stands for the unit normal vector (to Γ) which is outward to Ω.

The problem we deal with consists in finding a potential φ created by a (point-wise) source term $f (\cdot)$ . It is subject to a homogeneous Dirichlet boundary condition, for sake of simplicity. Different boundary conditions are possible and do not greatly alter the analysis. The purpose is to solve the following Poisson model $\begin{matrix} (1) & \begin{matrix} - div (a \nabla φ) = f in Ω, \\ φ = 0 on Γ . \end{matrix} \end{matrix}$ The conduction parameter $a (\cdot)$ is a given bounded non-degenerate function in the sense that $\begin{matrix} (2) & \forall x \in Ω, 0 < a_{*} ⩽ a (x) ⩽ a^{*} . \end{matrix}$ Unless explicitly mentioned, no additional smoothness is granted to $a (\cdot)$ other than it belongs to $L^{\infty} (Ω)$ .

The basic aim is to provide a relevant practical and computational way to handle problem (1), where the potential is created by point-wise Dirac measures $\begin{matrix} (3) & f (\cdot) = \sum_{s \in S} λ_{s} δ_{s} . \end{matrix}$ The set S of point-wise sources is a finite subset of Ω. The points in S are the different locations of as many internal Dirac sources; they are all at a given distance away from the boundary Γ. The ${(λ_{s})}_{s \in S}$ are their related intensities. The particularity and cause of difficulty here is that the conduction parameter $a (\cdot)$ is space varying, may have jumps and some of the Dirac sources are located at the interfaces of discontinuity. The target is to exhibit accurate informations on the singular behavior of the solution φ, by confronting it to the convolutional Green function, the fundamental solution of the Laplace operator.

A variety of applications can be found in the literature. We quote first from the introduction in [9] “Problems of this type occur in many applications from different areas, like in the mathematical modeling of electromagnetic fields [ 21 ]. Dirac [sources] can also be found on the right-hand side of adjoint equations in optimal control of elliptic problems with state constraints [ 11 ]. As further examples where such measures play an important role, we mention […] parameter identification problems with point-wise measurements [ 26 ].”

More examples come from the modeling of epilepsy seizures. They are associated with a diffusive abnormal electrical activity that originate from small cortex generators. Epilepsy foci are hence represented by point-wise Dirac current sources. The characteristic of the multi-layer brain is its piecewise continuous conductivity. At last, epileptic sources may be located at the cortex interfaces (see [19]). This is the very context the paper is involved with: diffusion equations, space dependent conductivities with jumps and Dirac point-wise sources. Studying the singular behavior of the solution is worth doing in its own interest and useful for direct simulations [29,31] or for the inverse Dirac sources detection [6,7,13,16]. The same mathematical problem arises in the study of the electro-cardiac potential in the heart [23].

The outline of the paper is as follows. Section 2 is devoted to the ‘non-regular’ variational formulation to problem (1). A duality argument sets the ‘natural’ variational problem where the Lebesgue regularity is the only requirement. The qualifier ‘natural’ is to contrast with the variational formulation in the weighted Sobolev spaces (see, eg, [4]). Given the wide range of admissible conductivities, they are only bounded (and bounded away from zero), the well posedness relies on a sharp continuity result, stated in the late fifties by Di Giorgi and Nash (see [15,25]). In Section 3, a key expansion of the solution φ is established, where the singular contribution is explicited by means of the fundamental solution to the Laplace operator. The correction function (φ minus the singular part) solves a regular elliptic variational problem – set in the Sobolev space $H^{1}$ . This has a welcome advantage in the practical grounds: many numerical methods can be used for the computation of that correction term. Some mild smoothness on the conductivity, limited to the vicinity of Dirac sources, are expected for the proof. Define ‘some limits’ of the conductivity at those points is the least one can demand. An analogous singularity splitting process is studied in Section 4, when the conductivity has jumps along straight-lines and the Dirac source is located at those discontinuity lines. Cut-off functions are then introduced around each Dirac source, so as to delimit the propagation of the singular extraction to small regions surrounding Dirac points. Section 5 is an extension of the configuration where the jumps occur along curved interfaces. In the last Section 6, some analytical examples are presented to discuss the assumptions used in the study.

Some notations. The Lebesgue space $L^{2} (Ω)$ of square integrable complex valued functions is endowed with the natural norm $‖ \cdot ‖_{L^{2} (Ω)}$ . We need also some Sobolev spaces, $H^{1} (Ω)$ involves all the functions that are in $L^{2} (Ω)$ as well as their partial derivatives. The subspace containing the functions in $H^{1} (Ω)$ that vanish on Γ is denoted by $H_{0}^{1} (Ω)$ . The set of the traces over Γ of all the functions of $H^{1} (Ω)$ is denoted $H^{1 / 2} (Γ)$ and $H^{- 1 / 2} (Γ)$ is its dual (see [2]). For a positive real number $ν < 1$ , the Hölder space is defined by $\begin{matrix} C^{0, ν} (Ω) = {χ \in C (Ω), | χ |_{C^{0, ν}} = sup_{x, y \in Ω} \frac{| χ (x) - χ (y) |}{| x - y |^{ν}} < \infty} . \end{matrix}$ It is a Banach space if endowed with the natural norm $‖ \cdot ‖_{C^{0, ν}} = ‖ \cdot ‖_{C} + | \cdot |_{C^{0, ν}}$ . The symbol [·] stands for the jump across a given boundary. We will also use the symbol ${(ν)}_{-}$ for an arbitrary real number $< ν$ , generally close to ν.

2. Variational framework and problem setting

The starting point is a functional Hilbertian setting where one can write a variational formulation of (1), when the data $f (\cdot)$ has a poor regularity. This is the case for point-wise Dirac sources. The option we follow comes from duality and the transposition of the variational formulation currently used, when $f (\cdot)$ belongs to $L^{2} (Ω)$ (see, eg, [10,27]). Assume first that $f \in L^{2} (Ω)$ , then the potential φ will be sought for in the Sobolev space $H_{0}^{1} (Ω)$ as the solution to $\begin{matrix} (4) & \int_{Ω} a \nabla φ \cdot \nabla ψ d x = \int_{Ω} f ψ d x, \forall ψ \in H_{0}^{1} (Ω) . \end{matrix}$ This problem is well posed and has only one solution, stable with respect to the data f. More regularity of φ can be expected according to the smoothness of the parameter $a (\cdot)$ . If it belongs to $W^{1, \infty} (Ω)$ then φ is in $H^{2} (Ω)$ (see [3,18]). Morrey’s Theorem shows that, for $d ⩽ 3$ , the solution is continuous (see [2]). As will be seen later, this stability is required for studying the problem with non-regular sources. Though, Morrey’s approach may fail to apply.

The possible lack of smoothness of the conductivity, when it is limited to $a \in L^{\infty} (Ω)$ , is cause of technical difficulty. Estimate (5) can no longer be obtained by means of Sobolev embeddings. It has to be directly proved, a pretty hard result, established in a series of technical papers, by Di-Giorgi first (see [15]), and then by Nash [25], Moser [24] and Stampacchia [27]. The following lemma can be found in [10, Theorem 9.34].

Lemma 2.1.
Assume $d ⩽ 3$ and $a (\cdot)$ satisfies ( 2 ). If $f \in L^{2} (Ω)$ then the solution φ to Problem ( 4 ) is Hölderian and there exists a real number $ν \in] 0, 1 [$ and a constant $C > 0$ such that $\begin{matrix} (5) & ‖ φ ‖_{C^{0, ν} (Ω)} ⩽ C ‖ f ‖_{L^{2} (Ω)} . \end{matrix}$ Both ν and C depend upon the parameter $a (\cdot)$ .

Now, we are well equipped to investigate the singular data f given in (3). The suggested variational framework to equation (1) is based on a duality argument. Indeed, the potential φ fails to be in $H^{1} (Ω)$ , and the space where it should be sought for has no better regularity than the Lebesgue space $L^{2} (Ω)$ . The variational formulation is processed as follows. Let $F (\cdot)$ be given in $L^{2} (Ω)$ , an arbitrary test function. Then, denote by $r_{F} (\cdot)$ the unique solution of the ‘regular formulation’ of Poisson’s problem (4), where $f (\cdot)$ is replaced by $F (\cdot)$ . The transposition method applied to (1) yields the alternative variational equation that consists in: finding φ in $L^{2} (Ω)$ such that $\begin{matrix} (6) & \int_{Ω} φ (x) F (x) d x = \sum_{s \in S} λ_{s} r_{F} (s), \forall F \in L^{2} (Ω) . \end{matrix}$ By Lemma 2.1, $r_{F} (\cdot)$ belongs to $C (\overline{Ω})$ and the duality defined in the right hand side makes sense. The following existence result holds
Lemma 2.2.
Assume $d ⩽ 3$ . Problem ( 6 ) has a unique solution $φ \in L^{2} (Ω)$ with $\begin{matrix} ‖ φ ‖_{L^{2} (Ω)} ⩽ C \sum_{s \in S} | λ_{s} | . \end{matrix}$ The constant C is independent of the cardinality of S, the number of the Dirac point-wise sources.
Proof.
The linear form $F \mapsto ⟨ f, r_{F} ⟩ = \sum_{s \in S} λ_{s} r_{F} (s)$ is continuous on $L^{2} (Ω)$ . Indeed, according to the stability (5) we have that $\begin{matrix} ⟨ f, r_{F} ⟩ ⩽ ‖ r_{F} ‖_{C (Ω)} \sum_{s \in S} | λ_{s} | ⩽ C ‖ F ‖_{L^{2} (Ω)} \sum_{s \in S} | λ_{s} | . \end{matrix}$ The constant C is the one in Lemma 2.1 and is then independent of the cardinality of S. The Riesz representation Theorem results in the existence and uniqueness. The proof is complete. □
Remark 2.1.
Let us consider the Hilbert space $\begin{matrix} V = {ψ \in H_{0}^{1} (Ω); - div (a \nabla ψ) \in L^{2} (Ω)}, \end{matrix}$ equipped with the graph norm $‖ \cdot ‖_{V} = ‖ div (a \nabla (\cdot)) ‖_{L^{2} (Ω)}$ and denote by $V^{}$ the dual space. The operator $A_{a} : ψ \mapsto - div (a \nabla ψ)$ determines an isomorphism from V into $L^{2} (Ω)$ . Its adjoint $A_{a}^{}$ is therefore an isomorphism from $L^{2} (Ω)$ into $V^{}$ . If the data f is given in $L^{2} (Ω)$ then the variational problem (4) solves the equation $A_{a} φ = f$ . Otherwise, when $f \in V^{}$ , the formulation (6) solves the equation $A_{a}^{} φ = f$ (see [18]). The point-wise Dirac source f we deal with here belongs to $V^{}$ .
Remark 2.2.
Owing to Morrey’s inequality (see [2]), we have that $f \in H^{{(- 1)}_{-}} (Ω)$ in $R^{2}$ , and $f \in H^{{(- 3 / 2)}_{-}} (Ω)$ in $R^{3}$ . By elliptic regularity, the solution φ lies in $H^{{(1)}_{-}} (Ω)$ in 2 dimensions, and is expected to belong to $H^{{(1 / 2)}_{-}} (Ω)$ in 3 dimensions (see [18]).
Remark 2.3.
Despite the singularities generated by Dirac sources, still working in some $H^{1}$ -space is possible and prompts authors concerned with to resort to weighted Sobolev spaces. Introducing an appropriate weight reduces the power of the singularities to blow up the integrals at the basis of the variational principle (see, e.g. [4]).

The variational problem (6) brings a functional framework to address Poisson’s equation with singular data. The indisputable mathematical interest of it cannot hide its weakness in the practical ground. When it comes to the numerical simulation of problem (1), formulation (6) can hardly be effective. Elaborating a tractable alternative approach is hence mandatory. The main point to keep in mind is that the new method shall enhance the feasibility and the efficiency of widespread approximation methods such as the finite elements.

Let us begin with a preliminary remark. Assume that $a (\cdot)$ is defined in the whole $R^{d}$ while still fulfilling (2). If $G_{a} ((\cdot), s)$ is the fundamental solution of the operator $A_{a} = - div (a \nabla (\cdot))$ (in $R^{d}$ ), then the potential φ solving (6) can be split into a singular contribution $K_{a} (\cdot)$ , the sum of the individual mono-polar potentials created by each Dirac source $(λ_{s} δ_{s})$ , and a regular correction term χ. This correction accounts for the boundary condition. Then the following decomposition holds $\begin{matrix} (7) & φ (\cdot) = \sum_{s \in S} λ_{s} G_{a} ((\cdot), s) + χ (\cdot) = K_{a} (\cdot) + χ (\cdot) . \end{matrix}$ The correction $χ (\cdot)$ lies in $H^{1} (Ω)$ , satisfies $χ_{| Γ} = - K_{a}$ and is solution to the Laplace equation (4) (with $f = 0$ ). This process is conditioned by the knowledge of $G_{a} ((\cdot), s)$ .

Unfortunately, it is scarcely accessible, because of the space-dependency of the conductivity $a (\cdot)$ . Some contributions have dealt with this expansion, called a splitting or subtraction process (see [7,8,14,30]). In all of them, the conductivity is chosen constant at the vicinity of the sources, which enables the use of the convolutional Green kernel $G (\cdot) = G_{1} (\cdot)$ , the fundamental solution of the Laplace equation(1
¹
Recall the explicit expressions of the Laplace-Green function $\begin{matrix} G (x) = \{\begin{matrix} - \frac{1}{2 π} log | x |, & if d = 2, \\ \frac{1}{ω_{d} (d - 2) | x |^{d - 2}}, & if d ⩾ 3 . \end{matrix} \end{matrix}$ $ω_{d}$ is the surface of the unit sphere in $R^{d}$ ( $ω_{3} = 4 π$ ).

).

An alternative to (7) is to find a way to use the Green kernel $G ((\cdot) - s)$ (instead of $G_{a} ((\cdot), s)$ ), for space-varying conductivities.

The challenge of the current work is then to derive a new expansion similar to ( 7 ), where the singular contribution is expressed in terms of $G ((\cdot) - s)$ while the correction part still belongs to $H^{1} (Ω)$ . This last requirement is the central criterion for our investigation.

Recalling the following valuable statements gives a preliminar idea to the study we undertake. One of these statements is proven in [22, Equation (7.5)]: for $d = 3$ , both $G_{a} ((\cdot), s)$ and $G ((\cdot) - s)$ are positive and there exists a constant $β = β (a, s, d) > 1$ depending on $a (\cdot)$ , $s$ and d, such that $\begin{matrix} 0 < \frac{1}{β} G ((\cdot) - s) ⩽ G_{a} ((\cdot), s) ⩽ β G ((\cdot) - s) . \end{matrix}$ In the other, for $d = 2$ , the bounds are changed. There exists two positive constants $β = β (a, s)$ and $η = η (a, s)$ such that (see [28, Theorem 6.7]) $\begin{matrix} | G_{a} ((\cdot), s) | ⩽ η + β | G ((\cdot) - s) | . \end{matrix}$ The singularity extraction we deal with, if successful, brings incidentally supplementary information on $G_{a} ((\cdot), s)$ . We explore and confront its singular behavior to that of $G ((\cdot) - s)$ . The specific question is to figure out under which assumption (on $a (\cdot)$ ) the function $G ((\cdot) - s)$ is able to capture the ‘full singular behavior’ ( $\notin H_{loc}^{1}$ ) of the fundamental solution $G_{a} ((\cdot), s)$ , in the sense that $G_{a} ((\cdot), s) - λ G ((\cdot) - s)$ belongs to $H_{loc}^{1}$ , for an appropriate number λ.
3. Singularities extraction

Developments intended here requires the continuity of the conductivity $a (\cdot)$ at points $s \in S$ . We set the numbers $a_{s} = a (s)$ . The following splitting will be a remedy to the unfeasibility of (7), $\begin{matrix} (8) & φ (\cdot) = \sum_{s \in S} \frac{λ_{s}}{a_{s}} G ((\cdot) - s) + χ (\cdot) = K_{a} (\cdot) + χ (\cdot) . \end{matrix}$ We need to tell whether the correction χ is solution to any tractable variational problem and if it is exempt from any singular behavior. Inserting this expansion in (1) implies that χ is subject to the Dirichlet condition $χ_{| Γ} = - K_{a}$ , and satisfies (formally so far) the regular variational formulation, $\begin{matrix} (9) & \int_{Ω} a \nabla χ \cdot \nabla ψ d x = - \sum_{s \in S} \frac{λ_{s}}{a_{s}} \int_{Ω} (a - a_{s}) \nabla G ((\cdot) - s) \cdot \nabla ψ d x, \forall ψ \in H_{0}^{1} (Ω) . \end{matrix}$ Regular/singular will qualify henceforth functions that do/do not belong to $H^{1} (Ω)$ . Thereupon, $χ (\cdot)$ is regular and $K_{a} (\cdot)$ is singular.

Additional regularity is needed for the conductivity $a (\cdot)$ , at least at the vicinity of the Dirac sources $s$ ( $\in S$ ), so to secure a mathematical sense to equation (9). Hence, we assume that $\begin{matrix} (10) & x \mapsto \frac{a (x) - a_{s}}{| x - s |^{d - 1}} \in L^{2} (Ω), \forall s \in S . \end{matrix}$

Remark 3.1.
Assumption (10) requires a mild smoothness on $a (\cdot)$ , at the vicinity of the Dirac sources. A sufficient condition on $a (\cdot)$ to meet (10) reads as $\begin{matrix} (11) & sup_{x \in Ω} \frac{| a (x) - a_{s} |}{| x - s |^{τ}} < \infty, \forall s \in S, \end{matrix}$ with $τ > (d / 2 - 1)$ . This is known as the Hölder point-wise condition on S, symbolized by $a \in C^{0, τ} (S)$ . We refer to [5] for an instructive discussion about this type of properties.

The existence of the regular correction χ is the result to start with, before assessing the expansion (8).
Lemma 3.1.
Under assumption ( 10 ), the variational problem ( 9 ) has a unique solution $χ \in H^{1} (Ω)$ .
Proof.
Assumption (10) tells that the function $\begin{matrix} x \mapsto \sum_{s \in S} \frac{λ_{s}}{a_{s}} (a (x) - a_{s}) \nabla G (x - s) \in L^{2} (Ω) . \end{matrix}$ Now, for the Dirichlet boundary condition, notice that neither of the sources is located on Γ. As a consequence $G ((\cdot) - s)$ restricted to the boundary Γ is regular and lies in $H^{1 / 2} (Γ)$ ; so does $K_{a} (\cdot)$ . This is enough to apply the Lax-Milgram Lemma to problem (9). The proof is complete. □

Checking out whether the function φ given in (8) is solution to the variational equation (6) is the next step.
Proposition 3.2.
Let $a (\cdot)$ satisfy ( 10 ) and $χ \in H^{1} (Ω)$ be the solution to ( 9 ). Then, the potential $φ \in L^{2} (Ω)$ in ( 8 ) is the solution to problem ( 6 ).
Proof.
After observing that $K_{a} \in L^{2} (Ω)$ ( $d ⩽ 3$ ), let F be a test function in $L^{2} (Ω)$ . We start from the singular/regular integrals splitting $\begin{matrix} \int_{Ω} φ F d x = \int_{Ω} K_{a} F d x + \int_{Ω} χ F d x . \end{matrix}$ Each term is handled apart. The regular term can be transformed into $\begin{matrix} (12) & \int_{Ω} χ F d x = \int_{Ω} χ (- div (a \nabla r_{F})) d x = \int_{Ω} a \nabla χ \cdot \nabla r_{F} d x - {⟨ (a \partial_{n}) r_{F}, χ ⟩}_{1 / 2, Γ} . \end{matrix}$ Considering the small ball $B_{ϵ} = B (s, ϵ)$ centered at $s$ , the singular integral can be put as $\begin{matrix} \int_{Ω} K_{a} F d x = \sum_{s \in S} \frac{λ_{s}}{a_{s}} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) F d x . \end{matrix}$ Applying Green’s integration, we notice that $\begin{matrix} \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) F d x = & \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) (- div (a \nabla r_{F})) d x . \\ = & \int_{Ω ∖ B_{ϵ}} a \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x - {⟨ (a \partial_{n}) r_{F}, G ((\cdot) - s) ⟩}_{1 / 2, Γ} \\ - {⟨ (a \partial_{n}) r_{F}, G ((\cdot) - s) ⟩}_{1 / 2, \partial B_{ϵ}} . \end{matrix}$ The unit normal $n$ in both boundary integrals is oriented outward $Ω ∖ B_{ϵ}$ , then inward $B_{ϵ}$ . After easy manipulations, we derive that $\begin{matrix} \int_{Ω ∖ B_{ϵ}} a \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x = & \int_{Ω ∖ B_{ϵ}} (a - a_{s}) \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x \\ + a_{s} \int_{Ω ∖ B_{ϵ}} \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x . \end{matrix}$ Another use of Green’s formula on the last integral, after observing that $Δ G ((\cdot) - s) = 0$ away from $s$ , together with ${(r_{F})}_{| Γ} = 0$ , we obtain that $\begin{matrix} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} a \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x = & \int_{Ω} (a - a_{s}) \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x \\ + a_{s} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} \partial_{n} G ((\cdot) - s) r_{F} d τ . \end{matrix}$ Sum up various contributions produces the formula $\begin{matrix} \int_{Ω} φ F d x = & \int_{Ω} a \nabla χ \cdot \nabla r_{F} d x - {⟨ (a \partial_{n}) r_{F}, χ ⟩}_{1 / 2, Γ} - {⟨ (a \partial_{n}) r_{F}, K ⟩}_{1 / 2, Γ} \\ + \sum_{s \in S} \frac{λ_{s}}{a_{s}} \int_{Ω} (a - a_{s}) \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x + \sum_{s \in S} λ_{s} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} \partial_{n} G ((\cdot) - s) r_{F} d τ \\ - \sum_{s \in S} \frac{λ_{s}}{a_{s}} lim_{ϵ \to 0} ⟨ (a \partial_{n}) r_{F}, G ((\cdot) - s) ⟩_{1 / 2, \partial B_{ϵ}} . \end{matrix}$ By construction, $χ (\cdot)$ is solution to (9), hence we get $\begin{matrix} \int_{Ω} φ F d x = \sum_{s \in S} λ_{s} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} \partial_{n} G ((\cdot) - s) r_{F} d τ - \sum_{s \in S} \frac{λ_{s}}{a_{s}} lim_{ϵ \to 0} {⟨ (a \partial_{n}) r_{F}, G ((\cdot) - s) ⟩}_{1 / 2, \partial B_{ϵ}} . \end{matrix}$ The limits are calculated only in three dimensions. They are readily carried out in two dimensions following the same arguments. We have that $\begin{matrix} \int_{Ω} φ F d x = \sum_{s \in S} λ_{s} lim_{ϵ \to 0} \frac{1}{4 π ϵ^{2}} \int_{\partial B_{ϵ}} r_{F} d τ - \sum_{s \in S} \frac{λ_{s}}{a_{s}} lim_{ϵ \to 0} \frac{1}{4 π ϵ} {⟨ (a \partial_{n}) r_{F}, 1 ⟩}_{1 / 2, \partial B_{ϵ}} . \end{matrix}$ The last integral may be handled as follows ( $n$ is the internal normal to $B_{ϵ}$ ) $\begin{matrix} {⟨ (a \partial_{n}) r_{F}, 1 ⟩}_{1 / 2, \partial B_{ϵ}} = - \int_{B_{ϵ}} div (a \nabla r_{F}) d x = \int_{B_{ϵ}} F d x . \end{matrix}$ This yields the bound $\begin{matrix} (13) & \frac{1}{4 π ϵ} | {⟨ (a \partial_{n}) r_{F}, 1 ⟩}_{1 / 2, \partial B_{ϵ}} | ⩽ \frac{1}{4 π ϵ} \sqrt{| B_{ϵ} |} ‖ F ‖_{L^{2} (Ω)} ⩽ \sqrt{ϵ} ‖ F ‖_{L^{2} (Ω)} . \end{matrix}$ The limit is hence zero when ϵ decays towards zero. We turn now to the first integral. For any small ϵ, the relevant change of variables $(τ, d τ) = (s + ϵ \hat{τ}, ϵ^{2} d \hat{τ})$ gives the formula $\begin{matrix} (14) & \frac{1}{4 π ϵ^{2}} \int_{\partial B_{ϵ}} r_{F} d τ = \frac{1}{4 π} \int_{| \hat{τ} | = 1} r_{F} (s + ϵ \hat{τ}) d \hat{τ} . \end{matrix}$ The continuity of $r_{F} (\cdot)$ on Ω allows to pass to the limit so that $\begin{matrix} lim_{ϵ \to 0} \frac{1}{4 π ϵ^{2}} \int_{\partial B_{ϵ}} r_{F} d τ = r_{F} (s) . \end{matrix}$ Thereupon, $φ (\cdot)$ is solution to (6). The proof is complete. □

More results can be established on the correction function $χ (\cdot)$ , owing to the Hölderian smoothness exposed in [27] and to the elliptic regularity stated in [18]. Lemma 3.3.
Let $d = 3$ and ω be strongly included in $Ω ∖ S$ . Under assumption ( 10 ), the correction $χ \in C^{0, ν} (ω)$ for some $ν \in] 0, 1 [$ . Moreover, if $a \in W^{1, \infty} (Ω)$ then $χ \in H^{{(3 / 2)}_{-}} (Ω)$ .
Proof.
Let ω be a subset such that $\overline{ω} \subset Ω ∖ S$ . Then, it is straightforward that, for all $q > d = 3$ , the function $\begin{matrix} (15) & x \mapsto (a (x) - a_{s}) \nabla G ((\cdot) - s) \in L^{q} (ω) . \end{matrix}$ Using Theorem 7.6 in [27] ends to the desired result. In addition, if $a \in W^{1, \infty} (Ω)$ then the function in (15) belongs to $H^{{(1 / 2)}_{-}} (Ω)$ . By elliptic regularity we get $χ \in H^{{(3 / 2)}_{-}} (Ω)$ (see [18]). The proof is complete. □
Lemma 3.4.
For $d = 2$ , if the conductivity $a (\cdot)$ fulfills assumption ( 11 ) with $τ > 0$ , then $χ \in C^{0, ν} (Ω)$ for some $ν \in] 0, 1 [$ .
Proof.
The assumption on $a (\cdot)$ makes the result (15) valid on the whole domain (ω is changed into Ω), for all $q \in] 2, 2 / (1 - τ) [$ . Hence, using again Theorem 7.6 in [27] implies that $χ \in C^{0, ν} (Ω)$ . The Hölderian property is global. The proof is complete. □
Remark 3.2.
A by-product of Proposition 3.2 is that the Green kernel $G_{a} ((\cdot), s)$ in (7) can be split as fillows $\begin{matrix} G_{a} ((\cdot), s) = \frac{1}{a (s)} G ((\cdot) - s) + χ (\cdot) . \end{matrix}$ The correction function $χ (\cdot)$ lies in $H_{loc}^{1} (R^{3})$ .

A significant difference between $G_{a} ((\cdot), s)$ and $G ((\cdot) - s)$ is: the first kernel is bivariate while the second is a univariate convolutional kernel.
Remark 3.3.
The authors in [8,30] mention that if $a \in W^{1, \infty} (Ω)$ and it is constant at the vicinity of the sources $s \in S$ then $χ \in H^{2} (Ω)$ . This comes directly from the fact that $(a (\cdot) - a_{s})$ vanishes at the vicinity of S. The singular behavior of $\nabla G ((\cdot) - s)$ near the source point $s \in S$ is fully cancelled, after multiplication by $(a (\cdot) - a_{s})$ . In our case, $(a (\cdot) - a_{s})$ is zero only at the singular points, but it is still space-varying around S. This downgrades the singular behavior of $\nabla G ((\cdot) - s)$ without entirely killing it.

4. Conductivities with jumps

An interesting class of space varying conductivities includes those having discontinuities around the sources. The jumps take place at some interfaces of the domain. Those interfaces split Ω into a number of subregions ${(Ω_{i})}_{i ⩽ i_{*}}$ , where $a (\cdot)$ is regular enough. An illustration is provided by the simple example depicted in Fig. 1. A mono-polar source is located at the interface of two or more subregions. Studying these configurations is different and the technical complexity is increased. To make the presentation legible, we develop the main mathematical issues for two dimensional domains. This really brings ease to the exposition and the graphical illustrations. Hints will be provided at the end for the generalization to the three dimensional case. Additional assumptions are also introduced, some do not restrict the generality and we will explain how to get rid of others.

We also suppose for a while that the interfaces are straight lines, as in Fig. 1. We will stick with this example with a single source, only for clarity.

Fig. 1.

The source-point $s$ is the vertex of three sectors. The symbols $ω_{i}, 1 ⩽ i ⩽ 3$ , are the measures of different angles. The interfaces are numbered counterclockwise ${(γ_{j})}_{1 ⩽ j ⩽ 3}$ , where the edge $γ_{1}$ , eg, is the interface to $Ω_{1}$ and $Ω_{2}$ .

4.1. Piecewise constant conductivities

To proceed with the removal of the singularity, we explore first piecewise constant parameters. They are denoted $α (\cdot)$ , instead of $a (\cdot)$ , for a reason to appear later on. In each subregion $Ω_{i}$ , the conductivity $α (\cdot)$ is constant, i.e. $α = α_{i} > 0$ , for all $i ⩽ i_{*}$ .

Then, the fundamental solution of $A_{α}$ is a modulation of the Green function, $\begin{matrix} G_{α} (\cdot) = \frac{1}{α_{*}} G ((\cdot) - s) = - \frac{1}{2 π α_{*}} ln | (\cdot) - s | . \end{matrix}$ We have set $α_{*} = \frac{1}{2 π} \sum_{1 ⩽ i ⩽ i_{*}} α_{i} ω_{i}$ .

The function $G_{α} (\cdot)$ is harmonic within any sector $Ω_{i}$ , away from the vertex $s$ . Moreover, along the interfaces, the function $G_{α} (\cdot)$ is continuous, i.e., ${[G_{α}]}_{| γ_{j}} = 0$ for all $j ⩽ j_{*}$ . Its normal derivative vanishes on $γ_{j}^{'} s$ ; this secures the flux conservation, that is ${([α \partial_{n} G_{α}] = [α] \partial_{n} G_{α})}_{| γ_{j}} = 0$ . As a result, $G_{α} (\cdot)$ satisfies in the point-wise sense $\begin{matrix} (16) & A_{α} G_{α} = - div (α \nabla G_{α}) = 0, in R^{2} ∖ {s} . \end{matrix}$ The coefficient $α_{*}$ is chosen so that $G_{α} (\cdot)$ is the fundamental solution of the operator $A_{α}$ .

4.2. Piecewise continuous conductivities

The extraction process has to be updated in the case where the coefficient $a (\cdot)$ is space-varying with jumps. Assume therefore that $a (\cdot)$ is continuous inside the closure of the regions ${(Ω_{i})}_{i ⩽ i_{*}}$ and have jumps across the interfaces ${(γ_{j})}_{j ⩽ j_{*}}$ . The preliminary step is to compute different sectorial limits at the vertex $s$ , $\begin{matrix} α_{i} = lim_{x \in Ω_{i} \to s} a (x), for all i ⩽ i_{*} . \end{matrix}$ Then, we define $α (\cdot)$ , the piecewise constant function, to take the value $α_{i}$ within $Ω_{i}$ . The solution φ to (6) can therefore be split as $\begin{matrix} (17) & φ (\cdot) = λ_{s} G_{α} (\cdot) + χ (\cdot) = K_{α} (\cdot) + χ (\cdot) . \end{matrix}$ The regular correction $χ \in H^{1} (Ω)$ is subject to the Dirichlet condition $χ_{| Γ} = - K_{α}$ , and is solution to $\begin{matrix} (18) & \int_{Ω} a \nabla χ \cdot \nabla ψ d x = - \int_{Ω} (a - α) \nabla K_{α} \cdot \nabla ψ d x, \forall ψ \in H_{0}^{1} (Ω) . \end{matrix}$ The following well posedness result is promptly needed.

Lemma 4.1.
Assume that $a \in W^{τ, \infty} (Ω_{i})$ , for all $i ⩽ i_{}$ , with* $τ > 0$ . The variational problem ( 18 ) has a unique solution $χ \in H^{1} (Ω)$ .
Proof.
Before engaging in the proof, we shall underline the fact that the regularity of $a (\cdot)$ is mild and local. It says that $a (\cdot)$ is continuous in ${\overline{Ω}}_{i}$ for all $i ⩽ i_{}$ , but jumps across ${(γ_{j})}_{j ⩽ j_{}}$ are allowed.

The continuity of the linear form in (9) is the key-point. The clue for the proof is that $a \in W^{τ, \infty} (Ω_{i})$ . Hence, the following bound holds, $\begin{matrix} (19) & | (a (x) - α_{i}) \nabla K_{α} (x) | ⩽ C | λ_{s} | | x - s |^{τ - 1}, \forall x \in Ω_{i} . \end{matrix}$ As a result $x \mapsto (a (x) - a_{s}) \nabla K_{α} (x)$ belongs to $L^{2} (Ω)$ . Lax-Milgram Lemma applies to problem (9). The proof is complete. □

As in Section 3, one has to prove that the function φ, constructed through the splitting (17), is the solution to the dual variational problem (6). Proposition 4.2.
Under the assumptions of Lemma 4.1 , the potential $φ (\cdot)$ in ( 17 ) is the solution to problem ( 6 ).
Proof.
Parts of the technical tools are similar to those used for Proposition 3.2. Only few amendments should be introduced.

Let $F (\cdot)$ be in $L^{2} (Ω)$ . We start as in Proposition 3.2. The regular contribution is identical to (12) while the singular part undergoes some changes. It reads as $\begin{matrix} \int_{Ω} K_{α} F d x = \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) F d x = \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) (- div (a \nabla r_{F})) d x . \end{matrix}$ Performing a by part integration and then plugging $α (\cdot)$ , we obtain $\begin{matrix} \int_{Ω ∖ B_{ϵ}} G ((\cdot) - s) F d x = & \int_{Ω ∖ B_{ϵ}} a \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x - \int_{Γ} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ \\ - \int_{\partial B_{ϵ}} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ \\ = & \int_{Ω ∖ B_{ϵ}} (a - α) \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x + \int_{Ω ∖ B_{ϵ}} α \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x \\ - \int_{Γ} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ - \int_{\partial B_{ϵ}} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ . \end{matrix}$ Using again Green’s formula for the last integral, and in view of the harmonicity of $G ((\cdot) - s)$ away from $s$ together with ${(r_{F})}_{| Γ} = 0$ , we come up with $\begin{matrix} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} α \nabla G ((\cdot) - s) \cdot \nabla r_{F} d x = lim_{ϵ \to 0} \int_{\partial B_{ϵ}} (α \partial_{n}) G ((\cdot) - s) r_{F} d τ . \end{matrix}$ Thereafter, summing up different contributions produces $\begin{matrix} \int_{Ω} φ F d x = & \int_{Ω} a \nabla χ \cdot \nabla r_{F} d x - \int_{Γ} χ (a \partial_{n}) r_{F} d τ - \int_{Γ} K_{α} (a \partial_{n}) r_{F} d τ \\ + \int_{Ω} (a - α) \nabla K_{α} \cdot \nabla r_{F} d x + \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} (α \partial_{n}) G ((\cdot) - s) r_{F} d τ \\ - \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ . \end{matrix}$ The function $χ (\cdot)$ fulfills (18), thus we derive that $\begin{matrix} \int_{Ω} φ F d x & = \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} (α \partial_{n}) G ((\cdot) - s) r_{F} d τ - \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ . \\ = \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \frac{1}{ϵ} \int_{\partial B_{ϵ}} α r_{F} d τ - \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} ln ϵ \int_{\partial B_{ϵ}} (a \partial_{n}) r_{F} d τ . \end{matrix}$ The last integral is delt with as in (13), and the limit is zero when ϵ tends toward zero. We turn then to the first integral. Proceeding by change of variables as in (14), we infer for any small ϵ, $\begin{matrix} \frac{1}{ϵ} \int_{\partial B_{ϵ}} α r_{F} d τ = \int_{| \hat{τ} | = 1} α r_{F} (s + ϵ \hat{τ}) d \hat{τ} . \end{matrix}$ Notice that the function $α (\cdot)$ is constant within each sector delimiting $Ω_{i}$ , at the vicinity of the vertex $s$ ; it is by no means dependent on ϵ. This yields in particular that $α (s + ϵ \hat{τ}) = α_{i}$ , within the arc ${s + ϵ \hat{τ}, | \hat{τ} | = 1} \cap Ω_{i}$ . The smoothness of $r_{F} (\cdot)$ allows to pass to the limit so that $\begin{matrix} \frac{1}{α_{}} lim_{ϵ \to 0} \frac{1}{ϵ} \int_{\partial B_{ϵ}} α r_{F} d τ & = \frac{1}{α_{}} lim_{ϵ \to 0} \sum_{1 ⩽ i ⩽ i_{}} α_{i} \int_{{s + ϵ \hat{τ}, | \hat{τ} | = 1} \cap Ω_{i}} r_{F} (s + ϵ \hat{τ}) d \hat{τ} \\ = \frac{1}{α_{}} r_{F} (s) \sum_{1 ⩽ i ⩽ i_{}} α_{i} ω_{i} = r_{F} (s) . \end{matrix}$ The proof is complete. □
Remark 4.1.
A corollary of Proposition 4.2 is that the Green function $G_{a} ((\cdot), s)$ can be split as $\begin{matrix} G_{a} ((\cdot), s) = G_{α} (\cdot) + χ (\cdot) = \frac{1}{α_{}} G ((\cdot) - s) + χ (\cdot) . \end{matrix}$ The correction function $χ \in H_{loc}^{1} (R^{2})$ . Moreover, proceeding as in Lemma 3.4, the bound (19) shows that $(a (\cdot) - α (\cdot)) \nabla K_{α} (\cdot)$ belongs actually to $L_{loc}^{q} (R^{2})$ for some $q > 2$ , it follows that $χ (\cdot)$ is Hölderian, since it belongs to $C_{loc}^{0, ν} (R^{2})$ . The Laplace-Green function captures the whole singular behavior of $G_{a} ((\cdot), s)$ .

4.3. Localized singularity extraction

We aim to control the process of the singularities extraction so to restrict it to a small region around the point-source $s$ . The truncation consists in multiplying the Green kernel by a cut-off (function) and delimit its effects to the proximity of $s$ . This induces a less complexity in the practical (numerical) treatment.

The direct consequence of the localization is that the corrector $χ (\cdot)$ , determined by (18), must be revised. To do so, let $ϱ > 0$ be fixed and consider a regular non-negative function $ξ_{ϱ} \in C_{c}^{1} ([0, \infty [)$ , supported in $[0, 2 ϱ]$ and is equal to unity on $[0, ϱ]$ . We have $\begin{matrix} 0 ⩽ ξ_{ϱ} ⩽ 1, {ξ_{ϱ}}_{| [0, ϱ]} = 1, {ξ_{ϱ}}_{| [2 ϱ, \infty [} = 0 . \end{matrix}$ We set $ξ_{s} = ξ_{ϱ} (| (\cdot) - s |)$ . The fact that the cut-off function is radial is critical. Then, expansion (17) is transformed into $\begin{matrix} (20) & φ (\cdot) = λ_{s} G_{α} (\cdot) ξ_{s} (\cdot) + χ (\cdot) = K_{α} (\cdot) + χ (\cdot) . \end{matrix}$ In case of multiple sources, the localization is intended to prevent them from interacting with each other. The cut-off functions with appropriate shrinked supports ensure their separation. Likewise, selecting ϱ lower than the distance of $s$ to the boundary excludes interactions between the singularities of the solution with Γ.

The regular correction $χ \in H^{1} (Ω)$ is finally subject to the homogeneous Dirichlet condition $χ_{| Γ} = 0$ , and satisfies the variational problem $\begin{array}{l} \int_{Ω} a \nabla χ \cdot \nabla ψ d x = - \int_{Ω} (a - α) \nabla K_{α} \cdot \nabla ψ d x + \int_{Ω} α \nabla ξ_{s} \cdot (\nabla K_{α} ψ - K_{α} \nabla ψ) d x, \\ (21) & \forall ψ \in H_{0}^{1} (Ω) . \end{array}$ The extra integral term (with respect to (18)) in the linear form is well defined. In fact, $ξ_{ϱ} (\cdot)$ being constant close to zero, the gradient $\nabla ξ_{s} (\cdot)$ vanishes in the ball $B (s, ϱ)$ . As a result, the singularity contained in $G_{α} (\cdot)$ is fully canceled. Problem (23) has thus a unique solution in $H^{1} (Ω)$ . Moreover, we have the following result,

Proposition 4.3.
Assume that $a \in W^{τ, \infty} (Ω_{i})$ , for all $i ⩽ i_{}$ , with* $τ > 0$ . The potential φ in ( 20 ), with χ that solves ( 21 ), is the solution to problem ( 6 ).
Proof.
Handling of the singular part in problem (6) needs to be revisited because of the action of the cut-off function. We have to start from (after observing that $χ_{| Γ} = 0$ ) $\begin{matrix} \int_{Ω} (K_{α} + χ) F d x = λ_{s} lim_{ϵ \to 0} \int_{Ω ∖ B_{ϵ}} ξ_{s} G_{α} (- div (a \nabla r_{F})) d x + \int_{Ω} a \nabla χ \cdot \nabla r_{F} d x . \end{matrix}$ Green’s formula provides $\begin{matrix} \int_{Ω ∖ B_{ϵ}} ξ_{s} G_{α} F d x = \int_{Ω ∖ B_{ϵ}} a \nabla (ξ_{s} G_{α}) \cdot \nabla r_{F} d x - \int_{\partial B_{ϵ}} G_{α} (a \partial_{n}) r_{F} d τ . \end{matrix}$ The integral term on the natural boundary vanishes, because of the localizing effect, i.e., $ξ_{s} = 0$ on Γ, while $ξ_{s} = 1$ on $\partial B_{ϵ}$ , for ϵ small enough ( $⩽ ϱ$ , the cut-off parameter). Replicating the process with the insertion of the piecewise function $α (\cdot)$ we are left with the modified term $\begin{matrix} \int_{Ω ∖ B_{ϵ}} α \nabla (ξ_{s} G_{α}) \cdot \nabla r_{F} d x = \int_{Ω ∖ B_{ϵ}} α G_{α} \nabla ξ_{s} \cdot \nabla r_{F} d x + \int_{Ω ∖ B_{ϵ}} α ξ_{s} \nabla G_{α} \cdot \nabla r_{F} d x . \end{matrix}$ Operating again with Green’s formula for the last integral, we derive that $\begin{matrix} \int_{Ω ∖ B_{ϵ}} α ξ_{s} \nabla G_{α} \cdot \nabla r_{F} d x & = \int_{Ω ∖ B_{ϵ}} - div (α ξ_{s} \nabla G_{α}) r_{F} d x + \int_{\partial B_{ϵ}} (α \partial_{n}) G r_{F} d τ \\ = \int_{Ω ∖ B_{ϵ}} - α \nabla ξ_{s} \cdot \nabla G_{α} r_{F} d x + \int_{\partial B_{ϵ}} (α \partial_{n}) G_{α} r_{F} d τ . \end{matrix}$ Observe that the integral on $- ξ_{s} div (α \nabla G_{α}) = 0$ , away from $s$ . This is because the $γ_{j}$ ’s, the jump lines of $α (\cdot)$ , are radial. Putting together different terms yields $\begin{matrix} \int_{Ω} φ F d x = & \int_{Ω} a \nabla χ \cdot \nabla r_{F} d x + \int_{Ω} (a - α) \nabla K_{α} \cdot \nabla r_{F} d x \\ - \int_{Ω ∖ B_{ϵ}} α \nabla ξ_{s} \cdot (\nabla K_{α} r_{F} - K_{α} \nabla r_{F}) d x + λ_{s} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} (α \partial_{n}) G_{α} r_{F} d τ \\ - λ_{s} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} G_{α} (a \partial_{n}) r_{F} d τ . \end{matrix}$ Accounting for (21), provides that $\begin{matrix} \int_{Ω} φ F d x = \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} (α \partial_{n}) G ((\cdot) - s) r_{F} d τ - \frac{λ_{s}}{α_{}} lim_{ϵ \to 0} \int_{\partial B_{ϵ}} G ((\cdot) - s) (a \partial_{n}) r_{F} d τ . \end{matrix}$ The proof can be achieved as in Proposition 4.2. □
Remark 4.2.
The equivalent configuration in three dimensions is when the subdomains ${(Ω_{i})}_{i ⩽ i_{}}$ are conical sectors with the source location $s$ as the common apex. The internal boundary $γ_{i} = \partial Ω_{i} \cap Ω$ (of $Ω_{i}$ ) is hence a developable conical nappe. Each line contained in $γ_{i}$ passing par the apex $s$ is a generatrix. The union of all the generatrices gives the nappe $γ_{i}$ . Moreover, the directional vector of any generatrix is a radius vector orthogonal to the normal vector $n_{i}$ to $γ_{i}$ . As a consequence, the normal derivative of any radial function is zero. Such is the case for the Green function $G ((\cdot) - s)$ . Finally, the symbol $ω_{i}$ stands for the solid angle of $Ω_{i}$ at the vertex $s$ .

The fundamental solution is a modulation of the three dimensional Green function $\begin{matrix} G_{α} (\cdot) = \frac{1}{α_{}} G ((\cdot) - s) = \frac{1}{4 π α_{}} \frac{1}{| (\cdot) - s |}, \end{matrix}$ where $α_{} = \frac{1}{4 π} \sum_{1 ⩽ i ⩽ i_{}} α_{i} ω_{i}$ . For reasons similar to those invoqued in two dimensions, the function $G_{α}$ is such that $\begin{matrix} - div (α \nabla G_{α}) = δ_{s}, in R^{3} . \end{matrix}$ Lemma 4.1 and Proposition 4.2, stated in two dimensions, extend as well to the three dimensions for a more regular conductivites that is $a \in W^{τ, \infty} (Ω_{i})$ , for all $i ⩽ i_{}$ , with $τ > 1 / 2$ . The proofs still operate following exactly the same lines.

5. Curved interfaces

Removing the assumption on the linear interfaces where the jumps occur on $a (\cdot)$ is the subject. Instead of straight-lines, they are parameterized curves, that is $\begin{matrix} δ_{j} (t) = (ζ_{j} (t), η_{j} (t)), \forall t \in I_{j} = [0, {(t_{j})}_{*}], δ_{j} (0) = s . \end{matrix}$ The vertex $s$ is therefore a point of concurrency of all the interfaces. Each curve $δ_{j} (\cdot)$ is smooth enough so that it has a tangent line at the vertex $s$ , denoted ${(γ_{j})}_{j ⩽ j_{*}}$ . We suppose that $s$ is a regular point to each curve, that is $δ_{j}^{'} (0) \neq 0$ for all $j ⩽ j_{*}$ . It is the direction vector of the tangent. The configuration that will be focussed on here is when the tangent vectors at point $s$ are all different. We refer to the domain Ω on the left of Fig. 2. All these assumptions are only to avoid tedious technicalities. Different configurations can be handled after necessary adaptations are added. To be complete with the geometrical features, if $Ω_{i}$ is delimited by the curves $δ_{j}$ and $δ_{j^{'}}$ , we denote $Ω_{i}^{'}$ the sector defined by the tangent lines $γ_{j}$ and $γ_{j^{'}}$ . The measure of internal angle of that sector is then $ω_{i}$ (see Fig. 2).

Fig. 2.

Dashed lines are the tangent lines to the curved interfaces $δ_{j}$ ’s, at the vertex $s$ .

The geometry being fixed, we turn now to more functional analysis tasks. We try to exhibit the singularity around the Dirac point source $s$ . Starting from the piecewise continuity of $a (\cdot)$ , we compute the limit $α_{i} = {lim}_{x \in Ω_{i} \to s} a (x)$ , for all $i ⩽ i_{*}$ . Next, we define the piecewise constant $α (\cdot)$ that equal $α_{i}$ in the region $Ω_{i}^{'}$ (rather than $Ω_{i}$ ). We therefore suggest to analyze the following splitting of the solution to (4) $\begin{matrix} (22) & φ (\cdot) = λ_{s} G_{α} (\cdot) ξ_{s} (\cdot) + χ (\cdot) = K_{α} (\cdot) + χ (\cdot) . \end{matrix}$ Here, the cut-off function is used as explained in Section 4.3. The correction $χ \in H^{1} (Ω)$ fulfills the condition $χ_{| Γ} = 0$ , and is solution to the variational problem $\begin{array}{l} \int_{Ω} a \nabla χ \cdot \nabla ψ d x = - \int_{Ω} (a - α) \nabla K_{α} \cdot \nabla ψ d x + \int_{Ω} α \nabla ξ_{s} \cdot (\nabla K_{α} ψ - K_{α} \nabla ψ) d x, \\ (23) & \forall ψ \in H_{0}^{1} (Ω) . \end{array}$ The main fact to point out is that $a (\cdot)$ is discontinuous at the curved interfaces, the $δ_{j}$ ’s, while $α (\cdot)$ jumps at the lines $γ_{j}$ ’s.

Theorem 5.1.

Assume that $a_{| Ω_{i}} \in W^{τ, \infty} (Ω_{i})$ , for all $i ⩽ i_{*}$ , and $δ_{j}^{'} \in W^{τ, \infty} (I_{j})$ , for all $j ⩽ j_{*}$ , with $τ > 0$ . The potential φ in ( 22 ), with χ that solves ( 23 ), is the solution to problem ( 6 ).

Proof.

The poof is similar to Proposition 4.3. The only point we need to assess is the validity of assumption (10), and therefore to show that $\begin{matrix} x \mapsto (a (x) - α (x)) \nabla K_{α} (x) = \frac{λ_{s}}{α_{*}} (a (x) - α (x)) \nabla G (x - s) \in L^{2} (Ω) . \end{matrix}$ This is carried out locally in each curved sector $Ω_{i}$ .

Let us focus on $Ω_{1}$ , represented in Fig. 3. Recall that $Ω_{1}^{'}$ is the straight sector defined by both tangent lines $γ_{3}$ and $γ_{1}$ . For the sake of simplification, we use a local coordinates system $(t, s)$ around the vertex $s$ , where the t-axis is the tangent line $γ_{3}$ . At least at the vicinity of $s$ , the new parameterization of the boundaries of $Ω_{1}$ is therefore $δ_{3} = (ζ_{3} (t) = t, η_{3} (t))$ , with $η_{3}^{'} (0) = 0$ and $δ_{1} = (ζ_{1} (t), η_{1} (t))$ , with $η_{1}^{'} (0) \neq 0$ . We begin by observing that as in (19), the following result holds $\begin{matrix} x \mapsto | (a (x) - α_{1}) \nabla G (x - s) | = \frac{| a (x) - α_{1} |}{| x - s |} \in L^{2} (Ω_{1} \cap Ω_{1}^{'}) . \end{matrix}$ It remains to check out the same integrability in $Ω_{1}^{″} = Ω_{1} ∖ Ω_{1}^{'}$ , the dashed region in Fig. 3 delimited by $δ_{3}$ and $γ_{3}$ . We start by $\begin{matrix} \int_{Ω_{1}^{″}} {(a (x) - α (x))}^{2} | \nabla G (x - s) |^{2} d t d s & ⩽ ‖ a ‖_{L^{\infty}}^{2} \int_{0}^{r_{*}} \int_{η_{3} (t)}^{0} \frac{1}{t^{2} + s^{2}} d s d t \\ = ‖ a ‖_{L^{\infty}}^{2} \int_{0}^{r_{*}} \frac{1}{t} arctan (\frac{| η_{3} (t) |}{t}) d t \\ ⩽ ‖ a ‖_{L^{\infty}}^{2} \int_{0}^{r_{*}} \frac{| η_{3} (t) |}{t^{2}} d t . \end{matrix}$ Here, the radius $r_{*}$ is related to the support of the cut-off function $ξ_{s} (\cdot)$ . Be aware also that we did not try to benefit from the fact that $α (x) = α_{1}$ , because a different configuration is possible. Think for example of the curve $δ_{3} (\cdot)$ located above the tangent line $γ_{3} (\cdot)$ , the function $α (\cdot)$ will not be equal to $α_{1}$ .

Fig. 3.

Curved sector of $Ω_{i}$ and local coordinates.

Accounting for the facts that $η_{3}^{'} (0) = 0$ and $η_{3}^{'} (\cdot) \in W^{τ, \infty} (I_{3})$ , with $τ > 0$ , we derive that $| η_{1} (t) | ⩽ t^{1 + τ}$ . Back to the integral $\begin{matrix} \int_{Ω_{1}^{″}} {(a (x) - α (x))}^{2} | \nabla G (x - s) |^{2} d t d s ⩽ C \int_{0}^{r_{*}} \frac{1}{t^{1 - τ}} d t < \infty . \end{matrix}$ This completes the proof. □

Remark 5.1.

Once again, the smoothness requirement is on the behavior of $a (\cdot)$ at the vicinity of the source point $s$ . So to make the assumption more accurate, we only need that $\begin{matrix} | a_{| Ω_{i}} (x) - a_{| Ω_{i}} (s) | = | a_{| Ω_{i}} (x) - α_{i} | ⩽ C | x - s |^{τ}, \forall x \in Ω_{i} . \end{matrix}$ A little more than the continuity around $s$ is demanded.

The same comment holds for the tangent to the interface $δ_{j}$ . The necessary condition is then $\begin{matrix} | δ_{j}^{'} (t) - δ_{j}^{'} (0) | ⩽ C t^{τ}, \forall t \in I_{j} . \end{matrix}$ We need thus a little more than the tangency of $δ_{j}$ to the line $γ_{j}$ .

Remark 5.2.

An alternative proof is to straighten the curved interfaces by using an appropriate coordinates system. The conductivity is therefore changed. The new one is dependent not only on the old one $a (\cdot)$ but also on the derivatives of the curves parameterizations ${(δ_{j}^{'})}_{j ⩽ j_{*}}$ . This explains why it is necessary for $δ_{j}^{'} (\cdot)$ to belong to $W^{τ, \infty} (I_{j})$ , for all $j ⩽ j_{*}$ .

6. Analytical examples

Some analytical examples are exposed so to assess the validity of Proposition 3.2, and check out how essential the required smoothness on the conductivity can be to secure the singularity extraction. The geometrical domain is invariably the unit ball in $R^{d}$ , $d = 2, 3$ , and the conductivities are all radial, for the feasibility of the calculations. For the sake of simplicity, we limit the investigation to a single Dirac source with a unit mass, centered at the origin.

6.1. Exact solutions

The symmetry context suggests that the solution is radial and depends only on $r = | x |$ , that is $φ (x) = φ (r)$ . Solving, in spherical coordinates, the Laplace equation away from the singularity yields that $\begin{matrix} (24) & φ (x) = φ (r) = - \int_{r}^{1} \frac{β}{t^{d - 1} a (t)} d t, r \in] 0, 1] . \end{matrix}$ The constant β is given by $β = - \frac{1}{2 π}$ if $d = 2$ and equal to $β = \frac{1}{4 π}$ when $d = 3$ . Notice that if $a = 1$ , we get the expected following Green functions $\begin{matrix} φ (x) = - \frac{1}{2 π} log | x |, if d = 2, φ (x) = \frac{1}{4 π} (\frac{1}{| x |} - 1), if d = 3 . \end{matrix}$ In the range of examples we discuss, the conductivities are all highly smooth away from the location of the Dirac source. According to the expectations, the solutions are therefore analytic except at the origin, where a singularity arises. The nature of that singular behavior is related to the conductivity.

Let α, ϵ be a couple of real numbers such that $0 < ϵ < α$ . In the first example, the conductivity is $\begin{matrix} \frac{1}{a (x)} = α + ϵ | x |^{λ}, | x | ⩽ 1 . \end{matrix}$ If the exponent $λ > 0$ , then assumption (2) is valid. Calculating the integral (24) in $R^{2}$ yields $\begin{matrix} (25) & φ (x) = \frac{1}{2 π} \int_{| x |}^{1} \frac{α + ϵ t^{λ}}{t} d t = - \frac{α}{2 π} ln | x | + \frac{ϵ}{2 π λ} (1 - | x |^{λ}) = \frac{1}{a (0)} G (x) + χ (x) . \end{matrix}$ The function $χ (\cdot)$ is continuous (see Lemma 3.4) and belongs to $H^{1} (Ω)$ . The same solution in $R^{3}$ reads as $\begin{matrix} φ (x) & = - \frac{1}{4 π} \int_{| x |}^{1} \frac{α + ϵ t^{λ}}{t^{2}} d t = \frac{α}{4 π} (\frac{1}{| x |} - 1) - \frac{ϵ}{4 π (λ - 1)} (\frac{1}{| x |^{1 - λ}} - 1) \\ = \frac{1}{a (0)} G (x) + χ (x) . \end{matrix}$ For the correction function $χ (\cdot)$ to belong to $H^{1} (Ω)$ , the parameter λ needs to be larger than $1 / 2$ , in agreement with assumption (10) and Lemma 3.1. Indeed, only then the conductivity $a (\cdot)$ satisfies the Hölder condition at the origin, with the exponent $λ > 1 / 2$ . Otherwise, if $λ ⩽ 1 / 2$ , the correction function χ lies in $H^{{(1 / 2 + λ)}_{-}} (Ω)$ . Lemma 3.1 cannot be substantially improved.

The next examples help figuring out how the limit failure at the origin of the conductivity are impacted on the solution. Assume hence that $\begin{matrix} a (x) = \frac{1}{α + ϵ sin | x |^{- 1}}, | x | ⩽ 1 . \end{matrix}$ This parameter fulfills (2). However, it is highly oscillatory at $x = 0$ , and the limit does not exist. This case does not frame with any of our assumptions.

Computations in the $R^{3}$ turn out to be more straightforward, we start from there. Applying in (24), the variable change $ρ = r^{- 1}$ , entails $\begin{matrix} φ (x) & = \frac{1}{4 π} \int_{1}^{\frac{1}{| x |}} (α + ϵ sin (ρ)) d ρ = \frac{α}{4 π | x |} - \frac{α}{4 π} - \frac{ϵ}{4 π} (cos \frac{1}{| x |} - cos 1) \\ = \frac{1}{4 π | x |} (α - ϵ | x | cos \frac{1}{| x |})) - \frac{1}{4 π} (α - ϵ cos 1) = η (x) G (x) + χ (x) . \end{matrix}$ The same variable change $ρ = r^{- 1}$ in $R^{2}$ results in $\begin{matrix} φ (x) & = - \frac{1}{2 π} \int_{1}^{\frac{1}{| x |}} \frac{1}{ρ} (α + ϵ sin (ρ)) d ρ = - \frac{α}{2 π} ln | x | - \frac{ϵ}{2 π} (Si (\frac{1}{| x |}) - Si (1) \\ = - \frac{1}{2 π} ln | x | (α + \frac{ϵ}{ln | x |} Si (\frac{1}{| x |})) + \frac{ϵ}{2 π} Si (1)) \\ = η (x) G (x) + χ (x) . \end{matrix}$ The symbol Si $(\cdot)$ is for the sine integral function (see [1]).

In both splittings the modulation function $η (\cdot)$ is continuous at the source $s = 0$ . Indeed, we have ${lim}_{x \to 0} η (x) = α$ . Connecting this limit with the behavior of the conductivity at the source location is not possible. Indeed, because of the oscillations, $a (\cdot)$ has no limit at $s = 0$ . Besides, be aware that $(η (\cdot) - α) G (\cdot)$ does not belong to $H^{1} (Ω)$ , in either cases. The results stated in Proposition 3.2 are not valid here.

A different behavior surfaces when the singular part is not determined. This can be observed in three dimensions when the conductivity is $\begin{matrix} a (x) = \frac{1}{α + ϵ sin (ln | x |)}, | x | ⩽ 1 . \end{matrix}$ Operating the variable changing $r = ln (t)$ ends to $\begin{matrix} φ (x) & = \frac{α}{4 π | x |} - \frac{α}{4 π} + \frac{ϵ}{4 π} \int_{| x |}^{1} \frac{1}{t^{2}} sin (ln (t)) d t = \frac{α}{4 π | x |} - \frac{α}{4 π} + \frac{ϵ}{4 π} \int_{ln | x |}^{0} e^{- r} sin r d r \\ = [α + \frac{ϵ}{\sqrt{2}} sin (ln | x | + \frac{π}{4})] G (x) - \frac{1}{4 π} (α + \frac{ϵ}{2}) = η (x) G (x) + χ (x) . \end{matrix}$ The singular contribution to the expansion has no determination, because of the oscillations. No limit exists on the multiplier $η (\cdot)$ . However, we have that $η \in H^{{(3 / 2)}_{-}} (Ω)$ which a better regularity than for the previous three dimensional example.

6.2. Finite element solutions

We realize some basic simulations using the finite element method to solve problem (6) through formula (8). The aim is to get a first insight concerning the only uncertainty related to the right hand side (because of its particular form). The computations are realized in Freefem $_{+ +}$ , the finite element code developed by F. Hecht (see [20]).

The two-dimensional problem whose solution is given in (25) is the one we deal with. In fact, only the regular corrector $χ (\cdot)$ is simulated using linear finite elements. The computed counterpart $χ_{h} (\cdot)$ is therefore solution to the finite element discretization of problem (9). Notice that the ‘discrete’ potential $φ_{h} (\cdot)$ is obtained after adding the exact singular contribution $K (\cdot) (= α G (\cdot))$ to $χ_{h} (\cdot)$ . This induces the identity $(φ - φ_{h}) = (χ - χ_{h})$ . The gap $(χ - χ_{h})$ is therefore evaluated in the $H^{1}$ – and $L^{2}$ –norms and the convergence curves are depicted in Fig. 4.

Fig. 4.

Convergence curves for $χ (\cdot)$ (then for $φ (\cdot)$ ), with $λ = 0.1$ (left) and $λ = 0.5$ (right).

That the convergence rate is dependent on the smoothness of the solution to approximate ( $χ (\cdot)$ here) is widely known. Looking closely at the corrector $χ (\cdot)$ shows that it belongs to $H^{{(1 + λ)}_{-}}$ . The $H^{1}$ -error is then expected to decay like $C h^{λ_{-}}$ (see [12]). The slopes of the linear regressions of the convergence curves are evaluated to $(0.19, 0.52)$ for the exponents $λ = (0.1, 0.5)$ . The observations are quite close to the predictions.

To step further, we run computations where the potential $φ (\cdot)$ is crudely approximated by the finite element solution $φ_{h, F} (\cdot)$ to problem (4) – the subscript $_{F}$ is for Freefem. An automatic mesh generation, a metric-based anisotropic mesh adaptation, is activated. The refining process relies on the computation of the Hessian of $φ_{h, F} (\cdot)$ . The corrector $χ (\cdot)$ in expansion (8) is approximated by $χ_{h} (\cdot)$ with the same adapted mesh, and as described above $φ_{h} (\cdot)$ is reconstructed by adding the singular contibution $K (\cdot)$ . Table 1 provides the variation of the gaps $e_{h, F} = (φ - φ_{h, F})$ and $e_{h} = (φ - φ_{h})$ , during the adaptation procedure. The errors are computed in $(L^{2}, H^{1})$ -norms. Notice that, the solution φ is not in $H^{1}$ , the error is therefore evaluated by considering the interpolation of φ instead. The computational results show a clear improvement brought about by the special extraction process of the singularity. In particular, the $H^{1}$ -norm of the errors (last row of both tables, for instance) reveal a noticeable difference between both approaches in favor of the singularity extraction.

Table 1

The $(L^{2}, H^{1})$ -errors

$N v$	$e_{h, F} \times 10$	$e_{h} \times 10$
100	$(0.25, 2.1)$	$(0.068, 0.43)$
536	$(0.027, 1.1)$	$(0.0093, 0.2)$
1755	$(0.031, 0.56)$	$(0.009, 0.06)$

$N v$	$e_{h, F} \times 10$	$e_{h} \times 10$
100	$(0.22, 2.0)$	$(0.02, 0.13)$
525	$(0.036, 2.7)$	$(0.0048, 0.035)$
1343	$(0.028, 0.52)$	$(0.005, 0.008)$

Left: $λ = 0.1$ . Right: $λ = 0.5$ . $N v$ is the number of mesh vertices.

7. Conclusion and forthcoming work

Splitting the solution to elliptic problems (1) into a singular/regular contributions, for space varying conductivities, is useful in the practical grounds, not to mention its mathematical interest. The singular contribution is calculated by an explicit formula while the computation of the regular correction reduces to solve the same elliptic problem with regular data. The analysis carried out here lays a path to multiple extensions. One is when the source is supported by a curve in $R^{3}$ . Preliminary developments can be found in [17]. Our belief is that our approach is liable to improve the results proven there. This work is currently under way and some substantial progress is already obtained. Another question to explore is for higher dimensions $d ⩾ 4$ . The hilbertian framework is not valid anymore and the Banach spaces $L^{p}$ and $W^{1, p}$ , for $p < 2$ , are used instead of $L^{2}$ and $H^{1}$ . A different direction to investigate is involved in the time dependent diffusion equation, where the results by J. B. Nash [25] play an important role. These researches are planned in the coming future, in both theoretical and numerical folders.

Footnotes

Acknowledgements

To the anonymous reviewers, we would like to say “Un GRAND MERCI!”. Your insightful comments and suggestions have had an important impact on the readability of the paper. We are indebted to you for your valuable time.

The authors are thankful to Prof. Henda El Fekih for many fruitful discussions. Working with her is always a pleasure.

Research by Eya Bejaoui is made during her master’s degree and is supported by a SFRI fellowship from UTC-Sorbonne Université (April–August 2022) – SFRI = Structuration de la Formation par la Recherche dans les Initiatives d’excellence.

References

Abramowitz and

I.A.

Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1964.

R.A.

Adams and

J.J.

Fournier, Sobolev Spaces, Pure and Applied Mathematics, Vol. 140, Academic Press, New-York, London, 2003.

Agmon,

Douglis and

Nirenberg, Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions. I, Comm. Pure Appl. Math. 12 (1959), 623–727. doi:10.1002/cpa.3160120405.

J.P.

Agnelli,

E.M.

Garau and

Pedro, A posteriori error estimates for elliptic problems with Dirac measure terms in weighted spaces, ESAIM M 2AN(48) (2014), 1557–1581. doi:10.1051/m2an/2014010.

Andersson, Characterization of pointwise Hölder regularity, Applied and Computational Harmonic Analysis 4 (1997), 429–443. doi:10.1006/acha.1997.0219.

Azizollahi,

Darbas,

Diallo,

A.E.

Badia and

Lohrengel, EEG in neo-nates: Forward modeling and sensitivity analysis with respect to variations of the conductivity, Mathematical Biosciences and Engineering 15 (2018), 905–932. doi:10.3934/mbe.2018041.

Baratchart,

Ben Abda,

Ben Hassen and

Leblond, Recovery of pointwise sources or small inclusions in 2d domains and rational approximation, Inverse Problems 21 (2004), 51–74. doi:10.1088/0266-5611/21/1/005.

Beltrachini, The analytical subtraction approach for solving the forward problem in eeg, Journal of neural engineering 16 (2019), 056029. doi:10.1088/1741-2552/ab2694.

Bertoluzza,

Decoene,

Lacouture and

Martin, Local error estimates of the finite element method for an elliptic problem with a Dirac source term, Numerical Methods for Partial Differential Equations 34 (2017), 97–120. doi:10.1002/num.22186.

10.

Brezis, Functional Analysis, Sobolev Spaces and Partial Differential Equations, Universitext, Springer, New York, 2010.

11.

Casas, Control of an elliptic problem with pointwise state constraints, SIAM J. Control 24 (1986), 1309–1318. doi:10.1137/0324078.

12.

P.G.

Ciarlet, Basic error estimates for elliptic problems, in: Handbook of Numerical Analysis, Vol. II,

P.-G.

Ciarlet and

J.-L.

Lions, eds, Handb. Numer. Anal., II, North-Holland, Amsterdam, 1991, pp. 17–351.

13.

Darbas,

M.M.

Diallo,

El Badia and

Lohrengel, An inverse dipole source problem in inhomogeneous media: Application to the eeg source localization in neonates, Journal of Inverse and Ill-posed Problems 27 (2019), 255–281. doi:10.1515/jiip-2017-0120.

14.

Darbas and

Lohrengel, Review on mathematical modelling of electroencephalography (EEG), Jahresbericht der Deutschen Mathematiker-Vereinigung 121 (2019), 3–39. doi:10.1365/s13291-018-0183-z.

15.

De Giorgi, Sulla differenziabilitá e l’analiticitá delle estremali degli integrali multipli regolari. (Italian), Mem. Accad. Sci. Torino 3 (1957), 25–43.

16.

El Badia and

Ha-Duong, An inverse source problem in potential analysis, Inverse Problems 16 (2000), 651–663. doi:10.1088/0266-5611/16/3/308.

17.

I.G.

Gjerde,

Kumar,

J.M.

Nordbotten and

Wohlmuth, Splitting method for elliptic equations with line sources, ESAIM M 2AN(53) (2019), 1715–1739. doi:10.1051/m2an/2019027.

18.

Grisvard, Elliptic problems in nonsmooth domains, Society for Industrial and Applied Mathematics (2011).

19.

Hämäläinen,

Hari,

R.J.

Ilmoniemi,

Knuutila and

O.V.

Lounasmaa, Magnetoencephalography–theory, instrumentation, and applications to noninvasive studies of the working human brain, Rev. Mod. Phys. 65 (1993), 413–497. doi:10.1103/RevModPhys.65.413.

20.

Hecht, New development in freefem++, J. Numer. Math. 20 (2012), 251–265.

21.

J.D.

Jackson, Classical Electrodynamics, John Wiley & Sons, Inc., New York–London–Sydney, 1975.

22.

Littman,

Stampacchia and

Weinberger, Regular points for elliptic equations with discontinuous coefficients, Annali Della Scuola Normale Superiore Di Pisa-classe Di Scienze 1–2 (1963), 43–77.

23.

Malmivuo and

Plonsey, Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields, Oxford University Press, New York, Oxford, 1995.

24.

Moser, A new proof of de Giorgi’s theorem concerning the regularity problem for elliptic differential equations, Comm. Pure Appl. Math. 13 (1960), 457–468. doi:10.1002/cpa.3160130308.

25.

Nash, Continuity of solutions of parabolic and elliptic equations, American Journal of Mathematics 80 (1958), 931–954. doi:10.2307/2372841.

26.

Rannacher and

Vexler, A priori error estimates for the finite element discretization of elliptic parameter identification problems with pointwise measurements, SIAM J. Control Optim. 44 (2005), 1844–1863. doi:10.1137/040611100.

27.

Stampacchia, Equations elliptiques du second ordre à coefficients discontinus, Séminaire Jean Leray 3 (1963–1964), 1–77.

28.

J.L.

Taylor,

Kim and

R.M.

Brown, The Green function for elliptic systems in two dimensions, Communications in Partial Differential Equations 38 (2013), 1574–1600. doi:10.1080/03605302.2013.814668.

29.

C.H.

Wolters, The finite element method in EEG/MEG source analysis, SIAM News 40 (2007), 165–177.

30.

C.H.

Wolters,

Köstler,

Möller,

Härdtlein,

Grasedyck and

Hackbusch, Numerical mathematics of the subtraction method for the modeling of a current dipole in eeg source reconstruction using finite element head models, SIAM Journal on Scientific Computing 30 (2008), 24–45. doi:10.1137/060659053.

31.

C.H.

Wolters,

Kuhn,

Anwander and

Reitzinger, A parallel algebraic multigrid solver for finite element method based source localization in the human brain, Computing and Visualization in Science 5 (2002), 165–177. doi:10.1007/s00791-002-0098-0.