Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 26, pp. 1–15.

ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu

BOUNDARY-DOMAIN INTEGRAL EQUATIONS FOR

DIRICHLET DIFFUSION PROBLEMS WITH

NON-SMOOTH COEFFICIENT

CARLOS FRESNEDA-PORTILLO, ZENEBE W. WOLDEMICHEAL

Abstract. We obtain a system of boundary-domain integral equations (BDIE)
equivalent to the Dirichlet problem for the diffusion equation in non-homoge-

neous media. We use an extended version of the boundary integral method for

PDEs with variable coefficients for which a parametrix is required. We gener-
alize existing results for this family of parametrices considering a non-smooth

variable coefficient in the PDE and source term in Hs−2(Ω), 1/2 < s < 3/2
defined on a Lipschitz domain. The main results concern the equivalence be-

tween the original BVP and the corresponding BDIE system, as well as the

well-posedness of the BDIE system.

1. Introduction

The popularity of the boundary integral equation method (BIE), also called
the potential method, is owed to the reduction of dimension of a boundary value
problem (BVP) with constant coefficients and homogeneous right hand side de-
fined on a domain of Rn. By applying the BIE method, one can reformulate the
original BVP in terms of an equivalent integral equation defined exclusively on the
boundary of the domain. This method has already been extensively studied for
many boundary value problems, for instance: Laplace, Helmholtz, Stokes, Lamé in
[13, 14, 24]. This approach requires an explicit formula for the fundamental solution
of the PDE operator in the BVP which is not always available when the BVP has
variable coefficients [6, 21].

The overcome this issue, one can construct a parametrix (Levi function) [6, 21]
for the PDE operator and use it to derive an equivalent system of Boundary-Domain
Integral Equations following a similar approach as for the BIE method. However,
the reduction of dimension no longer applies as volume integrals will appear in
the new formulation as a result of the remainder term. This is also the case for
non-homogeneous problems with constant coefficients [13, Chapter 1 and 2].

To preserve the reduction of dimension, one can use the radial integration method
(RIM) which allows to transform volume integrals into boundary only integrals [11].
This method has been successfully implemented to solve boundary-domain integral

2020 Mathematics Subject Classification. 35J25, 31B10, 45K05, 45A05.
Key words and phrases. Non-smooth coefficients; parametrix; Lipschitz domain;

diffusion equation; boundary-domain integral equations; minimal wave speed.
©2022. This work is licensed under a CC BY 4.0 license.

Submitted May 31, 2021. Published March 4, 2022.

1



2 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

equations derived from BVPs with variable coefficients [1, 3] and has the ability to
remove singularities appearing in the domain integrals.

Recent developments on numerical approximation of the solution of BDIEs show
that there are effective and fast algorithms able to compute the solution. For
example: the collocation method [22, 23] which, although leads to fully populated
matrices, can be further enhanced by using hierarchical matrix compression and
adaptive methods as shown in [12] to reduce the computational cost. Localised
approaches to reduce the matrix dimension and storage have also been developed
[15] which lead to sparse matrices.

Moreover, reformulating the original BVP in the Boundary Domain Integral
Equation form can be beneficial, for instance, in inverse problems with variable
coefficients [5].

In this article, we consider the diffusion equation in non homogeneous media, i.e.

Au(x) :=

3∑
i=1

∂

∂xi

(
a(x)

∂u(x)

∂xi

)
, (1.1)

A possible parametrix is given by

P y(x, y) = P (x, y; a(y)) =
−1

4πa(y)|x− y|
,

where the superscript indicates P y(x, y) is a function of the variable coefficient
depending on y. This parametrix has been extensively studied [6, 7, 18] and more
references therein. However, a parametrix for a given partial differential operator
is not necessarily unique. For instance, another possible parametrix for the same
operator A is given by

P x(x, y) = P (x, y; a(x)) =
−1

4πa(x)|x− y|
.

In this case, the parametrix depends on the variable coefficient a(x). This parametrix
was introduced in [19] for the mixed problem for the operator A in Lipschitz 3D
domains with smooth coefficient in [19] and for the Dirichlet problem with smooth
coefficient in Lipschitz domains and f ∈ H−1(Ω) in [10]. In this article, we gener-
alize these results for the Dirichlet problem with non-smooth coefficient.

The study of new families of parametrices is helpful at the time of constructing
parametrices for systems of PDEs. For instance, for the Stokes system. In this
case, the fundamental solution for the pressure does not present any relationship
with the viscosity coefficient, whereas the parametrix for the pressure depends on
two variable viscosity coefficients: one depending on y and another depending on
x, see also [20].

However, most of the numerical methods to solve BDIEs aforementioned are
tested for the Dirichlet problem [12, 22, 23, 1] with smooth coefficient. However,
there are applications in Science and Engineering of elliptic PDEs of this type with
non-smooth coefficient, see for example [2, 9]. In order to compare the performance
of the parametrices P x(x, y) and P y(x, y), one needs first to prove the equivalence
between the original Dirichlet BVP and the system of BDIEs as well as the unique-
ness of the solution (well-posedness) of the system of BDIEs, which corresponds to
another of the various purposes of this article.



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 3

2. Preliminaries

2.1. Function spaces. Let Ω = Ω+ be a bounded simply connected open Lipschitz
domain, Ω− := R3 \ Ω̄+ the complementary (unbounded) domain. We assume that
their common boundary ∂Ω is a compact and simply connected surface.

In what follows D(Ω) := C∞0 (Ω) denotes the space of test functions with compact
support, and D∗(Ω) denotes its topological dual space, i.e. the space of Schwartz
distributions or generalized functions. Furthermore, we use the Bessel potential
spaces on the domain Ω and its boundary ∂Ω denoted respectively by Hs(Ω) =
{r

Ω
g : g ∈ Hs(R3)} and Hs(∂Ω), where s ∈ R (see e.g. [14, 24] for more details).

We recall that the spaces Hs coincide with the Sobolev-Slobodetski spaces W s,2

for any non-negative s [4].

