Open AccessArticle

Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory

Maxence Bigerelle

^1,*

Hakim Naceur

^1,† and

Alain Iost

^2,†

LAMIH UMR CNRS 8201, Université de Valenciennes, Valenciennes 59313, France

Laboratoire MSMP, Arts et Métiers ParisTech, Lille 59046, France

Author to whom correspondence should be addressed.

^†

These authors contributed equally to this work.

Entropy 2016, 18(4), 155; https://doi.org/10.3390/e18040155

Submission received: 26 October 2015 / Revised: 15 March 2016 / Accepted: 12 April 2016 / Published: 21 April 2016

(This article belongs to the Special Issue Entropy Generation in Thermal Systems and Processes 2015)

Download

Browse Figures

Versions Notes

Abstract

In a previous investigation (Bigerelle and Iost, 2004), the authors have proposed a physical interpretation of the instability λ = Δt/Δx² > 1/2 of the parabolic partial differential equations when solved by finite differences. However, our results were obtained using integration techniques based on erf functions meaning that no statistical fluctuation was introduced in the mathematical background. In this paper, we showed that the diffusive system can be divided into sub-systems onto which a Brownian motion is applied. Monte Carlo simulations are carried out to reproduce the macroscopic diffusive system. It is shown that the amount of information characterized by the compression ratio of information of the system is pertinent to quantify the entropy of the system according to some concepts introduced by the authors (Bigerelle and Iost, 2007). Thanks to this mesoscopic discretization, it is proved that information on each sub-cell of the diffusion map decreases with time before the unstable equality λ = 1/2 and increases after this threshold involving an increase in negentropy, i.e., a decrease in entropy contrarily to the second principle of thermodynamics.

Keywords:

instability; entropy; parabolic partial differential equations; Monte Carlo simulations; data compression; information theory

PACS:

89.70.Cf; 02.70.Bf; 05.10.Ln; 05.40.Jc

Graphical Abstract

1. Introduction

Numerous physical phenomena can be modeled by partial differential equations (PDE) [1] and there are many numerical methods for solving PDE [2,3]. The Finite Element Method (FEM) and the Finite Volume Method (FVM) are particularly relevant to solve PDE. The FEM or FVM are mainly useful to solve PDE in complex geometries as well as to optimize the computation times by using different mesh sizes (or different discretization times). Among all these, the Finite Difference Method (FDM) stands out as being universally appropriate and straightforward for solving both linear and nonlinear problems, particularly for uniform mesh sizes and time increments. Stability criteria are analytically easily expressed for uniform mesh sizes and time increments. As the aim of this paper consists in analyzing criteria stability of PDE via Information theory, we deal exclusively with the FDM [4]. To solve these PDE numerically, FDM consists in expanding all variables of the PDE in a Taylor series on a grid at different points of the domain and to limit this expansion to the first derivatives. Then these finite expansions are introduced into the PDE problems and finally the algebraic system of equations is solved using adequate numerical algorithms. The main problems raised in this method are that PDE become discretized and all differential elements do not tend to 0 but to a fixed value that depends on the number of discretized points. However, discretization involves local truncation errors and all rounding errors introduced during the computation due to the finite representation number might grow during the iterative process and leads to numerical instabilities. A high number of mathematical tools and theorems can be used to study the stability of the system, which we shall briefly summarize. However, no physical interpretation in term of Information theory has ever been proposed so far. The fundamental question can be summed up as follows: “Is there a link between the stability threshold and the quantity of information of the diffusive system at a given scale and how to quantify this amount of information?” Brillouin [5] proved that information could be considered as the negative of the system’s entropy called “Negentropy”. Entropy measures the lack of information; it gives the total amount of information missing in the ultramicroscopic structure of a system. In another scientific field, the Information theory can be a powerful tool to analyze data compression [6,7]. The more informative, the lower the ratio of data compression is. The discretized system will be integrally described if the size of the system after compression is minimal, regardless of the width of the data compression of all possible algorithms. Then the size of the program becomes a mathematical measure of negentropy. With regard to the numerous mathematical publications related to PDE stability, we propose to study the stability of parabolic PDE. Moreover, an important category of physical problems can easily be modeled at a microscopic scale and characterized by the Information theory concept via data compression [8].

2. Diffusion Equations

2.1. Parabolic Differential Equations

The Parabolic PDE is given by the following mathematical definition:

\frac{\partial u}{\partial t} - \sum_{i, j}^{n} \frac{\partial}{\partial x_{i}} a_{i, j} \frac{\partial u}{\partial x_{j}} = f (x, t)

(1)

where

(x, t) \in Ω \times

IR+,

Ω

an open set of IRn.

These equations characterize a high number of transport phenomena such as heat transfer, viscous fluid mechanics, transport of atoms-vacancies, Ohm law, contagious epidemics, etc. While remaining quite general, we shall treat in this investigation the simple mono-dimensional Fick Equation [9] that reduces Equation (1) to:

\partial C (x, t) / \partial t = D \partial^{2} C (x, t) / \partial x^{2}

(2)

where

C (x, t)

is the concentration of particles at position

x

after a time

t

and

D

is the diffusion coefficient. If we impose the initial condition,

C (0, 0) = C_{0} δ_{0}

(explain what is

δ_{0}

?) and

Ω = (- \infty, + \infty)

, then the well-known solution is obtained:

