# The Method of Lines for Ternary Diffusion Problems

Authors: Henryk Leszczyński; Milena Matusik
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The method of lines (MOL) for diffusion equations with Neumann boundary conditions is considered. These equations are transformed by a discretization in space variables into systems of ordinary differential equations. The proposed ODEs satisfy the mass conservation law. The stability of solutions of these ODEs with respect to discreteL
2 norms and discrete W
1
,
∞ norms is investigated. Numerical examples confirm the parabolic behaviour of this model and very regular dynamics.


---

## Body

## 1. Introduction

Diffusion is one of several transport phenomena that occur in nature. Many physical processes in metallurgy are controlled by diffusion, for example, oxidation and sintering. There are two ways to introduce the notion of diffusion: an atomistic approach and continuum approach. The central role in the diffusion theory of multicomponent systems plays the Kirkendall effect [1], i.e., the discovery that drift in solids can be generated by diffusion. The explanation of the effect and an original method of its mathematical analysis was proposed by Darken [2]. The Kirekndall effect describes what happens when two solids diffuse into each other at different rates. This subject is treated in many textbooks on solid-state and publications [3–7]. In [8, 9] the authors deal with stability of Kirkendall planes. The peculiarities of the deformation in one-dimensional binary diffusion couple with a moving interface are described in [10].Our ternary system arises from the Fick second law. If the overall molar densityc
=
c
(
t
,
x
) and component molar ratio N
i
=
N
i
(
t
,
x
) are defined by
(1)
c
=
∑
i
=
1
3
‍
c
i
,
(
c
i is the molar concentration of the ith component in the mixture, [mole/dm3])
(2)
N
i
=
c
i
c
,
where this quatity is a ratio of ith component atoms to averall number of atoms. Then we have
(3)
∂
N
i
∂
t
=
D
i
∂
2
N
i
∂
x
2
-
∂
(
N
i
ν
D
)
∂
x
,
i
=
1,2
,
3
,
where D
i denotes the intrinsic diffusivity of the ith component and ν
D is the drift velocity. In the paper N
i, i
=
1,2
,
3 will be denoted by u, v, and w, respectively.The method proposed by Darken [2], was extended by Danielewski et al. in [11, 12]. Their generalization of Darken’s method describes the interdiffusion process in bounded mixtures with constant concentrations and with variable diffusivities of the elements.We are interested in establishing an approximation method of solutions to the ternary system mentioned above by solutions of associated systems of ordinary differential equations. These systems of ordinary differential equations are obtained by using a discretization in spatial variables.From the abundant literature concerning the numerical method of lines (MOL) for classical PDEs we mention the monographs [13–15]. MOL for time-dependent one-dimensional systems of parabolic partial differential equations is proposed in [16]. Numerical examples for problems such as Burger’s equation or nonlinear diffusion equation with nonlinear boundary conditions were presented in [16]. A posteriori error estimates for evolution problems (general nonlinear parabolic problem with a strongly monotonic elliptic operator, a linear nonstationary convection-diffusion problem, and linear second-order hyperbolic problem) solved by the method of lines are derived in [17].The aim of the paper is to construct a method of lines for diffusion equations in three-component system with Neumann boundary conditions. We prove stability of approximate solutions with respect to discreteL
2 norms and discrete W
1
,
∞ norms.Our research is concerned with the situation where the initial data and solutions are at least of classC
2. This regularity is assumed throughout the paper. It is observed that even very irregular data in a short time give highly smooth solutions (concentrations of the components). Moreover, these concentrations are mixed, which means that they are inside the interval (
0,1
) for t
>
0. Therefore, we assume that the initial data are of class C
2. This observation also justifies our assumptions that the Lipschitz constants and second derivatives of solutions are bounded by some constant independent of the discretization parameter. Under these assumptions we prove that solutions of ODEs satisfy the maximum principle and remain inside the interval (
0,1
). We also show convergence and stability of the approximate solutions. Our results confirm the parabolic nature of the equations in question. The existence and uniqueness of solutions of the interdiffusion problem in W
1,2 space was shown by the variational methods in [12]. In [18] authors constructed an example of an interdiffusion equation that does not fulfill a peculiar type of parabolicity, meant in the sense of a decreasing L
2 norm. In Section 4 we show the convergence of MOL in discrete W
1
,
∞ norms, while in Section 3 we attempt to show convergence in discrete L
2 norms. It turns out that we have stability in L
2 norms for certain values of the diffusion coefficients D
1, D
2, D
3, which highlights the parabolic nature of the equation.In Section5 we present numerical examples where we consider Darken’s trajectory besides solutions. Since the right-hand sides of ODEs are polynomial, we attempt in the appendix to estimate the radius of convergence of Taylor series solutions. It turns out that this radius is proportional to h
2. This means that Cauchy-Kovalevskaya theorem has limited (restricted) applicability to our system. Because R-K methods are close to Taylor expansions, our observation points out that one should be careful with these methods.


## 2. Formulation of the Problem

