Application of Shannon Entropy to Reaction–Diffusion Problems Using the Stochastic Finite Difference Method

TL;DR

This paper explores how Shannon entropy can be used to study diffusion processes, offering new insights into mixing and particle distribution.

Contribution

The paper introduces Shannon entropy as a novel metric for analyzing diffusion dynamics and mixing efficiency.

Findings

Shannon entropy correlates with diffusion kinetics in numerical simulations.

The method provides a robust way to describe diffusion-driven mixing.

Applications span engineering, environmental science, and biophysics.

Abstract

In this study, we introduce Shannon entropy as a key metric for assessing concentration variability in diffusion processes. Shannon entropy quantifies the uncertainty or disorder in the spatial distribution of diffusing particles, providing a novel perspective on diffusion dynamics. This proposed approach enables a more comprehensive characterization of mixing efficiency, equilibrium states, and transient diffusion behavior. Numerical simulations performed using the finite difference method in the MAPLE 2025 symbolic computing environment illustrate how entropy evolution correlates with diffusion kinetics. The computational model used in this study is based on a previously developed framework from our earlier research, ensuring consistency and validation of the results. The findings suggest that Shannon entropy can serve as a robust descriptor of diffusion-driven mixing, with potential applications in engineering, environmental science, and biophysics.

Figures

Equations

Full text

1. Introduction

The finite difference method (FDM) [1,2] is a fundamental numerical technique used in computational engineering to approximate the solutions of differential equations. It has been widely applied in fields such as heat transfer [3], electromagnetics [4], and geomechanics, where partial differential equations (PDEs) often arise. The classical FDM approach replaces continuous derivatives with discrete finite difference approximation, offering a straightforward and efficient tool for solving boundary value problems on structured grids [5]. In recent years, research has increasingly focused on extending the FDM to stochastic settings, where system parameters or boundary conditions exhibit randomness. This leads to implementation and application of the stochastic finite difference method (SFDM), which enables the analysis of uncertainty propagation in systems modeled by PDEs. A prominent approach within this context is the generalized stochastic perturbation technique [6,7], which involves Taylor series expansions of the stochastic input variables. This method has previously been applied in the context of the finite element method (FEM) [6], and the boundary element method (BEM) [7], and here we adapt it for the FDM with novel capabilities: (1) arbitrary-order stochastic expansions, (2) compatibility with general probability distributions [8], and (3) symbolic computation of hierarchical equations using the MAPLE 2025 system. The stochastic perturbation method is frequently contrasted with the Karhunen-Loeve and polynomial chaos expansion [9] or the Monte-Carlo simulation approach [10,11]. In this study, we apply the SFDM to a reaction–diffusion problem—a class of models that describe how chemical substances diffuse and react over time and space. Reaction–diffusion equations close to the advection-diffusion ones [12] are central in various scientific disciplines, including biophysics [13], chemical engineering, and ecology, due to their ability to capture both spatial transport and reaction dynamics. Our implementation focuses on one-dimensional steady and unsteady cases with a single uncertain parameter, although the methodology is general and can be extended to multi-dimensional and nonlinear problems.

A key contribution of this work is the introduction of Shannon (information) entropy as a metric for quantifying uncertainty in concentration profiles over time [14]; such type of the entropy has been used also in structural damage identification [15]. While classical diffusion analysis often relies solely on concentration gradients or variances, entropy provides a complementary perspective, capturing disorder and mixing efficiency in a probabilistic sense [16,17]. By tracking the temporal evolution of entropy in stochastic diffusion processes, we gain insight into system behavior, equilibrium tendencies, and the influence of input variability [18,19]. Although entropy-based methods have been used in theoretical studies of diffusion [20,21,22], their integration into computational frameworks for uncertainty quantification is still relatively rare [23,24]. Let us note that fractional diffusion and fractional order derivatives became recently very popular in mathematical and numerical modeling [25,26].