By H̃s(Ω) we denote the closure of D(Ω) in Hs(R3). Furthermore, for Lipschitz
domains, the following characterization applies

H̃s(Ω) := {g ∈ Hs(R3) : supp g ⊂ Ω}.

Sometimes, we will also require the space

H̃s
•(Ω) := rΩH̃

s(Ω) := {rΩg : g ∈ H̃s(Ω)} ⊂ Hs(Ω),

where r
Ω

denotes the restriction operator on Ω.
Moreover, we shall introduce the fundamental concepts and notation of the

Hölder spaces which will be useful for the treatment of non-smooth coefficients.
Let θ be a real number with 0 < θ ≤ 1. Then a function f : Ω→ R is said to be

Hölder continuous with exponent θ if its Hölder coefficient |f |0,θ, defined as

|f |0,θ := sup
x 6=y; x,y∈Ω

|f(x)− f(y)|
‖x− y‖θ

,

is finite. The set of all the Hölder continuous functions with exponent θ is denoted
by C0,θ(Ω). This space is a locally convex topological vector space endowed with
the seminorm given by the Hölder coefficient

With this definition in mind, we can generalize the definition of the Hölder space
for functions in Cm,θ(Ω) as the subset of Cm(Ω) whose derivatives of order m are
Hölder continuous with exponent θ on the closure of Ω. The space Cm,θ(Ω) is a
Banach space endowed with the norm

‖f‖m,θ := ‖f‖Cm(Ω) + max
α=m

|Dαf |0,θ

where α ranges over multindices.
Following the notation from [18, Definition 2.5], given µ a non negative real

number, we define the space Cµ+(Ω)

Cµ+(Ω) :=


L∞(Ω) µ = 0

Cµ−1,1(Ω) µ ∈ Z+

Cm,θ+ε(Ω) where ε > 0 when µ = m+ θ, m ∈ Z+, θ ∈ (0, 1).

2.2. Partial differential equation. Let us consider the scalar elliptic partial dif-
ferential equation

Au(x) :=

3∑
i=1

∂

∂xi

(
a(x)

∂u(x)

∂xi

)
= f(x), x ∈ Ω, (2.1)



4 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

where the coefficient a(x) ∈ C |s−1|
+ (Ω) with s ∈ (1/2, 3/2), is bounded and satisfies

the assumption

0 < a0 ≤ a(x) ≤ a1,

for some constants a0, a1 ∈ R+. Additionally, this equation can be written in weak
form as

〈Au, v〉Ω := −EΩ(u, v), ∀v ∈ D(Ω), (2.2)

where u ∈ Hs(Ω) and the bilinear functional EΩ is given by

EΩ(u, v) := 〈a∇u,∇v〉Ω :=

3∑
i=1

〈a∂iu, ∂iv〉Ω.

Indeed, because of the deseness of D(Ω) in H̃2−s(Ω), the operator A from (2.2) de-

fines a continuous mapping A : Hs(Ω) → Hs−2(Ω) = [H̃2−s(Ω)]∗, (cf. [16, Section
3.1])

〈Au, v〉Ω := −EΩ(u, v) ∀u ∈ Hs(Ω), v ∈ H̃2−s(Ω).

Let us introduce the operator, Ǎ : Hs(Ω)→ H̃s−2(Ω) = [H2−s(Ω)]∗, defined by

〈Ǎu, v〉Ω := −ĚΩ(u, v) = 〈∇ · E̊Ω(a∇u), v〉Ω, ∀u ∈ Hs(Ω), v ∈ H2−s(Ω), (2.3)

where E̊Ω denotes the extension operator of functions defined in Ω, by zero outside
Ω in R3. Alternatively, the operator Ǎ can also be defined as

Ǎu := ∇ · E̊ΩrΩ(a∇u), (2.4)

where rΩ is the restriction operator to Ω. Let us note the operator Ǎ is bounded
as well as a result of the boundedness of A, E̊Ω and rΩ.

2.3. Traces and conormal derivatives. The Trace Theorem [8, Lemma 3.7]
states that if a distribution u belongs to Hs(Ω), with 1/2 < s < 3/2, then

γ+u ∈ Hs− 1
2 (∂Ω), where γ+ := γ+

∂Ω is the trace operator on ∂Ω from Ω. Let

γ−1 : Hs− 1
2 (∂Ω)→ Hs(R3) denote a (non-unique) continuous right inverse to the

trace operator γ+, i.e., γ+
∂Ωγ

−1w = w for any w ∈ Hs− 1
2 (∂Ω).

For u ∈ Hs(Ω), s > 3/2, and a ∈ C(Ω), we can define on ∂Ω the co-normal
derivative operator T±, in the classical sense as

T±x u :=

3∑
i=1

a(x)γ±
( ∂u
∂xi

)
n±i (x),

where n+(x) is the exterior unit normal vector directed outwards from the interior
domain Ω at a point x ∈ ∂Ω. Similarly, n−(x) is the unit normal vector directed
inwards to the interior domain Ω at a point x ∈ ∂Ω. Sometimes, we will also
use the notation T±x u or T±y u to emphasise with respect to which variable we
are differentiating. When the variable of differentiation is obvious or is a dummy
variable, we will simply use the notation T±u.

However, the classical co-normal derivative operator is generally not well defined
if u ∈ Hs(Ω), 1

2 < s < 3
2 , see [18, Appendix A]. To overcome this issue, one can use

the generalized co-normal derivative instead; see for instance [16, Definition 3.1].

Definition 2.1. Let 1/2 < s < 3/2, u ∈ Hs(Ω), a ∈ C |s−1|
+ (Ω) and rΩAu = rΩf̃

for some f̃ ∈ H̃s−2(Ω). Furthermore, let (γ−1)∗ : H−s(Rn) → H−s+
1
2 (∂Ω) be

the operator satisfying 〈(γ−1)∗ϕ,w〉∂Ω := 〈ϕ, γ−1w〉Rn for any w ∈ Hs− 1
2 (∂Ω) and



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 5

for any ϕ ∈ H−s(Rn). Then, the generalized co-normal derivative T+(f̃ , u) ∈
Hs− 3

