Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 84, pp. 1–32.

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

EXISTENCE OF A SOLUTION AND ITS NUMERICAL

APPROXIMATION FOR A STRONGLY NONLINEAR COUPLED

SYSTEM IN ANISOTROPIC ORLICZ-SOBOLEV SPACES

FRANCISCO ORTEGÓN GALLEGO, HAKIMA OUYAHYA, MOHAMED RHOUDAF

Abstract. We study the existence of a capacity solution for a nonlinear el-

liptic coupled system in anisotropic Orlicz-Sobolev spaces. The unknowns are
the temperature inside a semiconductor material, and the electric potential.

This system may be considered as a generalization of the steady-state ther-

mistor problem. The numerical solution is also analyzed by means of the least
squares method in combination with a conjugate gradient technique.

1. Introduction

This work concerns a generalization of the steady-state thermistor problem. It
consists of two coupled nonlinear elliptic equations governing the temperature, u,
and the electric potential, ϕ, inside a semiconductor device, namely,

−A(u) = ρ(u)|∇ϕ|2 in Ω

div(ρ(u)∇ϕ) = 0 in Ω,

ϕ = ϕ0 on ∂Ω,

u = 0 on ∂Ω,

(1.1)

where Ω ⊂ Rd (the thermistor geometry) is a bounded domain, d ≥ 2 is an integer,
and the operator A, given by

A(u) =

d∑
i=1

∂i
(
ai(x, u, ∂iu)

)
, ∂i =

∂

∂xi
,

which is assumed to be of the Leray-Lions type on certain Orlicz-Sobolev spaces. For
each i = 1, . . . , d, the function ai(x, s, ζ) : Ω×R×R 7→ R is a Carathéodory function,
that is, measurable with respect to x in Ω for all (s, ζ) ∈ R2, and continuous with
respect to (s, ζ) for almost every x ∈ Ω. The vector function a = (a1, . . . , ad)
satisfies certain monotonicity and coercivity conditions in the anisotropic Orlicz-
Sobolev space

W 1LM(Ω) = {u ∈ LM0
(Ω) : ∂iu ∈ LMi

(Ω), i = 1, . . . , d},

2020 Mathematics Subject Classification. 35J70, 35J66, 46E30, 65N22.
Key words and phrases. Nonlinear elliptic equations; capacity solution; least squares method;

anisotropic Orlicz-Sobolev spaces; conjugate gradient algorithm.
©2022. This work is licensed under a CC BY 4.0 license.

Submitted April 21, 2022. Published December 21, 2022.

1



2 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

where M = (M1, . . . ,Md), M1, . . . ,Md are N -functions which, in general, do not
satisfy the ∆2-condition, and M0 = min1≤i≤dMi. Also, ρ ∈ C (R) ∩ L∞(R) stands
for the temperature dependent electric conductivity and ϕ0 ∈ H1(Ω) ∩ L∞(Ω) is
given.

System (1.1) is a mathematical model which generalizes the so-called thermistor
problem [2, 6, 10, 13]. In most practical cases, one has ρ(s) > 0 for all s ∈ R and
ρ(s) → 0 as s → +∞. In particular, the equation for ϕ is nonuniformly elliptic,
and consequently, no a priori estimates for ∇ϕ will be available so that ϕ may not
belong to a Sobolev space. This means that the search for weak solutions to (1.1)
is not well-suited in this setting.

To deal with this difficulty, we consider the function Φ = ρ(u)|∇ϕ|2 as a whole
and then show that it belongs to L2(Ω)d. A new formulation of system (1.1) will
lead us to the introduction of the notion of capacity solution.

The concept of capacity solution was first introduced by Xu in [17] in the analysis
of a modified version of the evolution thermistor problem. He also applied this
concept to more general settings where weaker assumptions [18] or mixed boundary
conditions [19] are considered. Later on, it has also been used by other authors
in different situations [6, 12, 14]. For instance, in [6] the authors analyzed the
existence of a capacity solution for the evolution thermistor problem in W 1,p for
p ≥ 2. Moussa et al. [12] studied this system in isotropic Orlicz-Sobolev spaces,
whereas Talbi et al.
citeHajar considered the anisotropic case with polynomial growth with respect to
the variable ζ.

The goals of this paper are twofold. First, we analyze the existence of a capacity
solution to (1.1) for arbitrary and different growths of the functions ai(x, s, ζ),
i = 1, . . . , d (in particular, some or all of the functions Mi, i = 1, . . . , d, do
not need to satisfy the ∆2-condition). For instance, in d = 2, we may have
a1(x, s, ζ) = |ζ|p−2ζ and a2(x, s, ζ) = 2β exp(βζ2)ζ where the corresponding N -
functions are given by M1(ζ) = |ζ|p/p and M2(ζ) = exp(βζ2) − 1, respectively; in
this case, M1 satisfies the ∆2-condition whereas M2 does not. Secondly, we de-
scribe a numerical algorithm for the approximation of the solution to problem (1.1)
based on the least squares method combined with a conjugate gradient technique
[4, 5]. Though some numerical simulations have led to good results by using this
algorithm, the numerical resolution of (1.1) for arbitrary functions ai remains a
challenge; this is related to the machine precision. Indeed, since this algorithm
generates a sequence of approximate solutions for the temperature, say (umh ), by

means of descent directions, (zmh ), that is um+1
h = umh − λmzmh , where λm > 0 is an

optimal value, it may occurs that λm becomes smaller than the machine precision
which may result in an underflow situation. Consequently, λm is taken to be zero
inside the machine, and the algorithm would produce um+1

h = umh after that value
of m.

Notice that the numerical resolution of system (1.1) is an important issue. These
numerical simulations may yield a very useful information while designing a ther-
mistor for certain specific purposes. In this sense, knowing the steady-state tem-
perature distribution inside the thermistor for a given geometry Ω and potential
ϕ0 is crucial, although this is only a first step in this analysis (we also want to
know the whole history of both, temperature and potential, from a known initial



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 3

temperature up to being close enough to the steady-state by solving the evolution
problem).

This article is organized as follows. In Section 2 we introduce some notation,
concepts and functional spaces together with certain technical results that will be
needed along this paper. Section 3 states the assumptions on data and introduces
the concept of a capacity solution to (1.1) in the framework given in the previous
section. In Section 4 we present an existence theorem along with its proof. Finally,
Section 5 is devoted to the description of a numerical algorithm for the approximate
solution to problem (1.1) for certain choices of the functions ai, 1 ≤ i ≤ d, including
some numerical results obtained by the implementation of this algorithm.

2. Preliminaries

We begin by recalling some definitions and properties of Orlicz spaces [1, 11] and
then we introduce the anisotropic Orlicz-Sobolev spaces.

2.1. N-functions. The basic concept in an Orlicz normed space is that of N -
function.

Definition 2.1. A function M : R → R is called an N -function if it fulfills the
following conditions:

(i) M is convex in R: M(λs1 + (1 − λ)s2) ≤ λM(s1) + (1 − λ)M(s2), for all
s1, s2 ∈ R and for all λ ∈ [0, 1].

(ii) M is an even function: M(s) = M(−s) for all s ∈ R.
(iii) M(0) = 0 and M(s) > 0 for all s ∈ R.

(iv) M(s)
s → 0 as s→ 0 and M(s)

s → +∞ as s→ +∞.

An N -function M is said to satisfy the ∆2-condition for all s ∈ R if, for some
k > 0,

M(2s) ≤ kM(s) or all s ∈ R.

We say that M satisfies the ∆2-condition for s large if there exist s0 ≥ 0 and k > 0
such that

M(2s) ≤ kM(s) for all s ≥ s0.

An equivalent definition [11] of an N -function is a function M that admits the
representation

M(s) =

∫ |s|
0

m(σ) dσ,

where m : R+ → R+ is a non-decreasing and right-continuous function, m(s) > 0 for
all s > 0 and m(s)→ +∞ as s→ +∞. For an N -function M , the complementary
or conjugate is defined by

M̄(s) =

∫ |s|
0

m̄(σ) dσ,

where m̄ : R+ → R+ is given by m̄(t) = sup{s : m(s) ≤ t}.
We have the Young’s inequality

|ts| ≤M(t) + M̄(s) for all t, s ∈ R.



4 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

Let Ω be an open set in Rd and d ∈ N. The Orlicz class LM (Ω) (resp. the Orlicz
space LM (Ω)) is defined as the set of (equivalence classes of) real-valued Lebesgue
measurable functions u in Ω such that∫

Ω

M(u(x)) dx < +∞ (resp.

∫
Ω

M
(u(x)

λ

)
dx < +∞ for some λ > 0).

Notice that LM (Ω) is a Banach space under the so-called Luxemburg norm

‖u‖M = inf
{
λ > 0 :

∫
Ω

M
(u(x)

λ

)
dx ≤ 1

}
,

and LM (Ω) is a convex subset of LM (Ω). Indeed, LM (Ω) is the linear hull of
LM (Ω). The closure in LM (Ω) of the set of bounded measurable functions with
compact support in Ω̄ is denoted by EM (Ω). The equality EM (Ω) = LM (Ω) holds if
and only if M satisfies the ∆2-condition, for all s or for s large according to whether
Ω has infinite measure or not. The dual of EM (Ω) can be identified with LM̄ (Ω)
by means of the duality pairing

∫
Ω
u(x)v(x) dx, and the dual norm on LM̄ (Ω) is

equivalent to ‖ · ‖M̄ .
In LM (Ω) we define the Orlicz norm

‖u‖(M) = sup

∫
Ω

u(x)v(x) dx (2.1)

where the supremum is taken over all v ∈ EM̄(Ω) such that ‖v‖M̄ ≤ 1. It turns out

that the norms ‖ · ‖M and ‖ · ‖(M) are equivalent. In fact, it can be shown that

‖u‖M ≤ ‖u‖(M) ≤ 2‖u‖M for all u ∈ LM (Ω). (2.2)

Also, the Hölder inequality holds∫
Ω

|u(x)v(x)|dx ≤ ‖u‖M‖v‖(M̄) for all u ∈ LM (Ω) and v ∈ LM̄ (Ω),

and by (2.2)∫
Ω

|u(x)v(x)|dx ≤ 2‖u‖M‖v‖M̄ for all u ∈ LM (Ω) and v ∈ LM̄ (Ω).