Our work addresses this gap by embedding entropy computation directly into the SFDM algorithm. We present numerical experiments that demonstrate how entropy evolves in time and space, reflecting the interplay between diffusion kinetics and stochastic variability. The results are verified using symbolic differentiation and parametric studies, with particular attention to the role of the input uncertainty [27]. In summary, this study offers a computational framework that combines stochastic finite difference modeling with entropy-based uncertainty analysis. This approach enables a more comprehensive understanding of diffusion-driven processes and provides a foundation for future extensions to complex, real-world systems with uncertain inputs.

2. Deterministic Model

We begin by formulating a deterministic boundary value problem describing the steady-state concentration profile C(x), where x∈R, within a classical reaction–diffusion framework. This formulation applies to a linear system governed by reaction–diffusion dynamics over a one-dimensional spatial domain, as discussed in the literature [28,29]. This scenario assumes a planar region in which molecular concentration exhibits spatial variability, while the boundary conditions and model parameters remain fixed and deterministic. A conceptual representation of the concentration field within this domain is provided in Figure 1—although presented schematically, the analysis and numerical model refer strictly to a one-dimensional domain.

This problem is described by a system of differential equations:

$$[eqn]$$

for which a closed-form analytical solution is available and given by $[eqn]$ , where e denotes Euler’s number, x denotes the current coordinate, and l stands for the channel length. This closed-form expression serves as a benchmark for validating the numerical accuracy of the proposed finite difference scheme, particularly in the context of preparing the computational framework for entropy-based uncertainty quantification. Alternatively, the same problem can be approached numerically using the finite difference method (FDM). In this approach, the spatial domain [0,l] is discretized into a finite number of the uniform subintervals of length Δ, where the nodal positions are defined as $[eqn]$ , for $[eqn]$ . The approximations of concentration values at these discrete locations are denoted by $[eqn]$ . Spatial derivatives appearing in the governing equations are approximated using central difference schemes, as detailed in [5]:

$$[eqn]$$ $$[eqn]$$

The further numerical solution needs a definition of the imposed boundary conditions, namely

$$[eqn]$$ $$[eqn]$$

for x_n+_2 being a fictitious node out of the initial domain. Then, we solve for $[eqn]$ , $[eqn]$ from the following linear equations system:

$$[eqn]$$

which can be represented as the matrix equation

$$[eqn]$$

with tri-diagonal coefficients:

$$[eqn]$$

As is known, the numerical solution to this problem by the FDM has quadratic convergence since the discretization errors for both first- and second-order derivatives are equal to $[eqn]$ . The solution to the non-linear boundary value problem is equivalent to the reaction–diffusion problem [28]:

$$[eqn]$$

with slightly modified boundary conditions introduced to test the sensitivity of the stochastic solution

$$[eqn]$$

One may then proceed with the discretization $[eqn]$ in quite a similar way until one obtains

$$[eqn]$$

As we will see in the next section, we can also provide the multipoint versions of the finite difference method for the needs of the stochastic computational solution or apply some error analysis techniques to determine the accuracy of determination of the particular probabilistic moments [30].

Figure 2 presents the numerical validation of the finite difference method (FDM) applied to the reaction–diffusion system. The computed solution is compared with the reference deterministic profile, demonstrating good consistency across the spatial domain. The observed agreement confirms the numerical stability and accuracy of the implemented scheme, providing a reliable foundation for the subsequent entropy-based uncertainty quantification. Ensuring fidelity at this stage is essential, as even minor numerical artifacts could distort entropy metrics in stochastic extensions.

3. A Stochastic Process Approach

The theory of stochastic processes provides a mathematical framework for describing systems that evolve over time and space according to probabilistic laws. This formalism is extensively applied in diverse fields of science and engineering to model complex, inherently uncertain phenomena [31]. In the context of reaction–diffusion problems, we now extend the formulation from a random variable $[eqn]$ , which induces a spatially random field $[eqn]$ , to a time-dependent random variable $[eqn]$ , leading to a spatiotemporal stochastic field $[eqn]$ . Let us consider a probability space $[eqn]$ , and a discrete time domain $[eqn]$ . A function $[eqn]$ is defined as a stochastic process if the following condition holds:

$$[eqn]$$

which means that for each fixed time $[eqn]$ , the function $[eqn]$ is a random variable on $[eqn]$ . In this framework, the concentration field $[eqn]$ can be expressed as a function of space, probability, and time:

$$[eqn]$$

Analogously to classical random variables, stochastic processes are characterized by statistical descriptors such as the expectation, variance, and higher-order moments [31]:

• expected value

$$[eqn]$$

variance

$$[eqn]$$

autocovariance

$$[eqn]$$

autocorrelation

$$[eqn]$$

autocorrelation coefficient

$$[eqn]$$

We consider a specific example of a stochastic process given by a power-type time series of the form

$$[eqn]$$

where $[eqn]$ is a Gaussian (normally distributed) random variable. This formulation introduces randomness through the multiplicative coefficient, while the time dependence follows a deterministic power-law structure. It should be noted that the use of a polynomial time series (19) does not limit the presented method solely to the analysis of polynomial entropy.

Based on this definition, the key central probabilistic moments of the process—namely, the expectation, variance, and the kth central moment—can be derived analytically:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The computational algorithm including the time series in the SFDM analysis using the generalized stochastic perturbation technique described before for a single random variable is displayed using the flowchart in Table 1.

The core concept, illustrated schematically in Table 1, involves discretizing the temporal domain using a uniform time step Δt. At each discrete time point, the first two probabilistic moments (i.e., the mean and variance) of the input stochastic time series describing the channel height are computed. Subsequently, the boundary value problem is solved iteratively using the response function method (RFM) integrated with the stochastic finite difference method. In this approach, the stochastic response function is approximated by a low-order polynomial expansion in terms of the input random variable, and the coefficients of this expansion are determined using symbolic or numerical techniques. This allows efficient estimation of output probabilistic moments by analytical integration over the input probability density function. The RFM-based framework enables a parametric analysis of the system behavior, capturing the dependence of the diffusion process on both time and the variability of uncertain parameters, while significantly reducing the computational cost compared to traditional Monte Carlo simulations.

All the computations are completed with iterative determination of the Shannon entropy [14,18,19], which is presented as a function of time during the diffusion process. An explanation of the final point in this algorithm is delivered for a continuous random variable b having the probability density function $[eqn]$ discretized by a set of subintervals of the same length δ; this formulation directly follows Monte Carlo simulation-based histograms. The following representation of $[eqn]$ can be made in each ith subinterval [24]:

$$[eqn]$$

Further, let us introduce the additional quantized random variable $[eqn]$ defined as

$$[eqn]$$

whose probability equals

$$[eqn]$$

So, probabilistic entropy according to the Shannon definition corresponding to a partition of the histogram of the density function $[eqn]$ into n equal subsets can be calculated as

$$[eqn]$$

Symbol h refers here to the h-theorem proposed by Boltzmann and is further used to traditionally denote probabilistic entropy under consideration.

4. Reaction–Diffusion with Uncertainty

We now consider the following boundary value problem defined on a stochastic domain:

$$[eqn]$$

where $[eqn]$ is a truncated Gaussian random variable. The concentration $[eqn]$ is modeled as a polynomial time series whose coefficients are also truncated Gaussian variables with prescribed first and second probabilistic moments. To solve this stochastic boundary value problem, the generalized stochastic perturbation method is employed. In this framework, the random variable b is characterized by a known probability density function $[eqn]$ , and the stochastic response function $[eqn]$ is analyzed in terms of its central probabilistic moments. The mth central moment is defined as

$$[eqn]$$