We consider the system of diffusion equations:(4)
∂
u
∂
t
=
D
1
∂
2
u
∂
x
2
-
∂
∂
x
(
u
ν
D
)
,
∂
v
∂
t
=
D
2
∂
2
v
∂
x
2
-
∂
∂
x
(
v
ν
D
)
,
∂
w
∂
t
=
D
3
∂
2
w
∂
x
2
-
∂
∂
x
(
w
ν
D
)
,
on [
0
,
T
]
×
[
-
L
,
L
]
⊂
ℝ
2 with the initial boundary conditions
(5)
u
(
0
,
x
)
=
u
0
(
x
)
,
v
(
0
,
x
)
=
v
0
(
x
)
,
w
(
0
,
x
)
=
w
0
(
x
)
for
x
∈
[
-
L
,
L
]
∂
u
∂
x
=
0
,
∂
v
∂
x
=
0
,
∂
w
∂
x
=
0
for
x
=
±
L
,
where the diffusion coefficients D
1, D
2, D
3 are different and positive, say D
1
>
D
2
>
D
3
>
0, u
0
,
v
0
,
w
0
∈
C
2
(
[
-
L
,
L
]
,
ℝ
+
), and u
0
+
v
0
+
w
0
=
1 for x
∈
[
-
L
,
L
], and
(6)
ν
D
=
D
1
∂
u
∂
x
+
D
2
∂
v
∂
x
+
D
3
∂
w
∂
x
,
is the drift velocity. These data imply u
+
v
+
w
≡
1.We formulate the method of lines corresponding to (4) on a regular mesh with step h such that N
h
=
L. The difference operators δ
(
2
), δ
+, δ
- are defined in the following way:
(7)
δ
+
u
(
i
)
=
u
(
i
+
1
)
-
u
(
i
)
h
,
δ
-
u
(
i
)
=
u
(
i
)
-
u
(
i
-
1
)
h
,
δ
(
2
)
u
(
i
)
=
u
(
i
+
1
)
-
2
u
(
i
)
+
u
(
i
-
1
)
h
2
=
δ
+
δ
-
u
(
i
)
.Denote byv
(
j
) the discrete drift velocity v
(
j
)
=
D
1
δ
+
u
(
j
)
+
D
2
δ
+
v
(
j
)
+
D
3
δ
+
w
(
j
). We will consider the following ODE system:
(8)
d
d
t
u
(
j
)
=
D
1
δ
(
2
)
u
(
j
)
-
u
(
j
)
δ
-
v
(
j
)
-
1
2
δ
+
u
(
j
)
·
v
(
j
)
-
1
2
δ
-
u
(
j
)
·
v
(
j
-
1
)
,
d
d
t
v
(
j
)
=
D
2
δ
(
2
)
v
(
j
)
-
v
(
j
)
δ
-
v
(
j
)
-
1
2
δ
+
v
(
j
)
·
v
(
j
)
-
1
2
δ
-
v
(
j
)
·
v
(
j
-
1
)
,
d
d
t
w
(
j
)
=
D
3
δ
(
2
)
w
(
j
)
-
w
(
j
)
δ
-
v
(
j
)
-
1
2
δ
+
w
(
j
)
·
v
(
j
)
-
1
2
δ
-
w
(
j
)
·
v
(
j
-
1
)
,
for j
=
-
N
,
…
,
N with the initial conditions
(9)
u
(
j
)
(
0
)
=
u
0
(
h
j
)
≥
0
,
v
(
j
)
(
0
)
=
v
0
(
h
j
)
≥
0
,
w
(
j
)
(
0
)
=
w
0
(
h
j
)
≥
0
,
j
=
-
N
,
…
,
N
.
The discrete Neumann boundary conditions
(10)
u
(
N
+
1
)
=
u
(
N
-
1
)
,
v
(
N
+
1
)
=
v
(
N
-
1
)
,
w
(
N
+
1
)
=
w
(
N
-
1
)
,
u
(
-
N
-
1
)
=
u
(
-
N
+
1
)
,
v
(
-
N
-
1
)
=
v
(
-
N
+
1
)
,
w
(
-
N
-
1
)
=
w
(
-
N
+
1
)
can be regarded as a convenient definition of auxiliary quantities: u
(
N
+
1
), v
(
N
+
1
), w
(
N
+
1
), u
(
-
N
-
1
), v
(
-
N
-
1
), w
(
-
N
-
1
).Remark 1.
Local existence and uniqueness of solutions of (8) with the above conditions follows from the general theory of ODE systems.Suppose that we have two constantsM
0, M, independent of h, satisfying the inequality M
≥
M
0
>
0. This assumption is valid throughout the paper. M
0 will be the Lipschitz constant for the initial functions u
0, v
0, w
0. M will play the same role for u, v, w
(11)
𝔛
0
=
{
(
u
0
,
v
0
,
w
0
)
∈
ℝ
3
(
2
N
+
1
)
:
u
0
(
j
)
,
v
0
(
j
)
,
w
0
(
j
)
≥
0
,
u
0
(
j
)
+
v
0
(
j
)
+
w
0
(
j
)
=
1
for
j
=
-
N
,
…
,
N
|
u
0
(
i
)
-
u
0
(
j
)
|
≤
M
0
|
i
-
j
|
h
,
|
v
0
(
i
)
-
v
0
(
j
)
|
≤
M
0
|
i
-
j
|
h
|
w
0
(
i
)
-
w
0
(
j
)
|
≤
M
0
|
i
-
j
|
h
}
,
𝔛
=
{
(
u
,
v
,
w
)
:
[
0
,
T
]
⟶
ℝ
3
(
2
N
+
1
)
:
u
(
j
)
,
v
(
j
)
,
w
(
j
)
≥
0
,
u
(
j
)
+
v
(
j
)
+
w
(
j
)
=
1
for
j
=
-
N
,
…
,
N
|
u
(
i
)
-
u
(
j
)
|
≤
M
|
i
-
j
|
h
,
|
v
(
i
)
-
v
(
j
)
|
≤
M
|
i
-
j
|
h
|
w
(
i
)
-
w
(
j
)
|
≤
M
|
i
-
j
|
h
}
.
We will refer to these Lipschitz-type inequalities in the definitions of 𝔛
0 and 𝔛 as M
0-Lipschitz and M-Lipschitz, respectively. It is assumed that approximate solutions are Lipschitz continuous with respect to spatial variable. The Lipschitz continuity of approximate solutions can be proved rigorously; however, we omit the proof for brevity. To justify our proceeding, consider solutions of ternary systems. They are bounded together with their first and second derivatives.If the initial data belong to𝔛
0 the equations for w can be eliminated from (8). The following lemma allows us to express (8) as the system of 2
(
2
N
+
1
) ordinary differential equations.Lemma 2.
Suppose that(
u
,
v
,
w
) is the solution of (8) with the initial condition (
u
0
,
v
0
,
w
0
)
∈
𝔛
0. Assume that u, v, w are M-Lipschitz. Then (
u
,
v
,
w
)
∈
𝔛.Proof.
First, we prove thatu
(
j
)
+
v
(
j
)
+
w
(
j
)
=
1 for j
=
-
N
,
…
,
N. Let g
=
u
+
v
+
w. We have
(12)
d
d
t
(
g
(
j
)
)
=
[
1
-
g
(
j
)
]
δ
-
v
(
j
)
-
1
2
δ
+
(
g
(
j
)
)
·
v
(
j
)
-
1
2
δ
-
(
g
(
j
)
)
·
v
(
j
-
1
)
,
for j
=
-
N
,
…
,
N with the initial condition g
(
j
)
(
0
)
=
1 and the discrete Neumann boundary conditions. It follows from the general ODE theory that there exists exactly one solution g
(
j
)
≡
1.
We next show thatu
(
j
)
,
v
(
j
)
,
w
(
j
)
≥
0 on [
0
,
T
] for j
=
-
N
,
…
,
N. Suppose first that u
(
j
)
(
0
)
,
v
(
j
)
(
0
)
,
w
(
j
)
(
0
)
>
0. We claim that u
(
j
)
,
v
(
j
)
,
w
(
j
)
>
0 on [
0
,
T
] for j
=
-
N
,
…
,
N. To obtain a contradiction, we suppose that there is t
~
>
0 such that u
(
j
)
(
t
)
,
v
(
j
)
(
t
)
,
w
(
j
)
(
t
)
>
0 for t
<
t
~ and j
=
-
N
,
…
,
N and there is i
∈
{
-
N
,
…
,
N
} such that u
(
i
)
(
t
~
)
=
0. We have u
(
j
)
(
t
~
)
,
v
(
j
)
(
t
~
)
,
w
(
j
)
(
t
~
)
≥
0 for j
=
-
N
,
…
,
N. Then
(13)
d
d
t
u
(
i
)
(
t
~
)
≤
0
,
δ
+
u
(
i
)
(
t
~
)
≥
0
,
δ
-
u
(
i
)
(
t
~
)
≤
0
.
We exploit (8) at the point t
~ to get the following inequality:
(14)
0
≥
d
d
t
u
(
i
)
(
t
~
)
=
D
1
δ
(
2
)
u
(
i
)
(
t
~
)
-
1
2
δ
+
u
(
i
)
(
t
~
)
v
(
i
)
-
1
2
δ
-
u
(
i
)
(
t
~
)
v
(
i
-
1
)
=
δ
+
u
(
i
)
(
t
~
)
[
D
1
h
-
v
(
i
)
]
-
δ
-
u
(
i
)
(
t
~
)
[
D
1
h
+
v
(
i
-
1
)
]
.
The expressions in the brackets are positive for small h. Hence
(15)
δ
+
u
(
i
)
(
t
~
)
=
0
,
δ
-
u
(
i
)
(
t
~
)
=
0
.
We see at once that u
(
i
)
(
t
~
)
=
0 for i
=
-
N
,
…
,
N. It follows from the general ODE theory [19, page 108] that u
(
t
)
=
0 for t
∈
[
0
,
T
], contrary to strict initial inequalities. This argument also settles the case of v
(
i
)
(
t
~
)
=
0 or w
(
i
)
(
t
~
)
=
0. We now turn to the case u
(
0
,
·
)
,
v
(
0
,
·
)
,
w
(
0
,
·
)
≥
0. It follows from the theorem on the continuous dependence on initial data [19, page 145] that u
(
j
)
,
v
(
j
)
,
w
(
j
)
≥
0  for  j
=
-
N
,
…
,
N.Denote byh
∑
j
′ the trapezoidal sum. More precisely,
(16)
h
∑
j
‍
′
u
(
j
)
=
h
∑
j
=
-
N
N
-
1
‍
u
(
j
)
+
u
(
j
+
1
)
2
.
This quantity corresponds to the total mass of u. The MOL system is conservative. It preserves the mass of u, v.Lemma 3.
Suppose thatu, v, w are solutions of (8) such that (
u
0
,
v
0
,
w
0
)
∈
𝔛
0. Then
(17)
h
∑
j
‍
′
u
(
j
)
=
c
o
n
s
t
,
h
∑
j
‍
′
v
(
j
)
=
c
o
n
s
t
.Proof.
Note thatδ
+
u
(
j
)
=
δ
-
u
(
j
+
1
) for j
=
-
N
,
…
N. Hence
(18)
h
∑
j
‍
′
D
1
δ
(
2
)
u
(
j
)
=
h
∑
j
‍
′
[
D
1
δ
+
u
(
j
)
-
D
1
δ
+
u
(
j
-
1
)
]
=
0
,
h
∑
j
‍
′
(
-
u
(
j
)
δ
-
v
(
j
)
-
δ
+
u
(
j
)
2
·
v
(
j
)
-
δ
-
u
(
j
)
2
·
v
(
j
-
1
)
)
=
h
D
1
′
2
∑
j
‍
′
(
-
u
(
j
)
δ
+
u
(
j
)
+
u
(
j
)
δ
+
u
(
j
-
1
)
-
u
(
j
+
1
)
δ
+
u
(
j
)
+
u
(
j
-
1
)
δ
+
u
(
j
-
1
)
)
+
h
D
2
′
2
∑
j
‍
′
(
-
u
(
j
)
δ
+
v
(
j
)
+
u
(
j
)
δ
+
v
(
j
-
1
)
-
u
(
j
+
1
)
δ
+
v
(
j
)
+
u
(
j
-
1
)
δ
+
v
(
j
-
1
)
)
=
0
,
where D
1
′
=
D
1
-
D
3 and D
2
′
=
D
2
-
D
3. It follows from these equations that
(19)
d
d
t
(
h
∑
j
‍
′
u
(
j
)
)
=
h
∑
j
‍
′
d
d
t
u
(
j
)
=
0
.