C (x, t) = \frac{C_{0}}{2 \sqrt{π D t}} \exp (- \frac{x^{2}}{4 D t})

(3)

This equation is in fact a Gaussian probability density function (if we substitute in Equation (3)

C_{0}

by Equation (1) with zero mean and

σ (t) = \sqrt{2 D t}

standard deviation). Using Taylor’s expansion (explicit in time and central in space), we obtain:

C_{i}^{j + 1} = C_{i}^{j} + D Δ t / {(Δ x)}^{2} [(C_{i - 1}^{j} - 2 C_{i}^{j} + C_{i + 1}^{j})]

(4)

where

C_{i}^{j}

represents the concentration at point

x = i Δ x

and time

t = j Δ t

Usually linear parabolic problems are solved via implicit or Crank–Nicolson’s schemes in order to avoid paying the price of the time step restriction. Unfortunately, implicit schemes (or mixed explicit/implicit schemes) involve that concentrations at time

t + Δ t

(partially or totally) depend on the gradients at time

t + Δ t

. In 1905, Einstein proposed a physical explanation of

D

in Equation (2) by postulating that the concentration at time

t + Δ t

only depends on the concentration at time t. Then the explicit scheme will be retained because it is physically causal (i.e., evolution of the system only depends on previous time). It can be noticed that microscopic simulations of diffusion mechanisms are only explicit methods (Monte Carlo methods, Cellular automata, etc.).

2.2. Stability Study

In this Equation (4) scheme, concentration at time

t + Δ t

is calculated by means of the concentration at time

t

and then the explicit scheme becomes a dynamical system. As a result, it will be possible that this scheme becomes unstable for a set of values

D Δ t / {(Δ x)}^{2}

. The main idea behind stability is that a numerical process should limit the amplification of all components of the initial condition including rounding errors, when applied exactly. The basic analysis consider the growth of perturbations in initial data or the growth of errors introduced at mesh points at a given time level. This mathematical philosophy seems to us to be near the theory of chaos [10] and consequently we think that a duality exists between a mathematical interpretation and its physical representation. There are three common methods to investigate stability: the Von Neumann method, the matrix method, and the energy method. Using one of these three methods, it can be shown that the explicit scheme will be stable if [4]:

λ = D Δ t / {(Δ x)}^{2} \leq 1 / 2

(5)

2.3. Monte Carlo Simulation

The solution of the Cell-PDE could also be modeled using a Brownian motion. The discretization of the diffusion system takes place as follows. A grid of size

r_{x} \times r_{y}

r_{x} = r_{y} = 256

, is constructed. A 2D Monte Carlo (MC) simulation was used rather than a 1D one to visualize the complexity of the discretized diffusion front. In the middle of this grid, each elementary cell is affected to the value 1 on a length

Δ x = r_{x} / 10

(good approximation, result not shown) that represents the entity (we shall call “particle”), which will diffuse; elsewhere, the value 0 is given (no particle). This system constitutes the initial state at time zero (original time). Figure 1 represents this initial state (

t = 0

) where black cells represent the value “1” (particle) and white ones the value “0” (no particle).

Then, for each

p

particle,

n_{s}

jumps are processed on the left or on the right with a probability of 1/2 (more precisely, a jump occurs by choosing at random one particle from all particles). This constitutes a Monte-Carlo Step (MCS). From the microscopic theory of diffusion [9], the diffusion coefficient

D

can be expressed by:

D = a^{2} Γ / 2

(6)

where

Γ

is the particle jump frequency of length

a

from one micro cell to one adjacent micro cell.

The stability condition then can be expressed by including Equation (5) into Equation (6):

\frac{Δ t}{{(Δ x)}^{2}} \leq \frac{1}{a^{2} Γ}

(7)

Then, the number of jumps

n_{s}

over time

Δ t

is given by

n_{s} = Γ Δ t

(8)

Without any lack of generality, we shall consider

D = 1

. The elementary cell size is

a = 1

(pixel unit). Then the stability criterion becomes:

\sqrt{n_{s}} \leq Δ x = r_{x} / 10

(9)

Figure 1 represents the snapshots of the evolution of the diffusion system at different times,

n_{s} = 0, 10, 300, 655, 1000

MCS. The particular value

n_{s}^{s c} = Δ x^{2} = {(255 / 10)}^{2} = 655

corresponds to the stability–instability threshold. Firstly, we shall verify the convergence rate of the Monte Carlo simulation to the Gaussian solution. Figure 1 shows also the concentration curves for different evolution times. The theoretical Gaussian curves given by Equation (3) are plotted for each concentration map. As can be observed, convergence to the Gaussian approximation is well suited if

n_{s} > 300

. Similar to these analyses, we shall develop and induce the Gaussian Probability Density Function approximation for the case

n_{s} = n_{s}^{s c}

; it becomes clear from the Monte-Carlo Simulation that the Gaussian PDF is an adequate assumption to model the diffusion process.

3. Quantitative Description of the System in the Area of the Information Theory

3.1. Mathematical Formalism

The system X of diffusion states is an ordered set whose dimension

\dim (X) = r_{x} r_{y} = 256^{2}

. This system contains sub-systems of the same dimension represented by

x_{i}

with length equal to

Δ x = r_{x} / 10

and then

\dim (x_{i}) = r_{x} r_{y} / 10 = 256^{2} / 10

Let us note

T

, the nonlinear bijective algebra that transforms the system