In particular, if Ω has finite measure, Hölder’s inequality yields the continuous
inclusion LM (Ω) ⊂ L1(Ω).

An important inequality in LM (Ω) is the following:∫
Ω

M(u(x)) dx ≤ ‖u‖(M) for all u ∈ LM (Ω) such that ‖u‖(M) ≤ 1, (2.3)

wherefrom we readily deduce∫
Ω

M
( u(x)

‖u‖(M)

)
dx ≤ 1 for all u ∈ LM (Ω) \ {0}. (2.4)

Definition 2.2. We say that (un) ⊂ LM (Ω) converges to u ∈ LM (Ω) for the
modular convergence in LM (Ω) if, for some λ > 0, one has∫

Ω

M
(un(x)− u(x)

λ

)
dx→ 0 as n→∞.

Modular convergence is weaker than the convergence in the norm of LM (Ω).
However, it is enough to our purposes. The next result tells us that the modular
convergence in LM implies the convergence in the weak-∗ topology σ(LM , LM̄ ).



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 5

Lemma 2.3 ([3, 7]). Let (un) ⊂ LM (Ω), u ∈ LM (Ω) and v ∈ LM̄ (Ω) such that
un → u with respect to the modular convergence. Then

(1) unv → uv strongly in L1(Ω). In particular,
∫

Ω
unv →

∫
Ω
uv.

(2) Furthermore, if (vn) ⊂ LM̄ (Ω) is such that vn → v with respect to the
modular convergence, then unvn → uv strongly in L1(Ω).

2.2. Anisotropic Orlicz-Sobolev spaces. Let Ω be an open subset of Rd, and
Mi be an N -function for each i = 1, . . . , d. We write M = (M1, . . . ,Md), M̄ =
(M̄1, . . . , M̄d). The anisotropic Orlicz space LM(Ω) (respectively, EM(Ω)) is de-
fined by

LM(Ω) =

d∏
i=1

LMi
(Ω)

(
respectively, EM(Ω) =

d∏
i=1

EMi
(Ω)
)
,

endowed with the norm

‖u‖ =

d∑
i=1

‖ui‖Mi , (2.5)

To introduce the anisotropic Orlicz-Sobolev spaces it will be interesting to define
the function

M0(s) = min
1≤i≤d

Mi(s) . (2.6)

Remark 2.4. It is easy to check that:

(i) The function M0 is an N -function.
(ii) The embedding LMi

(Ω) ↪→ LM0
(Ω) is continuous for each i ∈ {1, . . . , d}.

The anisotropic Orlicz-Sobolev spaces are defined by

W 1LM(Ω) = {u ∈ LM0(Ω) : ∂iu ∈ LMi(Ω), i = 1, . . . , d},
W 1EM(Ω) = {u ∈ EM0

(Ω) : ∂iu ∈ EMi
(Ω), i = 1, . . . , d},

which are Banach spaces under the norm

‖u‖1,M = ‖u‖M0
+

d∑
i=1

‖∂iu‖Mi
. (2.7)

Both spaces, W 1LM(Ω) and W 1EM(Ω), can be identified as subspaces of the

product space Π = LM0
(Ω) × LM(Ω). Then, the predual space of Π, Π̂, is

Π̂ = EM̄0
(Ω)×EM̄(Ω). We will use the weak-∗ topology σ(Π, Π̂). Let D(Ω) be the

space of functions in C∞(Ω) with compact support in Ω. The space W 1
0EM(Ω)

is defined as the (norm) closure of the space D(Ω) in W 1EM(Ω), and the space

W 1
0LM(Ω) as the σ(Π, Π̂)-closure of D(Ω) in W 1LM(Ω).

Lemma 2.5 ([12]). Let Ω be a bounded and open set in Rd. Assume that mi(t) ≥ t
for all t ≥ 0 and all i = 1, . . . , d. Then the following continuous embeddings hold
for i = 1, . . . , d:

LMi(Ω) ↪→ L2(Ω) ↪→ LM̄i
(Ω).

In particular, W 1
0LM(Ω) ↪→ H1

0 (Ω) and H−1(Ω) ↪→W−1LM̄(Ω).

Remark 2.6. Assume that, for each i = 1, . . . , d, one has mi(t) ≥ t for all t ≥ 0.
Then ∫

Ω

v2 dx ≤ 2

∫
Ω

Mi(v) dx, for all v ∈ LMi
(Ω). (2.8)



6 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

Theorem 2.7 ([16]). Let Ω ⊂ Rd be an open and bounded set with locally Lipschitz
boundary. Then the embedding W 1LM (Ω) ↪→ EM (Ω) is compact. Furthermore,
the compact imbedding W 1

0LM (Ω) ↪→ EM (Ω) holds without the locally Lipschitz
boundary assumption.

Corollary 2.8. Let Ω be an open and bounded set in Rd and M0 the N -function
defined in (2.6). Then, the embedding W 1

0LM(Ω) ↪→ EM0
(Ω) is compact.

Poincaré’s inequality in W 1
0LM(Ω) also holds.

Lemma 2.9 ([7]). Let Ω ⊂ Rd be an open and bounded set. Then, there exist
constants κ0 and κ1 = κ1(Ω) such that∫

Ω

M0(u) dx ≤ κ0

d∑
i=1

∫
Ω

Mi(κ1∂ui) dx for all u ∈W 1
0LM(Ω).

Corollary 2.10. The seminorm u ∈ W 1LM(Ω) 7→
∑d
i=1 ‖∂iu‖Mi

is a norm in
W 1

0LM(Ω) and it is equivalent to the norm ‖ · ‖1,M given in 2.7.

Since the elements of the space W 1
0LM(Ω) have been defined as the weak-∗

limit of convergent sequences in D(Ω), the following result states that, for certain
domains Ω, D(Ω) is ‘dense’ in W 1

0LM(Ω) with respect to the modular convergence
as well.

Definition 2.11. A bounded domain Ω ⊂ Rd is said to satisfy the segment prop-
erty, if there exist a locally finite open covering {Ui} of ∂Ω and corresponding
vectors {yi} ⊂ Rd such that for all x ∈ Ω̄∩Ui and any µ ∈ (0, 1) one has x+µyi ∈ Ω.

Lemma 2.12. Let Ω ⊂ Rd be an open and bounded set satisfying the segment
property and u ∈ W 1

0LM(Ω). Then there exists a sequence (un) ⊂ D(Ω) such that
un → u with respect to the modular convergence in W 1LM(Ω); that is, there exists
λ > 0 such that∫

Ω

M0

(
(un − u)/λ

)
+

d∑
i=1

∫
Ω

Mi

(
(∂iun − ∂iu)/λ

)
→ 0 as n→∞.

The proof of the above lemma is a straightforward adaptation of [8, Theorem 4]
for isotropic Orlicz-Sobolev spaces.

Finally, we introduce the following dual spaces

W−1LM̄(Ω) =
{
f ∈ D ′(Ω) : f =

d∑
i=1

∂ifi with fi ∈ LM̄i
(Ω), for all i, 1 ≤ i ≤ d}

W−1EM̄(Ω) = {f ∈ D ′(Ω) : f =

d∑
i=1

∂ifi with fi ∈ EM̄i
(Ω), for all i, 1 ≤ i ≤ d}

These spaces are equipped by their usual quotient norms.

3. Essential assumptions and main result

From this point on we will assume that Ω ⊂ Rd is an open and bounded set
satisfying the segment property. We now state the assumptions on the differential
operator in divergence form given by A : W 1

0LM(Ω)→W−1LM̄(Ω)

A(u) =

d∑
i=1

∂i
(
ai(x, u, ∂iu)

)
.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 7

(A1) For each i = 1, . . . , d the function ai : Ω × R × R → R, ai = ai(x, s, ζ) is a
Carathéodory function (measurable in x for all (s, ζ) and continuous with
respect to (s, ζ) almost everywhere (a.e.) in Ω).

(A2) There exist N -functions Mi, i = 1, . . . , d, a function a0 ∈ EM̄i
and positive

constants k and γ such that a.e. in Ω, for all (s, ζ) ∈ R2 and for i = 1, . . . , d.

|ai(x, s, ζ)| ≤ γ
(
a0(x) + M̄−1

i (Mi(k|ζ|))
)
. (3.1)

(A3) For a.e. in Ω, for all s ∈ R, and all ζ, ζ′ ∈ Rd with ζ 6= ζ′, ζ = (ζ1, . . . , ζd),
ζ′ = (ζ ′1, . . . , ζ

′
d) we have

d∑
i=1

[ai(x, s, ζi)− ai(x, s, ζ ′i)](ζi − ζ ′i) > 0 (3.2)

(A4) There exists λ0 > 0 such that a.e. in Ω, for all s ∈ R, and all (ζ1, . . . , ζd) ∈
Rd we have

d∑
i=1

ai(x, s, ζi)ζi ≥ α
d∑
i=1

Mi

( ζi
λ0

)
(3.3)

(A5) Let B ⊂ R be a bounded set. Then for a.e. in Ω, for all s ∈ B, and
ζ′ = (ζ ′1, . . . , ζ

′
d) ∈ Rd

d∑
i=1

[ai(x, s, ζi)− ai(x, s, ζ ′i)](ζi − ζ ′i)→ +∞ (3.4)

as |ζ| → +∞, ζ = (ζ1, . . . , ζd) ∈ Rd, uniformly in s.
(A6) For all i = 1, . . . , d, we have

ai(x, s, 0) = 0 (3.5)

(A7) ρ ∈ C (R) and there exists ρ̄ ∈ R such that for all s ∈ R,

0 < ρ(s) ≤ ρ̄ (3.6)

(A8) ϕ0 ∈ H1(Ω) ∩ L∞(Ω).

Now we are ready to state the definition of a capacity solution to problem (1.1).

Definition 3.1. A triplet (u, ϕ,Φ) is called a capacity solution to problem (1.1) if
it satisfies the following conditions:

(A9) u ∈ W 1
0LM(Ω), ai(·, u, ∂iu) ∈ LM̄i

(Ω) for 1 ≤ i ≤ d, ϕ ∈ L∞(Ω) and

Φ ∈ L2(Ω)d.
(A10) (u, ϕ,Φ) verifies the system of differential equations in the sense of distri-

butions

−
d∑
i=1

∂i
(
ai(x, u, ∂iu)

)
= div(ϕΦ) in Ω,

div Φ = 0 in Ω.