## 3. Stability with respect to DiscreteL
2 Norms

The stability of the method of lines with respect to discreteL
2 norm is established by the following theorem.Theorem 4.
Suppose that(1)
(
u
,
v
,
w
)
∈
𝔛 is the unique solution of system (8) with the initial condition (
u
0
,
v
0
,
w
0
)
∈
𝔛
0,
(2)
(
u
¯
,
v
¯
,
w
¯
)
∈
𝔛 is an approximate solution of (8) with the initial condition (
u
¯
0
,
v
¯
0
,
w
¯
0
)
∈
𝔛
0 and with perturbations of the right-hand sides denoted by ξ
1
(
j
), ξ
2
(
j
), ξ
3
(
j
) for j
=
-
N
,
…
,
N,
(3)
consider(20)
-
D
2
2
+
14
D
1
D
2
-
D
1
2
>
0
,
-
D
3
2
+
14
D
1
D
3
-
D
1
2
>
0
,
-
D
3
2
+
14
D
2
D
3
-
D
2
2
>
0
.
Then we have(21)
A
t
≤
C
1
[
A
0
+
∥
ξ
∥
]
exp
⁡
⁡
(
K
1
t
)
,
where K
1
,
C
1
>
0 are constants independent of h and
(22)
∥
ξ
∥
=
max
⁡
0
≤
t
≤
T
h
∑
j
‍
′
[
(
ξ
1
(
j
)
(
t
)
)
2
+
(
ξ
2
(
j
)
(
t
)
)
2
+
(
ξ
3
(
j
)
(
t
)
)
2
]
,
A
t
=
h
∑
j
‍
′
[
(
u
¯
(
j
)
(
t
)
-
u
(
j
)
(
t
)
)
2
+
(
v
¯
(
j
)
(
t
)
-
v
(
j
)
(
t
)
)
2
+
(
w
¯
(
j
)
(
t
)
-
w
(
j
)
(
t
)
)
2
]
.Proof.
Denote byΔ
u, Δ
v, Δ
w the differences between the approximate and exact solutions of (8).
The proof is based on the following observation:(23)
d
d
t
A
t
=
2
h
∑
j
‍
′
[
Δ
u
(
j
)
d
d
t
(
u
¯
(
j
)
-
u
(
j
)
)
+
Δ
v
(
j
)
d
d
t
(
v
¯
(
j
)
-
v
(
j
)
)
+
Δ
w
(
j
)
d
d
t
(
w
¯
(
j
)
-
w
(
j
)
)
]
=
h
∑
j
‍
′
2
Δ
u
(
j
)
{
1
2
D
1
δ
(
2
)
Δ
u
(
j
)
-
u
¯
(
j
)
δ
-
v
¯
(
j
)
+
u
(
j
)
δ
-
v
(
j
)
11111111111
-
1
2
δ
+
u
¯
(
j
)
·
v
¯
(
j
)
+
1
2
δ
+
u
(
j
)
·
v
(
j
)
11111111111
-
1
2
δ
-
u
¯
(
j
)
·
v
¯
(
j
-
1
)
+
1
2
δ
-
u
(
j
)
·
v
(
j
-
1
)
}
+
h
∑
j
‍
′
2
Δ
v
(
j
)
{
d
d
t
(
v
¯
(
j
)
-
v
(
j
)
)
}
+
h
∑
j
‍
′
2
Δ
w
(
j
)
{
d
d
t
(
w
¯
(
j
)
-
w
(
j
)
)
}
+
h
∑
j
‍
′
{
2
Δ
u
(
j
)
ξ
1
(
j
)
+
2
Δ
v
(
j
)
ξ
2
(
j
)
+
2
Δ
w
(
j
)
ξ
3
(
j
)
}
.
The terms (
d
/
d
t
)
(
v
¯
(
j
)
-
v
(
j
)
) and (
d
/
d
t
)
(
w
¯
(
j
)
-
w
(
j
)
) are calculated similarly as (
d
/
d
t
)
(
u
¯
(
j
)
-
u
(
j
)
). Set
(24)
Z
1
=
δ
+
Δ
u
(
j
)
(
t
)
,
Z
2
=
δ
+
Δ
v
(
j
)
(
t
)
,
Z
3
=
δ
+
Δ
w
(
j
)
(
t
)
.
We apply the summation by parts and obtain
(25)
d
d
t
A
t
=
h
∑
j
=
-
N
N
-
1
{
-
2
D
1
(
Z
1
)
2
-
2
D
2
(
Z
2
)
2
-
2
D
3
(
Z
3
)
2
+
Z
1
(
Δ
u
(
j
+
1
)
+
Δ
u
(
j
)
)
·
v
¯
(
j
)
+
Z
2
(
Δ
v
(
j
+
1
)
+
Δ
v
(
j
)
)
·
v
¯
(
j
)
+
Z
3
(
Δ
w
(
j
+
1
)
+
Δ
w
(
j
)
)
·
v
¯
(
j
)
+
D
1
(
Z
1
)
2
(
u
(
j
+
1
)
+
u
(
j
)
)
+
D
2
Z
1
Z
2
(
u
(
j
+
1
)
+
u
(
j
)
)
+
D
3
Z
1
Z
3
(
u
(
j
+
1
)
+
u
(
j
)
)
+
D
1
Z
1
Z
2
(
v
(
j
+
1
)
+
v
(
j
)
)
+
D
2
(
Z
2
)
2
(
v
(
j
+
1
)
+
v
(
j
)
)
+
D
3
Z
2
Z
3
(
v
(
j
+
1
)
+
v
(
j
)
)
+
D
1
Z
1
Z
3
(
w
(
j
+
1
)
+
w
(
j
)
)
+
D
2
Z
2
Z
3
(
w
(
j
+
1
)
+
w
(
j
)
)
+
D
3
(
Z
3
)
2
(
w
(
j
+
1
)
+
w
(
j
)
)
}
+
h
∑
j
‍
′
{
2
Δ
u
(
j
)
ξ
1
(
j
)
+
2
Δ
v
(
j
)
ξ
2
(
j
)
+
2
Δ
w
(
j
)
ξ
3
(
j
)
}
.
From Lemma3 we have Z
3
=
-
Z
1
-
Z
2. Set
(26)
2
A
(
i
)
=
(
u
(
i
+
1
)
+
u
(
i
)
)
,
2
B
(
i
)
=
(
v
(
i
+
1
)
+
v
(
i
)
)
,
2
C
(
i
)
=
(
w
(
i
+
1
)
+
w
(
i
)
)
.
Hence the sign of (
d
/
d
t
)
A
t depends on the quadratic form T
=
[
T
i
j
]
i
,
j
=
1,2:
(27)
T
11
=
-
2
(
D
1
+
D
3
)
+
2
(
D
1
-
D
3
)
(
A
(
i
)
-
C
(
i
)
)
,
T
12
=
T
21
=
-
2
D
3
+
D
2
′
(
A
(
i
)
-
C
(
i
)
)
+
D
1
′
(
B
(
i
)
-
C
(
i
)
)
,
T
22
=
-
2
(
D
2
+
D
3
)
+
2
D
2
′
(
B
(
i
)
-
C
(
i
)
)
,
which is negative definite for D
1, D
2, D
3 satisfying condition (20) of Theorem 4. We have
(28)
d
d
t
A
t
≤
-
K
h
∑
j
=
-
N
N
-
1
‍
(
Z
1
2
+
Z
2
2
+
Z
3
2
)
+
h
∑
j
‍
′
(
2
Δ
u
(
j
)
ξ
1
(
j
)
+
2
Δ
v
(
j
)
ξ
2
(
j
)
+
2
Δ
w
(
j
)
ξ
3
(
j
)
)
+
h
∑
j
=
-
N
N
-
1
‍
{
(
D
1
δ
+
u
¯
(
j
)
+
D
2
δ
+
v
¯
(
j
)
+
D
3
δ
+
w
¯
(
j
)
)
Z
1
(
Δ
u
(
j
+
1
)
+
Δ
u
(
j
)
)
1111111
×
(
D
1
δ
+
u
¯
(
j
)
+
D
2
δ
+
v
¯
(
j
)
+
D
3
δ
+
w
¯
(
j
)
)
1111111
+
Z
2
(
Δ
v
(
j
+
1
)
+
Δ
v
(
j
)
)
1111111
×
(
D
1
δ
+
u
¯
(
j
)
+
D
2
δ
+
v
¯
(
j
)
+
D
3
δ
+
w
¯
(
j
)
)
1111111
+
Z
3
(
Δ
w
(
j
+
1
)
+
Δ
w
(
j
)
)
1111111
×
(
D
1
δ
+
u
¯
(
j
)
+
D
2
δ
+
v
¯
(
j
)
+
D
3
δ
+
w
¯
(
j
)
)
}
=
I
1
+
I
2
+
I
3
.
It follows from the inequality
(29)
|
a
b
|
≤
1
2
(
K
a
2
+
b
2
K
)
for
K
>
0
that
(30)
|
I
2
|
≤
h
K
∑
j
‍
′
(
Δ
u
(
j
)
)
2
+
(
Δ
v
(
j
)
)
2
+
(
Δ
w
(
j
)
)
2
+
K
∥
ξ
∥
,
|
I
3
|
≤
h
M
1
(
D
1
+
D
2
+
D
3
)
×
∑
j
=
-
N
N
-
1
‍
(
Z
1
Δ
u
(
j
+
1
)
+
Z
2
Δ
v
(
j
+
1
)
+
Z
3
Δ
w
(
j
+
1
)
)
+
h
M
1
(
D
1
+
D
2
+
D
3
)
×
∑
j
=
-
N
N
-
1
‍
(
Z
1
Δ
u
(
j
)
+
Z
2
Δ
v
(
j
)
+
Z
3
Δ
w
(
j
)
)
≤
K
h
∑
j
=
-
N
N
-
1
‍
(
Z
1
2
+
Z
2
2
+
Z
3
2
)
+
h
M
1
2
(
D
1
+
D
2
+
D
3
)
2
K
×
∑
j
‍
′
(
(
Δ
u
(
j
)
)
2
+
(
Δ
v
(
j
)
)
2
+
(
Δ
w
(
j
)
)
2
)
.
Hence
(31)
d
d
t
A
t
≤
-
K
h
∑
j
=
-
N
N
-
1
‍
(
Z
1
2
+
Z
2
2
+
Z
3
2
)
+
K
h
∑
j
=
-
N
N
-
1
‍
(
Z
1
2
+
Z
2
2
+
Z
3
2
)
+
(
1
+
M
1
2
(
D
1
+
D
2
+
D
3
)
2
)
K
A
t
+
K
∥
ξ
∥
.
We conclude from this differential inequality that
(32)
A
t
≤
C
1
(
A
0
+
∥
ξ
∥
)
exp
⁡
⁡
(
(
1
+
M
1
2
(
D
1
+
D
2
+
D
3
)
2
)
K
t
)
,
where
(33)
C
1
=
K
2
(
1
+
M
1
2
(
D
1
+
D
2
+
D
3
)
2
)
.
This finishes the proof.