x

into a system

y

with

y = T (x)

. Let

{A}

be an algebraic set,

T_{\min}

is said

{A} - maxi

contractile in

X

if there exists one algebra noted

T_{\min} \in {A}

such that:

\forall x \subset X, \forall T \in {A}, \exists T_{\min} \in {A}, \dim (T_{\min} (x)) \leq \dim (T (y))

(10)

T_{o p t}

is said

{\infty} - maxi

contractile in

X

T_{o p t}

{A} - maxi

contractile in

X

where

{A}

denoted all the sets of all possible algebras defined by arithmetics. Regrettably, according to the Gödel incompleteness theorem, the proposition

≪ T_{o p t}

is said

{\infty} - maxi

contractile in

(E_{S}, Ψ, G) ≫

is undecidable [11] where

(E_{S}, Ψ, G)

are, respectively, a vector space, the sets of all possible subspaces of E_S, G a relation from E_S to E_S. However, a set of

{A}

algebra can be built and then

≪ {A} - maxi

contractile in

X ≫

becomes decidable.

3.2. The Different Classes of Algebra

Run Length Encoding (RLE) is often used for data compression [12]. Each cell of the Monte-Carlo system is encoded by one bit, representing the existence or the absence of a particle. The efficiency of the RLE algorithm comes from the fact that selecting a cell at random, the probability that cells in the vicinity possess the same state is high under the hypothesis that physical system is not a pure random one (a pure random system possess the maximal statistical entropy). The algorithm counts the cells in rows on the grid looking for runs of cells having the same state. Encoding the number of runs rather than all individual states significantly reduces the initial size. Then we applied the Huffman algebra based on a list of the alphabetical symbols in decreasing probability order [13] on this system. This algebra allows us to better describe the structure of the diffusion that can be related to the probabilistic approach to this algebra (some statistical tools are closed to those used in statistical thermodynamics). This composition of both algebra denoted

H U F \circ R L E (X)

is well adapted to describe the power laws met in the diffusion process [8,14].

4. Results

4.1. Compression Data Analyses

In the Monte Carlo simulation, the diffusive system

X (t)

is composed of 10 sub-systems

x_{i} (t)

with

X (t) = \cup_{i = 1}^{10} x_{i} (t)

. At initial time

t = 0

, the cells of the sub-system

x_{5} (t)

are all equal to 1 and all the other sub-system cells are equal to 0. We shall analyze by the

H U F \circ R L E

algebra the system

X

and the adjacent cells of

x_{5}

, i.e.,

x_{4}

(or

x_{6})

that are the cells that “received” the diffusion particles from the source located in

x_{5}

. To integrate the stochastic aspect of the diffusion, 1000 simulations are carried out. At

t = 0

, we get

\dim (H U F \circ R L E (X)) = 711

and

\dim (H U F \circ R L E (x_{4})) = 488

and these probability density functions are Dirac distributions (Figure 2). As it can be shown,

H U F \circ R L E

algebra is contractile on both

X

and

x_{4}

systems because:

\dim (H U F \circ R L E (X (0))) < \dim (X (0)) = 65536 and \dim (H U F \circ R L E (x_{4} (0))) < \dim (x_{4} (0)) = 6554

