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

Abstract: In a previous investigation (Bigerelle and Iost, 2004), the authors have proposed a physical interpretation of the instability λ = ∆t/∆x2 > 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.


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].

Parabolic Differential Equations
The Parabolic PDE is given by the following mathematical definition: where px, tq P Ωˆ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: BC px, tq {Bt " DB 2 C px, tq {Bx 2  (2) where C px, tq 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 p0, 0q " C 0 δ 0 (explain what is δ 0 ?) and Ω " p´8, `8q, then the well-known solution is obtained: 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 σ ptq " ?2Dt standard deviation).Using Taylor's expansion (explicit in time and central in space), we obtain: where C j i 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 Entropy 2016, 18, 155 3 of 16 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.).

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{ p∆xq 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]: (5)

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 ˆry , 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: 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): Then, the number of jumps n s over time ∆ t is given by Without any lack of generality, we shall consider D " 1.The elementary cell size is a " 1 (pixel unit).Then the stability criterion becomes: ?   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 sc s " ∆x 2 " p255{10q 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 sc s ; it becomes clear from the Monte-Carlo Simulation that the Gaussian PDF is an adequate assumption to model the diffusion process.

Mathematical Formalism
The system X of diffusion states is an ordered set whose dimension dimpXq " 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 dimpx i q " 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 " Tpxq.Let tAu be an algebraic set, T min is said tAu -maxi contractile in X if there exists one algebra noted T min P tAu such that: @x Ă X, @T P tAu , DT min P tAu , dimpT min pxqq ď dimpTpyqq (10) T opt is said t8u -maxi contractile in X if T opt is tAu -maxi contractile in X where tAu denoted all the sets of all possible algebras defined by arithmetics.Regrettably, according to the Gödel incompleteness theorem, the proposition !T opt is said t8u -maxi contractile in pE S , Ψ, Gq " is undecidable [11] where pE S , Ψ, Gq 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 tAu algebra can be built and then !tAu -maxi contractile in X " becomes decidable.

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 HUF ˝RLEpXq is well adapted to describe the power laws met in the diffusion process [8,14].

Compression Data Analyses
In the Monte Carlo simulation, the diffusive system Xptq is composed of 10 sub-systems x i ptq with Xptq " x i ptq.At initial time t " 0, the cells of the sub-system x 5 ptq are all equal to 1 and all the other sub-system cells are equal to 0. We shall analyze by the HUF ˝RLE algebra the system X and the adjacent cells of x 5 , i.e., x 4 (or x 6 q 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 dimpHUF ˝RLEpXqq " 711 and dimpHUF ˝RLEpx 4 qq " 488 and these probability density functions are Dirac distributions (Figure 2).As it can be shown, HUF ˝RLE algebra is contractile on both X and x 4 systems because: dimpHUF ˝RLEpXp0qqq ă dimpXp0qq " 65536 and dimpHUF ˝RLEpx 4 p0qqq ă dimpx 4 p0qq " 6554 Then, during the diffusion process (t ą 0), dimpHUF ˝RLEpXptqqq and dimpHUF ˝RLEpx 4 ptqqq 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 px, tq 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, dimpHUF ˝RLEpXptqqq 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.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( ( ( ))) HUF RLE 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.for the grid size of 255 2 cells (see Figure 1).Points are means of 40 simulations.

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 4 x .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 max 30 650 . This leads to an important remark.Indeed, statistically: (255 / 10) 655 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 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 t  that is the stability criterion (Figure 3).When the diffusion time is higher than the stability threshold, then 4 dim( ( ( ))) HUF RLE x 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)., , x x x of the discretized space grid 1 x   on which mesoscopic simulations are processed up to the instability criterion.Gaussian curves G1, G2 and G3 with standard deviation of ) are centered on each cell.

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 max HUF˝RLE " 650 ˘30 .This leads to an important remark.Indeed, statistically: 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 dimpHUF ˝RLEpx 4 p0qqq 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 dimpHUF ˝RLEpx 4 p0qqq until it reaches t max HUF˝RLE that is the stability criterion (Figure 3).When the diffusion time is higher than the stability threshold, then dimpHUF ˝RLEpx 4 ptqqq 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).for the grid size of 255 2 cells (see Figure 1).Points are means of 40 simulations.

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 4 x .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 max 30 650 . This leads to an important remark.Indeed, statistically: (255 / 10) 655 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 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 t  that is the stability criterion (Figure 3).When the diffusion time is higher than the stability threshold, then 4 dim( ( ( ))) HUF RLE x 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)., , x x x of the discretized space grid 1 x   on which mesoscopic simulations are processed up to the instability criterion.Gaussian curves G1, G2 and G3 with standard deviation of ) are centered on each cell.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 " ?2Dt, 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: From Equation (2), B 2 C i,j px, t, ∆xq {Bx 2 9BC i,j px, t, ∆xq {B t and it can be concluded that 50% of dx micro intervals (dx !∆x) of x 2 intervals get a positive temporal gradient (BC 1,0 px, t, ∆xq {Bt ą 0 and BC 3,0 px, t, ∆xq {Bt ą 0), and 50% a negative gradient (BC 1,0 px, t, ∆xq {Bt ă 0 and BC 3,0 px, t, ∆xq {Bt ă 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: B dimpHUF ˝RLEpx 4 ptqqq{Bt ă 0 if 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 !tHUF ˝RLEu -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: Figure 5 represents the normalized dimension obtained by dimpHUF ˝RLEpx 4 ptqqq{dimpx 4 ptqq 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 max HUF˝RLE psq 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].

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: 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: 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 in Equation ( 7), one gets and the local stability criterion becomes: Finally, the global stability criterion becomes: A Gaussian function will be used to simulate the variation of diffusion coefficient:

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: 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: 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 ∆apxq is introduced to the unitary jump length during the Brownian simulation process and the mean displacement of each particle becomes 1 `∆apxq.By posing apxq " 1 `∆apxq in Equation ( 7), one gets Dpxq " papxqq 2 and the local stability criterion becomes: Finally, the global stability criterion becomes: A Gaussian function will be used to simulate the variation of diffusion coefficient: ∆apx, c, σq " e 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).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 1024 / 8 128 x    .At t = 0, the diffusion cell is located the middle of the map (at t = 0, xi,j = 1 if i belongs to i4 = [448,576] otherwise xi,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 2 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 20 pixels, result not shown).According to Equation ( 21), for very small standard deviations, one gets: and for high standard deviations, one gets: 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( ( ( ))) / dim( ( ))  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 @j P r1, 1024s if i belongs to i 4 = [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 2 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 20 pixels, result not shown).According to Equation (21), for very small standard deviations, one gets: lim σÑ0 pn s px, c, σqq ď ∆x 2 " 128 2 " 16384, @x, @c, (24) and for high standard deviations, one gets: In these two cases, diffusion becomes linear and stability condition does not depend on x position.
Let us now analyze the two compression ratio, dimpHUF ˝RLEpi n ptqqq{dimpi n ptqq, of the two adjacent cells (i 3 , i 5 ) at i 4 (i 3 = [320,448], i 5 = [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 sc s px, c, σq that depend on the values of the standard deviations.0   , t = 50 0   , t = 16,384 0   , t = 28,377 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.
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.From our postulate, the theoretical stability criterion is equal to n sc s px, c, σq " p∆x{p1 `∆apx, c, σqqq 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.
given by Equation ( 21) for both cells i4 and i5 corresponding to maximal values of Figure 8.

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: 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: 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, . The stability criterion is given by: 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 2 , 256 2 and 512 2 pixels.Maximal values of the dimension are extracted (Figure 11).HUF˝RLE versus the theoretical stability criterion n sc s pxq " p∆x{p1 `∆apxqqq 2 given by Equation ( 21) for both cells i 4 and i 5 corresponding to maximal values of Figure 8.

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: In our case, the diffusion coefficient depends on the concentration C px, tq. Lee [17] proposed diffusion coefficients relating to the uptake of excess calcium by calcium chloride [18].The diffusion coefficient D takes the form: 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 " t´1, 0, 1, 2, 3, 4, 5u.The stability criterion is given by: λ " max pD pC px, tqqq ∆t 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 2 , 256 2 and 512 2 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.
Values of the theoretical stability threshold agree with simulated values (Figure 12), meaning that our method can be applied to nonlinear PDE.Values of the theoretical stability threshold agree with simulated values (Figure 12), meaning that our method can be applied to nonlinear PDE.Values of the theoretical stability threshold agree with simulated values (Figure 12), meaning that our method can be applied to nonlinear PDE.

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).

Discussion
To summarize, the mathematical criterion of stability 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].

Conclusions
In this paper, we showed that a relation exists between the well-known stability criterion 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 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].