The fundamental idea behind the stochastic perturbation approach is inspired by classical Taylor expansions. Both input variables and state functions are approximated using truncated Taylor series expanded about their mean values. The expansion is formulated in terms of a small perturbation parameter $[eqn]$ . For instance, the nth-order truncated expansion of a stochastic concentration function $[eqn]$ can be written as [27]

$$[eqn]$$

where

$$[eqn]$$

is the first variation of b about its expected value and, similarly,

$$[eqn]$$

is the second variation of b around its expected value, where nth-order variation can be expressed accordingly. Traditionally, the stochastic perturbation approach to all the physical problems is entered by the respective perturbed equations of the zeroth, first, and successively higher orders being a modification of the initial deterministic formulation. There holds for $[eqn]$ :

• the zeroth order partial differential equation

$$[eqn]$$

• the nth-order partial differential equation

$$[eqn]$$

Having solved those equations for $[eqn]$ and its higher orders, respectively, (specifically its partial derivatives with respect to random input within all discrete points of the grid), we derive the expressions for the expected values and the other moments. Since we propose the finite difference scheme to solve these equations numerically, we need to expand additionally Equation (7) to obtain the stochastic finite difference method scheme. There holds for $[eqn]$ :

• the zeroth equation

$$[eqn]$$

• the linear recursive equation for each next perturbation order k

$$[eqn]$$

to be solved all consecutively to determine the probabilistic characteristics of the random output. In order to calculate the expected values and higher-order probabilistic moments of $[eqn]$ , the same Taylor expansion is employed for the definitions of its probabilistic moments like

$$[eqn]$$

If there is a high random dispersion in the input random variable and the symmetric probability density function is chosen, then the generalized expansion simplifies to

$$[eqn]$$

where $[eqn]$ denotes the 2kth order probabilistic moment of the variable $[eqn]$ . When the probability density function is defined as the Gaussian with the standard deviation σ(b), we additionally obtain

$$[eqn]$$

Using such an extension of the random input, a desired efficiency of the expected values can be achieved by the appropriate choice of the perturbation parameter and maximum order corresponding to the particular input probability density function type, probabilistic moments interrelations, acceptable error of the computations, etc. This choice can reasonably be made by the comparative studies with the Monte Carlo simulations or theoretical results obtained by the direct (i.e., symbolic) integration. Consequently, the mth-order probabilistic moment of the system response is evaluated based on the polynomial representation constructed using the response function method. In this approach, the system response is approximated by a response function—typically a low-degree polynomial in the input random variable—whose coefficients are obtained through symbolic or numerical procedures. This method enables analytical integration with respect to the probability density function of the input parameter, yielding closed-form expressions for expectations, variances, and higher-order moments:

$$[eqn]$$

The computational effort remains comparable to the deterministic case and avoids the high cost of large-scale Monte Carlo simulations, which require a substantial number of samples and postprocessing steps such as regression via the least squares method. Let us also note that the expansion provided in Equation (39) may be completed automatically in computer algebra systems for the given natural order.

5. Numerical Illustration

The algorithm proposed has been implemented in the computer algebra environment MAPLE 2025. The length of the channel $[eqn]$ has been adopted as the input Gaussian parameter in this problem, where its initial expected value equals $[eqn]$ , while the coefficient of variation α(l) takes 0.1. The time series is represented as the linear polynomial, such that $[eqn]$ , where $[eqn]$ . Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 contain the results of the determination of probabilistic moments and coefficients for the concentration ratio $[eqn]$ using the stochastic perturbation-based finite difference method. We have in turn: the expectations (Figure 3), the coefficients of variation (Figure 4), the Shannon entropy (Figure 5) and skewness (Figure 6), and finally, kurtosis (Figure 7). We discretize the entire computational domain, $[eqn]$ , into 20 equidistant subintervals along its length.