Then, during the diffusion process (

t > 0

\dim (H U F \circ R L E (X (t)))

and

\dim (H U F \circ R L E (x_{4} (t)))

follow Gaussian densities meaning that the stochastic variation of the Monte Carlo is detected by the RLE and Huffman algebra. If a diffusion system is finite, the concentration

C (x, t)

will become a random variable that will follow a Gaussian law under a pure Brownian motion assumption [15]. The measure of the reduced space dimension follows the same density probability function meaning that an isomorphism exists between space dimension and the stochastic aspect of diffusion. As the diffusion time increases,

\dim (H U F \circ R L E (X (t)))

increases (Figure 3): in fact, during the diffusion process, the total entropy increases, which will decrease negentropy. The system tends to disorder and information can be less and less reduced.

4.2. The Relation between Stability Criteria and the Dimension of Reduced Space

Now we shall analyze the evolution of the program size for the sub-space

x_{4}

. Figure 3 represents the evolution of the program size versus diffusion time.

First the dimension increases to reach a maximum and then decreases after this threshold. The diffusion time corresponding to this maximum is equal to

t_{H U F \circ R L E}^{\max} = 650_{\pm 30}

. This leads to an important remark. Indeed, statistically:

t_{H U F \circ R L E}^{\max} = n_{s}^{s c} = Δ x^{2} = {(255 / 10)}^{2} = 655

(11)

In the mono-dimensional diffusion problem, a PDE case gets two adjacent PDE cells that will diffuse on this cell. We shall consider particles coming from the adjacent cell. At the original time, no particle from the adjacent cell is present on the PDE cell and then the description of the cells can be summarized by an algorithmic formalism to “n rep 0” and then

\dim (H U F \circ R L E (x_{4} (0)))

is minimal. Then these particles will diffuse on the cell with respect to time. Then entropy will increase thanks to the diffusion and the amount of information will increase therefore increasing

\dim (H U F \circ R L E (x_{4} (0)))

until it reaches

t_{H U F \circ R L E}^{\max}

that is the stability criterion (Figure 3). When the diffusion time is higher than the stability threshold, then

\dim (H U F \circ R L E (x_{4} (t)))

decreases. This means that there exists a link between the PDE stability criterion and the amount of information characterized by this original tool. We shall now give an explanation of the decrease in the dimension after the threshold criterion.

Now let us analyze very precisely this critical density function (Figure 4).

Let us now consider the

x_{2}

intervals. Concentration is the summation of all the Gaussian curves on these intervals and the concentration on the

x_{2}

cell only results from the flux from cells

x_{1}

and

x_{3}

. For the condition

Δ x = \sqrt{2 D t}

, the inflexion points of the Gaussians

G_{1}

and

G_{3}

are in the middle of the interval

x_{2}

. Then the following expressions can be stated:

\partial^{2} C_{1, 0} (x, Δ x^{2} / 2 D, Δ x) / \partial x^{2} < 0, \partial^{2} C_{3, 0} (x, Δ x^{2} / 2 D, Δ x) / \partial x^{2} > 0 \forall x \in [x_{2} - Δ x / 2, x_{2}]

(12)

\partial^{2} C_{1, 0} (x_{2}, Δ x^{2} / 2 D, Δ x) / \partial x^{2} = 0, \partial^{2} C_{3, 0} (x_{2}, Δ x^{2} / 2 D, Δ x) / \partial x^{2} = 0

(13)

\partial^{2} C_{1, 0} (x, Δ x^{2} / 2 D, Δ x) / \partial x^{2} > 0, \partial^{2} C_{3, 0} (x, Δ x^{2} / 2 D, Δ x) / \partial x^{2} < 0 \forall x \in [x_{2}, x_{2} + Δ x / 2]

(14)

From Equation (2),

\partial^{2} C_{i, j} (x, t, Δ x) / \partial x^{2} \propto \partial C_{i, j} (x, t, Δ x) / \partial t

and it can be concluded that 50% of

d x

micro intervals (

d x ≪ Δ x

) of

x_{2}

intervals get a positive temporal gradient (

\partial C_{1, 0} (x, t, Δ x) / \partial t > 0

and

\partial C_{3, 0} (x, t, Δ x) / \partial t > 0

), and 50% a negative gradient (

\partial C_{1, 0} (x, t, Δ x) / \partial t < 0

and

\partial C_{3, 0} (x, t, Δ x) / \partial t < 0

). After the stability threshold, more than 50% of micro cells get a negative gradient. This clearly means that particles issued from adjacent cells

x_{1}

and

x_{3}

will decrease in keeping with time, which will decrease the configuration entropy of the system

x_{2}

contrarily to the second principle of thermodynamics for a diffusive system. The mathematical criterion of stability of the explicit scheme is then physically related to the production of entropy characterized by our algebra operator: if the scheme is unstable, then the production of entropy on each discretized cell is negative and the size of the program will decrease.

We eventually get:

\partial \dim (H U F \circ R L E (x_{4} (t))) / \partial t > 0 if \sqrt{2 D t} < Δ x, t > 0

(15)

\partial \dim (H U F \circ R L E (x_{4} (t))) / \partial t = 0 if \sqrt{2 D t} = Δ x

(16)

\partial \dim (H U F \circ R L E (x_{4} (t))) / \partial t < 0 if \sqrt{2 D t} > Δ x

(17)

To illustrate these partial differential equations, different system sizes are simulated using the Monte Carlo diffusion process (from s = 250 to 1000). Then, according to Equation (11), and if assertion

≪ {H U F \circ R L E} - maxi

contractile in

X ≫

is true, the relation between the threshold that depends on the system size (because the system is always decomposed into 10 sub-systems with parts of equal length) becomes:

t_{H U F \circ R L E}^{\max} (s) = n_{s}^{s c} (s) = Δ x^{2} (s) = {(s / 10)}^{2}

(18)

Figure 5 represents the normalized dimension obtained by

\dim (H U F \circ R L E (x_{4} (t))) / \dim (x_{4} (t))

for the different system sizes. The maximal values for each system size are computed and compared to Equation (18).

As can be observed, data match to the theoretical equation meaning that whatever the system size, the time of maximal entropy

t_{H U F \circ R L E}^{\max} (s)

is obtained so that it corresponds to the stability criterion. Thanks to statistical thermodynamics, it was shown that stability threshold corresponding to a violation of the second principle of the thermodynamics that confirms all results shown by the Information theory tools [16].

5. Two Examples of PDE

5.1. Non Constant Diffusion Coefficient

The problems treated by the information theory were investigated with a constant diffusion coefficient. In this case, the diffusion coefficient does not depend on the position. Thus, we will treat the following problem:

\frac{\partial}{\partial t} C (x, t) = \frac{\partial}{\partial x} D (x) \frac{\partial}{\partial x} C (x, t)

(19)

For a real positive diffusion parameter D(x) that depends on the x-position, the classical explicit first order in time and second order in space method in the one-dimensional case is stable if:

λ = \frac{\max (D (x)) Δ t}{{(Δ x)}^{2}} \leq 1 / 2

(20)

The x-dependence of diffusion coefficient D can be due to different variations of material properties like crystal structure, local molar concentration, vacancy gradient, etc., or external conditions like residual stress, temperature, electrical field, chemical potential, etc.

A positive perturbation

Δ a (x)

is introduced to the unitary jump length during the Brownian simulation process and the mean displacement of each particle becomes

1 + Δ a (x)

. By posing

a (x) = 1 + Δ a (x)

in Equation (7), one gets

D (x) = {(a (x))}^{2}

and the local stability criterion becomes:

\sqrt{n_{s} (x)} \leq (Δ x / (1 + Δ a (x)))

(21)

Finally, the global stability criterion becomes:

\sqrt{n_{s}} \leq \min_{x} (Δ x / (1 + Δ a (x)))

(22)

A Gaussian function will be used to simulate the variation of diffusion coefficient:

Δ a (x, c, σ) = e^{\frac{1}{2} {(\frac{x - c}{σ})}^{2}}

(23)

To assess nonlinearity, the size of the system is increased and resolution map will be extended to N = 1024 (rather than N = 256 for the previous simulations). The Gaussian is centered at 3/4 of the diffusion map and then c = 3 × 1024/4 = 768 (Figure 6).

The mesh size is equal to

Δ x = 1024 / 8 = 128

. At t = 0, the diffusion cell is located the middle of the map (at t = 0, x_i,j = 1

\forall j \in [1, 1024]

if i belongs to i₄ = [448,576] otherwise x_i,j = 0). The diffusion only occurs in the x direction. Figure 7 represents snapshots of the diffusion process for two standard deviations. As can be observed, for the standard deviation of 11² pixels, a dissymmetry of the concentration appears and diffusion increases on the right of the diffusion map. Contrarily, diffusion map is symmetrical for diffusion map obtained with small standard deviation (and for high standard deviation, i.e., greater than 2²⁰ pixels, result not shown). According to Equation (21), for very small standard deviations, one gets:

\lim_{σ \to 0} (n_{s} (x, c, σ)) \leq Δ x^{2} = 128^{2} = 16384, \forall x, \forall c,

(24)

and for high standard deviations, one gets:

\lim_{σ \to \infty} (n_{s} (x, c, σ)) \leq Δ x^{2} / 4 = 4096, \forall x, \forall c

(25)

In these two cases, diffusion becomes linear and stability condition does not depend on x position.

Let us now analyze the two compression ratio,

\dim (H U F \circ R L E (i_{n} (t))) / \dim (i_{n} (t))

, of the two adjacent cells (i₃, i₅) at i₄ (i₃ = [320,448], i₅ = [576,704]) for different values of standard deviations (Figure 8). As can be observed, all curves present maximal values of compression leading to a MCS time

n_{s}^{s c} (x, c, σ)

that depend on the values of the standard deviations.

From our postulate, the theoretical stability criterion is equal to

n_{s}^{s c} (x, c, σ) = {(Δ x / (1 + Δ a (x, c, σ)))}^{2}

. Values of the theoretical stability threshold agree with simulated values (Figure 9), meaning that our method can be applied to linear PDE with non-constant diffusion coefficient.

5.2. Nonlinear PDE

The problems treated by the information theory were investigated on a linear PDE. A vast array of complex phenomena of motion, reaction, diffusion, equilibrium, conservation, and more lead to nonlinear PDE. A PDE is said to be nonlinear if the relations between the unknown functions and their partial derivatives involved in the equation are nonlinear. Thus, we will treat the following nonlinear problem:

\frac{\partial}{\partial t} C (x, t) = \frac{\partial}{\partial x} D (C (x, t)) \frac{\partial}{\partial x} C (x, t)

(26)

In our case, the diffusion coefficient depends on the concentration

C (x, t)

. Lee [17] proposed diffusion coefficients relating to the uptake of excess calcium by calcium chloride [18]. The diffusion coefficient D takes the form:

D (C (x, t)) = \frac{D_{0}}{1 + a C (x, t)}

(27)

where a is a coefficient that represents sur-diffusion processes (a < 0) or a sub-diffusion ones (a > 0) or a classical diffusion ones (a = 0). In our simulations,

a = {- 1, 0, 1, 2, 3, 4, 5}

. The stability criterion is given by:

λ = \frac{\max (D (C (x, t))) Δ t}{{(Δ x)}^{2}} \leq 1 / 2

(28)

The Monte Carlo simulation is based on an explicit scheme. Jump length at time t + dt is computed from concentration at time t (Figure 10). Too validate our hypothesis, 50 simulations are performed with three sizes of the system: 128², 256² and 512² pixels. Maximal values of the dimension are extracted (Figure 11).

Values of the theoretical stability threshold agree with simulated values (Figure 12), meaning that our method can be applied to nonlinear PDE.

6. Discussion

To summarize, the mathematical criterion of stability

λ = Δ t / Δ x^{2} > 1 / 2

of the explicit scheme leads to an unstable scheme due to the decrease of information versus time characterized by the reduction size of the system on each discretized cell of the whole system. In the stability vision of Von Neumann, a finite difference scheme is stable if the errors made at one time step of the calculation do not let the errors increase for the following time steps. Therefore, in our approach, this means that, for a fixed time of observation of the diffusion, an increase of the size of the sub-cells of the grid will irrevocably lead to a critical size from which discretized information will not be enough to contain information on concentration gradient described by the partial differential equation. In fact, the physical interpretation of instability of the PDE was never introduced in a high number of papers, which treat mathematical aspects of instability via information formalisms. Mishra has introduced the problems of violation of local thermodynamic laws in the transport equations via Information theory. He discretizes the transport equation via a simple centered finite difference scheme. The time derivative is replaced with a forward difference and the spatial derivative with a central difference. In this scheme of discretization, the central scheme leads to a growth of energy at every time step and is unstable whatever the size of both discretization of time and space [19]. Mishra explains why the solutions computed with the central scheme blow up. After all, the central scheme seems a reasonable approximation of the transport equation. The physical explanation can be deduced from the following argument: the exact solution moves to the right with entities speed (Figure 13).

Unfortunately, for the scheme described in this publication that we have physically proof to be physical relevant in term of the physics causality via Einstein theory of the Brownian motion [15], no justification of the conditional stabilities was found in the bibliography via Information theory. However, some results using entropy concept in PDE stability conclude [21] that many linear Partial Differential Equations or Integral equations with non-constant coefficients satisfy some entropy dissipation properties [22]. Finally, explicit time discretization leads to entropy production [23] and more specifically, explicit scheme with Lax–Friedrichs flux is entropy stable under conditions given by Equation (5) [24].

7. Conclusions

In this paper, we showed that a relation exists between the well-known stability criterion

λ = D Δ t / {(Δ x)}^{2} \leq 1 / 2

applied on the explicit numerical scheme

C_{x}^{t + Δ t} = C_{i}^{t} + λ (C_{x - Δ x}^{t} - 2 C_{x}^{t} + C_{x + Δ x}^{t})

to guarantee the convergence of the 1D non-stationary PDE

\partial C (x, t) / \partial t = D \partial^{2} C (x, t) / \partial x^{2}

to the solution and the quantity of information contained at a mesoscopic scale under the size of discretization

Δ x

. Implicit and semi implicit schemes (such as the unconditionally stable scheme of Crank–Nicolson method) lead to a violation of the principle of causality: the cause must precede its effect according to all inertial observers meaning that the cause and its effect are separated by a time-like interval, and the effect belongs to the future of its cause. Jaroszkiewicz pointed out the violation of causality in the fields of discrete mechanics by analyzing the implicit and explicit Euler Schemes and deduced that the implicit scheme involves that the flow of information occurs in a temporal loop. Jaroszkiewicz concluded his argument with the following question: “we should ask, given this violation of causal implication, how could we ever use implicit Euler scheme in practical calculation” [25]. Under Causality scheme assumptions, it was proved that instability occurs during iterative process if and only if information contained inside the cell

Δ x

decreases with time, thus violating the second principle of the thermodynamics. As microscopic simulations of diffusion mechanisms are only explicit methods (Monte Carlo methods, Cellular automata, etc.), the criterion of sub-cell information seems to be a relevant tool to assess that a stable explicit scheme meets both the causality and the second principle of the thermodynamics.

Author Contributions

Maxence Bigerelle found the concept of negentropy production in explicit numerical scheme. Alain Iost applied this to material science. Hakim Naceur formulated the explicit scheme. All authors have read and approved the final manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

Dautray, R.; Lions, J.-L. Analyse Mathématique et Calcul Numérique pour les Sciences et les Techniques. Volume 1: Modèles Physiques; Masson: Paris, France, 1987. (In French) [Google Scholar]
Dautray, R.; Lions, J.-L. Analyse Mathématique et Calcul Numérique pour les Sciences et les Techniques. Volume 7: Evolution, Fourier, Laplace; Masson: Paris, France, 1988. (In French) [Google Scholar]
Dautray, R.; Lions, J.-L. Analyse Mathématique et Calcul Numérique pour les Sciences et les Techniques. Tome 9: Evolution, Numérique, Transport; Masson: Paris, France, 1997. (In French) [Google Scholar]
Forsythe, G.E.; Wasow, W.R. Finite Difference Methods for Partial Differential Equations; Wiley: New York, NY, USA, 1960. [Google Scholar]
Brillouin, L. Science and Information Theory; Academic Press: New York, NY, USA, 1956. [Google Scholar]
Delahaye, J.P. Information, Complexité Hasard; Hermès: Paris, France, 1994. (In French) [Google Scholar]
Chaitin, G.J. Algorithmic Information Theory; Cambridge University Press: Cambridge, UK, 1887. [Google Scholar]
Bigerelle, M.; Iost, A. Relations entre l’entropie physique, le codage de l’information et l’énergie de simulation. Can. J. Phys. 2007, 85, 1381–1394. (In French) [Google Scholar] [CrossRef]
Adda, Y.; Philibert, J. La Diffusion Dans les Solides, Tome II; Presses Universitaires de France: Paris, France, 1966. [Google Scholar]
Peitgen, H.O.; Jürgens, H.; Saupe, D. Chaos and Fractals: New Frontiers of Science; Springer: New York, NY, USA, 1992. [Google Scholar]
Feferman, J.W. Kurt Gödel Collected Works. Volume 1: Publications, 1929–1936; Oxford University Press: Oxford, UK, 1986. [Google Scholar]
Salomon, D. Data Compression; Springer: New York, NY, USA, 1998. [Google Scholar]
Huffman, D. A method for the construction of minimum redundancy codes. Proc. IRE 1952, 40, 1098–1101. [Google Scholar] [CrossRef]
Bigerelle, M.; Iost, A. Relationship between information reduction and the equilibrium state description. Comput. Mater. Sci. 2002, 24, 133–138. [Google Scholar] [CrossRef]
Einstein, A. Investigation on the Theory of Brownian Motion; Dover Publications: New York, NY, USA, 1956. [Google Scholar]
Bigerelle, M.; Iost, A. Physical interpretation of the numerical instabilities in diffusion equations via statistical thermodynamics. Int. J. Nonlinear Sci. Numer. Simul. 2004, 5, 121–134. [Google Scholar] [CrossRef]
Lee, C.F. On the solution of some diffusion equations with concentration-dependent diffusion coefficients—II. IMA J. Appl. Math. 1972, 10, 129–133. [Google Scholar] [CrossRef]
Wagner, C. Diffusion processes during the uptake of excess calcium by calcium fluoride. J. Phys. Chem. Solids 1968, 29, 1925–1930. [Google Scholar] [CrossRef]
Mishra, S.; Veerappa Gowda, G.D. Optimal entropy solutions for conservation laws with discontinuous flux-functions. J. Hyperbolic Differ. Equ. 2005, 2, 783–837. [Google Scholar]
Mishra, S. Numerical Methods for Conservation Laws and Related Equations. Available online: https://www2.math.ethz.ch/education/bachelor/lectures/fs2013/math/nhdgl/numcl_notes_HOMEPAGE.pdf (accessed on 14 April 2016).
Tadmor, E. Entropy stability theory for difference approximations of nonlinear conservation laws and related time-dependent problems. Acta Numer. 2003, 12, 451–512. [Google Scholar] [CrossRef]
Carrillo, J.A.; Jüngel, A.; Markowich, P.A.; Toscani, G.; Unterreiter, A. Entropy dissipation methods for degenerate parabolic problems and generalized sobolev inequalities. Monatshefte für Mathematik 2001, 133, 1–82. [Google Scholar] [CrossRef]
Lefloch, P.G.; Mercier, J.M.; Rohde, C. Fully discrete, entropy conservative schemes of arbitrary order. SIAM J. Numer. Anal. 2002, 40, 1968–1992. [Google Scholar] [CrossRef]
Noble, P.; Vila, J.P. Stability theory for difference approximations of some dispersive shallow water equations and application to thin film flows. 2013; arXiv:1304.3805. [Google Scholar]
Jaroszkiewicz, G. Principles of Discrete Time Mechanics; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]