Conclusions
In this paper, we showed that a relation exists between the well-known stability criterion λ " D∆t{ p∆xq 2 ď 1{2 applied on the explicit numerical scheme C t`∆t x " C t i `λ ´Ct x´∆x ´2C t x `Ct x`∆x to guarantee the convergence of the 1D non-stationary PDE BC px, tq {Bt " DB 2 C px, tq {Bx 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.

)Figure 1 .
Figure 1.Snapshot of the Monte Carlo simulation for a grid size of 255 2 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 1 .
Figure 1.Snapshot of the Monte Carlo simulation for a grid size of 255 2 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 .Figure 2 .
Figure2.Histograms of the program sizes for the system X (a-e) and the sub system 4 x (f-j) corresponding to the Monte Carlo simulation for a grid size of 255 2 at different MCS times (0, 10, 300, Figure2.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 2 at different MCS times (0, 10, 300, 655 (stability criterion), and 1000).(a) Mean and standard deviation are computed and Gaussian density is plotted.

2 (
)) / dim( ( )) HUF RLE x t x t  versus the diffusion time (in MCS) for different grid sizes of cells (250, 300, …, 950, 1000).Values of the time max HUF RLE t  are plotted (graph on the right corner) versus the theoretical stability criterion

Figure 6
Figure 6.Values of i3, i5) at i4 (i3 = [320,448], i5 = [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 ( , , ) sc s n x c  that depend on the values of the standard deviations.

Figure 9 .
Figure 9. Plot of MCS max HUF RLE t 

Figure 9 .
Figure 9. Plot of MCS t maxHUF˝RLE versus the theoretical stability criterion n sc s pxq " p∆x{p1 `∆apxqqq 2 given by Equation (21) for both cells i 4 and i 5 corresponding to maximal values of Figure8.

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

Figure 10 .Figure 10 .
Figure 10.Snapshots of the nonlinear Monte Carlo simulation for a grid size of 129 2 pixels at five Monte Carlo Steps (MCS) times with different values of a given by D pC px, tqq " D 0 { p1 `aC px, tqq and snapshots corresponding to stability criterion and their associated Monte Carlo Steps.
versus the diffusion time (in MCS) for different values of diffusion coefficient

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

Figure 13 .
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.

Figure 13 .
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.
Figure 12.Plot of MCS t max HUF˝RLE versus the theoretical stability criterion for both cells i 4 and i 5 corresponding to maximal values of Figure 11 (Size = 257 2 ) and Size = (129 2 , 512 2 ).