## 4. Stability with respect to DiscreteW
1
,
∞ Norms

We define the norm∥
·
∥
1
,
∞ in ℝ
2
N
+
1 as follows:
(34)
∥
u
∥
1
,
∞
=
∥
u
∥
∞
+
∥
δ
+
u
∥
∞
,
where ∥
u
∥
∞
=
max
⁡
-
N
≤
j
≤
N
⁡
|
u
(
j
)
| (with the discrete Neumann convention).Theorem 5.
Suppose that(1)
(
u
,
v
,
w
)
∈
𝔛 is the exact solution of (8) with the initial condition (
u
0
,
v
0
,
w
0
)
∈
𝔛
0,
(2)
(
u
¯
,
v
¯
,
w
¯
)
∈
𝔛 is an approximate solution of (8) with the initial condition (
u
¯
0
,
v
¯
0
,
w
¯
0
)
∈
𝔛
0 and with perturbations of the right-hand sides denoted by ξ
1
(
j
), ξ
2
(
j
), ξ
3
(
j
) for j
=
-
N
,
…
,
N.
Then we have(35)
∥
Δ
u
(
t
)
∥
1
,
∞
+
∥
Δ
v
(
t
)
∥
1
,
∞
≤
K
~
∥
ξ
∥
∞
t
1
/
2
,
where K
~ is some constant independent of h.Proof.
It is easily seen thatΔ
u satisfies differential equation
(36)
d
d
t
Δ
u
(
i
)
=
[
-
1
2
δ
-
Δ
u
(
i
)
(
D
1
′
δ
-
u
¯
(
i
)
+
D
2
′
δ
-
v
¯
(
i
)
)
D
1
δ
(
2
)
Δ
u
(
i
)
-
u
(
i
)
(
D
1
′
δ
(
2
)
Δ
u
(
i
)
+
D
2
′
δ
(
2
)
Δ
v
(
i
)
)
-
Δ
u
(
i
)
(
D
1
′
δ
(
2
)
u
¯
(
i
)
+
D
2
′
δ
(
2
)
v
¯
(
i
)
)
-
1
2
δ
+
u
(
i
)
(
D
1
′
δ
+
Δ
u
(
i
)
+
D
2
′
δ
+
Δ
v
(
i
)
)
-
1
2
δ
+
Δ
u
(
i
)
(
D
1
′
δ
+
u
¯
(
i
)
+
D
2
′
δ
+
v
¯
(
i
)
)
-
1
2
δ
-
u
(
i
)
(
D
1
′
δ
-
Δ
u
(
i
)
+
D
2
′
δ
-
Δ
v
(
i
)
)
-
1
2
δ
-
Δ
u
(
i
)
(
D
1
′
δ
-
u
¯
(
i
)
+
D
2
′
δ
-
v
¯
(
i
)
)
]
+
ξ
1
(
i
)
(
t
)
.
Δ
v satisfies a similar differential equation.
Suppose that we have Green functionsG
1
,
i
,
j, G
2
,
i
,
j corresponding to the differential-difference operator:
(37)
[
d
d
t
-
D
1
δ
(
2
)
+
D
1
′
u
(
i
)
δ
(
2
)
D
2
′
u
(
i
)
δ
(
2
)
D
1
′
v
(
i
)
δ
(
2
)
d
d
t
-
D
2
δ
(
2
)
+
D
2
′
v
(
i
)
δ
(
2
)
]
such that
(38)
h
∑
j
‍
G
k
,
i
,
j
(
t
,
s
)
≤
C
,
h
∑
j
‍
δ
+
G
k
,
i
,
j
(
t
,
s
)
≤
C
t
-
s
,
where C
>
0 and k
=
1,2 (see [20]). Using the above Green functions we express the differences Δ
u
(
i
) and Δ
v
(
i
) as follows:
(39)
[
Δ
u
(
i
)
(
t
)
Δ
v
(
i
)
(
t
)
]
=
∫
0
t
‍
h
∑
j
‍
[
G
1
,
i
,
j
(
t
,
s
)
P
1
,
j
(
s
)
G
2
,
i
,
j
(
t
,
s
)
P
2
,
j
(
s
)
]
d
s
,
where P
1
,
j
(
t
), P
2
,
j
(
t
) depend on Δ
u
(
j
), Δ
v
(
j
), δ
+
Δ
u
(
j
), δ
+
Δ
u
(
j
-
1
), δ
+
Δ
v
(
j
), δ
+
Δ
v
(
j
-
1
) and satisfy the estimates
(40)
∥
P
i
(
t
)
∥
∞
≤
C
0
(
∥
Δ
u
(
t
)
∥
∞
+
∥
δ
+
Δ
u
(
t
)
∥
∞
)
+
C
0
(
∥
Δ
v
(
t
)
∥
∞
+
∥
δ
+
Δ
v
(
t
)
∥
∞
)
+
∥
ξ
i
(
t
)
∥
∞
,
i
=
1,2
,
where C
0 depends on M
1, D
1
′, D
2
′. Hence
(41)
∥
Δ
u
(
t
)
∥
∞
≤
∫
0
t
‍
C
∥
P
1
(
s
)
∥
∞
d
s
≤
∫
0
t
‍
C
t
-
s
∥
P
1
(
s
)
∥
∞
d
s
≤
∫
0
t
‍
C
C
0
t
-
s
×
(
∥
Δ
u
(
s
)
∥
1
,
∞
+
∥
Δ
v
(
s
)
∥
1
,
∞
)
+
∥
ξ
∥
∞
d
s
,
∥
δ
+
Δ
u
(
t
)
∥
∞
≤
∫
0
t
‍
C
t
-
s
∥
P
1
(
s
)
∥
∞
d
s
≤
∫
0
t
‍
C
C
0
t
-
s
×
(
∥
Δ
u
(
s
)
∥
1
,
∞
+
∥
Δ
v
(
s
)
∥
1
,
∞
)
+
∥
ξ
∥
∞
d
s
,
where
(42)
∥
ξ
∥
∞
=
max
⁡
0
≤
t
≤
T
{
∥
ξ
1
(
t
)
∥
∞
+
∥
ξ
2
(
t
)
∥
∞
}
.
The same inequalities are satisfied by Δ
v
(
t
), δ
+
Δ
v
(
t
). We get
(43)
∥
Δ
u
(
t
)
∥
∞
+
∥
δ
+
Δ
u
(
t
)
∥
∞
+
∥
Δ
v
(
t
)
∥
∞
+
∥
δ
+
Δ
v
(
t
)
∥
∞
≤
∫
0
t
‍
4
C
C
0
t
-
s
(
∥
Δ
u
(
s
)
∥
1
,
∞
+
∥
Δ
v
(
s
)
∥
1
,
∞
+
∥
ξ
∥
∞
)
.
It follows from Lemma A.1 that
(44)
∥
Δ
u
(
t
)
∥
1
,
∞
+
∥
Δ
v
(
t
)
∥
1
,
∞
≤
K
~
∥
ξ
∥
∞
t
1
/
2
,
where K
~ depends on C, C
0.