2 (∂Ω) is defined by

〈T+(f̃ , u), w〉
∂Ω

:= 〈f̃ , γ−1w〉Ω + ĚΩ(u, γ−1w) = 〈f̃ − Ǎu, γ−1w〉Ω, (2.5)

for all w ∈ H 3
2−s(∂Ω), or alternatively by T+(f̃ , u) := (γ−1)∗(f̃ − Ǎu).

The key property of the generalized co-normal derivative T+(f̃ ;u) is its inde-
pendence from the (non-unique) choice of the operator γ−1, see [16, Theorem 5.3].
Moreover, the first Green identity holds in the following form for u ∈ Hs(Ω) with

rΩAu = rΩf̃ for some f̃ ∈ H̃s−2(Ω),

〈T+(f̃ , u), γ+v〉
∂Ω

= 〈f̃ , v〉Ω + ĚΩ(u, v) = 〈f̃ − Ǎu, v〉Ω, ∀v ∈ H2−s(Ω). (2.6)

Unfortunately, the generalized co-normal derivative operator is non-linear with
respect to u for a given fixed f̃ . Nevertheless, it is possible to gain linearity if the
domain of the generalized co-normal derivative is constrained to an appropriate
subspace of Hs(Ω), as the one given below.

Definition 2.2. Let s ∈ R, and A∗ : Hs(Ω) → D∗(Ω) be a linear operator. For
t ∈ R, we introduce the space

Hs,t(Ω;A∗) := {g : g ∈ Hs(Ω) : A∗g ∈ H̃t
•(Ω)},

endowed with the norm

‖g‖Hs,t(Ω;A∗) :=
(
‖g‖2Hs(Ω) + ‖A∗g‖2H̃t

•(Ω)

)1/2

.

The operator resulting from constraining the domain to a Sobolev space of the
type defined above is the canonical co-normal derivative [16, Definition 6.5].

Definition 2.3. For u ∈ Hs,−1/2(Ω;A) and a ∈ C
|s−1|
+ (Ω), 1/2 < s < 3/2, we

define the canonical co-normal derivative T+u ∈ Hs− 3
2 (∂Ω) as

〈T+u,w〉
∂Ω

:= 〈Ãu, γ−1w〉Ω + ĚΩ(u, γ−1w) = 〈Ãu− Ǎu, γ−1w〉Ω (2.7)

for all w ∈ H 3
2−s(∂Ω); that is, T+u := (γ−1)∗(Ãu− Ǎu).

By [16, Theorem 3.9] and [17, Theorem 6.6], the canonical co-normal derivative
T+u is independent of the (non-unique) choice of the operator γ−1, the operator

T+ : Hs,−1/2(Ω;A)→ Hs− 3
2 (∂Ω) is continuous, and the first Green identity holds

in the form

〈T+u, γ+v〉
∂Ω

:= 〈Ãu, v〉Ω + ĚΩ(u, v) = 〈Ãu− Ǎu, v〉Ω ∀v ∈ H2−s(Ω). (2.8)

The operator T+ : Hs,t(Ω;A) → Hs− 3
2 (∂Ω) in Definition 2.5 is continuous for

t ≥ − 1
2 . The canonical co-normal derivative is defined by the function u and the

operator A only, and thus does not depend separately on the right hand side f̃ (i.e.
its behavior on the boundary), unlike the generalized co-normal derivative defined
in (2.6). Additionally, the canonical co-normal derivative operator T+ is linear in
u.

For further reference on the generalized and canonical co-normal derivative, we
refer the reader to [16, 17].



6 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

2.4. Green identities. Let 1/2 < s < 3/2 and a ∈ C |s−1|
+ (Ω). If u ∈ Hs,−1/2(Ω;A),

then Definitions (2.5) and (2.7) imply that the generalized co-normal derivative for

arbitrary extension f̃ ∈ H̃s−2(Ω) of the distributions rΩAu can be expressed as

〈T+(f̃ , u), w〉
∂Ω

:= 〈T+u,w〉
∂Ω

+ 〈f̃ − Ãu, γ−1w〉Ω ∀w ∈ H 3
2−s(Ω).

If u ∈ Hs(Ω) and v ∈ H2−s,− 1
2 (Ω;A), then swapping over the roles of u and v in

(2.8), we obtain the first Green identity for v,

ĚΩ(u, v) + 〈u, Ãv〉Ω = 〈T+v, γ+u〉
∂Ω
. (2.9)

Furthermore, rΩAu = rΩf̃ with f̃ ∈ H̃s−2(Ω), then subtracting (2.9) from (2.6),
the following second Green identity is obtained,

〈f̃ , v〉Ω − 〈u, Ãv〉Ω = 〈T+(f̃ , u), γ+v〉
∂Ω
− 〈T+v, γ+u〉

∂Ω
. (2.10)

When u ∈ Hs,−1/2(Ω;A) and v ∈ H2−s,− 1
2 (Ω;A), we arrive at the familiar form of

the second Green identity,

〈v, Ãu〉Ω − 〈u, Ãv〉Ω = 〈T+u, γ+v〉
∂Ω
− 〈T+v, γ+u〉

∂Ω
.

2.5. Boundary value problem. For 1/2 < s < 3/2, we aim to derive boundary-
domain integral equation systems for the following Dirichlet boundary value prob-
lem. Find a function u ∈ Hs(Ω) satisfying

Au = f in Ω, (2.11a)

γ+u = ϕ0 on ∂Ω, (2.11b)

where ϕ0 ∈ Hs− 1
2 (∂Ω) and f ∈ Hs−2(Ω).

The following assertion is well known for s = 1 and follows from the first Green
identity and the Lax Milgram Theorem, see [18, Theorem 5.1].

Theorem 2.4. Let 1/2 < s < 3/2 and a ∈ C
|s−1|
+ (Ω). The Dirichlet problem

(2.11) has at most one solution in Hs(Ω).

3. Parametrices and remainders