Figure 1. Snapshot of the Monte Carlo simulation for a grid size of 255² at different Monte Carlo Steps (MCS) times (0, 10, 300, 655 (stability criterion), and 1000) and the concentration profiles under Gaussian approximation and given by the Monte Carlo method.

Figure 2. Histograms of the program sizes for the system

X

(a–e) and the sub system

x_{4}

(f–j) corresponding to the Monte Carlo simulation for a grid size of 255² at different MCS times (0, 10, 300, 655 (stability criterion), and 1000). (a) Mean and standard deviation are computed and Gaussian density is plotted.

Figure 2. Histograms of the program sizes for the system

X

(a–e) and the sub system

x_{4}

Figure 3.

\dim (H U F \circ R L E (X (t)))

and

\dim (H U F \circ R L E (x_{4} (t)))

versus the diffusion time (in MCS) for the grid size of 255² cells (see Figure 1). Points are means of 40 simulations.

Figure 3.

\dim (H U F \circ R L E (X (t)))

and

\dim (H U F \circ R L E (x_{4} (t)))

versus the diffusion time (in MCS) for the grid size of 255² cells (see Figure 1). Points are means of 40 simulations.

Figure 4. Three adjacent cells

x_{1}, x_{2}, x_{3}

of the discretized space grid

Δ x = 1