## 5. Numerical Simulations

Numerical Example
1. Set L
=
1, D
1
=
0.2, D
2
=
0.15, D
3
=
0.08 and the initial distributions u, v are almost step functions, regularized near zero, whereas w starts from a droplet:
(45)
u
(
0
,
x
)
=
Θ
(
100
x
)
-
Φ
(
100
x
)
+
1
(
.
01,1
]
(
x
)
,
Θ
(
x
)
=
(
-
1
4
x
3
+
3
4
x
+
9
20
)
1
[
-
1,1
]
(
x
)
,
v
(
0
,
x
)
=
Θ
(
-
100
x
)
-
Φ
(
100
x
)
+
1
[
-
1
,
.
01
)
(
x
)
,
w
(
0
,
x
)
=
2
·
Φ
(
100
x
)
,
Φ
(
x
)
=
(
1
10
x
4
-
1
5
x
2
+
1
10
)
1
[
-
1,1
]
(
x
)
,
for x
∈
[
-
1,1
] where 1 stands for a characteristic function. We use a second-order Runge-Kutta method with time step h
t
=
2
×
1
0
-
6. Here D
1, D
2, D
3 satisfy condition (20) of Theorem 4. Note that (
u
0
,
v
0
,
w
0
)
∈
𝔛
0 with M
=
M
0
=
30. In Figure 1(a) we give experimental values of the solution (
u
,
v
,
w
) of (8) for t
=
0.01 and h
=
.
002 with the initial condition (45) (a). Since h
=
.
002 we divide [
-
1,1
] into 1000 intervals. We put the values of u, v, w at t
=
0. Note that u
0 is increasing on the interval [
-
1,1
]. Furthermore v
0 is decreasing on the interval [
-
1,1
]. The even function w
0 (“marker”) is increasing on the interval [
-
1,0
] and decreasing on the interval [
0,1
]. We consider the differential equations describing the marker position (
d
z
/
d
t
)
(
t
)
=
ν
D
(
t
,
z
(
t
)
), where ν
D is the Darken velocity. We use a second-order Runge-Kutta method. The approximate position of the marker denoted by z (Figure 2(b)) forms a trajectory in (
u
,
v
)-plane shown in Figure 2(a). This object traces the evolution of concentrations of diffunding substances u, v.Initial values (a) and experimental values ofu, v, and w for t
=
0.01, h
=
.
002 (b).
(a)
(b)Evolution of concentrationsu, v (a) and experimental position of marker z
=
z
(
t
) for 0
≤
t
≤
0.01 (b).
(a)
(b)Numerical Example
2.   Set L
=
1, D
1
=
0.2, D
2
=
0.15, D
3
=
0.08 and the initial distribution of w is constant:
(46)
u
(
0
,
x
)
=
Ψ
(
100
x
)
+
5
11
1
B
(
x
)
+
10
11
1
(
.
01,1
]
(
x
)
,
v
(
0
,
x
)
=
-
Ψ
(
100
x
)
+
5
11
1
B
(
x
)
+
10
11
1
[
-
1
,
-
.
01
)
(
x
)
,
w
(
0
,
x
)
=
1
11
,
Ψ
(
x
)
=
(
30
176
x
5
-
50
88
x
3
+
150
176
x
)
1
[
-
1,1
]
(
x
)
,
for x
∈
[
-
1,1
] where 1
B is the characteristic function of B
:
=
[
-
.
01
,
.
01
]. Note that (
u
0
,
v
0
,
w
0
)
∈
𝔛
0 with M
=
M
0
=
30. In Figure 3(a) we give initial values of the solution (
u
,
v
,
w
) of (8) with the initial condition (46) for h
t
=
2
×
1
0
-
6 and h
=
.
002. Starting with a constant distribution w
0
=
1
/
11 (Figure 3(a)) we give experimental values of u, v, w for t
=
0.01 and h
=
.
002 (Figure 3(b)). It is seen that diffusion causes a perturbation for w. This observation is coherent with engineers’ practice.Experimental values ofu, v, and w for t
=
0.01, h
=
.
002.
(a)
(b)Our numerical experiments confirm the theory and stability analysis. It is seen that the proposed numerical method is stable, mass-conservative, and dissipative. Its solutions do not leave the interval[
0,1
]. Since our stability analysis relies on the M-Lipschitz condition for u, v, w, we show (numerically) that the Lipschitz constants for u, v, w decrease in time. We focus on numerical Example 1. Denote these time-dependent constants by M
u
(
t
), M
v
(
t
), M
w
(
t
), respectively, and show them in Figure 4: M
u (a) and M
v (b). Considering the same numerical example we demonstrate convergence in Figure 5 where numerical solutions for different steps h
=
.
008, h
=
.
004, and h
=
.
002 are very close to each other. We take the concentrations u (a) and v (b).Experimental values of Lipschitz constantsM
u
(
t
) and M
v
(
t
) for h
=
.
002, h
=
.
004, and h
=
.
008.
(a)
(b)Experimental values ofu and v for h
=
.
002, h
=
.
004, and h
=
.
008.
(a)
(b)



---

*Source: 517285-2014-04-28.xml*