One of the purposes of this article is to obtain a system of integral equations that
is equivalent to (2.11). To obtain such a system, one needs first to obtain an integral
representation of the solution in terms of surface and volume potentials. In the case
of constant coefficients for the Laplace equation, one substitutes the fundamental
solution in the second Green identity, see for example [13, 24]. However, when the
PDE has variable coefficients, although the fundamental solution might exist, it
may not always be available explicitly as required for numerical approximation of
the solution. One possible way around this is to introduce a parametrix.

Definition 3.1. A distribution P is said to be a parametrix for a given differential
operator B with remainder R if

BP = δ +R.

where δ(·) is the Dirac distribution.

For a given operator B, the parametrix might not be unique. For example, the
parametrix

P y(x, y) =
1

a(y)
P∆(x− y), x, y ∈ R3,



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 7

was employed in [15, 6] for the operator A defined in (2.1), where

P∆(x− y) =
−1

4π|x− y|
is the fundamental solution of the Laplace operator. The remainder corresponding
to the parametrix P y is

Ry(x, y) =

3∑
i=1

1

a(y)

∂a(x)

∂xi

∂

∂xi
P∆(x− y), x, y ∈ R3.

In this article, for the same operator A defined in (2.1), we use another parametrix,

P (x, y) := P x(x, y) =
1

a(x)
P∆(x− y), x, y ∈ R3, (3.1)

which leads to the corresponding remainder

R(x, y) = Rx(x, y) = −
3∑
i=1

∂

∂xi

( 1

a(x)

∂a(x)

∂xi
P∆(x, y)

)
= −

3∑
i=1

∂

∂xi

(∂ ln a(x)

∂xi
P∆(x, y)

)
, x, y ∈ R3.

(3.2)

Note that both the remainders Rx and Ry are weakly singular, i.e.,

Rx(x, y), Ry(x, y) ∈ O(|x− y|−2).

4. Domain and boundary integral operators

After replacing the parametrix in the Green identities, we obtain an integral
representation formula of the solution in terms of potential-type operators whose
kernel is somehow related to the parametrix. Therefore, these operators receive
the name of surface or volume parametrix based potential operators depending on
whether the integral operator is defined exclusively on the boundary or within the
domain.

In this section, we introduce the parametrix-based volume and surface potentials
and analyze their mapping properties.

4.1. Volume parametrix-based potentials. Let ρ ∈ D(Ω), the volume poten-
tial and the remainder potential operator, corresponding to parametrix (3.1) and
remainder (3.2) are defined as

Pρ(y) :=

∫
R3

P (x, y)ρ(x) dx, y ∈ R3,

Pρ(y) :=

∫
Ω

P (x, y)ρ(x) dx, y ∈ Ω,

Rρ(y) =

∫
Ω

R(x, y)ρ(x) dx, y ∈ Ω.

Since the parametrix and the remainder are related to the fundamental solution of
the Laplace equation, as shown in (3.1)-(3.2), then the volume parametrix-based
potentials also preserve a similar relation with their analogous volume potentials of
the constant coefficient case, i.e. based on the fundamental solution, see [20].

Pρ = P4

(ρ
a

)
, (4.1)



8 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

Pρ = P∆

(ρ
a

)
, (4.2)

Rρ = ∇ · [P∆(ρ∇ ln a)]− P∆(ρ∆ ln a). (4.3)

Taking into account relations (4.1)-(4.3) and that the variable coefficient is con-
tinuous and bounded, one can prove the following mapping properties which result
from applying the analogous mapping properties of the volume newtonian potential
for the constant coefficient case given in [18, Lemma 3.1] using a similar argument
as in [18, Theorem 3.2].

Theorem 4.1. Let Ω be a bounded Lipschitz domain in R3, the following operators
are continuous:

µP : Hs(R3)→ Hs+2(R3), s ∈ R, a ∈ C |s+2|
+ (R3), ∀µ ∈ D(R3);

P : H̃s(Ω)→ Hs+2(Ω), s ∈ R, a ∈ C |s+2|
+ (Ω);

: Hs(Ω)→ Hs+2(Ω), −1

2
< s <

1

2
, a ∈ Cs+2

+ (Ω);

R : H̃s(Ω)→ Hs+1(Ω), s ∈ R, a ∈ C |s+1|
+ (Ω);

: Hs(Ω)→ Hs+1(Ω), −1

2
< s <

1

2
, a ∈ Cs+1

+ (Ω);

: Hs(Ω)→ Hs,−1/2(Ω;A),
1

2
< s <

3

2
, a ∈ C3/2

+ (Ω).

R : Hs(Ω)→ Hs(Ω),
1

2
< s <

3

2
, a ∈ Cs+(Ω);

γ+R : Hs(Ω)→ Hs− 1
2 (∂Ω),

1

2
< s <

3

2
, a ∈ Cs+(Ω).

The following corollary follows directly from applying the Rellich Theorem on
compact embeddings for Sobolev spaces [14, Theorem 3.27] and the Trace Theorem
to the mapping properties given in the previous theorem for the R operator.

Corollary 4.2. Let Ω be a bounded Lipschitz domain in R3. Then, the following
operators are compact.

R : Hs(Ω)→ Hs(Ω),
1

2
< s <

3

2
, a ∈ Cs+(Ω);

γ+R : Hs(Ω)→ Hs− 1
2 (∂Ω),

1

2
< s <

3

2
, a ∈ Cs+(Ω).

4.2. Surface parametrix-based potentials. The parametrix-based single layer
and double layer surface potentials are defined for y ∈ R3 : y /∈ ∂Ω, as

V ρ(y) := −
∫
∂Ω

P (x, y)ρ(x) dS(x), Wρ(y) := −
∫
∂Ω

T+
x P (x, y)ρ(x) dS(x).

By relation (3.1), the operators V and W can also be written in terms of the
surface potentials associated with the Laplace operator,

V ρ = V∆

(ρ
a

)
, (4.4)

Wρ = W∆ρ− V∆

(
ρ
∂ ln a

∂n

)
. (4.5)

Similar to the previous section, one can derive the mapping properties of the single
layer and double layer parametrix-based surface potential operators taking into



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 9