on which mesoscopic simulations are processed up to the instability criterion. Gaussian curves G₁, G₂ and G₃ with standard deviation of

Δ x = \sqrt{2 D t} = 1

(

D = 1

t = 1 / 2

, i.e.,

λ = 1 / 2

) are centered on each cell.

Figure 4. Three adjacent cells

x_{1}, x_{2}, x_{3}

of the discretized space grid

Δ x = 1

on which mesoscopic simulations are processed up to the instability criterion. Gaussian curves G₁, G₂ and G₃ with standard deviation of

Δ x = \sqrt{2 D t} = 1

(

D = 1

t = 1 / 2

, i.e.,

λ = 1 / 2

) are centered on each cell.

Figure 5.

\dim (H U F \circ R L E (x_{4} (t))) / \dim (x_{4} (t))

versus the diffusion time (in MCS) for different grid sizes of cells (250, 300, …, 950, 1000). Values of the time

t_{H U F \circ R L E}^{\max}

are plotted (graph on the right corner) versus the theoretical stability criterion

n_{s}^{s c} = {(s / 10)}^{2}

given by Equation (18).

Figure 5.

\dim (H U F \circ R L E (x_{4} (t))) / \dim (x_{4} (t))

versus the diffusion time (in MCS) for different grid sizes of cells (250, 300, …, 950, 1000). Values of the time