Expected values remain the same for x = 0 independently of the time of a process and of the input coefficient of variation, and are equal to 1. On the opposite side of a computational domain, the concentration’s expectation decreases linearly from the initial unit value to 2, also independently of the input coefficient of random dispersion. Spatial variations of the expected value of the concentration are highly nonlinear along a computational domain showing definite convexity independent of a time. The standard deviation of this concentration coefficient is highly dependent upon the input coefficient of the variation—the higher $[eqn]$ , the larger the resulting $[eqn]$ . They have zero values at any time for $[eqn]$ , where we apply the deterministic boundary condition, and then they increase moderately to reach the maximum at the opposite edge of a domain. They also increase almost linearly together with time in this particular process.

The Shannon entropy, shown in Figure 5, decreases over time, indicating a progressive reduction in system uncertainty and convergence towards equilibrium. Simultaneously, the coefficients of variation of concentration, shown in Figure 4, exhibit significant fluctuations in early diffusion stages, followed by stabilization, suggesting that initial randomness dissipates as the system evolves. The presented approach offers a refined perspective on diffusion dynamics, with potential applications in reliability assessment, engineering simulations, and environmental modeling. The resulting diagram of the expected concentration profile, along with the corresponding standard deviation (Figure 4b), shows the output coefficients of variation. It is remarkable that for any combination of the parameters of time $[eqn]$ and length $[eqn]$ , the output coefficient of variation is significantly higher than the input one. The resulting ratio of output versus input coefficients increases together with an input random dispersion level. Finally, one may notice that the output coefficient $[eqn]$ increases nonlinearly together with l, and almost linearly with respect to the time $[eqn]$ .

The third central probabilistic moments given in Figure 6 are positive everywhere and they also significantly increase together with an input coefficient $[eqn]$ . Absolute values of the third moment are very close to 0, while skewness (in Figure 6a) belongs to the interval $[eqn]$ , which gives an almost symmetrical probability density function, especially for smaller values of the input $[eqn]$ . Both quantities are equal to 0 for $[eqn]$ , where the deterministic boundary conditions are applied and, of course, they increase together with l and with time $[eqn]$ . All the aforementioned surfaces representing parametric variability of the probabilistic moments are very smooth without any local oscillations and singularities, that may result from higher-order partial derivatives.

Finally, we study the fourth central probabilistic moments of the concentration as well as the kurtosis (Figure 7) for three different values of $[eqn]$ . Fourth moments appear to be positive, although even smaller than the third ones, while the resulting kurtosis keeps close to 3, especially for $[eqn]$ The surfaces obtained for larger values of the input $[eqn]$ are not as regular and smooth as those of the lower probabilistic moments and the coefficients. Their values both increase remarkably, together with the time of a process, but are almost independent of the parameter $[eqn]$ . A comparison of Figure 4, Figure 6 and Figure 7 shows that higher probabilistic moments become less and less dependent upon the domain length $[eqn]$ .

However, the most important conclusion here is that for smaller random dispersion of the coefficients in the time series representation of $[eqn]$ the final concentration factor keeps close to the Gaussian distribution. Its distance to this PDF increases together with the value of the input coefficient of variation $[eqn]$ .

6. Concluding Remarks

I.As documented above, Shannon entropy may serve as some additional uncertainty quantification and propagation measure close to the coefficients of variation not only in steady-state but also in non-stationary problems with random coefficients. It should be noted that the analysis of coefficients of variation provides precise information about the nature of the process at a given time point, enabling the assessment of the relative dispersion of random values at specific spatiotemporal locations. In contrast, Shannon entropy can serve as a complementary measure for uncertainty quantification and propagation, offering a more global insight into the process dynamics. Unlike the local characteristics described by coefficients of variation, entropy captures the trend of uncertainty propagation over time and space, both in steady-state and non-stationary problems involving random coefficients. As such, it represents a valuable tool for the comprehensive analysis of complex processes in which the distribution of uncertainty evolves significantly over the course of the system’s development. Although the determination of Shannon entropy demands the application of the Monte Carlo simulation, it gives accumulated knowledge about uncertainty propagation, which in probabilistic moments-based analysis follows the graphs of expectations and second-order statistics. It has also been shown that an application of the symbolic algebra environment enables a solution of the generalized stochastic-perturbation-based discrete equations’ system as well as for an introduction of the perturbation parameter ε into the final solution expansions and the resulting probabilistic moments. Thanks to the very extensive linear algebra tools in MAPLE 2025 (as well as in other symbolic computation systems), we can efficiently implement the stochastic perturbation scheme using the response function method [27]. This approach relies on constructing low-order polynomial approximations of the system response with respect to uncertain input parameters and enables analytical computation of output probabilistic moments. As a result, we obtain a compact and flexible framework that is compatible with both symmetric and non-symmetric coefficient matrices, and capable of handling nonlinear effects due to input uncertainty.II.The numerical analysis presented here demonstrates that Shannon entropy systematically decreases over time, indicating a progressive reduction in uncertainty and increased spatial uniformity of the concentration field. The highest entropy values are observed near the boundary, suggesting localized disorder at the system’s edge. Simultaneously, the coefficient of variation stabilizes rapidly, with initially elevated values diminishing significantly as the system evolves. This convergence highlights the decreasing influence of stochastic fluctuations and the growing regularity of the solution. The consistency between entropy and variance-based indicators confirms that Shannon entropy serves as a reliable measure of uncertainty and information content in stochastic diffusion–reaction processes. These findings support the utility of entropy-based diagnostics in identifying spatial zones of increased randomness and guiding the analysis of complex transport phenomena. Numerical results presented above also indicate that the system evolves towards a steady state, where concentration differences are minimized, and both entropy and the coefficient of variation decrease over time. The initially high entropy and variance highlight the presence of dynamic mixing processes and stochastic fluctuations, which gradually diminish, leading to a more deterministic concentration distribution. These findings confirm that the stochastic modeling approach effectively captures the initial uncertainty within the system and its subsequent reduction as diffusion progresses, providing a robust framework for analyzing reaction–diffusion dynamics in uncertain environments.III.The generalized stochastic finite difference method seems to be the efficient probabilistic computational tool to model the reaction–diffusion problems with random parameters and time series and should be further developed for unsteady problems, and also for inhomogeneous domains [32] in 2D or 3D [29]. As was demonstrated above, this method can be applied to problems with physical (mechanical or chemical) parameters having some stochastic boundary waviness, although essentially more complicated situations with multiple randomness sources may also be modeled. It has been shown that the stochastic methodology described above may also be applied directly to the implementations of the FDM on irregular grids [33] or combined with some meshless techniques [34]. Further comparative studies against the Monte Carlo simulation programs may give a direct answer concerning the precision of the perturbation-based SFDM in problems that have some input parameters defined as time series with random coefficients.IV.Although the proposed stochastic finite difference method (SFDM) shows promising results for linear and stationary reaction–diffusion problems, its extension to more realistic, nonlinear, and time-dependent systems remains a significant challenge. In such cases, several numerical difficulties must be addressed, i.e., nonlinearities often introduce stiffness and sensitivity to perturbations, which can severely constrain the choice of discretization parameters. Similarly, time-dependent problems impose additional stability conditions on the numerical scheme, resulting in high computational costs for long-time simulations. To overcome these challenges, future versions of the SFDM should incorporate implicit integration techniques, which are known to offer superior stability properties for stiff systems. Furthermore, adaptive time-stepping strategies could be employed to dynamically control the temporal resolution depending on the rate of change in the system. These enhancements, while increasing computational complexity, would allow the method to handle more general classes of reaction–diffusion systems, including those with strong nonlinearity and spatial heterogeneity.

Work is currently underway to investigate such extensions, with the aim of broadening the applicability of SFDM beyond textbook examples. (See Table 2 and Table 3).