account relations (4.4)-(4.5) along with the analogous mapping properties of their
counterpart operators of the constant coefficient case given in [18, Theorem 3.3]
using a similar argument as in [18, Theorem 3.5].

Theorem 4.3. Let Ω be a bounded Lipschitz domain.

(i) The following operators are continuous if 1/2 ≤ s ≤ 3/2:

µV : Hs− 3
2 (∂Ω)→ Hs(R3), a ∈ Cs+(R3), ∀µ ∈ D(R3);

µW : Hs− 1
2 (∂Ω)→ Hs(Ω), a ∈ Cs+(Ω);

µrΩ−W : Hs− 1
2 (∂Ω)→ Hs(Ω−), a ∈ Cs+(Ω−), ∀µ ∈ D(R3).

(ii) The following operators are continuous if 1/2 < s ≤ 3/2 and a ∈ C3/2
+ (Ω):

rΩV : Hs− 3
2 (∂Ω)→ Hs,−1/2(Ω;A);

µrΩ−V : Hs− 3
2 (∂Ω)→ H

s,−1/2
loc (Ω−;A), ∀µ ∈ D(R3);

rΩW : Hs− 1
2 (∂Ω)→ Hs,−1/2(Ω;A);

µrΩ−W : Hs− 1
2 (∂Ω)→ Hs,−1/2(Ω−;A), ∀µ ∈ D(R3);

(iii) The following operators are continuous if 1/2 < s < 3/2,

γ±V : Hs− 3
2 (∂Ω)→ Hs− 1

2 (∂Ω), a ∈ Cs+(Ω±);

γ±W : Hs− 1
2 (∂Ω)→ Hs− 1

2 (∂Ω), a ∈ Cs+(Ω±).

The next theorem follows from the well-known jump properties of the boundary
integral operators for the Laplace operator and relations (4.4)-(4.5), see also [14,
Theorem 6.11] or [18, Theorem 3.3 (iii)].

Theorem 4.4. Let ∂Ω be a compact Lipschitz boundary, 1/2 < s < 3/2, ϕ ∈
Hs− 1

2 (∂Ω) and ψ ∈ Hs− 3
2 (∂Ω). Then

γ+V ψ − γ−V ψ = 0, γ+Wϕ− γ−Wϕ = −ϕ, if a ∈ Cs+(R3). (4.6)

The previous jump properties lead to the definition of two new boundary integral
operators related with the traces of the single layer and double layer operators, also
called direct values of the single and double layer potentials.

V := γ+V, W :=
1

2
(γ+W + γ−W ).

The mapping properties of V and W are given by the following result which is a
direct consequence of applying the Trace Theorem to the mapping properties of the
operator V and W given in Theorem 4.3.

Corollary 4.5. Let ∂Ω be a compact Lipschitz boundary, 1/2 < s < 3/2. The
following operators are continuous:

V : Hs− 3
2 (∂Ω)→ Hs− 1

2 (∂Ω), a ∈ Cs+(Ω±); (4.7)

W : Hs− 1
2 (∂Ω)→ Hs− 1

2 (∂Ω), a ∈ Cs+(Ω±). (4.8)

Employing definitions (4.7) and (4.8), the jump properties (4.6) can be re-written
as

γ±V ψ = Vψ, γ±Wϕ = ∓1

2
ϕ+Wϕ, if a ∈ Cs+(R3). (4.9)



10 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

By Corollary 4.5 and relations (4.4) and (4.5), the operators V and W can be
expressed in terms of the volume and surface potentials and operators associated
with the Laplace operator, as follows [20],

Vρ = V∆

(ρ
a

)
,

Wρ =W∆ρ− V∆

(
ρ
∂ ln a

∂n

)
.

A key result for proving uniqueness of the solution of the integral equation systems
that we will derive in the upcoming sections is related to the invertibility of the
single layer potential. This result is well known (see [14, Theorem 8.12 and Corollary
8.13] or [18, Lemma 4.8 (i)].) for the constant coefficient case, and can easily be
extrapolated to the variable coefficient case thanks to the relation (4.4).

Lemma 4.6. Let Ω be a bounded simply-connected Lipschitz domain, let a ∈ Cs+(Ω)
and let 1/2 < s < 3/2.

(1) The operator V : Hs(∂Ω)→ Hs+1(∂Ω) is an isomorphism.

(2) If Ψ∗ ∈ Hs− 3
2 (∂Ω) and

rΩVΨ∗(y) = 0, y ∈ Ω

then Ψ∗(y) = 0.

4.3. Integral representation of the solution of the BVP. In this section, we
apply the usual procedures of the boundary-domain integral equation method [6, 10]
for obtaining an integral equation system derived from the integral representation
formula which results from substituting the parametrix (3.1) in the first or second
Green identities.

Let u ∈ Hs(Ω), 1/2 < s < 3/2, a ∈ Cs+(Ω), and let us substitute v(x) =
P (x, y), in the first Green identity (2.10). Then we obtain the following integral
representation formula, also called generalized third Green identity, for the solution
of (2.11),

u+Ru+Wγ+u = PǍu in Ω, (4.10)

where the term in the right hand side, considering the definitions of the operator
Ǎ given in (2.3) and (2.4), is given by

PǍu(y) := 〈Ǎu, P (., y)〉Ω = −Ě(u, P (., y)) = −〈Ẽs−1
Ω (a∇u),∇P (., y)〉Ω, y ∈ Ω.

When Au = r
Ω
f̃ in Ω, where f̃ ∈ H̃s−2(Ω), the integral representation formula

(4.10) can be reformulated as

u+Ru− V T+(f̃ , u) +Wγ+u = P f̃ in Ω. (4.11)

Since we aim to obtain a system of integral equations equivalent to the Dirichlet
problem (2.11), it will be necessary to obtain an integral representation formula for
the trace of the solution. Then, by applying the trace operator to both sides of
the identity (4.11) and taking into account the jump properties given in (4.9), we
obtain

1

2
γ+u+ γ+Ru− VT+(f̃ , u) +Wγ+u = γ+P f̃ , on ∂Ω. (4.12)



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 11