t_{H U F \circ R L E}^{\max}

are plotted (graph on the right corner) versus the theoretical stability criterion

n_{s}^{s c} = {(s / 10)}^{2}

given by Equation (18).

Figure 6. Values of

Δ a (x, c, σ) = \exp 0.5 {(x - c)}^{2} / σ^{2}

versus x for different standard deviations σ. These functions allow producing nonlinear PDEs involving different stability criteria

\sqrt{n_{s} (x, c, σ)} \leq Δ x / (1 + Δ a (x, c, σ))

along the x position for a given pair

(c, σ)

Figure 6. Values of

Δ a (x, c, σ) = \exp 0.5 {(x - c)}^{2} / σ^{2}

versus x for different standard deviations σ. These functions allow producing nonlinear PDEs involving different stability criteria

\sqrt{n_{s} (x, c, σ)} \leq Δ x / (1 + Δ a (x, c, σ))

along the x position for a given pair

(c, σ)

Figure 7. Snapshots of the Monte Carlo simulation for a grid size of 1024² pixels at different Monte Carlo Steps (MCS) times (50, stability criterion, 28,377). (a)

σ = 0

, t = 50; (b)

σ = 0

, t = 16,384; (c)

σ = 0

, t = 28,377; (d)

σ = 11^{2}