(A11) For every S ∈ C 1
0 (R) (that is, S ∈ C 1(R) and has compact support in R),

one has S(u)ϕ− S(0)ϕ0 ∈ H1
0 (Ω), and

S(u)Φ = ρ(u)[∇(S(u)ϕ)− ϕ∇S(u)].



8 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

Remark 3.2. Notice that the concept of capacity solution involves a third com-
ponent, namely, Φ ∈ L2(Ω)d. Where this vector field comes from? The original
problem only has two unknowns so that there must be a relationship between Φ and
(u, ϕ). Indeed, this is true and the relationship is given by the condition expressed
in (A11). In fact, assume that u ∈ L∞(Ω) and take S ∈ C 1

0 (R) such that S ≡ 1
in the interval [−‖u‖∞, ‖u‖∞]. Then, we deduce that ϕ ∈ H1(Ω), Φ = ρ(u)∇ϕ
and div(ϕΦ) = ρ(u)|∇u|2. Consequently, bounded capacity solutions are weak so-
lutions. In the general case where u 6∈ L∞(Ω) the expression given in (A11) allows
us to define the gradient of ϕ pointwise almost everywhere from the identity

Φχ{|u|<K} = ρ(u)∇ϕχ{|u|<K} for any K > 0.

4. An existence result

In this section we establish the main result of this article.

Theorem 4.1. Under the assumptions (A1)–(A8), system (1.1) admits a capacity
solution.

To prove the main result, we will need to show the existence of a weak solution
to a similar problem but under a more restrictive assumption, namely,

ρ ∈ C (R) and there exist ρ1 and ρ2 ∈ R such that

0 < ρ1 ≤ ρ(s) ≤ ρ2, for all s ∈ R.
(4.1)

Theorem 4.2. Assume (A1)–(A6), (A8) and (4.1) hold. Then, there exists a
weak solution (u, ϕ) to (1.1), that is u ∈ W 1

0LM(Ω), ai(·, u, ∂iu) ∈ LM̄i
(Ω), for all

1 ≤ i ≤ d, ϕ− ϕ0 ∈ H1
0 (Ω) ∩ L∞(Ω), and

d∑
i=1

∫
Ω

ai(x, u, ∂iu)∂iφ =

∫
Ω

ρ(u)|∇ϕ|2φ, for all φ ∈W 1
0LM(Ω),∫

Ω

ρ(u)∇ϕ∇ψ = 0, for all ψ ∈ H1
0 (Ω).

Proof. To prove the existence of a weak solution, Schauder’s fixed point theorem
will be applied together with a result on the existence and uniqueness of a weak
solution to an elliptic equation.

For a function ω ∈ EM0
(Ω) we consider the elliptic problem

div(ρ(ω)∇ϕ) = 0 in Ω,

ϕ = ϕ0 on ∂Ω.
(4.2)

Thanks to Lax-Milgram’s theorem, (4.2) has a unique solution ϕ ∈ H1(Ω). In this
case, from the maximum principle we also have ϕ ∈ L∞(Ω) and

‖ϕ‖L∞(Ω) ≤ ‖ϕ0‖L∞(Ω). (4.3)

Using ϕ− ϕ0 ∈ H1
0 (Ω) as a test function in (4.2) we obtain∫

Ω

ρ(ω)|∇ϕ|2 =

∫
Ω

ρ(ω)∇ϕ∇ϕ0,

hence

ρ1

∫
Ω

|∇ϕ|2 dx ≤
∫

Ω

ρ(ω)|∇ϕ‖∇ϕ0|dx ≤ ρ2

∫
Ω

|∇ϕ||∇ϕ0|dx.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 9

By the Cauchy-Schwarz inequality, we obtain∫
Ω

|∇ϕ|2 dx ≤ C(ρ1, ρ2, ϕ0) = C. (4.4)

Thanks to the elliptic equation for ϕ, the term ρ(ω)|∇ϕ|2 also belongs to the space
H−1(Ω). Indeed, let φ ∈ D(Ω) and take ψ = φϕ as a test function in (4.2). We
have ∫

Ω

ρ(ω)∇ϕ∇(φϕ) dx = 0,

that is∫
Ω

ρ(ω)|∇ϕ|2φdx = −
∫

Ω

ρ(ω)ϕ∇ϕ∇φ dx = 〈div(ρ(ω)ϕ∇ϕ), φ〉D′(Ω),D(Ω) .

This means that

ρ(ω)|∇ϕ|2 = div(ρ(ω)ϕ∇ϕ) in D ′(Ω). (4.5)

Since ρ(ω)ϕ∇ϕ ∈ L2(Ω)N we finally deduce the regularity

div(ρ(ω)ϕ∇ϕ) ∈ H−1(Ω).

Now we consider the elliptic problem

−
d∑
i=1

∂i
(
ai(x, ω, ∂iu)

)
= div(ρ(ω)ϕ∇ϕ) in Ω,

u = 0 on ∂Ω.

(4.6)

The variational formulation of this elliptic equation is

u ∈W 1
0LM(Ω), ai(·, ω, ∂iu) ∈ LM̄i

(Ω), for i = 1, . . . , d,

d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iφ = −
∫

Ω

ρ(ω)ϕ∇ϕ∇φ, for all φ ∈W 1
0LM(Ω).

(4.7)

We have div(ρ(ω)ϕ∇ϕ) ∈ H−1(Ω) ↪→ W−1LM̄(Ω) because of Lemma 2.5. The
proof of the existence and uniqueness of solution to (4.7) is a straightforward ap-
plication of the result given in [15].

Now, we show that ∂iu/λ0 ∈ LMi
(Ω) for i = 1, . . . , d, where λ0 is the constant

appearing in (3.3), and the estimates

d∑
i=1

∫
Ω

Mi(∂iu/λ0) dx ≤ C(ϕ0, α, λ0, ρ2) = C0, (4.8)

d∑
i=1

‖ai(·, ω, ∂iu)‖M̄i
≤ C1. (4.9)

Indeed, let λ > 0 such that |∂iu|/λ ∈ LMi(Ω) for i = 1, . . . , d. Since ϕ ∈ H1(Ω) ⊂
W 1LM̄(Ω), there exists µ > 0 such that 1

µρ2‖ϕ0‖L∞(Ω)∂iϕ ∈ LM̄i
(Ω) for i =

1, . . . , d. By taking φ = u as a test function in (4.7), from (3.3), (3.5), (4.1), (4.3)
and Young’s inequality, we obtain

α

λµ

d∑
i=1

∫
Ω

Mi(∂iu/λ0) dx



10 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

≤ 1

λµ

∫
Ω

d∑
i=1

ai(x, ω, ∂iu)∂iudx

≤
d∑
i=1

∫
Ω

M̄i(ρ2‖ϕ0‖L∞(Ω)∂iϕ/µ)dx+

d∑
i=1

∫
Ω

Mi(∂iu/λ)dx <∞

This shows that ∂iu/λ0 ∈ LMi
(Ω) for i = 1, . . . , d. To obtain (4.8), by using

Young’s inequality and the estimate (2.8), we deduce

α

d∑
i=1

∫
Ω

Mi(∂iu/λ0) dx ≤
d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iudx

≤
d∑
i=1

(λ2
0

α
ρ2

2‖ϕ0‖2L∞(Ω)

∫
Ω

|∇ϕ|2 +
α

4

∫
Ω

|∂iu/λ0|2
)

≤ C∗(ϕ0, α, λ0, ρ2) +
α

2

d∑
i=1

∫
Ω

Mi(∂iu/λ0)

and thus (4.8) holds for C0 = 2C∗(ϕ0, α, λ0, ρ2)/α. To obtain (4.9), first notice
that from the previous two inequalities we obtain

d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iudx ≤ αC0. (4.10)

Then, because of (3.2), for any φ ∈W 1
0EM(Ω) such that

∑d
i=1 ‖∂iφ‖(Mi) = 1/(k+1)

we have

0 ≤
d∑
i=1

∫
Ω

(
ai(x, ω, ∂iu)− ai(x, ω, ∂iφ)

)
(∂iu− ∂iφ) dx.

Owing to (4.10) and Young’s inequality, we deduce

d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iφ dx

≤
d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iudx−
d∑
i=1

∫
Ω

ai(x, ω, ∂iφ)(∂iu− ∂iφ) dx

≤ αC0 +

d∑
i=1

∫
Ω

|ai(x, ω, ∂iφ)∂iu|dx+

d∑
i=1

∫
Ω

ai(x, ω, ∂iφ)∂iφ dx

≤ αC0 + 2γλ0

d∑
i=1

∫
Ω

[
M̄i

(ai(x, ω, ∂iφ)

2γ

)
+Mi(∂iu/λ0)

]
dx

+ 2γ

d∑
i=1

∫
Ω

[
M̄i

(ai(x, ω, ∂iφ)

2γ

)
+Mi(|∂iφ|)

]
dx,

where γ is the constant appearing in (3.1). Since

d∑
i=1

M̄i

(ai(x, ω, ∂iφ)

2γ

)
≤ 1

2

d∑
i=1

(
M̄i(a0(x)) +Mi(k|∂iφ|)

)
a.e. in Ω,



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 11

using (2.3) we obtain

d∑
i=1

∫
Ω

M̄i

(ai(x, ω, ∂iφ)

2γ

)
dx ≤ 1

2

d∑
i=1

∫
Ω

M̄i(a0(x)) dx+
1

2
= C2.

Notice that C2 does not depend on ω. Therefore, gathering all these estimates, we

deduce that for all φ ∈W 1
0EM(Ω) such that

∑d
i=1 ‖∂iφ‖(Mi) = 1/(k + 1) we have

d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iφ dx ≤ C1,

from which, by considering the dual norm on LM̄i
(Ω), for each i = 1, . . . , d, we

obtain the estimate (4.9).
Now we introduce the operator G : ω ∈ EM0(Ω)→ G(ω) = u ∈W 1

0LM(Ω), with
u being the unique solution to (4.7). Our strategy is to show that G satisfies the
conditions of Schauder’s fixed point theorem. From Corollary 2.8, W 1

0LM(Ω) ↪→
EM0

(Ω) with compact embedding. Consequently, G maps EM0
(Ω) into itself and,