5. Boundary-domain integral equation system

We aim to obtain a system of integral equations equivalent to Dirichlet BVP
(2.11). The integral equations will be written in terms of the surface and vol-
ume parametrix-based potentials which are defined on the boundary and in the
domain respectively. Hence, this system of integral equations receives the name of
Boundary-Domain Integral Equation System (BDIES). The BDIES is obtained fol-
lowing the boundary-domain integral equations method which extends the Bound-
ary Integral Equation method for the variable coefficient case, see e.g., [19, 18, 6]
for more details.

Let a ∈ Cs+(Ω) with s ∈ ( 1
2 ,

3
2 ). Furthermore, let f̃ ∈ H̃s−2(Ω) be an extension

of f ∈ Hs−2(Ω) (i.e., f = rΩf̃), which always exists [16, Lemma 2.15 and Theo-
rem 2.16]. Then, to obtain the boundary domain integral equations of the system
from the integral representation formulas (4.12)-(6.1), we introduce the unknown

co-normal derivative ψ := T+(f̃ , u) ∈ Hs− 3
2 (∂Ω) which shall be regarded as for-

mally segregated from u. On the other hand, let us recall that the trace datum for
the Dirichlet problem is known. Hence, we replace in (4.12)-(6.1) γ+u = ϕ0 and we
move all these terms to the right hand side. We obtain the following system con-
sisting of two equations with two unknown functions (u, ψ) ∈ Hs(Ω)×Hs− 3

2 (∂Ω),

u+Ru− V ψ = F0 in Ω, (5.1a)

γ+Ru− Vψ = γ+F0 − ϕ0 on ∂Ω, (5.1b)

where

F0 = P f̃ −Wϕ0. (5.2)

Note that for ϕ0 ∈ Hs− 1
2 (∂Ω), we have the inclusion F0 ∈ Hs(Ω) if f̃ ∈ H̃s−2(Ω)

due to the mapping properties of the surface and volume potentials.
The system (5.1), given by (5.1a)-(5.1b) can be written in matrix notation as

A1U = F1,

where U represents the vector containing the unknowns of the system,

U = (u, ψ)> ∈ Hs(Ω)×Hs− 3
2 (∂Ω),

the right-hand side vector is

F1 := [F0, γ
+F0 − ϕ0]> ∈ Hs(Ω)×Hs− 1

2 (∂Ω),

and the matrix operator A1 is defined by

A1 =

[
I +R −V
γ+R −V

]
.

6. Main results

Let us prove that the Dirichlet boundary value problem (2.11) in Ω is equivalent
to the system of the Boundary Domain Integral Equations (5.1a), (5.1b). Before
proving the main equivalence theorem, let us prove an instrumental lemma which
will simplify the proof of the main result.

The following statement is similar to [18, Lemma 4.2] in which a different
parametrix was used. It extends [10, Lemma 4.1], where the corresponding as-

sertion was proven for f̃ ∈ L2(Ω), s = 1, a ∈ C∞(Ω).



12 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

Lemma 6.1. Let 1/2 < s < 3/2 and a ∈ Cs+(Ω). Let u ∈ Hs(Ω), Ψ ∈ Hs− 3
2 (∂Ω),

Φ ∈ Hs− 1
2 (∂Ω), and f̃ ∈ H̃s−2(Ω) satisfy

u+Ru− VΨ +WΦ = P f̃ in Ω. (6.1)

Then

Au = rΩf̃ in Ω,

V (Ψ− T+(f̃ , u))−W (Φ− γ+u) = 0 in Ω.

Proof. We subtract (6.1) from the first Green identity (4.10) to obtain

V (Ψ)−W (Φ− γ+u) = P(Ǎu− f̃). (6.2)

Then (6.2) can be rewritten in terms of the constant coefficient volume and surface
potentials by applying the relations (4.2), (4.4), and (4.5), as

V∆

(Ψ

a

)
−W∆(Φ− γ+u) + V∆

(∂ ln a

∂n
(Φ− γ+u)

)
= P∆

( Ǎu− f̃
a

)
.

Then, we apply the Laplace operator to both sides of the previous equation taking
into account that the terms with V∆ and W∆ are harmonic functions and hence
vanish. Furthermore, the Laplace operator applied to the newtonian potential gives

∆P∆

( Ǎu− f̃
a

)
=
Ǎu− f̃

a
.

Since the coefficient a(x) is strictly positive,

Ǎu− f̃ = 0 in Ω. (6.3)

This implies that rΩǍu = Au = rΩf̃ , which completes the proof of part (i).

On the other hand, part (i) implies that f̃ is an extension of the distribution

Au ∈ Hs−2(Ω) to H̃s−2(Ω). Therefore, by applying the first Green identity (2.6)
with v = P , we obtain (see [18, Lemma 4.2])

P(ǍΩu− f̃) = V T+(f̃ , u) in Ω. (6.4)

Finally, substituting (6.3) and (6.4) into (6.2), we arrive at the relation in part
(ii). �

The above lemma will simplify considerably the proof of the Equivalence The-
orem, one of the main results of this paper. In particular, the lemma will help us
showing that a solution of the BDIES (5.1) is a solution of the BVP.

Theorem 6.2. Let 1/2 < s < 3/2 and a ∈ Cs+(Ω). Let ϕ0 ∈ Hs− 1
2 (∂Ω), f ∈

Hs−2(Ω) and f̃ ∈ H̃s−2(Ω) is such that rΩf̃ = f .

(i) If a function u ∈ Hs(Ω) solves the Dirichlet BVP (2.11), then the pair

(u, ψ)> ∈ Hs(Ω)×Hs− 3
2 (∂Ω) where

ψ = T+(f̃ , u), on ∂Ω, (6.5)

solves the BDIE system (5.1).

(ii) If a pair (u, ψ)> ∈ Hs(Ω)×Hs− 3
2 (∂Ω) solves the BDIE system (5.1) then

u solves the BVP and the functions ψ satisfy (6.5).



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 13

Proof. (i) Let u ∈ Hs(Ω) be a solution of the boundary value problem (2.11) and