, t = 50; (e)

σ = 11^{2}

, t = 9938; (f)

σ = 11^{2}

, t = 28,377.

Figure 7. Snapshots of the Monte Carlo simulation for a grid size of 1024² pixels at different Monte Carlo Steps (MCS) times (50, stability criterion, 28,377). (a)

σ = 0

, t = 50; (b)

σ = 0

, t = 16,384; (c)

σ = 0

, t = 28,377; (d)

σ = 11^{2}

, t = 50; (e)

σ = 11^{2}

, t = 9938; (f)

σ = 11^{2}

, t = 28,377.

Figure 8.

\dim (H U F \circ R L E (i_{5} (t))) / \dim (i_{5} (t))

versus the diffusion time (in MCS) for different values of standard deviations.

Figure 8.

\dim (H U F \circ R L E (i_{5} (t))) / \dim (i_{5} (t))

versus the diffusion time (in MCS) for different values of standard deviations.

Figure 9. Plot of MCS

t_{H U F \circ R L E}^{\max}

versus the theoretical stability criterion

n_{s}^{s c} (x) = {(Δ x / (1 + Δ a (x)))}^{2}

given by Equation (21) for both cells i₄ and i₅ corresponding to maximal values of Figure 8.

Figure 9. Plot of MCS

t_{H U F \circ R L E}^{\max}

versus the theoretical stability criterion

n_{s}^{s c} (x) = {(Δ x / (1 + Δ a (x)))}^{2}

given by Equation (21) for both cells i₄ and i₅ corresponding to maximal values of Figure 8.

Figure 10. Snapshots of the nonlinear Monte Carlo simulation for a grid size of 129² pixels at five Monte Carlo Steps (MCS) times with different values of a given by

D (C (x, t)) = D_{0} / (1 + a C (x, t))

and snapshots corresponding to stability criterion and their associated Monte Carlo Steps.

Figure 10. Snapshots of the nonlinear Monte Carlo simulation for a grid size of 129² pixels at five Monte Carlo Steps (MCS) times with different values of a given by

D (C (x, t)) = D_{0} / (1 + a C (x, t))

and snapshots corresponding to stability criterion and their associated Monte Carlo Steps.

Figure 11.

\dim (H U F \circ R L E (i_{5} (t))) / \dim (i_{5} (t))

versus the diffusion time (in MCS) for different values of diffusion coefficient

D (C (x, t)) = D_{0} / (1 + a C (x, t))

Figure 11.

\dim (H U F \circ R L E (i_{5} (t))) / \dim (i_{5} (t))

versus the diffusion time (in MCS) for different values of diffusion coefficient

D (C (x, t)) = D_{0} / (1 + a C (x, t))

Figure 12. Plot of MCS

t_{H U F \circ R L E}^{\max}

versus the theoretical stability criterion for both cells i₄ and i₅ corresponding to maximal values of Figure 11 (Size = 257²) and Size = (129², 512²).

Figure 12. Plot of MCS

t_{H U F \circ R L E}^{\max}

versus the theoretical stability criterion for both cells i₄ and i₅ corresponding to maximal values of Figure 11 (Size = 257²) and Size = (129², 512²).

Figure 13. The central scheme used by Mishra [20] to illustrate physical interpretation of the unconditional instability of the simple centered finite difference scheme. Continuous arrows (green) indicate numerical propagation and dashed ones arrows (magenta) physical propagation.

© 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Bigerelle, M.; Naceur, H.; Iost, A. Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory. Entropy 2016, 18, 155. https://doi.org/10.3390/e18040155

AMA Style

Bigerelle M, Naceur H, Iost A. Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory. Entropy. 2016; 18(4):155. https://doi.org/10.3390/e18040155

Chicago/Turabian Style

Bigerelle, Maxence, Hakim Naceur, and Alain Iost. 2016. "Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory" Entropy 18, no. 4: 155. https://doi.org/10.3390/e18040155

APA Style

Bigerelle, M., Naceur, H., & Iost, A. (2016). Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory. Entropy, 18(4), 155. https://doi.org/10.3390/e18040155

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Analyses of the Instabilities in the Discretized Diffusion Equations via Information Theory

Abstract

1. Introduction

2. Diffusion Equations

2.1. Parabolic Differential Equations

2.2. Stability Study

2.3. Monte Carlo Simulation

3. Quantitative Description of the System in the Area of the Information Theory

3.1. Mathematical Formalism

3.2. The Different Classes of Algebra

4. Results

4.1. Compression Data Analyses

4.2. The Relation between Stability Criteria and the Dimension of Reduced Space

5. Two Examples of PDE

5.1. Non Constant Diffusion Coefficient

5.2. Nonlinear PDE

6. Discussion

7. Conclusions

Author Contributions

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI