International Journal of Analysis and Applications

Volume 17, Number 6 (2019), 1034-1051

URL: https://doi.org/10.28924/2291-8639

DOI: 10.28924/2291-8639-17-2019-1034

NUMERICAL QUENCHING FOR HEAT EQUATIONS WITH COUPLED NONLINEAR

BOUNDARY FLUX

KOUAMÉ BÉRANGER EDJA1,∗, KOFFI N’GUESSAN2, BROU JEAN-CLAUDE KOUA3

AND KIDJEGBO AUGUSTIN TOURÉ1

1Institut National Polytechnique Houphouët-Boigny Yamoussoukro, BP 2444, Côte d’Ivoire

2UFR SED, Université Alassane Ouattara , 01 BP V 18 Bouaké 01, Côte d’Ivoire

3UFR Mathématique et Informatique, Université Félix Houphouët Boigny, Côte d’Ivoire

∗Corresponding author: kouame.edja@inphb.ci

Abstract. In this paper, we study a numerical approximation of the following problem ut = uxx, vt = vxx,

0 < x < 1, 0 < t < T ; ux(0, t) = u
−m(0, t) + v−p(0, t), vx(0, t) = u

−q(0, t) + v−n(0, t) and ux(1, t) =

vx(1, t) = 0, 0 < t < T, where m, p, q and n are parameters. We prove that the solution of a semidiscrete

form of above problem quenches in a finite time only at first node of the mesh. We show that the time

derivative of the solution blows up at quenching node. Some conditions under which the non-simultaneous

or simultaneous quenching occurs for the solution of the semidiscrete problem are obtained. We establish

the convergence of the quenching time. Finally, some numerical results to illustrate our analysis are given.

1. Introduction

In this paper, we study the behavior of a semidiscrete approximation of the following heat equations

involving nonlinear boundary flux conditions :

ut(x,t) = uxx(x,t), vt(x,t) = vxx(x,t), (x,t) ∈ (0, 1) × (0,T), (1.1)

Received 2019-08-08; accepted 2019-09-23; published 2019-11-01.

2010 Mathematics Subject Classification. 65M06, 65M12, 35K05, 35K55.

Key words and phrases. Numerical quenching; non-simultaneous; heat equation; nonlinear boundary.

c©2019 Authors retain the copyrights

of their papers, and all open access articles are distributed under the terms of the Creative Commons Attribution License.

1034

https://doi.org/10.28924/2291-8639
https://doi.org/10.28924/2291-8639-17-2019-1034


Int. J. Anal. Appl. 17 (6) (2019) 1035

ux(0, t) = u
−m(0, t) + v−p(0, t), vx(0, t) = u

−q(0, t) + v−n(0, t), t ∈ (0,T), (1.2)

ux(1, t) = 0, vx(1, t) = 0, t ∈ (0,T), (1.3)

u(x, 0) = u0(x), v(x, 0) = v0(x), x ∈ [0, 1], (1.4)

where m,n ≥ 0, p,q > 0, u0 and v0 are positive smooth functions satisfying the compatibility conditions

u′0(0) = u
−m
0 (0) + v

−p
0 (0), v

′
0(0) = u

−q
0 (0) + v

−n
0 (0), u

′
0(1) = 0, v

′
0(1) = 0, and u

′
0,v
′
0 ≥ 0 and u′′0,v′′0 < 0 on

(0, 1].

Here [0,T) is the maximal time interval such that

∀t ∈ [0,T), inf min
0≤x≤1

{u(x,t),v(x,t)} > 0.

We have

lim
t→T−

inf min
0≤x≤1

{u(x,t),v(x,t)} = 0+.

The time T can be finite or infinite. If T is finite, then we say that the solution (u,v) quenches in a finite

time and T is called the quenching time of (u,v). If T is infinite, then we affirm that the solution (u,v)

quenches globally.

Nonlinear parabolic systems like (1.1)-(1.4) come from chemical reactions, heat transfer, etc, where u and

v represent the temperatures of two different materials during heat propagation. The quenching phenomenon

of parabolic problems has been the issue of intensive study (see for example [3, 4, 8–10] and the references

cited therein), particulary the study of heat equations system with nonlinear boundary conditions has been

the subject of investigation of several authors in recent years (see [6, 7, 14, 15, 17] and the references cited

therein). In [7] the authors study this problem, they prove that the solution (u,v) quenches in finite time T

and the quenching occurs only at the boundary x = 0 for 0 < u0,v0 ≤ 1. They show that

• if p < n + 1, there exist initial data such that the non-simultaneous quenching occurs ;

• if q ≤ n(m+1)
n+1

and p ≥ n + 1 (p ≤ m(n+1)
m+1

and q ≥ m + 1), the non-simultaneous quenching occurs

for any positive initial data ;

• if q ≥ m+1, p ≥ n+1, any quenching must be simultaneous and obtain of results on non-simultaneous

quenching rate.

Moreover, if quenching is simultaneous they found the quenching rate, which depends on the parameter in

the flux associated to the other component of the initial data.

To the best of our knowledge, no studies have been performed on the numerical approximation of equations

(1.1)-(1.4). In this paper, we investigate in the numerical study using a semidiscrete form of (1.1)-(1.4),



Int. J. Anal. Appl. 17 (6) (2019) 1036

especially in study of simultaneous and non-simultaneous quenching. For that, we consider a uniform mesh

on the interval [0, 1]

xi = (i− 1)h, i = 1, . . . ,I, h = 1/(I − 1),

Uh(t) = (U1(t), . . . ,UI(t))
T , Vh(t) = (V1(t), . . . ,VI(t))

T , where Ui(t) and Vi(t) are the values of the numerical

approximation of u and v at the nodes xi at time t. We also denote ϕ1,i and ϕ2,i, respectively, the values

of the numerical approximation of u0 and v0 at the nodes xi. By the finite difference method we obtain the

following system of ODEs whose the solution is (Uh,Vh) :

U′i(t) = δ
2Ui(t) − bi

(
U−mi (t) + V

−p
i (t)

)
, i = 1, . . . ,I, t ∈ (0,Th), (1.5)

V ′i (t) = δ
2Vi(t) − bi

(
U
−q
i (t) + V

−n
i (t)

)
, i = 1, . . . ,I, t ∈ (0,Th), (1.6)

Ui(0) = ϕ1,i Vi(0) = ϕ2,i, i = 1, . . . ,I, (1.7)

where

0 < ϕ1,i ≤ M, 0 < ϕ2,i ≤ N, i = 1, . . . ,I,

δ2Ui(t) =
Ui−1(t) − 2Ui(t) + Ui+1(t)

h2
, 2 ≤ i ≤ I − 1, t ∈ (0,Th),

δ2U1(t) =
2U2(t) − 2U1(t)

h2
, δ2UI(t) =

2UI−1(t) − 2UI(t)
h2

, t ∈ (0,Th),

b1 =
2

h
, and bi = 0, i = 2, . . . ,I.

Here [0,Th) is the maximal time interval such that

∀t ∈ [0,Th), inf min
1≤i≤I

{Ui(t),Vi(t)} > 0.

We have

lim
t→T−

h

inf min
1≤i≤I

{Ui(t),Vi(t)} = 0+.

The time Th can be finite or infinite. If Th is finite, then we say that the solution (Uh,Vh) quenches in a

finite time and Th is called the semidiscrete quenching time of (Uh,Vh). If Th is infinite, then we affirm that

the solution (Uh,Vh) quenches globally.

We show that our semidiscrete scheme reproduces well the conditions for the quenching, quenching set

or simultaneous and non-simultaneous quenching of system (1.1)-(1.4). By following, it is also proved that

when quenching occurs, the semidiscrete quenching time converges to the theoretical one when the mesh

size goes to zero and we give a result on numerical non-simultaneous quenching rate. For previous work on

numerical approximations of heat equations with non-linear boundary conditions we refer to [1,2,5,11–13,16]

and the references cited therein. The rest of the paper is organized as follows : in the next section, we give



Int. J. Anal. Appl. 17 (6) (2019) 1037

some properties concerning our semidiscrete scheme. In Section 3, under some conditions, we prove that

the solution of the semidiscrete scheme (1.5)-(1.7) quenches in a finite time, we give a result on numerical

quenching set. We also show that the time derivative of the solution blows up at quenching node. In Section

4 a criterion to identify simultaneous and non-simultaneous quenching is proposed. In Section 5, we show

the convergence of the semidiscrete scheme and the convergence of the quenching times to the theoretical

one when the mesh size goes to zero. Finally, in the last section, we give some numerical results to illustrate

our analysis.

2. Properties of the semidiscrete scheme

In this section, we give some auxiliary results for the problem (1.5)-(1.7).

Definition 2.1. We say that (Uh,Vh) ∈
(
C1([0,Th), R

I)
)2

is a lower solution of (1.5)-(1.7) if

Ui
′(t) ≤ δ2Ui(t) − bi(Ui−m(t) + Vi−p(t)), i = 1, . . . ,I, t ∈ (0,Th),

V ′i (t) ≤ δ
2Vi(t) − bi(Ui−q(t) + Vi−n(t)), i = 1, . . . ,I, t ∈ (0,Th),

0 < Ui(0) ≤ ϕ1,i, 0 < Vi(0) ≤ ϕ2,i, i = 1, . . . ,I,

where (Uh,Vh) is the solution of (1.5)-(1.7). On the other hand, we say that

(Uh,Vh) ∈
(
C1([0,Th), R

I)
)2

is an upper solution of (1.5)-(1.7) if these inequalities are reversed.

The following lemma is a discrete form of the maximum principle.

Lemma 2.1. Let eh, ch, αh, βh ∈ (C0([0,Th), RI) and Uh, Vh ∈ C1([0,Th), RI) such that

U′i(t) − δ
2Ui(t) + ei(t)Ui(t) + ci(t)Vi(t) ≥ 0, i = 1 . . . ,I, t ∈ (0,Th),

V ′i (t) −δ
2Vi(t) + αi(t)Ui(t) + βi(t)Vi(t) ≥ 0, i = 1 . . . ,I, t ∈ (0,Th),

Ui(0) ≥ 0, Vi(0) ≥ 0, i = 1 . . . ,I.

Then we have

Ui(t) ≥ 0, Vi(t) ≥ 0, i = 1 . . . ,I, t ∈ (0,Th).

Proof. Let T0 < Th and let (Zh(t),Wh(t)) = (e
λtUh(t),e

λtVh(t)) where λ is a real. We find that

(Zh(t),Wh(t)) satisfies the following inequalities :

Z′i(t) − δ
2Zi(t) + (ei(t) −λ)Zi(t) + ci(t)Wi(t) ≥ 0, i = 1 . . . ,I, t ∈ (0,Th), (2.1)

W ′i (t) − δ
2Wi(t) + αi(t)Zi(t) + (βi(t) −λ)Wi(t) ≥ 0, i = 1 . . . ,I, t ∈ (0,Th), (2.2)



Int. J. Anal. Appl. 17 (6) (2019) 1038

Zi(0) ≥ 0, Wi(0) ≥ 0, i = 1 . . . ,I. (2.3)

Set m = min

{
min

1≤i≤I,t∈[0,T0]
Zi(t), min

1≤i≤I,t∈[0,T0]
Wi(t)

}
. Since for i ∈ {0, . . . ,I}, Zi(t) and Wi(t) are

continuous functions on a compact, we can assume that m = Zi0 (ti0 ) for a certain i0 ∈{0, . . . ,I}.

Assume m < 0.

Taking λ negative such that

ei0 (ti0 ) −λ > 0 and βi0 (ti0 ) −λ > 0.

If ti0 = 0, then Zi0 (0) < 0, which contradicts (2.3), hence ti0 6= 0 ;

if 1 ≤ i0 ≤ I, we have

Z′i0 (ti0 ) = limk→0

Zi0 (ti0 ) −Zi0 (ti0 −k)
k

≤ 0.

Moreover by a straightforward computation we get

Z′i0 (ti0 ) − δ
2Zi0 (ti0 ) + (ei0 (ti0 ) −λ) Zi0 (ti0 ) + ci0 (ti0 )Wi0 (ti0 ) < 0,

but these inequalities contradict (2.1) and the proof is completed. �

Lemma 2.2. Let (Uh,Vh) and (Uh,Vh) be lower and upper solutions of (1.5)-(1.7) respectively such that,

(Uh(0),Vh(0)) ≤ (Uh(0),Vh(0)) then

(Uh(t),Vh(t)) ≤ (Uh(t),Vh(t)).

Proof. Let us define (Zh(t),Wh(t)) = (Uh(t),Vh(t)) − (Uh(t),Vh(t)). We obtain

Z′i(t) − δ
2Zi(t) −mbi(µi(t))−m−1Zi(t) −pbi(νi(t))−p−1Wi(t) ≥ 0, i = 1, . . . ,I (2.4)

W ′i (t) − δ
2Wi(t) −qbi(µi(t))−q−1Zi(t) −nbi(νi(t))−n−1Wi(t) ≥ 0, i = 1, . . . ,I (2.5)

Zi(0) ≥ 0, Wi(0) ≥ 0, i = 1, . . . ,I (2.6)

where µi(t), νi(t) lie, respectively, between Ui(t) and Ui(t), and between Vi(t) and Vi(t), for i ∈{1, . . . ,I}.

We can rewrite (2.4)-(2.5) as

Z′i(t) −δ
2Zi(t) + ei(t)Zi(t) + ci(t)Wi(t) ≥ 0, i = 1, . . . ,I, t ∈ (0,Th),

W ′i (t) −δ
2Wi(t) + αi(t)Zi(t) + βi(t)Wi(t) ≥ 0, i = 1, . . . ,I, t ∈ (0,Th),

where ei(t) = −mbi(µi(t))−m−1, ci(t) = −pbi(νi(t))−p−1, αi(t) = −qbi(µi(t))−q−1, and βi(t) =

−nbi(νi(t))−n−1 i = 1, . . . ,I, ∀t ∈ (0,Th). According to Lemma 2.1, Zi(t) ≥ 0, Wi(t) ≥ 0, for i = 1, . . . ,I,

∀t ∈ (0,Th) and the proof is completed. �



Int. J. Anal. Appl. 17 (6) (2019) 1039

The next lemma gives the properties of the semidiscrete solution.

Lemma 2.3. Let (Uh,Vh) ∈
(
C1([0,Th), R

I)
)2

be the solution of (1.5)-(1.7) with an initial data (ϕ1,h,ϕ2,h)

upper solution such that 0 < ϕ1,i < ϕ1,i+1 ≤ M and 0 < ϕ2,i < ϕ2,i+1 ≤ N for i = 1, . . . ,I − 1. Then we

have

(i) 0 < Ui(t) ≤ ϕ1,i ≤ M and 0 < Vi(t) ≤ ϕ2,i ≤ N, for i = 1, . . . ,I, t ∈ [0,Th);

(ii) (Ui+1(t),Vi+1(t)) > (Ui(t),Vi(t)), i = 1, . . . ,I − 1, t ∈ (0,Th) ;

(iii) (U′i(t),V
′
i (t)) ≤ 0, i = 1, . . . ,I, t ∈ (0,Th).

Proof. (i) Since (ϕ1,h,ϕ2,h) is an upper solution of (1.5)-(1.7), by the Lemma 2.1 and 2.2 we have

0 < Ui(t) ≤ M and 0 < Vi(t) ≤ N, for i = 1, . . . ,I, t ∈ [0,Th).

(ii) We argue by contradiction. Assume, that t0 the first t > 0, such that (Ki,Li)(t) = (Ui+1−Ui,Vi+1−

Vi)(t) > 0, for 1 ≤ i ≤ I − 1, but min{Ki0 (t0),Li0 (t0)} = 0 for a certain i0 ∈{1, ...,I − 1}. Assume

that Ki0 (t0) = Ui0+1(t0) − Ui0 (t0) = 0. Without lost of generality, we can suppose that i0 is the

smallest integer which satisfies the above equality. Therefore, by simple computation, (Kh,Lh)

verifies

K′h(t) = −A
′Kh(t) + B

′U−mh (t) + B
′V
−p
h (t),

L′h(t) = −A
′Lh(t) + B

′U
−q
h (t) + B

′V −nh (t),

where

A′ =
1

h2




3 −1 0 . . . 0

−1 2 −1
. . .

...

0
. . .

. . .
. . . 0

...
. . . −1 2 −1

0 . . . 0 −1 3



, B′ =




2
h

0 0 . . . 0

0 0
. . .

...

...
. . .

. . .
. . .

...

...
. . . 0 0

0 . . . . . . 0 0



.

On the one hand

K′i0 (t0) = lim�→0

Ki0 (t0) −Ki0 (t0 − �)
�

≤ 0,

and, on the other hand

−
I∑
j=1

a′i0,jKj(t) + b
′
i0
Umi0 (t0) + b

′
i0
V
p
i0

(t0) > 0.

Thus we have a contradiction, hence we obtain the desired result.



Int. J. Anal. Appl. 17 (6) (2019) 1040

(iii) Denote Fi(t) = Ui(t) −Ui(t + ε) and Gi(t) = Vi(t) −Vi(t + ε), for i = 1, . . . ,I, using (i) we obtain

Fi(0) ≥ 0, Gi(0) ≥ 0 for i = 1, . . . ,I. It is not hard to see that

F ′i (t) = δ
2Fi(t) + mbi(ξi(t))

−m−1Fi(t) + pbi(ηi(t))
−p−1Gi(t) ≥ 0,

G′i(t) = δ
2Gi(t) + qbi(ξi(t))

−q−1Fi(t) + nbi(ηi(t))
−n−1Gi(t) ≥ 0,

where ξi(t), ηi(t) lie, respectively, between Ui(t + ε) and Ui(t) and between Vi(t + ε) and Vi(t). From

Lemma 2.1 we get

Fi(t) ≥ 0 and Gi(t) ≥ 0 for i = 1, . . . ,I, t ∈ (0,Th).

This fact implies the desired result.

�

3. Quenching and blow-up

Let (Uh,Vh) be the solution of (1.5)-(1.7) with 0 < ϕ1,i ≤ M, 0 < ϕ1,i ≤ N for i = 1, . . . ,I. Inspired

by [4, 7] we prove that (Uh,Vh) quenches in a finite time and (U
′
h,V

′
h) blows up at quenching node.

Theorem 3.1. The solution (Uh,Vh) of (1.5)-(1.7) quenches in a finite time with the only quenching node

i =1.

Proof. Integrating (1.5) in time we find

Ui(t) −Ui(0) =
∫ t

0

δ2Ui(τ) + bi(U
−m
i (τ) + V

−p
i (τ))dτ

summing up the above inequality we get

I∑
i=1

hUi(t) =

I∑
i=1

hUi(0) +

∫ t
0

UI−1(τ) −UI(τ)
h

+
U2(τ) −U1(τ)

h
− 2(U−m1 (τ) + V

−p
1 (τ))dτ.

From (1.5) we have

h

2
UI(t) −

h

2
UI(0) =

∫ t
0

UI−1(τ) −UI(τ)
h

dτ, and

h

2
U1(t) −

h

2
U1(0) =

∫ t
0

U2(τ) −U1(τ)
h

− (U−m1 (τ) + V
−p
1 (τ))dτ.

Thus

h

2
UI(t) +

I−1∑
i=2

hUi(t) +
h

2
U1(t) =

h

2
UI(0) +

I−1∑
i=2

hUi(0) +
h

2
U1(0) −

∫ t
0

U−m1 (τ) + V
−p
1 (τ)dτ,

therefore

h

2
UI(t) +

I−1∑
i=2

hUi(t) +
h

2
U1(t) ≤ M − (M−m + N−p)t.



Int. J. Anal. Appl. 17 (6) (2019) 1041

Proceeding as before, we find that

h

2
VI(t) +

I−1∑
i=2

hVi(t) +
h

2
V1(t) ≤ N − (M−q + N−n)t,

which yield a contradiction because Uh and Vh are positive for all times. Then there exists 0 < Th < ∞ such

that

lim
t→T−

min{U1(t),V1(t)} = 0+.

To show i = 1 is the unique quenching node. In everything that follows i ∈ {1, . . . ,I − 1} and t ∈ (0,Th).

Set g(Ui(t)) = U
−m
i (t), f(Vi(t)) = V

−p
i (t), d(Ui(t)) = U

−q
i (t), j(Vi(t)) = V

−n
i (t), and

Zi(t) =
Ui+1(t) −Ui(t)

h
−φi(g(Ui(t)) + f(Vi(t))) (3.1)

Wi(t) =
Vi+1(t) −Vi(t)

h
−φi(d(Ui(t)) + j(Vi(t))) (3.2)

where φi, δ
2φi ≥ 0, δ+φi ≤ 0, φI = 0, φ1 = 1, φi(g(Ui(0))+f(Vi(0))) ≤ δ+Ui(0) and φi(d(Ui(0))+j(Vi(0))) ≤

δ+Vi(0).

By means of Taylor expansions we have

δ2(φik(Ji(t))) = φik
′(Ji(t))δ

2Ji(t) + k(Ji(t))δ
2φi + k

′(Ji(t))δ
+φiδ

+Ji(t) + k
′(Ji(t))δ

−φiδ
−Ji(t)

+φi
(δ+Ji(t))

2

2
k′′(ρi(t)) + φi

(δ−Ji(t))
2

2
(k′′(λi(t)), i = 2, . . . ,I − 1,

δ2(φ1k(J1(t))) = φ1k
′(J1(t))δ

2J1(t) + k(J1(t))δ
2φ1 + 2k

′(J1(t))δ
+φ1δ

+J1(t)

+φ1(δ
+J1(t))

2k′′(ρ1(t)).

If we use the fact that Ji, δ
+Ji(t) and δ

2Ji(t) are nonnegative and the hypothesis on φh, we arrive at

δ2(φik(Ji(t)) ≥ φik′(Ji(t))δ2Ji(t), i = 1, . . . ,I − 1. (3.3)

By using (3.3) we can get

Z′i(t) − δ
2Zi(t) ≥

bi
h

(g(Ui) + f(Vi)) + biφig
′(Ui) (g(Ui) + f(Vi)) + biφif

′(Ui) (d(Ui) + j(Vi)) .

The above inequalities implies that

Z′i(t) − δ
2Zi(t) + big

′(Ui(t))Zi(t) + bif
′(Vi(t))Wi(t) ≥ bi[

1

h
(g(Ui(t)) + f(Vi(t)))

+f′(Ui(t))(d(Ui(t)) + j(Vi(t))) + g
′(Ui(t))(g(Ui(t)) + f(Vi(t))].

We obtain

Z′i(t) − δ
2Zi(t) + big

′(Ui(t))Zi(t) + bif
′(Vi(t))Wi(t) ≥ 0,



Int. J. Anal. Appl. 17 (6) (2019) 1042

for the parameter h small enough. Thus we have

Z′i(t) −δ
2Zi(t) + big

′(Ui(t))Zi(t) + bif
′(Vi(t))Wi(t) ≥ 0,

W ′i (t) − δ
2Wi(t) + bid

′(Ui(t))Zi(t) + bij
′(Vi(t))Wi(t) ≥ 0,

Zi(0) ≥ 0,Wi(0) ≥ 0.

Using the Lemma 2.1 we have Zi(t) ≥ 0 and Wi(t) ≥ 0, for i = 1, . . . ,I − 1 and t ∈ (0,Th). This implies

that
Ui+1(t) −Ui(t)

h
≥ φi(g(Ui(t)) + f(Vi(t))) ≥

1

2

(
1

Mm
+

1

Np

)
for i = 1, . . . ,J, with φJ =

1

2
, where

J ∈{2, . . . ,I − 1}. Thus by summing we obtain

Ui(t) ≥ U1 +
(i− 1)h

2

(
1

Mm
+

1

Np

)
≥

(i− 1)h
2

(
1

Mm
+

1

Np

)
whenever i > 1.

The same happens for Vh. �

Theorem 3.2. If lim
t→T−

h

U1(t) = 0
+

(
lim
t→T−

h

V1(t) = 0
+

)
, then U′h(t) blows up (V

′
h(t) blows up).

Proof. Suppose U′h(t) is bounded. Then, there exists a negative constant M such that U
′
h(t) > M. We have

I−1∑
i=1

i∑
j=1

h2U′j(t) >

I−1∑
i=1

i∑
j=1

h2M.

I−1∑
i=1

i∑
j=1

h2M =

I−1∑
i=1

ih2M

=
M

2
+
hM

2
.

I−1∑
i=1

i∑
j=1

h2U′j(t) =

I−1∑
i=2


 i∑
j=2

h2U′j(t) + h
2U′1(t)


 + h2U′1(t)

From (1.5) we arrive at

I−1∑
i=1

i∑
j=1

h2U′j(t) = UI(t) −U1(t) −
(
V
−p
1 (t) + U

−m
1 (t)

)
+
h

2
U′1(t)

and Lemma 2.3 we arrive at

UI(t) −U1(t) −
(
V1(t)

−p + U1(t)
−m) > M

2
+ hM.

As t → T−h , the left-hand side tends to infinity while the right-side is finite. This contradiction shows that

U′h blows up. �



Int. J. Anal. Appl. 17 (6) (2019) 1043

4. Simultaneous vs. non-simultaneous quenching

In this Section we consider (Uh,Vh) the solution of (1.5)-(1.7) with h fixed, and we give some sufficient

conditions for the existence of simultaneous and non-simultaneous quenching.

Theorem 4.1. If Uh quenches and Vh does not quench in (1.5)-(1.7) then q < m + 1.

Proof. As Vh does not quench, by (1.5) there exists c > 0 such that

U′1(t) ≥−cU
−m
1 (t),

integrating this inequality from t to Th, we get

U1(t) ≤ C(Th − t)
1

m+1 , where C = ((m + 1)c)
1/(m+1)

. (4.1)

By using (4.1) and (1.6) , we obtain

V ′1 (t) ≤ δ
2V1(t) − b1

(
V −n1 (t) + C(Th − t)

− q
m+1

)
.

Thus V1(Th) ≤ C1 −C
∫Th

0
(Th − t)−

q
m+1 dt. We can see that this integral diverges if q ≥ m + 1, which is a

contradiction and the proof is completed. �

Corollary 4.1. Simultaneous quenching happens if q ≥ m + 1, p ≥ n + 1.

Lemma 4.1. Let (Uh,Vh) be the solution of (1.5)-(1.7). Assume that Uh quenches at time Th (Vh quenches

at time Th) and

δ2ϕ1,i − bi
(
ϕ−m1,i + ϕ

−p
2,i

)
+ c

(
ϕ−m1,i + ϕ

−p
2,i

)
≤ 0, (4.2)

δ2ϕ2,i − bi
(
ϕ
−q
1,i + ϕ

−n
2,i

)
+ c

(
ϕ
−q
1,i + ϕ

−n
2,i

)
≤ 0. (4.3)

Then there exists a positive constant C such that for t ∈ (0,Th)

U1(t)
m+1

C(m + 1)
≥ Th − t

(
V1(t)

n+1

C(n + 1)
≥ Th − t

)
,

U1(t) ≥ C(Th − t)
1

m+1

(
V1(t) ≥ C(Th − t)

1
n+1

)
. (4.4)

Proof. Set for i = 1, . . . ,I, t ∈ (0,Th),

Zi(t) = U
′
i(t) + c

(
U−mi (t) + V

−p
i (t)

)
and Wi(t) = V

′
i (t) + c

(
U
−q
i (t) + V

−n
i (t)

)
.

A straightforward calculation gives

Z′i(t) − δ
2Zi(t) + αi(t)Zi(t) + βi(t)Wi(t) ≤ 0, i = 1, . . . ,I, t ∈ (0,Th),

W ′i (t) − δ
2Wi(t) + ai(t)Zi(t) + bi(t)Wi(t) ≤ 0, i = 1, . . . ,I, t ∈ (0,Th),

Zi(0) ≤ 0, Wi(0) ≤ 0, i = 1, . . . ,I.



Int. J. Anal. Appl. 17 (6) (2019) 1044

By virtue of Lemma 2.1

Zi(t) ≤ 0, Wi(t) ≤ 0, i = 1, . . . ,I, t ∈ (0,Th).

Thus we get

U′i(t) ≤−cU
−m
i (t) and V

′
i (t) ≤−cV

−n
i (t), i = 1, . . . ,I, t ∈ (0,Th). (4.5)

Assume that Uh quenches (Vh quenches), integrating (4.5) from t to Th, we arrive at

U1(t)
m+1

C(m + 1)
≥ Th − t

(
V1(t)

n+1

C(n + 1)
≥ Th − t

)
,

which implies

U1(t) ≥ C(Th − t)
1

m+1

(
V1(t) ≥ C(Th − t)

1
n+1

)
.

�

Theorem 4.2. If p < n + 1, then there exist initial data such that Vh quenches but Uh doesn’t.

Proof. We argue by contradiction. Assuming that Uh and Vh quench simultaneously at time Th for any

initial data. We have∫ t
0

U′1(s) ds ≥
∫ Th

0

U′1(s) ds =
2

h2

∫ Th
0

U2(s) −U1(s) ds−
2

h

∫ Th
0

U−m1 (s) + V
−p
1 (s) ds

By using the Lemma 4.1, we obtain

U1(t) ≥ U1(0) +
2

h2

∫ Th
0

U2(s) − U1(s) ds−
2C

h

∫ Th
0

(Th −s)−
m

m+1 + (Th −s)−
p

n+1 ds

As p < n + 1 this integral is converged and

U1(t) ≥ C1 −C2T
1

m+1

h −C3T
n+1−p
n+1

h , with C1, C2, C3 > 0.

By summation of (1.6) we observe that

−
h

2
V ′1 (t) −

h

2
V ′I (t) −

I−1∑
i=2

hV ′i (t) = U
−q
1 (t) + V

−n
1 (t),

−
h

2
V ′1 (t) −

h

2
V ′I (t) −

I−1∑
i=2

hV ′i (t) ≥ U
−q
1 (0) + V

−n
1 (0) (4.6)

integrate (4.6) from 0 to Th, we can obtain

VI(0)
(
U
−q
1 (0) + V

−n
1 (0)

)−1
≥ Th,

then if Th is small enough (depending on Uh(0) and Vh(0)), U1(Th) ≥ c0 > 0. We have a contradiction with

the hypothesis that Uh quenches and the result is obtained as desired. �



Int. J. Anal. Appl. 17 (6) (2019) 1045

Theorem 4.3. If q ≤ n(m+1)
n+1

and p ≥ n + 1 (p ≤ m(n+1)
m+1

and q ≥ m + 1) then Uh (Vh) quenches alone

under any positive initial data.

Proof. Assume that there exists initial data such that Uh and Vh quench simultaneously at time Th. Without

lost of generality, we can suppose that this initial data satisfies (4.2)-(4.3). According to (1.6)

V ′1 (t) = δ
2V1(t) − b1(U

−q
1 (t) + V

−n
1 (t)),

V ′1 (t) ≥−b1(U
−q
1 (t) + V

−n
1 (t)),

V1(t) ≤ b1
∫ Th
t

U
−q
1 (s) + V

−n
1 (s) ds.

From Lemma 4.1 we know U1(t) ≥ C(Th−t)
1

m+1 , V1(t) ≥ C(Th−t)
1

n+1 , moreover q ≤ n(m+1)
n+1

. Hence there

exists C′ > 0 such that V1(t) ≤ C′(Th − t)
1

n+1 . Let us consider (1.5)

U′1(t) = δ
2U1(t) − b1(U−m1 (t) + V

−p
1 (t)),

U′1(t) ≤ δ
2U1(t) − b1V

−p
1 (t),

U′1(t) ≤ δ
2U1(t) − b1C′−p(Th − t)−

p
n+1 .

Integrating both sides from 0 to Th, we obtain

−U1(0) ≤ c1 − c2
∫ Th

0

(Th − t)−
p

n+1 dt.

We can see that the integral diverges if p ≥ n + 1, which is a contradiction. The result is obtained. �

Remark 4.1. Let (Uh,Vh) be the solution of (1.5)-(1.7) such that the initial data satisfies (4.2)-(4.3).

We can see of the Lemma 4.1 and the proof of Theorem 4.3 that if Uh (Vh) quenches at time Th, then

U1(t) ∼ (Th − t)
1

m+1

(
V1(t) ∼ (Th − t)

1
n+1

)
for t close enough to Th.

5. Convergence of the semidiscrete quenching time

In this section, we study the convergence of the semidiscrete quenching time. Now we will show that for

each fixed time interval [0,T] where (u,v) is defined, the solution (Uh,Vh) of (1.5)-(1.7) approximates (u,v)

when the mesh parameter h goes to zero. We denote

uh(t) = (u(x1, t), . . . ,u(xI, t))
T , vh(t) = (v(x1, t), . . . ,v(xI, t))

T ,

‖Uh(t)‖∞ = max
1≤i≤I

|Ui(t)|, ‖Uh(t)‖inf = min
1≤i≤I

|Ui(t)|.

Theorem 5.1. Assume that the problem (1.1)-(1.4) has solution

(u,v) ∈
(
C4,1 ([0, 1] × [0,T∗])

)2
and the initial data (ϕ1,h,ϕ2,h) at (1.5)-(1.7) verifies

‖ϕ1,h −uh(0)‖∞ = o(1), ‖ϕ2,h −vh(0)‖∞ = o(1) h → 0. (5.1)



Int. J. Anal. Appl. 17 (6) (2019) 1046

Then, for h small enough, the semidiscrete problem (1.5)-(1.7) has a unique solution (Uh,Vh) ∈(
C1
(
[0,T∗], RI

))2
such that

max
t∈[0,T∗]

‖Uh(t) −uh(t)‖∞ = O(‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + h2), as h → 0,

max
t∈[0,T∗]

‖Vh(t) −vh(t)‖∞ = O(‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + h2), as h → 0.

Proof. Let σ > 0 be such that

(‖u‖∞,‖v‖∞) < σ, t ∈ [0,T∗]. (5.2)

Then the problem (1.5)-(1.7) has for each h, a unique solution (Uh,Vh) ∈
(
C1
(
[0,T∗], RI

))2
. Let t(h) ≤ T∗

be the greatest value of t > 0 such that

max{‖Uh(t) −uh(t)‖∞,‖Vh(t) −vh(t)‖∞} < 1. (5.3)

The relation (5.1) implies t(h) > 0 for h small enough. Using the triangle inequality, we obtain

‖Uh(t)‖∞ ≤ 1 + σ and, ‖Vh(t)‖∞ ≤ 1 + σ for t ∈ (0, t(h)). (5.4)

Let (e1,h,e2,h)(t) = (Uh−uh,Vh−vh)(t), ∀t ∈ [0,T∗] be the discretization error. these error functions verify

e′1,i(t) = δ
2e1,i(t) + mbi(θi(t))

−m−1e1,i(t) + pbi(Θi(t))
−p−1e2,i(t) + O(h

2),

e′2,i(t) = δ
2e2,i(t) + qbi(θi(t))

−q−1e1,i(t) + nbi(Θi(t))
−n−1e2,i(t) + O(h

2),

where θi(t) and Θi(t) lie, respectively, between Ui(t) and u(xi, t), and between Vi(t) and v(xi, t), for i ∈

{1, . . . ,I}. Using (5.2) and (5.4), there exist K and L positive constants such that

e′1,i(t) ≤ δ
2e1,i(t) + biL|e1,i(t)| + biL|e2,i(t)| + Kh2,

e′2,i(t) ≤ δ
2e2,i(t) + biL|e1,i(t)| + biL|e2,i(t)| + Kh2

let (z,w) ∈
(
C4,1 ([0, 1], [0,T∗])

)2
be such that

z(x,t) =
(
‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + Qh2

)
e(M+2)t−(1−x)

2

and w = z, ∀(x,t) ∈ [0, 1] × [0,T∗],

with M, Q positive constants. We can prove by the Lemma 2.2 that

|e1,i(t)| < z(xi, t), |e2,i(t)| < w(xi, t), 1 ≤ i ≤ I, for t ∈ (0, t(h)).

We deduce that

‖Uh(t) −uh(t)‖∞ ≤
(
‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + Qh2

)
e(M+2)t,

‖Vh(t) −vh(t)‖∞ ≤
(
‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + Qh2

)
e(M+2)t,

for t ∈ (0, t(h)). Suppose that T∗ > t(h) from (5.3), we obtain

1 = ‖Uh(t(h)) −uh(t(h))‖∞ ≤
(
‖ϕ1,h −uh(0)‖∞ + ‖ϕ2,h −vh(0)‖∞ + Qh2

)
e(M+2)t.



Int. J. Anal. Appl. 17 (6) (2019) 1047

Since the term on the right hand side of the above inequality goes to zero as h tends to zero, we deduce that,

1 ≤ 0, which is impossible. Hence we have t(h) = T∗, and the proof is completed. �

Theorem 5.2. Let (u,v) ∈
(
C4,1([0, 1] × [0,T))

)2
be solution of (1.1)-(1.4) with quenches time T and the

initial data at (1.5)-(1.7) satisfies (4.2)-(4.3) and (5.1). Then the solution (Uh,Vh) of (1.5)-(1.7) quenches

in a finite time Th and we have

lim
h→0

Th = T.

Proof. From Theorem 3.1, (Uh,Vh) quenches in a finite time Th. Assume that Uh quenches.

Set ε > 0. There exists η > 0 such that

y1+m

C(m + 1)
≤
ε

2
, 0 ≤ y ≤ η. (5.5)

There exists a time T0 ∈ (T −
ε

2
,T) such that 0 < |u(xi, t)| ≤ η2 , for i = 1, . . . ,I, t ∈ [T0,T). Setting

T1 =
T0+T

2
, it is not hard to see that 0 < ‖u(xi, t)‖inf , for t ∈ [0,T1]. From Theorem 5.1, it follows that for

h sufficiently small

‖Uh(t) −uh(t)‖∞ ≤
η

2
.

Applying the triangle inequality, we get

‖Uh(T1)‖inf ≤‖Uh(T1) −uh(T1)‖∞ + ‖uh(T1)‖inf ≤ η.

Since Uh quenches, we can deduce from Lemma 4.1 and (5.5) that

|Th −T | ≤ |Th −T1| + |T1 −T | ≤
‖Uh(T1)‖1+minf
C(m + 1)

+
ε

2
≤ ε.

The case where Vh quenches is analogous. �

6. Numerical Experiments

In this section, we present some numerical approximations to the quenching time of (1.5)-(1.7) for the

initial data ϕ1,i = ϕ2,i = 1 +
4

Π
sin
(

Π
2

(i− 1)h
)

for i = 1, . . . ,I − 1, with different values of m, n, p and q.

We also consider the implicit scheme below

U
(k+1)
i −U

(k)
i

∆tkh
= δ2U

(k+1)
i − bi

((
U

(k)
i

)−m
+
(
V

(k)
i

)−p)
, 1 ≤ i ≤ I,

V
(k+1)
i −V

(k)
i

∆tkh
= δ2V

(k+1)
i − bi

((
U

(k)
i

)−q
+
(
V

(k)
i

)−n)
, 1 ≤ i ≤ I,

U
(0)
i = ϕ1,i, V

(0)
i = ϕ2,i 1 ≤ i ≤ I,

where k ≥ 0, ∆tkh = h
2 min

{
‖U(k)h ‖

m+1
inf ,‖U

(k)
h ‖

q+1
inf ,‖V

(k)
h ‖

p+1
inf ,‖V

(k)
h ‖

n+1
inf

}
.



Int. J. Anal. Appl. 17 (6) (2019) 1048

Definition 6.1. We say that the discrete solution (U
(k)
h ,V

(k)
h ) of the implicit scheme quenches in a finite

time if

lim
k→∞

inf{‖U(k)h ‖inf,‖V
(k)
h ‖inf} = 0 and the series

∑+∞
k=0 ∆t

k
h converges. The quantity t

k
h =

∑k−1
j=0 ∆t

j
h is

called the numerical quenching time of the solution (U
(k)
h ,V

(k)
h ) and Th =

∑+∞
k=0 ∆t

k
h is called the numerical

quenching time of the solution (Uh,Vh).

In Tables 1, 2 and 3, in rows, we present the numerical quenching times, the numbers of iterations and

the orders of the approximations corresponding to meshes of 16, 32, 64, 128, 256, 512, 1024. We take for

the numerical quenching time Th =
∑+∞
k=0 ∆t

k
h which is computed at the first time when ∆t

k
h = |t

k+1
h −t

k
h| ≤

10−16. The order(s) of the method is computed from

s =
log((T4h−T2h)/(T2h−Th))

log(2)
, where h = 1/(I − 1).

Table 1. Numerical quenching

times obtained with the implicit

Euler method for m = 0.5, p = 1,

q = 2, n = 0.5.

I Th k s

16 0.15390794 34896 -

32 0.14878519 48454 -

64 0.14737661 66676 1.86

128 0.14699264 92814 1.87

256 0.14688843 137205 1.88

512 0.14686025 238258 1.89

1024 0.14685267 545941 1.89

Table 2. Numerical quenching

times obtained with the explicit Eu-

ler method for m = 1, p = 2.5,

q = 0.5, n = 1.

I Th k s

16 0.13630655 168 -

32 0.13147195 484 -

64 0.13016571 1540 1.89

128 0.12981187 5384 1.88

256 0.12971578 20063 1.88

512 0.12968969 77515 1.88

1024 0.12968263 305050 1.88

Figure 1. On the left, no quenching of Uh and on the right, quenching of Vh for m = 0.5,

p = 1, q = 2, n = 0.5.



Int. J. Anal. Appl. 17 (6) (2019) 1049

Table 3. Numerical quenching times obtained with

the implicit Euler method for m = 0.3, p = 2, q = 2,

n = 0.3.

I Th k s

16 0.12862938 127 -

32 0.12271047 372 -

64 0.12106075 1213 1.84

128 0.12060846 4323 1.87

256 0.12048544 16309 1.88

512 0.12045217 63432 1.89

1024 0.12044321 250457 1.89

Figure 2. On the left, quenching of Uh and on the right, no quenching of Vh for m = 1,

p = 2.5, q = 0.5, n = 1.

Figure 3. On the left, quenching of Uh and on the right, quenching of Vh for m = 0.3,

p = 2, q = 2, n = 0.3.

Remark 6.1. The various tables of our numerical results show that there is a relationship between the

quenching time and flows on the boundaries. If we consider the problem (1.5)-(1.7) in the case where the

initial data ϕ2,i = ϕ1,i = 1 +
4

Π
sin
(

Π
2

(i− 1)h
)
, i = 1, . . . ,I, from figures 1-3, we observe respectively the

quenching of (Uh,Vh). From tables 1-3, we observe the convergence of quenching time Th of the solution of

(1.5)-(1.7), since the rate of convergence is near 2. This result does not surprise us because of the result



Int. J. Anal. Appl. 17 (6) (2019) 1050

established in the previous section. Moreover we can see that of the figure 1, Vh quenches while Uh doesn’t

when p < n + 1, of the figure 2, Uh quenches while Vh doesn’t when q ≤
n(m+1)
n+1

and p ≥ n + 1 and of the

figure 3, Uh and Vh quench simultaneously when p ≥ n + 1 and q ≥ m + 1. These numerical results are in

fact consistent with the Corollary 4.1, Theorem 4.2 and Theorem 4.3.

Conclusion

In this work, we proposed a semi-discrete scheme, based on finite difference method in space for system

of heat equations, coupled by nonlinear boundary flux. The stability and the convergence of semi-discrete

scheme are proved respectively. Under some conditions the semidiscrete scheme reproduces well the condi-

tions for the quenching, quenching set and simultaneous and non-simultaneous quenching. The analysis in

this paper can be extended to more general to some systems of nonlinear parabolic equations.

References

[1] K.A. Adou, K.A. Touré, A. Coulibaly, Numerical study of the blow-up time of positive solutions of semilinear heat equations,

Far East J. Appl. Math., 4 (2018), 91-308.

[2] T. K. Boni, Extinction for discretizations of some semilinear parabolic equations, C. R. Acad. Sci. Paris, 333 (2001), 79-800.

[3] T. K. Boni, H. Nachid, D. Nabongo, Quenching time of semilinear heat equations, Miskolc Math. Notes, 11 (2010), 27-41.

[4] C.Y. Chan, S.I. Yuen, Parabolic problems with nonlinear absorptions and releases at the boundaries, Appl. Math. comput.,

121 (2001), 203-209.

[5] K.B. Edja, K.A. Touré, B. J.-C. Koua, Numerical Blow-up for a Heat Equation with Nonlinear Boundary Conditions, J.

Math. Res., 10 (2018), 119-128.

[6] R. Ferreira, A. de Pablo, M.P. LLanos, J.D. Rossi, Incomplete quenching in a system of heat equations coupled at the

boundary, J. Math. Anal. Appl., 346 (2008), 1145-154.

[7] R.H. Ji, C.Y. Quc, L.D. Wang, Simultaneous and non-simultaneous quenching for coupled parabolic system, Appl. Anal.,

94 (2015), 233-250.

[8] H. Kawarada, On Solutions of Initial-Boundary Problem for ut = uxx +
1

1−u , Publ. RIMS, Kyoto Univ., 10 (1975), 729-736.

[9] H.A. Levine, The Quenching of Solutions of Linear Parabolic and Hyperbolic Equations with Nonlinear Boundary Condi-

tions, SIAM J. Math. Anal., 4 (1983), 1139-1153.

[10] H.A. Levine, J.T. Montgomer The quenching of solutions of some nonlinear parabolic equations, SIAM J. Math. Anal., 11

(1980), 842-847.

[11] Liang K.W, Lin AP, Tan RCE. Numerical Solution of Quenching Problems Using Mesh-Dependent Variable Temporal

Steps, Appl. Numer. Math., 57 (2007), 791-800.

[12] D. Nabongo, T.K. Boni, Quenching for semidiscretizations of a heat equation with a singular boundary condition, Asymp-

totic Anal., 59 (2008), 27-38.

[13] K.C. N’dri, K.A. Touré, G. Yoro, Numerical blow-up time for a parabolic equation with nonlinear boundary conditions,

Int. J. Numer. Methods Appl., 17 (2018), 141-160.

[14] H. Pei, Z. Li, Quenching for a parabolic system with general singular terms, J. Nonlinear Sci. Appl., 7 (2016), 1-10.

[15] B. Seluk, Quenching behavior of a semilinear reaction-diffusion system with singular boundary condition, Turk. J. Math.,

40 (2016), 166-180.



Int. J. Anal. Appl. 17 (6) (2019) 1051

[16] M. Taha, K. Toure, E. Mensah, Numerical approximation of the blow-up time for a semilinear parabolic equation with

nonlinear boundary equation, Far East J. Appl. Math., 60 (2012), 125-167.

[17] S.N. Zheng, X.F. Song, Quenching rates for heat equations with coupled nonlinear boundary flux, Sci. China Ser. A., 51

(2008), 1631-1643.


	1. Introduction
	2. Properties of the semidiscrete scheme
	3. Quenching and blow-up
	4. Simultaneous vs. non-simultaneous quenching
	5. Convergence of the semidiscrete quenching time
	6. Numerical Experiments
	Conclusion
	References