let ψ be defined by (6.5), evidently implying ψ ∈ Hs− 3
2 (∂Ω). Then, it immediately

follows from Theorem 6.1 and relations (6.1) and (4.12) that the pair (u, ψ) solves
BDIE system (5.1) with the right hand side F1 which completes the proof of item
(i).

(ii). Let now the pair (u, ψ)> ∈ Hs(Ω) × Hs− 3
2 (∂Ω) solve the BDIE system.

Taking the trace of the equation (5.1a) and subtract it from the equation (5.1b),
we obtain

γ+u = ϕ0, on ∂Ω. (6.6)

i.e. u satisfies the Dirichlet condition (2.11b).
Equation (5.1b) and Lemma 6.1 with Ψ = ψ,Φ = ϕ0 imply that u is a solution

of PDE (2.11a) and
VΨ∗ −WΦ∗ = 0 in Ω.

where Ψ∗ = ψ−T+(f̃ , u) and Φ∗ = φ−γ+u. From (6.6), Lemma 4.6 implies Ψ∗ = 0
which completes the proof of condition (6.6). Thus u obtained from solution of the
BDIE system (5.1) solves the Dirichlet problem. �

Theorem 6.3. Let Ω be a bounded simply-connected Lipschitz domain, 1/2 < s <
3/2, and σ = max{1, s}. Then the operator

A1 : Hs(Ω)×Hs− 3
2 (∂Ω)→ Hs(Ω)×Hs− 1

2 (∂Ω), if a ∈ Cσ+(Ω) (6.7)

satisfies:

(i) the operator A1 is bounded and Fredholm with zero index;
(ii) the operator A1 has a bounded inverse.

Proof. (i) The operator is continuous because the mapping properties of the oper-
ators involved in the matrix A1, see Theorem 4.1 and Theorem 4.3. To prove the
Fredholm property of operator (6.8), let us consider the operator

A1
0 : Hs(Ω)×Hs− 3

2 (∂Ω)→ Hs(Ω)×Hs− 1
2 (∂Ω), (6.8)

given by

A1
0 :=

[
I −V
0 −V

]
⇒ A1 −A1

0 =

[
R 0
γ+R 0

]
.

As a result of compactness properties of the operators R and γ+R (cf. Corollary
4.2), the operator A1

0 is a compact perturbation of operator (6.8), i.e. the operator
A1 −A1

0 is compact. The operator A1
0 is an upper triangular matrix operator with

scalar diagonal invertible operators I : Hs(Ω) → Hs(Ω) and V : Hs− 3
2 (∂Ω) →

Hs− 1
2 (∂Ω) where the invertibility of the operator V follows from Theorem 4.6.

This implies the invertibility of the operator A1
0. As a result, the operator (6.8) is

Fredholm with index zero.
To prove part (ii) which we consider two cases: (a) 1 ≤ s < 3/2 and (b) 1/2 <

s ≤ 1.
(a) Since 1 ≤ s < 3/2, we have σ = s. From the Equivalence Theorem (6.2)

(ii) and the uniqueness of the BVP provided by Theorem (2.4), we can deduce the
injectivity of the operator (6.7). Precisely, considering the homogeneous system
AU = 0, the zero right hand side F1 = 0 can be represented as in (5.2) in terms

of f̃ = 0 and ϕ0 = 0. By applying the Equivalence Theorem (6.2) (ii), then the
solution of this homogeneous system can be represented as U = (u, T+(0;u))T ,
where u is a solution of the Dirichlet problem (2.11) with right hand sides f = 0



14 C. FRESNEDA-PORTILLO, Z. W. WOLDEMICHEAL EJDE-2022/26

and ϕ0 = 0. However, the Dirichlet BVP has only the trivial solution u = 0,
because of Theorem (2.4). Hence, the only solution of the homogeneous BDIES
(5.1) is the trivial solution. Therefore, the Fredholm index of the operator (6.7),
provided by part (i), implies its invertibility for 1 ≤ s < 3/2.

(b) Since s ∈ (1/2, 1], we have σ = 1, and hence a ∈ C1
+(Ω). By Theorem (i), the

operator (6.7) is also Fredholm with zero index. Since the invertibility has already
been proved for the case s = 1, the result [18, Lemma 7.5] implies that the kernels
(null-spaces) of the operator (6.8) with s = 1 and the operator (6.8) for s ∈ ( 1

2 , 1)
consist of only the zero element for any s ∈ (1/2, 1], which implies that the operator
is invertible for all s ∈ (1/2, 1]. �

Remark 6.4. For a given function f ∈ Hs−2(Ω), its extension f̃ ∈ H̃s−2(Ω) is not
unique. Nevertheless, since the solution of a Dirichlet BVP (2.11) does not depend
on this extension, equivalence Theorem (6.2) (ii) implies that u in the solution of

BDIE systems A1 does not depend on the particular choice of extension f̃ , however,
ψ obviously does, see (6.5).

7. Conclusions

A new parametrix for the diffusion equation in non homogeneous media (with
variable coefficient) with Lipschitz boundary has been analyzed in this paper. Map-
ping properties of the corresponding parametrix-based surface and volume poten-
tials have been shown in corresponding Sobolev spaces.

A BDIE system, based on a new parametrix for the originalBVP has been ob-
tained. Equivalence of the BDIE system to the original BVP was proved in the
case when the PDE right-hand side function from Hs−2(Ω), 1

2 < s < 3
2 and the

Dirichlet data is from the space Hs− 3
2 (∂Ω). The invertibility of the matrix operator

defining the BDIE was proved in the corresponding Sobolev space.
Now, we have obtained an analogous system to the BDIE (5.1) of [6, 19] with a

new family of parametrices which is uniquely solvable. Hence, further investigation
about the numerical advantages of using one family of parametrices over another
will follow.

Analyzing BDIEs for different parametrices, i.e. depending on the variable coef-
ficient a(x) or a(y), is crucial in order to understand the analysis of BDIEs derived
with parametrices that depend on the variable coefficient a(x) and a(y) at the
same, as is the case for the Stokes system, see [19]. Therefore, applying a similar
approach, one could extend the results of this paper to elliptic and strongly elliptic
PDE systems.

References

[1] M. A. Al-Jawary, J. Ravnik, L. C. Wrobel, L. Škerget; Boundary element formulations for the

numerical solution of two-dimensional diffusion problems with variable coefficients, Comput.
Math. Appl. 64 (2012), no. 8, 2695–2711.

[2] C. Bernardi, R. Verfürth; Adaptive finite element methods for elliptic equations with non-

smooth coefficients, Numerische Mathematik 85, 2000.

[3] A. Beshley, R. Chapko, B. Tomas Johansson; An integral equation method for the numerical
solution of a Dirichlet problem for second-order elliptic equations with variable coefficients,

Journal of Engineering Mathematics 112 (2018), 63–73.
[4] H. Brezis; Functional analysis, Sobolev spaces and partial differential equations, Springer,

2011.



EJDE-2022/26 BOUNDARY-DOMAIN INTEGRAL EQUATIONS 15

[5] R. Chapko, B. Tomas Johansson; A boundary integral equation method for numerical solution

of parabolic and hyperbolic Cauchy problems, Appl. Numer. Math. 129 (2018), 104—119.

[6] O. Chkadua, S. E. Mikhailov, D. Natroshvil; Analysis of direct boundary-domain integral
equations for a mixed BVP with variable coefficient, i: Equivalence and invertibility, J.

Integral Equations Appl. 21 (2009), no. 4, 499–543.

[7] O. Chkadua, S. E. Mikhailov, D. Natroshvili; Analysis of direct segregated boundary-domain
integral equations for variable-coefficient mixed BVPs in exterior domains, Anal. Appl. (Sin-

gap.) 11 (2013), no. 4, 1350006, 33 pp.

[8] M. Costabel; Boundary integral operators on lipschitz domains: Elementary results, SIAM
J. Math. Anal. 19 (1988), no. 3, 613—626.

[9] E. Fernández-Cara, E. Zuazua; On the null controllability of the one-dimensional heat equa-

tion with bv coefficients, Computational and Applied Mathematics 21 (2002), no. 1, 167–190.
[10] C. Fresneda-Portillo, Z. Woldemicheal; A new family of boundary-domain integral equations

for the dirichlet problem of the diffusion equation in inhomogeneous media with H−1(ω)
source term on Lipschitz domains, Math. Methods Appl. Sci. 44 (2021), no. 12, 9817—9830.

[11] X. W. Gao; The radial integration method for evaluation of domain integrals with boundary-

only discretization, Engineering Analysis with Boundary Elements 26 (2002) No. 10, 905–916.
[12] R. Grzhibovskis, S. Mikhailov, Sergej Rjasanow; Numerics of boundary-domain integral and

integro-differential equations for BVP with variable coefficient in 3d, Comput. Mech. 51

(2013), no. 4, 495—503.
[13] G. C. Hsiao, W. L. Wendland; Boundary integral equations, vol. 164, Springer Berlin Heidel-

berg, 2008.

[14] W. McLean; Strongly elliptic systems and boundary integral equations, Cambridge university
press, 2000.

[15] S. E. Mikhailov; Localized boundary-domain integral formulations for problems with variable

coefficients, Engineering Analysis with Boundary Elements 26 (2002), 681–690.
[16] S. E. Mikhailov; Traces, extensions and co-normal derivatives for elliptic systems on lipschitz

domains, J. Math. Anal. Appl. 378 (2011), no. 1, 324-–342.
[17] Se. E. Mikhailov; Solution regularity and co-normal derivatives for elliptic systems with non-

smooth coefficients on Lipschitz domains, J. Math. Anal. Appl. 400 (2013), no. 1, 48-–67.

[18] S. E Mikhailov; Analysis of segregated boundary-domain integral equations for BVPs with
non-smooth coefficients on lipschitz domains, Bound. Value Probl. 2018, Paper No. 87, 52

pp.

[19] S. E. Mikhailov, C. Fresneda-Portillo; Analysis of boundary-domain integral equations based
on a new parametrix for the mixed diffusion BVP with variable coefficient in an interior

lipschitz domain, J. Integral Equations Appl. 32 (2020), no. 1, 59—75.

[20] S. E. Mikhailov, C. Fresneda-Portillo; Boundary-domain integral equations equivalent to an
exterior mixed BVP for the variable-viscosity compressible stokes PDEs, Commun. Pure

Appl. Anal. 20 (2021), no. 3, 1103—1133.

[21] A. Pomp; The boundary-domain integral method for elliptic systems: With application to
shells, Springer, 2006.

[22] J. Ravnik, L. Škerget; A gradient free integral equation for diffusion-convection equation with
variable coefficient and velocity, Eng. Anal. Bound. Elem. 37 (2013), no. 4, 683–690.

[23] J. Ravnik, L. Škerget; Integral equation formulation of an unsteady diffusion–convection

equation with variable coefficient and velocity, Comput. Math. Appl. 66 (2014), no. 12, 2477—
2488.

[24] O. Steinbach; Numerical approximation methods for elliptic boundary value problems,

Springer New York, 2008.

Carlos Fresneda Portillo

Department of Quantitative Methods, Universidad Loyola Andalućıa, Sevilla, 41704,

Spain
Email address: cfresneda@uloyola.es

Zenebe W. Woldemicheal
Department of Mathematics, Madda Walabu University, Bale Robe, P.O. Box 247, Ethiopia

Email address: zenebe.wogderesegn@aau.edu.et


	1. Introduction
	2. Preliminaries
	2.1. Function spaces
	2.2. Partial differential equation
	2.3. Traces and conormal derivatives
	2.4. Green identities
	2.5. Boundary value problem

	3. Parametrices and remainders
	4. Domain and boundary integral operators
	4.1. Volume parametrix-based potentials
	4.2. Surface parametrix-based potentials
	4.3. Integral representation of the solution of the BVP

	5. Boundary-domain integral equation system
	6. Main results
	7. Conclusions
	References