due to the estimates (4.8), G is a compact operator. Moreover, from Corollary 2.10
and (4.8) we have, for R > 0 large enough G(BR) ⊂ BR where BR = {v ∈ EM0(Ω) :
‖v‖Mi ≤ R, for i = 1, . . . , d}.

To complete the proof, it remains to show that G is a continuous operator.
Indeed, let (ωn) ⊂ BR such that ωn → ω strongly in EM0

(Ω) and consider the
corresponding functions to ωn, that is, un = G(ωn), ϕn, u = G(ω), and ϕ. Since the
injection EM0(Ω) ⊂ L2(Ω) is continuous, we can assume that for a subsequence, still
denoted in the same way, it is ωn → ω a.e. in Ω. Then, putting Fn = ρ(ωn)ϕn∇ϕn
and F = ρ(ω)ϕ∇ϕ, it is easy to check that Fn → F strongly in L2(Ω)d. On the
other hand, since (un) ⊂ W 1

0LM(Ω) is bounded in this space, from Corollary 2.8,
there exist a subsequence, still denoted in the same way, and a function U ∈ EM0

(Ω)
such that un → U strongly in EM0

(Ω) and a.e. in Ω. We need to prove that u = U .
To do so, we first need the following result.

Lemma 4.3. There exists a subsequence, still denoted in the same way, such that
∇un → ∇U a.e. in Ω.

This result is a special case of Proposition 4.4 below. For the proof of Lemma 4.3,
one may repeat the arguments given in the proof of Proposition 4.4.

As a consequence of Lemma 4.3, there exists a subsequence, still denoted in the
same way such that

ai(x, ωn, ∂iun)→ ai(x, ω, ∂iU) a.e. in Ω, for i = 1, . . . , d.

On the other hand, from the estimate (4.9), (ai(x, ωn, ∂iun)) ⊂ LM0
(Ω) is bounded

in this space and thus, there exist a subsequence, still denoted in the same way,
and functions Φi ∈ LM̄i

(Ω), i = 1, . . . , d, such that

ai(x, ωn, ∂iun)→ Φi weak−∗ in LM̄i
(Ω),

and therefore Φi = ai(x, ω, ∂iU) for i = 1, . . . , d.
By using all these convergences and passing to the limit in the equation of un,

that is,
d∑
i=1

∫
Ω

ai(x, ωn, ∂iun)∂iφ dx = −
∫

Ω

ρ(ωn)ϕn∇ϕn∇φ dx,



12 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

we obtain
d∑
i=1

∫
Ω

ai(x, ω, ∂iU)∂iφdx = −
∫

Ω

ρ(ω)ϕ∇ϕ∇φ dx,

Since u = G(ω), we also have

d∑
i=1

∫
Ω

ai(x, ω, ∂iu)∂iφ dx = −
∫

Ω

ρ(ω)ϕ∇ϕ∇φdx,

Thus,

d∑
i=1

∫
Ω

[ai(x, ω, ∂iu)− ai(x, ω, ∂iU)]∂iφ dx = 0 for all φ ∈W 1
0LM(Ω).

By taking φ = u−U and using (3.2) we deduce that u = U . In particular, it is the
whole sequence (un) that converges to u in EM0(Ω) and this shows the continuity
of the operator G. This completes the proof of Theorem 4.2. �

Proof of Theorem 4.1. The proof is divided into several steps. We start by
introducing a sequence of approximate problems and derive some a priori estimates
for the respective solutions to these approximate problems, then we show two in-
termediate results, namely, the strong convergence in L1(Ω), up to a subsequence,
of both ∂iun for i = 1, . . . , d and ϕn, where (un, ϕn) is a weak solution to the
n-th approximate problem of (1.1). The passing to the limit in the approximate
problems will yield the main result.

Step 1. Approximate problems and a priori estimates. We define the truncation
function at height K > 0, TK , by

TK(s) = s if |s| ≤ K, TK(s) = K
s

|s|
if |s| > K,

and we introduce the following regularization of the data, for every n ∈ N,

ρn(s) = ρ(s) +
1

n
, (4.11)

The n-th approximate problem is

−
d∑
i=1

∂i
(
ai(x, un, ∂iun)

)
= ρn(un)|∇ϕn|2 in Ω, (4.12)

div(ρn(un)∇ϕn) = 0 in Ω, (4.13)

un = 0, ϕn = ϕ0 on ∂Ω, (4.14)

In view of (3.6), we have that

n−1 ≤ ρn(s) ≤ ρ̄+ 1 = ρ3, for all s ∈ R. (4.15)

Thus, we can apply Theorem 4.2 to deduce the existence of a weak solution (un, ϕn)
to (4.12)-(4.14).

By the maximum principle we have

‖ϕn‖L∞(Ω) ≤ ‖ϕ0‖L∞(Ω), (4.16)

therefore, there exists a function ϕ ∈ L∞(Ω) and a subsequence, still denoted in
the same way, such that

ϕn → ϕ weakly-∗ in L∞(Ω). (4.17)



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 13

Multiplying (4.13) by ϕn − ϕ0 ∈ H1
0 (Ω) and integrating in Ω we obtain∫

Ω

ρn(un)|∇ϕn|2 =

∫
Ω

ρn(un)∇ϕn∇ϕ0 dx

≤
(∫

Ω

ρn(un)|∇ϕn|2
)1/2(∫

Ω

ρn(un)|∇ϕ0|2
)1/2

≤
(∫

Ω

ρn(un)|∇ϕn|2
)1/2(

ρ3

∫
Ω

|∇ϕ0|2
)1/2

,

hence ∫
Ω

ρn(un)|∇ϕn|2 dx ≤ C1, for all n ≥ 1, (4.18)

where C1 = C1(ρ3, ‖ϕ0‖H1(Ω)). Consequently, (ρn(un)∇ϕn) is bounded in L2(Ω).

Thus, there exists a function Φ ∈ L2(Ω)d and a subsequence, still denoted in the
same way, such that

ρn(un)∇ϕn → Φ weakly in L2(Ω)d. (4.19)

This weak limit function Φ ∈ L2(Ω)d is in fact the third component of the triplet
appearing in the Definition 3.1 of a capacity solution.

Taking un as a test function in (4.12), we obtain

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iun dx = −
d∑
i=1

∫
Ω

ρn(un)ϕn∂iϕn∂iun dx. (4.20)

Since un ∈ W 1
0LM(Ω) and ϕ ∈ H1(Ω) ⊂ W 1LM̄(Ω), there exist λn > 0 such that

for i = 1, . . . , d one has ∂iun/λn ∈ LMi
(Ω). By (2.8), (3.3), (3.5), (4.15), (4.16),

and Young’s inequality, we obtain

α

d∑
i=1

∫
Ω

Mi(∂iun/λ0) dx

≤
d∑
i=1

∫
Ω

ρn(un)|ϕn‖∇ϕn||∂iun|dx

≤
d∑
i=1

∫
Ω

√
ρ3‖ϕ0‖L∞(Ω)

√
ρn(un)|∇ϕn||∂iun|dx

≤
d∑
i=1

(ρ3‖ϕ0‖2L∞(Ω)λ
2
n

α

∫
Ω

ρn(un)|∇ϕn(un)|2 +
α

2

∫
Ω

Mi(∂iun/λn), dx
)

which implies that ∂iun/λ0 ∈ LMi(Ω). In particular, we may take λn = λ0 for all
n ≥ 1 and then

d∑
i=1

∫
Ω

Mi(∂iun/λ0) dx ≤ C, (4.21)

where C is a positive constant not depending on n. Thus, from Corollary 2.10
and (4.21), the sequence (un) is bounded in W 1

0LM(Ω), and since the embedding
W 1

0LM(Ω) ↪→ EM0
(Ω) is compact, there exist a subsequence of (un), still denoted

in the same way, and a function u ∈ EM0
(Ω) such that

un → u strongly in EM0
(Ω) and a.e. in Ω. (4.22)



14 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

On the other hand, since (∇un) ⊂ LM(Ω) is bounded and using (4.22), we also
have that, up to a subsequence,

∇un → ∇u weak-∗ in LM(Ω) (4.23)

Now let φ ∈ W 1
0EM(Ω) be such that

∑d
i=1 ‖∂iφ‖(Mi) = 1/(k + 1). In view of

the monotonicity of a = (a1, . . . , ad), we easily find that

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iφ

≤
d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iun −
d∑
i=1

∫
Ω

ai(x, un, ∂iφ)(∂iun − ∂iφ)

≤ C +

d∑
i=1

∫
Ω

|ai(x, un, ∂iφ)∂iun|+
d∑
i=1

∫
Ω

ai(x, un, ∂iφ)∂iφ,

(4.24)

The last integrals in (4.24) are bounded with respect to n. Indeed, for the first one,
owing to Young’s inequality

d∑
i=1

∫
Ω

|ai(x, un, ∂iφ)∂iun| ≤ 2γλ0

d∑
i=1

∫
Ω

[
M̄i

(ai(x, un, ∂iφ)

2γ

)
+Mi(∂iun/λ0)

]
,

by using (3.1) we have

2γ

d∑
i=1

M̄i

(ai(x, un, ∂iφ)

2γ

)
≤ γ

d∑
i=1

(
M̄i(a0(x)) +Mi(k∂iφ)

)
,

and thus
∑d
i=1

∫
Ω
|ai(x, un, ∂iφ)∂iun| ≤ C, for all n ≥ 1 and φ ∈ W 1

0EM(Ω) such

that
∑d
i=1 ‖∂iφ‖(Mi) = 1/(k + 1). In the same way, we can show that the second

integral in (4.24) is bounded. Gathering all these estimates, and using the dual
norm, it is easily deduced that, for i = 1, . . . , d,(

ai(x, un, ∂iun)
)

is bounded in LM̄i
(Ω). (4.25)

Thus, we have that, for a subsequence still denoted in the same way, there exists
δi ∈ LM̄i

(Ω), for each i = 1, . . . , d,, such that

ai(x, un, ∂iun)→ δi weak-∗ in LM̄i
(Ω). (4.26)

Step 2. Introduction of regularized sequences and the almost everywhere conver-
gence of the gradients. By using Lemma 2.12, there exists a sequence (vj) ⊂ D(Ω)
such that

(1) vj → u in W 1
0LM(Ω) for the modular convergence;

(2) vj → u and ∂ivj → ∂iu a.e. in Ω for i = 1, . . . , d.

We will establish the following result.

Proposition 4.4. Let (un, ϕn) be a solution of the approximate problem (4.12)-
(4.14). Then there exists a subsequence, still denoted in the same way, such that
for i = 1, . . . , d,

∂iun → ∂iu a.e. in Ω, (4.27)

as n tends to +∞.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 15

Proof. We denote by χjs and χs, respectively, the characteristic functions of the sets

Ωjs =
{
x ∈ Ω :

d∑
i=1

|∂iTK(vj)| ≤ s
}

and Ωs =
{
x ∈ Ω :

d∑
i=1

|∂iTK(u)| ≤ s
}
.

We denote by ε(n, j) and ε(n, j, s) any quantities such that

lim sup
j→∞

lim sup
n→∞

ε(n, j) = 0, lim sup
s→∞

lim sup
j→∞

lim sup
n→∞

ε(n, j, s) = 0.

For ν > 0 and i, j, n ≥ 1 using the test function ϕνn,j = Tν(un − TK(vj)) in (4.12),
we obtain

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iTν(un − TK(vj)) dx =

∫
Ω

ρn(un)|∇ϕn|2ϕνn,j dx. (4.28)

From (4.18) it follows that

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iTν(un − TK(vj)) dx ≤ C1ν. (4.29)

On the other hand
d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iTν(un − TK(vj)) dx

=

d∑
i=1

∫
{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂i(un − TK(vj)) dx

=

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂i(un − TK(vj)) dx

+

d∑
i=1

∫
{|un|≤K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂i(un − TK(vj)) dx

=

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj)) dx

+

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂iun dx

−
d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂iTK(vj) dx.

Then, using (3.3), we have

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iTν(un − TK(vj)) dx

≥
d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj))

−
d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂iTK(vj) dx.

(4.30)



16 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

Also, in the set {|un−TK(vj)| ≤ ν}, we have |un| ≤ |un−TK(vj)|+|TK(vj)| ≤ ν+K,
and thus, we can write

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, un, ∂iun)∂iTK(vj) dx

=

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, Tν+K(un), ∂iTν+K(un))∂iTK(vj) dx.

(4.31)

By (4.31), (4.30) becomes

d∑
i=1

∫
Ω

ai(x, un, ∂iun)∂iTν(un − TK(vj)) dx

≥
d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj))

−
d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, Tν+K(un), ∂iTν+K(un))∂iTK(vj).

(4.32)

We put

J1 =

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, Tν+K(un), ∂iTν+K(un))∂iTK(vj) dx.

Since
(
ai(x, TK+ν(un), ∂iTK+ν(un))

)
n

is bounded in LM̄i
(Ω), for each i = 1, . . . , d,

we obtain, for certain liK+ν ∈ LM̄i
(Ω), and up to a subsequence, that

ai(x, TK+ν(un), ∂iTK+ν(un))→ liK+ν weakly-∗ in LM̄i
(Ω).

Consequently,

d∑
i=1

∫
{|un|>K}∩{|un−TK(vj)|≤ν}

ai(x, Tν+K(un), Tν+K(un))∂iTK(vj) dx

→
d∑
i=1

∫
{|u|≥K}∩{|u−TK(vj)|≤ν}

liK+ν∂iTK(vj) dx

as n approaches infinity. Using Lemma 2.3, we obtain, as j tends to infinity, that

d∑
i=1

∫
{|u|≥K}∩{|u−TK(vj)|≤ν}

liK+ν∂iTK(vj)

→
d∑
i=1

∫
{|u|≥K}∩{|u−TK(u)|≤ν}

liK+ν∂iTK(u) = 0.

since ∂iTK(u) = 0 in the set {|u| ≥ K}. This implies

J1 = ε(n, j). (4.33)

Using (4.29) and (4.33) in (4.32), we obtain

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj)) dx

≤ C1ν + ε(n, j). (4.34)



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 17

On the other hand,

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj)) dx

=

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj)χ
s
j)

+

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(vj)χ
s
j − ∂iTK(vj))

= J2 + J3. (4.35)

The second integral, J3 tends to 0 as, first n, then j approach infinity. Indeed, since

ai(x, TK(un), ∂iTK(un)) ⇀ liK weakly-∗ in (LM̄i
(Ω))

and (
∂iTK(vj)χ

s
j − ∂iTK(vj)

)
χ{|TK(un)−TK(vj)|≤ν}

→
(
∂iTK(vj)χ

s
j − ∂iTK(vj)

)
χ{|TK(u)−TK(vj)|≤ν}

strongly in EM̄i
(Ω) as n→∞, for i = 1, . . . , d, it follows that

lim
n→∞

J3 =

d∑
i=1

∫
{|TK(u)−TK(vj)|≤ν}

liK · (∂iTK(vj)χ
s
j − ∂iTK(vj)) dx.

Finally letting j, then s, approach infinity, we deduce that

J3 = ε(n, j, s). (4.36)

Consequently, from (4.34), (4.35), and (4.36), we have

J2 =

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))(∂iTK(un)− ∂iTK(vj)χ
s
j)

≤ C1ν + ε(n, j, s). (4.37)

Let Sn be the non-negative expression

Sn =

d∑
i=1

(
ai(x, TK(un), ∂iTK(un))−ai(x, TK(un), ∂iTK(u))

)
(∂iTK(un)−∂iTK(u)),

and for each 0 < θ < 1, we write In,r =
∫

Ωr
Sθn dx. We have

In,r =

∫
Ωr

Sθnχ{|TK(un)−TK(vj)|≤ν} +

∫
Ωr

Sθnχ{|TK(un)−TK(vj)|>ν}. (4.38)

By using Hölder’s inequality, the second term of the right-side hand is bounded by(∫
Ωr

Sn dx
)θ
·
(∫

Ωr

χ{|TK(un)−TK(vj)|>ν} dx
)1−θ

.

Note that ∫
Ωr

Sn dx =

d∑
i=1

[ ∫
Ωr

ai(x, TK(un), ∂iTK(un))∂iTK(un) dx

−
∫

Ωr

ai(x, TK(un), ∂iTK(un))∂iTK(u) dx



18 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

+

∫
Ωr

ai(x, TK(un), ∂iTK(u))∂iTK(u) dx

−
∫

Ωr

ai(x, TK(un), ∂iTK(u))∂iTK(un) dx
]
.

Since, for each i = 1, . . . , d, the quantity
(
ai(x, TK(un), ∂iTK(un))

)
n

is bounded

in LM̄i
(Ω),

(
∂iTK(un)

)
n

is bounded in LMi
(Ω), and

(
ai(x, TK(un), ∂iTK(u))

)
n

is

bounded in EM̄i
(Ωr), it follows that (Sn) is bounded in L1(Ωr). Then there exists

a constant C3 > 0 such that∫
Ωr

Sθnχ{|TK(un)−TK(vj)|>ν} dx ≤ C3 meas{|TK(un)− TK(vj)| > ν}1−θ. (4.39)

Using again Hölder’s inequality, we deduce that∫
Ωr

Sθnχ{|TK(un)−TK(vj)|≤ν} dx

≤
(∫

Ωr

1 dx
)1−θ(∫

{|TK(un)−TK(vj)|≤ν}∩Ωr

Sn dx
)θ

≤ C4

(∫
{|TK(un)−TK(vj)|≤ν}∩Ωr

Sn dx
)θ
.

(4.40)

From (4.39) and (4.40), we obtain

In,r ≤ C3 meas{|TK(un)− TK(vj)| > ν}1−θ

+ C4

(∫
{|TK(un)−TK(vj)|≤ν}∩Ωr

Sn dx
)θ
.

(4.41)

Let s ≥ r > 0. We have∫
{|TK(un)−TK(vj)|≤ν}∩Ωr

Sn dx

≤
∫
{|TK(un)−TK(vj)|≤ν}∩Ωs

Sn dx

=

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}∩Ωs

[ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(u)χs)]

× [∂iTK(un)− ∂iTK(u)χs] dx

≤
d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

[ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(u)χs)]

× [∂iTK(un)− ∂iTK(u)χs] dx

=

d∑
i=1

[ ∫
{|TK(un)−TK(vj)|≤ν}

[ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(vj)χ
s
j)]

× [∂iTK(un)− ∂iTK(vj)χ
s
j ] dx

+

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un)) · [∂iTK(vj)χ
s
j − ∂iTK(u)χs] dx

+

∫
{|TK(un)−TK(vj)|≤ν}

[ai(x, TK(un), ∂iTK(vj)χ
s
j)− ai(x, TK(un), ∂iTK(u)χs)]



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 19

× ∂iTK(un) dx

−
∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(vj)χ
s
j) · ∂iTK(vj)χ

s
j dx

+

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(u)χs) · ∂iTK(u)χs dx
]

=

d∑
i=1

(
Ii,1 + Ii,2 + Ii,3 + Ii,4 + Ii,5

)
= I1 + I2 + I3 + I4 + I5.

where Ik =
∑d
i=1 Ii,k for 1 ≤ k ≤ 5. We will take the limit first in n, then in j and

s, as they tend to infinity in these last five integrals.

Estimate for I1. We rewrite

I1 =

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

[ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(vj)χ
s
j)]

× [∂iTK(un)− ∂iTK(vj)χ
s
j ]

=

d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(un))[∂iTK(un)− ∂iTK(vj)χ
s
j ]

−
d∑
i=1

∫
{|TK(un)−TK(vj)|≤ν}

ai(x, TK(un), ∂iTK(vj)χ
s
j)[∂iTK(un)− ∂iTK(vj)χ

s
j ]

= J2 − J4.

The estimate for J2 is given by (4.37). As J4 is concerned, we may repeat the same
arguments above as we did for J3 to obtain that J4 = ε(n, j, s). Therefore,

I1 ≤ Cν + ε(n, j, s) (4.42)

Estimate for I2. Since
(
ai(x, TK(un), ∂iTK(un))

)
n

converges to liK weakly-∗ in

LM̄i
(Ω) and

(
(∂iTK(vj)χ

s
j − ∂iTK(u)χs)χ{|TK(un)−TK(vj)|≤ν}

)
n

converges to

(∂iTK(vj)χ
s
j − ∂iTK(u)χs)χ{|TK(u)−TK(vj)|≤ν} strongly in EMi(Ω), we obtain

I2 =

d∑
i=1

∫
{|TK(u)−TK(vj)|≤ν}

liK(∂iTK(vj)χ
s
j − ∂iTK(u)χs) dx+ ε(n).

By letting now j →∞, and using Lemma 2.3, we obtain that

I2 = ε(n, j). (4.43)

Estimates for I3–I5. Using similar arguments as above yields

I3 = ε(n, j), (4.44)

I4 = −
d∑
i=1

∫
Ω

ai(x, TK(u), ∂iTK(u)χs)∂iTK(u)χs + ε(n, j, s), (4.45)

I5 =

d∑
i=1

∫
Ω

ai(x, TK(u), ∂iTK(u)χs)∂iTK(u)χs + ε(n, j, s). (4.46)



20 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

Gathering estimates (4.41)-(4.46), we obtain

In,r ≤ (Cν + ε(n, j, s))θ + C3 meas{|TK(un)− TK(vj)| > ν}1−θ. (4.47)

By taking, first the lim sup with respect to n→∞, then j →∞, s→∞ and ν → 0
yields

lim sup
n→∞

d∑
i=1

∫
Ωr

[
(ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(u)))

× (∂iTK(un)− ∂iTK(u))
]θ

dx = 0.

Consequently, there exist a subsequence, still denoted in the same way, and a
negligible subset Z ⊂ Ω such that for all x ∈ Ω \ Z one has

d∑
i=1

(
ai(x, TK(un), ∂iTK(un))− ai(x, TK(un), ∂iTK(u))

)
× (∂iTK(un)− ∂iTK(u))→ 0.

(4.48)

Let x ∈ Ω \ Z be fixed. Then, according to assumption (A5), the sequence(
∂iTK(un)(x)

)
n
⊂ R is bounded for i = 1, . . . , d. By extracting a convergent

subsequence, still denoted in the same way, for some ζ = (ζ1, . . . , ζd) ∈ Rd, we have
that ∂iTK(un)(x)→ ζi, for i = 1, . . . , d. Passing to the limit in (4.48) yields

d∑
i=1

(
ai(x, TK(u)(x), ζi)− ai(x, TK(u)(x), ∂iTK(u)(x))

)
·
(
ζi − ∂iTK(u)(x)

)
= 0,

which, according to (3.2), it is only possible if ζi = ∂iTK(u)(x) for i = 1, . . . , d.
Therefore, for any K > 0, we have deduced that, up to a subsequence, ∇TK(un)→
∇TK(u) a.e. in Ω. Since K is arbitrary, we finally obtained the desired result. This
completes the proof of Proposition 4.4. �

Remark 4.5. A straightforward consequence of Proposition 4.4 is that, owing to
(4.26), δi = ai(x, u, ∂iu) that is,

ai(x, un, ∂iun)→ ai(x, u, ∂iu) weakly-∗ in LM̄i
(Ω). (4.49)

Step 4. L1−convergence of (ϕn). In this step, we will show that ϕn → ϕ strongly
in L1(Ω) up to a subsequence. The strongly convergence of (ϕn) in L1(Ω) is based in
the next results which generalize that of González Montesinos and Ortegón Gallego
[13, Lemma 4] which, in its turn, generalize the original results due to Xu in [17].

Lemma 4.6 ([12]). Let Mi be an N -function for each i = 1, . . . , d which admits

the representation Mi(t) =
∫ |t|

0
mi(s) ds and such that s ≤ mi(s) for all s > 0 and

all i = 1, . . . , d. Let M0 be the N -function defined in (2.6). Let (un) be a bounded
sequence in W 1LM(Ω) such that un → u strongly in EM0

(Ω). Then there exists a
subsequence (un(j)) ⊂ (un) such that, for every ε > 0, there exists a constant value

M = M(ε) and a function ψ ∈W 1,1(Ω) satisfying the following properties:

0 ≤ ψ ≤ 1, (4.50)

‖ψ − 1‖L1(Ω) + ‖∇ψ‖L1(Ω) ≤ ε, (4.51)

|u|, |un(j)| ≤M on {ψ > 0} for all j ≥ 1. (4.52)



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 21

Lemma 4.7 ([13]). Let (un, ϕn) be a weak solution to the system (4.12)-(4.14),
u ∈ EM0(Ω) and ϕ ∈ L∞(Ω) the limit functions appearing, respectively, in (4.17)
and (4.22). Then, for any function S ∈ C 1

0 (R), there exists a subsequence, still
denoted in the same way, such that

S(un)ϕn ⇀ S(u)ϕ weakly in H1(Ω). (4.53)

Moreover, if 0 ≤ S ≤ 1, then there exists a constant C > 0, independent of S, such
that

lim sup
n→∞

∫
Ω

ρn(un)|∇[S(un)ϕn − S(u)ϕ]|2 ≤ C‖S′‖∞(1 + ‖S′‖∞). (4.54)

Lemma 4.8. There exists a subsequence (ϕn(j)) ⊂ (ϕn) such that

lim
j→∞

∫
Ω

|ϕn(j) − ϕ| = 0. (4.55)

The proof of this result is a straightforward adaptation to that of [12, Lemma 5.7].

Step 5. Passing to the limit. According to (4.17), (4.19), (4.23), and (4.25), it is
straightforward that condition (A9) of Definition 1 is fulfilled. The convergences
in Proposition 4.4 and Lemma 4.8 lead us to (A10), and to obtain (A11), using
Proposition 4.4 and Lemma 4.8 again with (4.53), it is sufficient to let j approach
infinity in the expression

S(un(j))ρn(j)(un(j))∇ϕn(j) = ρn(j)(un(j))[∇(S(un(j))ϕn(j))− ϕn(j)∇S(un(j))].

This completes the proof of Theorem 4.1.

Remark 4.9. Condition (3.4) is not necessary when the constant k appearing
in (3.1) is less than 1.

5. Numerical approximation

In practical situations, it is very interesting to know the behavior of the solution
(u, ϕ) of (1.1) for different choices of the functions ai(x, s, ζ), i = 1, . . . , d, and ρ(s),
not only from the quantitative standpoint but from the qualitative one as well. This
information can be then used in order to design a thermistor useful for a particular
task, for instance, for not letting pass the electric current (by self-destruction) in
the event of an unexpected voltage increase.

The numerical resolution of a problem like (1.1) faces certain challenges. One
reason is that this system is strongly coupled through nonlinear terms. Also, we do
not know whether or not the solution given by Theorem 4.1 is unique.

In this section, we describe a numerical algorithm based on the least squares
method to obtain a numerical approximation of a solution to a system like (1.1)



22 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

with certain mixed Dirichlet-natural boundary conditions in d = 1, 2 or 3, namely,

−
d∑
i=1

∂i
(
ai(x, u, ∂iu)

)
= ρ(u)|∇ϕ|2 in Ω,

div(ρ(u)∇ϕ) = 0 in Ω,

u = u0 on ΓD,

d∑
i=1

ai(x, u, ∂iu)ni = 0 on ΓN ,

ϕ = ϕ0 on ΓD,
∂ϕ

∂n
= 0 on ΓN ,

(5.1)

where Ω ⊂ Rd is a Lipschitz bounded domain with boundary ∂Ω = ΓD ∪ ΓN ,
ΓD ∩ ΓN = ∅, ΓN being an open set with respect to the induced topology of ∂Ω,∫

Γ0
> 0, u0 ≥ 0 is a constant value, ϕ0 ∈ H1(Ω) ∩ L∞(Ω), n = n(x) is the unit,

outer normal vector at x ∈ ΓN (a.e. in ΓN ), n = (n1, . . . , nd), and the functions
a1, . . . , ad and ρ satisfy the assumptions (A1)–(A7) and (A9)–(A11). A similar
analysis like the one developed in the previous sections yields to the existence of a
capacity solution to the system (5.1).

5.1. Least squares and conjugate gradient method. In [14] the authors im-
plemented a fixed-point iterative method to obtain the numerical approximation of
a system like (5.1) for d = 2 in which ΓN = ∅ and the elliptic operator of the first
equation is the anisotropic ~p-Laplacian, that is a1 = a1(∂1u) and a2 = a2(∂2u) are
given by

a1(ζ) = |ζ|p1−2ζ, a2(ζ) = |ζ|p2−2ζ, for all ζ ∈ R, (5.2)

for some p1, p2 ∈ R with p1, p2 ≥ 2. The numerical simulations described in [14]
based in the referred iterative method have shown good convergence properties for
values of p1 and p2 in the interval (2, 5). However, when we try to apply this
same technique when a1 or a2 have not a polynomial growth, the algorithm does
not converge. Instead, we have developed a technique based on the least squares
method [4, 5]. The application of this technique needs more regularity to both
a = (a1, . . . , ad) and ρ, namely,

(A12) for a.e. x ∈ Ω, all ζ ∈ R and all i = 1, . . . , d, the function s ∈ R→ ai(x, s, ζ)
is of class C 1.

(A13) for a.e. x ∈ Ω, all s ∈ R and all i = 1, . . . , d the function ζ ∈ R→ ai(x, s, ζ)
is of class C 1.

(A14) ρ ∈ C 1(R).

We apply the finite element method to approximate the functions in the spaces
involved in the solution of (5.1). We first consider polygonal/polyhedral approxima-

tions of Ω, ΓD and ΓN (all of them still denoted in the same way). Let Th = {Tj}NT
j=1

be a triangulation of the domain Ω, NT = card Th. We consider the discrete space

Vh = {vh ∈ C (Ω̄) : vh|Tj
∈P`(Tj), for j = 1, . . . , NT }, (5.3)

where ` ≥ 1 is an integer and P`(Tj) is the space of polynomial functions in Tj of
degree `. Notice that Vh ⊂ W 1,∞(Ω) ⊂ W 1LM(Ω). We consider the projection of
ϕ0 ∈ H1(Ω) onto Vh and we still denote this projection as ϕ0. We also consider the
subspace V0h ⊂ Vh defined as V0h = {vh ∈ Vh : vh|ΓD

= 0}.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 23

Now we introduce the functional J : u0 + V0h → R as

J(uh) =
1

2

∫
Ω

|∇ξh|2, (5.4)

where ξh = ξh(uh) is defined in the following way. First, we compute ϕh ∈ Vh as
the unique solution to the problem

Find ϕh ∈ Vh such that ϕ = ϕ0 on ΓD and∫
Ω

ρ(uh)∇ϕh∇φh = 0, for all φh ∈ Vh such that φh|ΓD
= 0.

(5.5)

Then, ξh ∈ Vh is the unique solution to the problem

Find ξh ∈ Vh such that ξh = 0 on ΓD and∫
Ω

∇ξh∇vh =

d∑
i=1

∫
Ω

ai(x, uh, ∂iuh)∂ivh −
∫

Ω

ρ(uh)|∇ϕh|2vh,

for all vh ∈ Vh such that vh|ΓD
= 0.

(5.6)

Clearly, J(uh) ≥ 0 for all uh ∈ u0 + V0h and J(uh) = 0 if and only if (uh, ϕh) is
the solution to the discrete problem

uh ∈ Vh, uh = u0 on ΓD, ϕh ∈ Vh, ϕ = ϕ0 on ΓD and

d∑
i=i

∫
Ω

ai(x, uh, ∂iuh)∂ivh =

∫
Ω

ρ(uh)|∇ϕh|2vh,

for all vh ∈ Vh such that vh|ΓD
= 0.∫

Ω

ρ(uh)∇ϕh∇φh = 0, for all φh ∈ Vh such that φh|ΓD
= 0.

(5.7)

Problem (5.7) is a nonlinear discrete version of the variational formulation of
problem (5.1). Thus, if (uh, ϕh) is a solution to (5.7) then uh is a global minimizer
of the functional J on u0 +V0h. To obtain an approximation of the solution to (5.7),
we generate a minimizing sequence (umh )m≥1 ⊂ u0 + V0h so that the non-negative
sequence (J(umh ))m≥1 ⊂ R is a decreasing sequence. We can do this by means of
the Polack-Ribière version of the conjugate gradient method [4, 5]. Denoting by
J ′(uh) ∈ V ′0h the derivative of J at uh ∈ u0 + V0h, V ′0h being the dual space of V0h,
this algorithm consists of the following four steps:

Step 1. Initialization. Given u0
h ∈ u0 + V0h, compute g0

h ∈ V0h as the solution to
the variational problem

g0
h ∈ Vh, g0

h = 0 on ΓD and∫
Ω

∇g0
h∇vh = 〈J ′(u0

h), vh〉, for all vh ∈ Vh such that vh|ΓD
= 0.

(5.8)

and put z0
h = g0

h.
Then, for m ≥ 0, assuming that umh , gmh and zmh , are already known, compute

um+1
h , gm+1

h and zm+1
h by:

Step 2. Descent.

um+1
h = umh − λmzmh , (5.9)



24 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

λm > 0 being the value where the function λ ∈ [0,∞) → J(umh − λzmh ) attains its
minimum, that is

J(umh − λmzmh ) ≤ J(umh − λzmh ), for all λ ∈ [0,∞). (5.10)

Step 3. Construction of the new descent direction. Let gm+1
h ∈ V0h be the solution

to the variational problem

gm+1
h ∈ Vh, gm+1

h = 0 on ΓD and∫
Ω

∇gm+1
h ∇vh = 〈J ′(um+1

h ), vh〉, for all vh ∈ Vh such that vh|ΓD
= 0.

(5.11)

Then, define the number γm ∈ R and the new descent direction zm+1
h as follows

γm =

∫
Ω
∇gm+1

h ∇(gm+1
h − gmh )∫

Ω
|∇gmh |2

, (5.12)

zm+1
h = gm+1

h + γmz
m
h , (5.13)

Step 4. Stopping test. If a certain stopping test is not satisfied, then increase m
by one, go to Step 2 and proceed.

For the generation of the sequence (gmh ) we first need to compute the gradient
of J at umh . To do so, let uh ∈ u0 + V0h and wh ∈ V0h, then from (5.4) we obtain

〈J ′(uh), wh〉 =

∫
Ω

∇〈ξ′h, wh〉∇ξh, (5.14)

where ξ′h ∈ L (V0h) is the derivative of ξh with respect to uh (here, L (V0h) stands
for the space of linear and continuous mappings from V0h into itself). Owing
to (5.6), and taking into account that ai = ai(x, s, ζ) a.e. in Ω, for all s ∈ R
and ζ ∈ R, we deduce the problem for 〈ξ′h, wh〉 ∈ V0h

〈ξ′h, wh〉 ∈ Vh such that 〈ξ′h, wh〉 = 0 on ΓD and∫
Ω

∇〈ξ′h, wh〉∇vh

=

d∑
i=1

[ ∫
Ω

∂sai(x, uh, ∂iuh)wh∂ivh +

∫
Ω

∂ζai(x, uh, ∂iuh)∂iwh∂ivh

]
−
∫

Ω

ρ′(uh)wh|∇ϕh|2vh − 2

∫
Ω

ρ(uh)∇〈ϕ′h, wh〉∇ϕhvh,

for all vh ∈ Vh such that vh|ΓD
= 0.

(5.15)

In our computations, we neglected the last term in the right side of (5.15). The
reason is twofold. First, this term has appeared to be of a smaller order in front
of the other terms in this equation and, secondly, the computation of 〈J ′(uh), wh〉
becomes much simpler.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 25

Thus, making vh = ξh in (5.15) yields

〈J ′(uh), wh〉

=

d∑
i=1

[ ∫
Ω

∂sai(x, uh, ∂iuh)wh∂iξh +

∫
Ω

∂ζai(x, uh, ∂iuh)∂iwh∂iξh

]
−
∫

Ω

ρ′(uh)wh|∇ϕh|2ξh

(5.16)

The expression (5.16) is then used to compute the sequence (gmh )m≥1 from the
solution of problems (5.8) and (5.11).

Remark 5.1. Notice that the description of the full algorithm described above only
involves functions lying in the discrete space Vh defined in (5.3). One may check
that each term in the above expressions makes sense from the assumptions (A1)–
(A11), (A12)–(A14), and from the fact that Vh ⊂ W 1,∞(Ω). If we try to develop
the same algorithm within the spaces W 1LM(Ω) and H1(Ω) then, in general, these
terms would not be well defined, and this algorithm would be meaningless.

Γ2
D

Γ1
D

ΓN

18 mm

20 mm

(a) Thermistor boundary and dimensions. (b) Thermistor triangulation.

Figure 1. Description of the thermistor boundary ∂Ω = Γ1
D ∪

Γ2
D ∪ ΓN and the mesh used in the finite element simulations.

5.2. Numerical simulations. The algorithm described in the previous section
has been implemented to compute the numerical approximation of a solution to
problem (5.1). To do so, we have used the Freefem++ software package ([9]). These
numerical simulations have been carried out in the 2D domain Ω of Figure 1(a) and
with the following data.

a1 = a1(ζ) = |ζ|p−2ζ, a2 = a2(ζ) = 2βζeβ|ζ|
2

,



26 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

for the particular choices p = 2.8, 3 and 3.2, and β = 10−7. The electric conduc-
tivity ρ is given by

ρ(s) = 10e−|s−30|/20, s ∈ R.
The domain Ω is a barril shape set in R2. In Figure 1(a) we may find the description
of the boundary and the actual dimensions of Ω. The Dirichlet boundary ΓD has
two connected components, namely, Γ1

D and Γ2
D, and the natural boundary ΓN is

the complement set of ΓD in ∂Ω.
As far as the boundary conditions are concerned, we have taken

u0 = 30, ϕ0 =

{
0 on Γ1

D,

V0 on Γ2
D,

where V0 is a constant voltage, namely, V0 = 10 V or V0 = 40 V.
Figure 1(b) shows the finite element triangulation of Ω from which the discrete

space Vh is defined as in (5.3). It consists of 6,382 triangles, and 3,257 vertices.

Remark 5.2. One may wonder why the mesh density in Figure 1(b) is higher
along the horizontal central line than elsewhere. The reason is that we first try the
algorithm described in the previous section for the case of the ~p-Laplacian operator
where the functions a1 and a2 are given by (5.2). We checked the convergence of
this algorithm for this case. However, the error function ξh presented a singularity
along this horizontal central line. This led us to consider a denser mesh along this
central line, then obtaining better results.

In the definition of the space Vh we have taken ` = 1. On the other hand, the
initial guess u0

h of the initialization stage (Step 1) have been u0
h = u0 in all cases.

0 0.2×10−4 0.4×10−4 0.6×10−4 0.8×10−4 1×10−4 1.2×10−4

20738

20748

20758

20768

20778

20788

λ

J
(u
m h

−
λ
z
m h
)

λ 7→ J(um
h − λz

m
h )

Figure 2. The function λ ∈ [0,+∞) 7→ J(umh − λzmh ) is plotted
around the optimal value λm for m = 40, p = 3.0 and V0 = 10.

Figure 2 shows the behavior of the function J along a descent direction at iter-
ation m = 40 for p = 3.0 and V0 = 10. This plot puts in evidence one of the issues
concerning the execution of this algorithm in this context: the range of values near
λm is very small versus the actual values of the error function J(umh − λzmh ). This
is not a particular case: it has happened in the six situations we have run, that is,
p ∈ {2.8, 3.0, 3.2} and V0 ∈ {10, 40}, and for each m ≥ 0.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 27

We have used two different stopping tests (see Step 4 above). The first one indi-
cates a maximum number of iterations, Nmax, so that if m > Nmax then the algo-
rithm stops and we keep uNmax

h and ϕNmax

h as an approximation of the corresponding
problem. The second stopping test compares the value of γm, given in (5.12), with
a very small value, say εmin > 0. In this case, if γm < εmin, the algorithm stops as
well, and we keep umh and ϕmh as an approximation of the corresponding problem.
This second test avoids underflow situations: whenever γm0

< εmin occurs, it would
yield um0+1

h = um0

h , and the sequence becomes stationary for m ≥ m0. We took
Nmax = 10, 000 and εmin = 2× 10−19.

0 227 769 2000 3029 3938 5000 6000 7000 8000 9000 10000

0.01

0.1

1

Iterations m of the conjugate gradient method

J
(u
m h
)/
J
(u

0 h
)

p = 2.8, V0 = 10 p = 3.0, V0 = 10 p = 3.2, V0 = 10

p = 2.8, V0 = 40 p = 3.0, V0 = 40 p = 3.2, V0 = 40

Figure 3. Descent of the normalized error functional
J(umh )/J(u0

h) along the iterations generated by the conju-
gate gradient method uh = umh .

Figure 3 plots the descent of J/J(u0
h) as a function of the iterations for the

six different cases we have considered. Only when p = 3.0 and p = 3.2, with
V0 = 10 in both cases, the maximum number of iterations Nmax was reached and
the convergence is very slow. In the other cases, though the convergence is much
faster, the execution was stopped at iterations m = 227 (p = 2.8, V0 = 40), m = 769
(p = 3.0, V0 = 40), m = 3, 029 (p = 2.8, V0 = 10) andm = 3, 938 (p = 3.2, V0 = 40),
respectively.

In Figure 4 six different iterations of the sequence (umh ) are shown for p = 3.0
and V0 = 10. Starting from u0

h = 30, these iterations seem to have an increasing
character. Notice the different scale in each plot.

Figures 5 and 6 show the numerical approximation (umh , ϕ
m
n ) obtained from the

conjugate gradient algorithm described in the previous section where m is the index
of the last computed iteration. Figure 5 plots the distribution of the temperature uh
in the six considered cases, whereas Figure 6 plots the corresponding potential ϕh.
Obviously, it is expected that the higher the voltage V0 the higher the maximum
temperature. For this reason, we think that the algorithm has underestimated
the computed temperature in the case p = 2.8, V0 = 40. We also remark that
for V0 = 10 the maximum temperature decreases with p, whereas for V0 = 40 is
just the contrary. This behavior may be interesting for the design of a particular



28 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

thermistor that should switch off the current passing through itself in the event of
an unexpected voltage increase. In this situation, the thermistor is used as a fuse
protecting certain circuit components which are much more expensive than a single
thermistor.

6. Conclusions

We have analyzed a nonlinear strongly coupled system of two partial differential
equations of elliptic type, the second equation not being uniformly elliptic. This
system is a generalization of the so-called thermistor problem in which the physical
unknowns are the temperature, u, and the electric potential, ϕ, in a semiconductor
device. The special anisotropic structure of the operator in the first equation leads
us to consider this analysis in the framework of the anisotropic Orlicz-Sobolev
spaces. On the other hand, since the second equation is not uniformly elliptic,
we have introduced the concept of capacity solution adapted to this situation, and
show an existence result of a capacity solution.

To obtain a numerical solution of this problem, we first consider a projection
of the original problem from a straightforward application of the finite element
method. This yields a discrete variational formulation in certain finite dimension
vector spaces.

We have described a least squares method for the numerical approximation of
this discrete variational formulation. The minimizing sequence is generated by
means of the Polack-Ribière version of the conjugate gradient method. We have
implemented this whole algorithm in a 2D domain by using the Freefem++ software
package and run some numerical simulations for different choices of the functions
a1(x, s, ζ) and a2(x, s, ζ), and for a given conductivity ρ(s). These numerical results
may provide the necessary information in order to design a specific thermistor in
an electric circuit.

Finally, the convergence of the minimizing sequence obtained from this algo-
rithm may be very slow. This could be improved by the introduction of some
preconditioning technique, which may be considered in future works.

Acknowledgements

This research was partially supported by Ministerio de Ciencia e Innovación
of the Spanish Government [grant PID2020-117201RB-C21] and by Consejeŕıa de
Transformación Económica, Industria, Conocimiento y Universidades of the Junta
de Andalućıa [grant UAL2020-FQM-B2046], in both cases with the participation
of the European Regional Development Fund (ERDF/FEDER).

References

[1] R. A. Adams; J. J. F. Fournier; Sobolev spaces, Elsevier, 2003.
[2] S. N. Antontsev, M. Chipot; The thermistor problem: Existence, smoothness uniqueness,

blowup, SIAM Journal on Mathematical Analysis, 25 (1994), 1128–1156.
[3] A. Elmahi, D. Meskine; Strongly nonlinear parabolic equations with natural growth terms in

Orlicz spaces, Nonlinear Analysis: Theory, Methods & Applications, 60 (2005), 1–35.

[4] R. Glowinski; Numerical Solution of Nonlinear Boundary Value Problems by Variational
Methods. Applications, North-Holland, Amsterdam, 1984, 1455–1508, Proceedings of the

International Congress of Mathematicians August 16-24, 1983, Warszawa.

[5] R. Glowinski; Numerical Methods for Nonlinear Variational Problems, Scientific Computa-
tion, Springer Berlin Heidelberg, 2008.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 29

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30

30.3

30.6

30.9

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30

32

34
35

37

(a) m = 10, ‖u‖∞ = 30.9095. (b) m = 500, ‖u‖∞ = 37.6472.

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

30
35

40
45

50

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
40

49
58

67

76

(c) m = 2, 000, ‖u‖∞ = 54.2487. (d) m = 5, 000, ‖u‖∞ = 75.4213.
0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
41

51

62

72

83

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

30
42

53
64

75

86

(e) m = 8, 000, ‖u‖∞ = 82.4945. (f) m = 10, 000, ‖u‖∞ = 85.3979.

Figure 4. Iterations m = 10, 500, 2,000, 8,000 and 10,000 of the
conjugate gradient method: temperature u for p = 3 and voltage
V0 = 10.



30 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
43

56
69

82

95

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
42

54
66

78

90

(a) p = 2.8, V0 = 10 V. (b) p = 2.8, V0 = 40 V.

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
42

53
64

75

86

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
44

58
72

86

100

(c) p = 3.0, V0 = 10 V. (d) p = 3.0, V0 = 40 V.

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

30
37

43
49

55

61

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

30
45

59

74

88

102

(e) p = 3.2, V0 = 10 V. (f) p = 3.2, V0 = 40 V.

Figure 5. Distribution of the temperature u for p = 2.8, 3.0 and
3.2 and voltages V0 = 10 (left) and V0 = 40 (right).

[6] M. T. González Montesinos, F. Ortegón Gallego; The evolution thermistor problem with
degenerate thermal conductivity, Communications on Pure & Applied Analysis, 1 (2002),

313–325.



EJDE-2022/84 EXISTENCE OF A SOLUTION AND ITS PPROXIMATION 31

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

0
2

4
6

8

10

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

0
8

16
24

32

40

(a) p = 2.8, V0 = 10 V. (b) p = 2.8, V0 = 40 V.

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

0
2

4
6

8

10

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

0
8

16
24

32

40

(c) p = 3.0, V0 = 10 V. (d) p = 3.0, V0 = 40 V.

0

0.0036

0.0072

0.0108

0.0144
0

0.004

0.008

0.012

0.016

0.02

0
2

4
6

8

10

0

0.0036

0.0072

0.0108

0.0144

0

0.004

0.008

0.012

0.016

0.02

0
8

16
24

32

40

(e) p = 3.2, V0 = 10 V. (f) p = 3.2, V0 = 40 V.

Figure 6. Potential ϕ for p = 2.8, 3.0 and 3.2 and voltages V0 =
10 (left) and V0 = 40 (right).

[7] J.-P. Gossez; Nonlinear elliptic boundary value problems for equations with rapidly (or slowly)

increasing coefficients, Trans. Amer. Math. Soc., 190 (1974), 163–205.



32 F. ORTEGÓN GALLEGO, H. OUYAHYA, M. RHOUDAF EJDE-2022/84

[8] J.-P. Gossez; Some approximation properties in Orlicz-Sobolev spaces, Studia Mathematica,

74 (1982), 17–24.

[9] F. Hecht; New development in Freefem++, Journal of Numerical Mathematics, 20 (2012),
251–266.

[10] J. R. S. Howison, M. Shillor; Stationary solutions to the thermistor problem, Journal of

Mathematical Analysis and Applications, 174 (1993), 573–588.
[11] M. A. Krasnosel’skĭı, Y. B. Rutickĭı; Convex functions and Orlicz spaces, Groningen, P.

Noordhoff, 1969.

[12] H. Moussa, F. Ortegón Gallego, M. Rhoudaf; Capacity solution to a nonlinear elliptic coupled
system in Orlicz-Sobolev spaces, Mediterranean Journal of Mathematics, 17 (2020), 1–28.

[13] F. Ortegón Gallego, M. T. González Montesinos; Existence of a capacity solution to a coupled

nonlinear parabolic–elliptic system, Communications on Pure & Applied Analysis, 6 (2007),
23–42.

[14] F. Ortegón Gallego, M. Rhoudaf, H. Talbi; Capacity solution and numerical approximation to
a nonlinear coupled elliptic system in anisotropic Sobolev spaces, Journal of Applied Analysis

& Computation, 12 (2022), 2184–2207.

[15] H. Ouyahya, B. Lahmi, R. Elarabi; Note on some nonlinear anisotropic elliptic equations,
Submitted.

[16] M. Tienari; A degree theory for a class of mappings of monotone type in Orlicz-Sobolev

spaces, Annales Academiae scientiarum Fennicae: Mathematica Fennica, (1994), 1–68.
[17] X. Xu; A strongly degenerate system involving an equation of parabolic type and an equation

of elliptic type, Communications in Partial Differential Equations, 18 (1993), 199–213.

[18] X. Xu; A degenerate Stefan-like problem with Joule’s heating, SIAM Journal on Mathemat-
ical Analysis, 23 (1992), 1417–1438.

[19] X. Xu, On the existence of bounded temperature in the thermistor problem with degeneracy,

Nonlinear Analysis: Theory, Methods & Applications, 42 (2000), 199–213.

Francisco Ortegón Gallego

Departamento de Matemáticas, Facultad de Ciencias del Mar y Ambientales, Universi-

dad de Cádiz, Campus del Ŕıo San Pedro, 11510 Puerto Real, Cádiz, Spain
Email address: francisco.ortegon@uca.es

Hakima Ouyahya, Équipe EDP et calcul scientifique, laboratoire de mathématiques

et leurs intéractions, Faculté des Sciences, Moulay Ismail University, Meknes, Mo-
rocco

Email address: hakima.ouyahya@edu.umi.ac.ma

Mohamed Rhoudaf

Équipe EDP et calcul scientifique, laboratoire de mathématiques et leurs intéractions,

Faculté des Sciences, Moulay Ismail University, Meknes, Morocco
Email address: m.rhoudaf@umi.ac.ma


	1. Introduction
	2. Preliminaries
	2.1. N-functions
	2.2. Anisotropic Orlicz-Sobolev spaces

	3. Essential assumptions and main result
	4. An existence result
	Proof of Theorem ??

	5. Numerical approximation
	5.1. Least squares and conjugate gradient method
	5.2. Numerical simulations

	6. Conclusions
	Acknowledgements
	References