Calculating how long it takes for a diffusion process to effectively reach steady state without computing the transient solution

TL;DR

This paper introduces a novel method to estimate the finite transition time for diffusion processes to reach near steady state without solving the transient equation explicitly, using moments and asymptotic analysis.

Contribution

It proposes a new approach based on asymptotic behavior of the cumulative distribution function to accurately compute finite transition times using moments.

Findings

The new approach provides highly accurate estimates of transition times.

Mean action time methods are less accurate for diffusion processes.

Accuracy improves with higher moments used in the new method.

Abstract

Mathematically, it takes an infinite amount of time for the transient solution of a diffusion equation to transition from initial to steady state. Calculating a \textit{finite} transition time, defined as the time required for the transient solution to transition to within a small prescribed tolerance of the steady state solution, is much more useful in practice. In this paper, we study estimates of finite transition times that avoid explicit calculation of the transient solution by using the property that the transition to steady state defines a cumulative distribution function when time is treated as a random variable. In total, three approaches are studied: (i) mean action time (ii) mean plus one standard deviation of action time and (iii) a new approach derived by approximating the large time asymptotic behaviour of the cumulative distribution function. The new approach leads to a simple formula for calculating the finite transition time that depends on the prescribed tolerance $\delta$ and the $(k-1)$th and $k$th moments ($k \geq 1$) of the distribution. Results comparing exact and approximate finite transition times lead to two key findings. Firstly, while the first two approaches are useful at characterising the time scale of the transition, they do not provide accurate estimates for diffusion processes. Secondly, the new approach allows one to calculate finite transition times accurate to effectively any number of significant digits, using only the moments, with the accuracy increasing as the index $k$ is increased.

Tables

Table 1. Table 1: Errors ( 38 ), corresponding to the three different global transition time estimates (defined in equations 26 , 27 , 28 and 36 ), for each of the three test cases. The value of ε s ( 3 ) superscript subscript 𝜀 𝑠 3 \varepsilon_{s}^{(3)} is calculated using different combinations of tolerance δ 𝛿 \delta and moment index k 𝑘 k .

Case	${\epsilon }_s^{\left(1\right)}$	${\epsilon }_s^{\left(2\right)}$	${\epsilon }_s^{\left(3\right)}$	${\epsilon }_s^{\left(3\right)}$	${\epsilon }_s^{\left(3\right)}$	${\epsilon }_s^{\left(3\right)}$
			[$k=2$, $\delta =0.02$]	[$k=2$, $\delta ={10}^{-3}$]	[$k=2$, $\delta ={10}^{-5}$]	[$k=5$, $\delta =0.02$]
A	0.3708	0.1354	0.0189	8.69e-04	7.64e-06	0.0200
B	0.3721	0.1356	0.0188	8.58e-04	7.43e-06	0.0200
C	0.3708	0.1354	0.0189	8.69e-04	7.64e-06	0.0200

Table 2. Table 2: Comparison of the exact global transition time t ^ s subscript ^ 𝑡 𝑠 \widehat{t}_{s} (defined in equation 2 ) and the global transition time estimate t ^ s ( 3 ) superscript subscript ^ 𝑡 𝑠 3 \widehat{t}_{s}^{(3)} (defined by equations 36 and 26 ) for each combinations of tolerance δ = 10 − 1 , 10 − 2 , 10 − 3 , 10 − 4 , 10 − 5 , 10 − 6 𝛿 superscript 10 1 superscript 10 2 superscript 10 3 superscript 10 4 superscript 10 5 superscript 10 6 \delta=10^{-1},10^{-2},10^{-3},10^{-4},10^{-5},10^{-6} and moment index k = 1 , 2 , 5 , 10 𝑘 1 2 5 10 k=1,2,5,10 .

			$\delta ={10}^{-1}$	$\delta ={10}^{-2}$	$\delta ={10}^{-3}$	$\delta ={10}^{-4}$	$\delta ={10}^{-5}$	$\delta ={10}^{-6}$
	${\hat{t}}_s$		1.0311	1.9643	2.8975	3.8307	4.7639	5.6971
	${\hat{t}}_s^{\left(3\right)}$	[$k=1$]	1.1513	2.3026	3.4539	4.6052	5.7565	6.9078
	${\hat{t}}_s^{\left(3\right)}$	[$k=2$]	1.0354	1.9948	2.9542	3.9136	4.8730	5.8324
	${\hat{t}}_s^{\left(3\right)}$	[$k=5$]	1.0311	1.9643	2.8975	3.8308	4.7640	5.6973
Case A	${\hat{t}}_s^{\left(3\right)}$	[$k=10$]	1.0311	1.9643	2.8975	3.8307	4.7639	5.6971
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=1$]	1.17e-01	1.72e-01	1.92e-01	2.02e-01	2.08e-01	2.12e-01
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=2$]	4.14e-03	1.55e-02	1.96e-02	2.16e-02	2.29e-02	2.38e-02
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=5$]	4.54e-05	2.63e-06	1.23e-05	2.05e-05	2.52e-05	2.84e-05
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=10$]	2.28e-09	5.79e-08	3.98e-07	5.45e-11	1.06e-10	2.39e-08
	${\hat{t}}_s$		31.0746	59.1707	87.2666	115.3624	143.4582	171.5541
	${\hat{t}}_s^{\left(3\right)}$	[$k=1$]	34.5967	69.1934	103.7901	138.3867	172.9834	207.5801
	${\hat{t}}_s^{\left(3\right)}$	[$k=2$]	31.1946	60.1603	89.1312	118.1046	147.0794	176.0552
	${\hat{t}}_s^{\left(3\right)}$	[$k=5$]	31.0689	59.1697	87.2706	115.3715	143.4724	171.5733
Case B	${\hat{t}}_s^{\left(3\right)}$	[$k=10$]	31.0749	59.1707	87.2665	115.3624	143.4582	171.5541
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=1$]	1.13e-01	1.69e-01	1.89e-01	2.00e-01	2.06e-01	2.10e-01
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=2$]	3.86e-03	1.67e-02	2.14e-02	2.38e-02	2.52e-02	2.62e-02
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=5$]	1.83e-04	1.63e-05	4.64e-05	7.87e-05	9.86e-05	1.12e-04
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=10$]	8.09e-06	1.55e-07	6.07e-08	1.26e-08	1.66e-08	3.63e-08
	${\hat{t}}_s$		0.6444	1.2277	1.8109	2.3942	2.9774	3.5607
	${\hat{t}}_s^{\left(3\right)}$	[$k=1$]	0.7196	1.4391	2.1587	2.8782	3.5978	4.3173
	${\hat{t}}_s^{\left(3\right)}$	[$k=2$]	0.6471	1.2467	1.8464	2.4460	3.0456	3.6453
	${\hat{t}}_s^{\left(3\right)}$	[$k=5$]	0.6444	1.2277	1.8110	2.3942	2.9775	3.5608
Case C	${\hat{t}}_s^{\left(3\right)}$	[$k=10$]	0.6444	1.2277	1.8109	2.3942	2.9774	3.5607
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=1$]	1.17e-01	1.72e-01	1.92e-01	2.02e-01	2.08e-01	2.12e-01
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=2$]	4.14e-03	1.55e-02	1.96e-02	2.16e-02	2.29e-02	2.38e-02
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=5$]	4.56e-05	2.78e-06	1.26e-05	2.08e-05	2.51e-05	2.83e-05
	$\left|{\hat{t}}_s-{\hat{t}}_s^{\left(3\right)}\right|/\left|{\hat{t}}_s\right|$	[$k=10$]	2.76e-07	2.02e-07	9.79e-08	3.20e-07	5.92e-08	7.36e-08

Table 3. Table 3: Computed values of the constants γ k subscript 𝛾 𝑘 \gamma_{k} and θ k subscript 𝜃 𝑘 \theta_{k} that appear in the global transition time estimate ( 41 ) for the generalised version of Case A.

$k$	${\gamma }_k$	${\theta }_k$	$\left|{\gamma }_k-4/{\pi }^2\right|$	$\left|{\theta }_k-4/\pi \right|$
2	0.4167	1.2000	1.14e-02	7.32e-02
4	0.4054	1.2712	1.60e-04	2.08e-03
6	0.4053	1.2732	2.03e-06	3.90e-05
8	0.4053	1.2732	2.51e-08	6.41e-07
10	0.4053	1.2732	3.10e-10	9.86e-09
12	0.4053	1.2732	3.83e-12	1.46e-10
14	0.4053	1.2732	4.72e-14	2.10e-12
16	0.4053	1.2732	5.83e-16	2.95e-14
18	0.4053	1.2732	7.20e-18	4.10e-16
20	0.4053	1.2732	8.89e-20	5.62e-18

Equations

$$\frac{u(x,t)-u_{\infty }(x)}{u_0(x)-u_{\infty }(x)}=\delta ,$$

$${\hat{t}}_s:=\underset{x\in \mathcal{R}}{\sup }{t}_s(x)$$

$$F(t;x):=1-\frac{u(x,t)-u_{\infty }(x)}{u_0(x)-u_{\infty }(x)},$$

$$\text{ts(x) is the value of t>0 satisfying F(t;x)=1−δ}.$$

$$f(t;x)=\frac{1}{u_{\infty }(x)-u_0(x)}\frac{\partial }{\partial t}\left[u(x,t)-u_{\infty }(x)\right],$$

$$M_k(x):={\int }_0^{\infty }{t}^{k}f(t;x)dt.$$

$$\frac{\partial u}{\partial t}=D\frac{{\partial }^{2}u}{\partial x^2},$$

$$u(x,0)=u_0(x),$$

$$a_{L}u(l_0,t)-b_L\frac{\partial u}{\partial x}(l_0,t)=c_L,$$

$$a_{R}u(l_m,t)+b_R\frac{\partial u}{\partial x}(l_m,t)=c_R,$$

$$D{u}_{\infty }^{\prime \prime }(x)=0,x\in \mathcal{L},$$

$$a_L{u}_{\infty }(l_0)-b_L{u}_{\infty }^{\prime }(l_0)=c_L,$$

$$a_R{u}_{\infty }(l_m)+b_R{u}_{\infty }^{\prime }(l_m)=c_R.$$

$${\int }_{l_0}^{l_m}{u}_{\infty }(x)dx={\int }_{l_0}^{l_m}{u}_0(x)dx,$$

$$M_k(x)=\frac{1}{h(x)}{\int }_0^{\infty }{t}^k\frac{\partial }{\partial t}\left[u(x,t)-u_{\infty }(x)\right]dt,$$

$$M_k(x)=\frac{k}{h(x)}{\int }_0^{\infty }{t}^{k-1}\left[u_{\infty }(x)-u(x,t)\right]dt.$$

$${\overline{M}}_k(x):=M_k(x)h(x)=k{\int }_0^{\infty }{t}^{k-1}\left[u_{\infty }(x)-u(x,t)\right]dt,$$

$${\overline{M}}_k^{\prime }(x)=k{\int }_0^{\infty }{t}^{k-1}\left[u_{\infty }^{\prime }(x)-\frac{\partial u}{\partial x}(x,t)\right]dt,$$

$${\overline{M}}_k^{\prime \prime }(x)=k{\int }_0^{\infty }{t}^{k-1}\left[u_{\infty }^{\prime \prime }(x)-\frac{{\partial }^{2}u}{\partial x^2}(x,t)\right]dt.$$

$${\overline{M}}_k^{\prime \prime }(x)=\frac{k}{D}{\int }_0^{\infty }{t}^{k-1}\frac{\partial }{\partial t}\left[u_{\infty }(x)-u(x,t)\right]dt.$$

$${\overline{M}}_k^{\prime \prime }(x)=-\frac{k}{D}{\overline{M}}_{k-1}(x),x\in \mathcal{L},$$

$$a_L{\overline{M}}_k(l_0)-b_L{\overline{M}}_k^{\prime }(l_0)=0,$$

$$a_R{\overline{M}}_k(l_m)+b_R{\overline{M}}_k^{\prime }(l_m)=0,$$

$${\overline{M}}_k(x)$$

$$G_k(x):=-\frac{k}{D}\int \int {\overline{M}}_{k-1}(x)dxdx.$$

$$\mathbf{Ac}=\mathbf{b},$$

$$\mathbf{A}$$

$$\mathbf{b}$$

$$G_k^{\prime }(x)=-\frac{k}{D}\int {\overline{M}}_{k-1}(x)dx.$$

$${\int }_{l_0}^{l_m}u(x,t)dx={\int }_{l_0}^{l_m}{u}_{\infty }(x)dx,$$

$${\int }_{l_0}^{l_m}{\overline{M}}_k(x)dx=0.$$

$$\mathbf{A}$$

$$\mathbf{b}$$

$$K_k(x)=\int G_k(x)dx,$$

$$M_k(x)$$

Full text

1 Introduction

A common question of practical interest that arises when modelling transport phenomena is how long does the process take? For example, how long does it take for a bar to heat up or how long does it take for an initial concentration field to distribute evenly? Mathematically, the correct answer is the highly impractical answer: it takes an infinite amount of time for the transient solution to transition from initial to steady state.

Finite (and hence practically useful) answers to the above question have been studied and proposed by many authors [9, 12, 10, 11, 24, 6, 16, 22, 21]. In this paper, we consider a finite transition time defined as the amount of time required for the transient solution to transition from its initial state to within a prescribed tolerance of its steady state [24]. To illustrate this definition, consider a transport process on an interval $\mathcal{L}:=(l_{0},l_{m})$ and let $u(x,t)$ be the transient solution at position $x\in\mathcal{L}$ and time $t>0$, $u_{0}(x)$ be the given initial solution and $u_{\infty}(x)$ be the steady state solution. We make the following assumptions:

1. (i)

$u_{0}(x)\neq u_{\infty}(x)$ for all $x\in\mathcal{L}$, 2. (ii)

For a given $x\in\mathcal{L}$, $u(x,t)$ monotonically increases to $u_{\infty}(x)$ (implying $u(x,t)<u_{\infty}(x)$) for all $t>0$ if $u_{0}(x)<u_{\infty}(x)$ or $u(x,t)$ monotonically decreases to $u_{\infty}(x)$ (implying $u(x,t)>u_{\infty}(x)$) for all $t>0$ if $u_{0}(x)>u_{\infty}(x)$, 3. (iii)

$u_{0}(x)$ is piecewise continuous on $\mathcal{L}$.

In this work, we distinguish between local and global transition times as follows. The local transition time provides a finite measure of the time required to transition from initial to steady state at position $x$. On the other hand, the global transition time is a finite measure of the time required for the entire transport process to effectively reach steady state, and hence provides a practical finite answer to the question of how long the process takes. Both of these quantities are defined below, where we remark that the assumption (ii) listed above ensures that the left-hand side of equation (1) is always positive.

{definition}

[Local transition time]

The local transition time, denoted by $t_{s}(x)$, provides the transition time as a function of position $x\in\mathcal{R}:=\{x\in\mathcal{L}\,|\,u_{0}(x)\neq u_{\infty}(x)\}$, and is defined as the value of $t>0$ satisfying:

[TABLE]

where $0<\delta\ll 1$ is a specified tolerance.

{definition}

[Global transition time]

The global transition time is defined as the supremum of the local finite transition time:

[TABLE]

where the supremum is taken over the domain of $t_{s}(x)$.

A common finite measure of the time required to reach steady state is the mean action time [24, 14] introduced by McNabb and Wake [22]. The mean action time takes a probabilistic approach to the problem by treating time $t$ as continuous random variable with support $[0,\infty)$ and computing its mean or expected value. This is achieved by utilising the observation that the function

[TABLE]

satisfies $F(0;x)=0$ and $\lim_{t\rightarrow\infty}F(t;x)=1$ and hence defines a cumulative distribution function of $t$, parametrised in terms of the position $x$. Using the cumulative distribution function, the local transition time (Definition 1) can now be reformulated as follows:

[TABLE]

Note that equation (4) is now equivalent to a classical inverse problem in probability theory: find $t_{s}(x)$ such that $P(t\geq t_{s}(x))=\delta$.

Differentiating (3) with respect to $t$ yields the probability density function [5, 24, 23]:

[TABLE]

which allows the $k$th moment, at position $x$, to be defined as follows:

[TABLE]

The mean action time at position $x$, denoted henceforth by $\text{MAT}(x)$, is then defined by the mean or first moment $M_{1}(x)$.

The attractiveness of the mean action time is that it is possible to obtain the explicit form of $\text{MAT}(x)$, and hence an estimate of the finite amount of time required for the transient solution to effectively reach steady state at position $x$, without explicitly calculating the transient solution $u(x,t)$ [6, 24, 23]. For this reason, the mean action time is a popular transition time estimate that has been used in several applications including freezing and thawing [16], morphogen formation in the formation of tissues and organs [1, 23] and groundwater modelling [24, 14, 15]. Another argument in favour of the mean action time is that it does not suffer from the subjectivity of choosing a threshold tolerance [6], as in Definition 1. While this is true, the question of how close the transient solution, evaluated at the supremum of the mean action time, is to the steady state still remains and involves some discretion. The opinion of this author is that the tolerance $\delta$ is useful and should be inferred by the physical problem under consideration. For example, in many problems involving comparison to experimental data, a suitable choice for $\delta$ is the error associated with the measurement device.

Since McNabb and Wake’s [22], notable contributions to the theory of mean action time, almost exclusively focussed on one-dimensional linear homogeneous problems, have appeared in several papers [16, 13, 6, 5, 24, 14, 23, 15]. For discussion on higher-dimensional, nonlinear and/or heterogeneous problems, which do not form the focus of this work, the reader is directed to three papers: [13, 6, 16].

Research on mean action time was revived in the early 2010s by both Berezhkovskii et al. [1], who proposed the definition of local accumulation time for a reaction-diffusion model related to the study of morphogen formation and Ellery et al. [6], who demonstrated that this definition was equivalent to McNabb’s mean action time. Over two papers, Ellery et al. [5, 6] derived the mean action time for a linear advection-diffusion-reaction process and presented a framework for computing the higher central moments (called moments of action by the authors), demonstrating how in each case exact expressions can be found without explicit calculation of the transient solution. The second moment of action, the variance of action time, was noted as being particularly insightful as a small value implies that the mean action time is a useful estimate of the time required to effectively reach steady state [6, 24]. Simpson et al. [24] applied and extended the theory to a linearised groundwater flow model, governed by a reaction-diffusion equation with constant source term, to study the time required for a transient response (such as an aquifer recharge/discharge process) to effectively reach steady state. Comparing the model predictions to a laboratory scale experiment, it was found that while the mean action time underestimates the time required to transition to steady state, significant improvement can be obtained by adding one standard deviation of action time. Jazaei et al. [14] extended the derivations for groundwater flow to time-dependent boundary conditions, where it was again argued, via visual observation of the experimental data, that the mean plus one standard deviation of action time was a good approximation for the time required for the system to effectively reach steady-state. Most recently, Simpson [23] and Jazaei et al. [15] have extended the theory to estimate the time required for the gradient or flux of the transient solution to approach their corresponding steady-state values.

In this paper, we study transition time estimates for a linear homogeneous diffusion equation with general time-independent boundary conditions. Two key contributions to the literature are presented, namely we:

1. (i)

demonstrate that, while both the mean action time and the mean plus one standard deviation of action time are useful at characterising the associated time-scale, they do not provide accurate estimates of the transition time for diffusion processes. 2. (ii)

propose a highly accurate alternative approach for estimating transition times based on using higher-order consecutive moments to approximate the large time asymptotic behaviour of the cumulative distribution function $F(t;x)$.

The second item is the main contribution of the work, with the novel approach providing a simple highly-accurate formula for calculating the local transition time (1) (and hence global transition time (2)) using only two consecutive moments $M_{k-1}(x)$ and $M_{k}(x)$ ($k\geq 1$). The result links, for the first time, the higher-order moments to calculation of finite transition times. This represents a significant breakthrough as previously it was believed that such calculation required the transient solution $u(x,t)$ [6, 24].

2 Diffusion model

The diffusion model considered in this paper is described below. Consider a diffusion process on the interval $\mathcal{L}:=[l_{0},l_{m}]$ governed by the linear homogeneous diffusion equation:

[TABLE]

In the equations listed above, $u(x,t)$ is the solution (e.g. temperature/concentration) at position $x$ and time $t$, $u_{0}(x)$ is the specified initial condition and $D>0$ is the constant diffusion coefficient. In the boundary conditions, $a_{L}$, $b_{L}$, $c_{L}$, $a_{R}$, $b_{R}$ and $c_{R}$ are constants satisfying: $a_{L}\geq 0$, $b_{L}\geq 0$, $a_{R}\geq 0$, $b_{R}\geq 0$, $a_{L}+b_{L}>0$ and $a_{R}+b_{R}>0$ [4, 2].

The corresponding steady-state solution of (7), denoted by $u_{\infty}(x)$, is the linear function satisfying the boundary value problem:

[TABLE]

For the Neumann problem ($u_{\infty}^{\prime}(l_{0})=u_{\infty}^{\prime}(l_{m})=0$), it is well-known [3] that the solution of (8) is unique only up to an additive constant and that including the additional constraint:

[TABLE]

gives the correct steady state solution of (7).

3 Computing the moments

In this section, we present an algorithm for calculating the first $q$ moments (6), namely $M_{k}(x)$ ($k=1,\ldots,q$), for the diffusion problem (7). The derivation follows closely the procedure taken by several authors [5, 7, 23] with the exception that we are interested in the raw moments as opposed to the central moments.

Combining equations (5) and (6) gives the following expression for the $k$th moment:

[TABLE]

for $k=0,1,\ldots$, where $h(x):=u_{\infty}(x)-u_{0}(x)$. Applying integration by parts and noting that $\lim\limits_{t\rightarrow\infty}t^{k}[u(x,t)-u_{\infty}(x)]=0$ [7] yields:

[TABLE]

For all integers $k>1$, each of the above moments can be calculated without requiring the transient solution $u(x,t)$ by deriving a boundary value problem satisfied by $M_{k}(x)$ as follows. Define:

[TABLE]

and consider the derivatives:

[TABLE]

Using equations (7a) and (8a) in equation (14) yields

[TABLE]

The above expressions lead to the following boundary-value problem for $\overline{M}_{k}(x)$:

[TABLE]

where the right-hand side of the differential equation (16a) is identified from (15) using $\overline{M}_{k-1}(x)=M_{k-1}(x)h(x)$ and equation (10) and the boundary conditions (16b) and (16c) are derived by utilising equations (12) and (13) and the boundary conditions for $u(x,t)$ and $u_{\infty}(x)$: (7c), (7d), (8b) and (8c).

In this paper, we consider only relatively simple initial condition functions $u_{0}(x)$, where the boundary value problem (16) can be solved analytically by integration. For more complicated initial conditions, a numerical method could easily be applied to solve (16). Integrating (16a) yields the general solution:

[TABLE]

where

[TABLE]

The constants $c_{k,1}$ and $c_{k,2}$ are determined by substituting the form of (17) into the boundary conditions (16b) and (16c) and solving the resulting linear system:

[TABLE]

where $\mathbf{c}=\left[c_{k,1},c_{k,2}\right]^{T}$ and

[TABLE]

with:

[TABLE]

We remark that the matrix $\mathbf{A}$ is the same as the one that appears when using a similar strategy to solve the boundary value problem (8) for the steady-state solution. For the Neumann problem ($\overline{M}_{k}^{\prime}(l_{0})=\overline{M}_{k}^{\prime}(l_{m})=0$), where the solution of the linear system (19) (and hence the solution of the boundary value problem (16)) is unique only up to an additive constant, we require an additional constraint on $\overline{M}_{k}(x)$ similar to equation (9). Note that equation (9) together with Neumann boundary conditions implies conservation, that is:

[TABLE]

for all $t\geq 0$. Integrating (12) from $x=l_{0}$ to $x=l_{m}$, reversing the order of integration in the resulting double integral and using (22) gives the required constraint on the solution of (16) for the Neumann problem:

[TABLE]

Incorporating this constraint into the linear system (19) yields the slightly modified form:

[TABLE]

where

[TABLE]

which has a unique solution provided that $G_{k}^{\prime}(l_{0})=G_{k}^{\prime}(l_{m})$. Once $c_{k,1}$ and $c_{k,2}$ are identified and $\overline{M}_{k}(x)$ is determined, the $k$th moment is computed as follows:

[TABLE]

Noting that $M_{0}(x)=1$ and hence $\overline{M}_{0}(x)=h(x)=u_{\infty}(x)-u_{0}(x)$, allows the $G_{k}(x)$ functions (18) and the moments to be calculated recursively, as outlined in Algorithm 3.

4 Transition time estimation

Using the moments, derived in the previous section, we now present three estimates of the local transition time (Definition 1), labelled $t_{s}^{(n)}(x)$ for $n=1,2,3$. For each estimate, the corresponding estimate of the global transition time is defined according to Definition 1 as:

[TABLE]

4.1 Low accuracy estimates using the first and second moments

4.1.1 Mean action time

The mean action time $\text{MAT}(x)$, which is often used as an estimate of the time required to reach steady state [22, 16, 24], defines the following estimate of the local transition time (Definition 1):

[TABLE]

Although $t_{s}^{(1)}(x)$ and the corresponding estimate of the global transition time $\widehat{t}_{s}^{(1)}$ do not depend on the tolerance $\delta$, it will be interesting to investigate how closely the transient solution $u(x,t)$, evaluated at $t=\widehat{t}_{s}^{(1)}$, is to the steady state solution.

4.1.2 Mean plus one standard deviation of action time

For most problems, it is probably unreasonable to expect that $\widehat{t}_{s}^{(1)}$ is an accurate estimate of the time required to reach steady state since its unlikely that the cumulative distribution function $F(t;x)$ evaluated at $t=\widehat{t}_{s}^{(1)}(x)$ is close to one111Or equivalently that $P(t\geq\widehat{t}_{s}^{(1)}(x))$ is small. for all $x\in\mathcal{R}$ since $t_{s}^{(1)}(x)$ is the mean of the distribution. Simpson et al. [24] suggest the mean action time plus one standard deviation of action time (square root of the variance of action time) as a way to improve estimation. Hence, we will also investigate the following estimate of the local transition time (Definition 1):

[TABLE]

which has been equivalently expressed in terms of the first and second moments. Analogously to $t_{s}^{(1)}$, equation (28) also does not depend on the tolerance $\delta$. We remark that the choice of one standard deviation is subjective and there is no reason why any number of standard deviations could not be used. However, we will not pursue this further as, unlike for the normal distribution, in general one cannot determine the probability that a random variable is greater than the sum of its mean and a given number of standard deviations.

4.2 High accuracy estimates using higher order moments

To estimate transition times, an accurate approximation of the cumulative distribution function $F(t;x)$ is required only for large $t$ (relative to the problem) where $u(x,t)\approx u_{\infty}(x)$. The exact solution of the diffusion problem (7) has the following functional form [26, 25]:

[TABLE]

where $\gamma_{n}$ and $\xi_{n}>0$ depend on the eigenvalues and eigenfunctions of the transient solution. Inserting (29) into (3) it follows that the cumulative distribution function $F(t;x)$ has the functional form:

[TABLE]

where $\zeta_{n}(x)=\gamma_{n}(x)/(u_{0}(x)-u_{\infty}(x))$. Assuming the $\xi$’s are arranged in ascending order (i.e., $\xi_{1}<\xi_{2}<\ldots$) and $\zeta_{1}(x)\neq 0$, then the cumulative distribution function satisfies the following asymptotic relation:

[TABLE]

The above analysis suggests an approximation to the large time behaviour of $F(t;x)$ should be sought in the form given below:

[TABLE]

where $\alpha(x)$ and $\beta(x)$ are as yet unspecified functions.

If we can devise a method for accurately computing values of $\alpha$ and $\beta$ that are close to $\zeta_{1}$ and $\xi_{1}$ then highly accurate transition time estimates can be calculated. The aim being, of course, to achieve this without explicit calculation of the transient solution $u(x,t)$ (i.e. without explicit computation of any eigenvalues or eigenfunctions) and using only the moments: $M_{k}(x)$ for $k=1,\ldots,q$. The probability density function corresponding to the cumulative distribution function (30) is given via differentiation:

[TABLE]

Using this form of $f(t;x)$ in equation (6), the $k$th moment is found to be:

[TABLE]

Since $\xi_{1}<\xi_{n}$ for all $n=2,3,\ldots$ it follows that $M_{k}(x)$ satisfies the asymptotic relation:

[TABLE]

This latter observation motivates the following pair of coupled equations satisfied by $\alpha(x)$ and $\beta(x)$, formulated by matching the $(k-1)$th and $k$th moments:

[TABLE]

Provided both $M_{k-1}(x)\neq 0$ and $M_{k}(x)\neq 0$, the above system of equations can be solved exactly to obtain the following explicit formulae:

[TABLE]

where we have included the subscript $k$ on $\alpha$ and $\beta$ to denote dependence on the $(k-1)$th and $k$th moments. We remark that $M_{k}(x)>0$ and hence $\alpha_{k}(x)>0$ and $\beta_{k}(x)>0$. The analysis above leads to the following approximation of the large time behaviour of the exact cumulative distribution function (30):

[TABLE]

which involves only the moments $M_{k-1}(x)$ and $M_{k}(x)$. Following definition (4), equating the right-hand side of (34) with $1-\delta$ and solving for $t$ yields the following local transition time estimate:

[TABLE]

where $\log$ is the natural (base $e$) logarithm. Note that $t_{s}^{(3)}(x)$ increases with decreasing tolerance $\delta$, which is consistent with the fact that it takes a longer amount of time for the transient solution to transition to within a smaller tolerance of its steady state. Moreover, we require $\delta\geq\alpha_{k}(x)$ for all $x\in\mathcal{R}$ to ensure the obvious physical constraint $t_{s}^{(3)}(x)\geq 0$ is satisfied. Inserting the expressions for $\alpha_{k}(x)$ and $\beta_{k}(x)$ into (35), we obtain the following simple formula for estimating the local transition time depending on the $(k-1)$th and $k$th moments:

[TABLE]

As a result of the asymptotic relation (33), we expect the accuracy of $t_{s}^{(3)}(x)$ to increase as $k$ increases, and this hypothesis is tested in the Section 5.

Assuming $\delta=10^{-p}$, where $p>0$, gives the following alternative form:

[TABLE]

which leads to the conclusion that for fixed $k$ and increasing $p$, the local transition time estimate $t_{s}^{(3)}(x)$ increases linearly with slope $\log(10)/\beta_{k}(x)$, e.g., the additional transition time required when decreasing $\delta$ from $10^{-(p-1)}$ to $10^{-p}$ is equal to the additional time when decreasing $\delta$ from $10^{-p}$ to $10^{-(p+1)}$.

Interestingly, $k=1$ gives $\alpha_{1}(x)=1$ and $\beta_{1}(x)=1/M_{1}(x)$ since $M_{0}(x)=1$. In this case, the approximation (34) simplifies to $F(t;x)\simeq 1-e^{-t/M_{1}(x)}$ for large $t$, which is nothing more than the cumulative distribution function of the exponential distribution with mean $M_{1}(x)$. For this simplest of cases, the estimate of the local transition time (36) reduces to the mean action time multiplied by a correction factor depending on the tolerance: $t_{s}^{(3)}(x):=M_{1}(x)\log(\delta^{-1})\equiv\text{MAT}(x)\log(\delta^{-1})$.

5 Results and Discussion

To investigate the accuracy of the transition time estimates presented and developed in sections 4.1 and 4.2, we consider the following three test cases:

Case A:

[TABLE]

Case B:

[TABLE]

Case C:

[TABLE]

where in each case $[l_{0},l_{m}]=[0,1]$. Together, the above problems test each of the three types of boundary conditions (Dirichlet, Neumann, Robin). Cases A and B resemble classical problems in heat conduction. For example, in Case A, a bar initially at temperature zero is suddenly heated at its left-boundary, how long does it take for the whole bar to heat up? On the other hand, Case C is a typical problem in mathematical biology [8], where a region initially fully occupied by cells ($u(x,0)=1$ for $0.25<x<0.75$) is left to diffuse.

Each of the test cases satisfy the requirement that $F(t;x)$ (3) defines a cumulative distribution function since, at each position $x$, the solution $u(x,t)$ is either non-decreasing for all $t>0$ if $u_{\infty}(x)>u_{0}(x)$ (as in Case A for all $0<x<1$ and Case C for $0<x<0.25$ and $0.75<x<1$) or non-increasing for all $t>0$ if $u_{\infty}(x)<u_{0}(x)$ (as in Case B for all $0<x<1$ and Case C for $0.25<x<0.75$). For Case C, if the region centered around $x=0$ that is fully occupied by cells, namely $0.25<x<0.75$, is shortened, however, the monotonicity property is violated.

Recall the three global transition time estimates presented in this paper, namely $\widehat{t}_{s}^{(n)}$ ($n=1,2,3$), the definitions of which are reiterated below:

• •

$\widehat{t}_{s}^{(1)}$ [Equations 27 and 26]

Based on using the mean or first moment of the probability distribution $f(t;x)$ (5) (or equivalently the mean action time) as an estimate of the local transition time, as described in section 4.1.1.

• •

$\widehat{t}_{s}^{(2)}$ [Equations 28 and 26]

Based on using the sum of the mean and standard deviation of the probability distribution $f(t;x)$ (5) (or equivalently the mean plus one standard deviation of action time) as an estimate of the local transition time, as described in section 4.1.2 and previously by Simpson et al. [24] and Jazaei et al. [14].

• •

$\widehat{t}_{s}^{(3)}$ [Equations 36 and 26]

Based on using the higher order moments of the probability distribution $f(t;x)$ (5) to approximate the large time asymptotic behaviour of the cumulative distribution function $F(t;x)$ (3) as described in section 4.2

For each global transition time estimate, the maximisation problem implied by equation (26) is solved by first noting that equation (26) is equivalent to $\widehat{t}_{s}^{(n)}=-\min_{x\in\mathcal{R}}[-t_{s}^{(n)}(x)]$ and then using MATLAB’s fminbnd function with option TolX = 1e-14 [18]. To calculate the moments we have implemented Algorithm 3 in MATLAB using the Symbolic Math Toolbox [20].

Recall that $\widehat{t}_{s}^{(3)}$ depends on two parameters that are free-to-choose: the prescribed tolerance $\delta$ and the index $k$ which specifies which two consecutive moments (i.e., $M_{k-1}(x)$ and $M_{k}(x)$) are utilised. Initially, we present results for $\delta=0.02$ and $k=2$, with the tolerance value chosen to ensure that the solution at the global transition time is visibly distinguishable, if only slightly, from the steady state $u_{\infty}(x)$. With $k=2$, the local transition time estimate $t_{s}^{(3)}(x)$ (36) simplifies to:

[TABLE]

which is comparable to $t_{s}^{(2)}(x)$ (28) in its simplicity and dependence on the first and second moments only.

Figure 1 plots the solution $u(x,t)$ of the diffusion problem (7) for each of the three test cases, depicting the transition from initial to steady state. In these plots, $u(x,t)$ is also given at each of the global transition times estimates, with the calculated values of $t=\widehat{t}_{s}^{(n)}$ for $n=1,2,3$ (rounded to four decimal places) tabulated in Figure 1. In addition to the visual comparison provided by Figure 1, in Table 1 we include errors:

[TABLE]

as a quantitive measure of how close the transient solution at each estimate is to the steady-state solution. The transient solution $u(x,\widehat{t}_{s}^{(n)})$ is evaluated by taking the first 50 terms in the classical eigenfunction expansion solution [25, 26], which is more than sufficient since the values of $\widehat{t}_{s}^{(n)}$ are relatively large. For $\widehat{t}_{s}^{(3)}$, we compute errors corresponding to three different tolerances $\delta=0.02,10^{-3},10^{-5}$. Noting the definition of the local transition time (Definition 1), $\varepsilon_{s}^{(n)}$ ideally should be close to the prescribed tolerance $\delta$.

The following observations can be drawn from Figure 1 and Table 1:

$\widehat{t}_{s}^{(1)}$ underestimates the time required to effectively reach steady state for all three test cases. The transition from initial to steady state is far from complete and this is confirmed by the large values of $\varepsilon_{s}^{(1)}$ in Table 1. 2. 2.

$\widehat{t}_{s}^{(2)}$ significantly improves on $\widehat{t}_{s}^{(1)}$ (as has been reported previously by Simpson et al. [24] and Jazaei et al. [14] for a groundwater modelling problem), however, in all three test cases it is clearly visible that the transition from initial to steady state is still not complete. 3. 3.

$\widehat{t}_{s}^{(3)}$ uses the the same moments as $\widehat{t}_{s}^{(2)}$ but produces a far superior estimate of the global transition time with the transient solutions at $t=\widehat{t}_{s}^{(3)}$ (see Figure 1) very close to steady state. For all three test cases, the accuracy of $\widehat{t}_{s}^{(3)}$ is quite remarkable. Using only the first and second moments, $\widehat{t}_{s}^{(3)}$ leads to errors $\varepsilon_{s}^{(3)}$ that are less than, and very close to, the prescribed tolerances of $\delta=0.02,10^{-3},10^{-5}$ (Table 1). Taking higher consecutive moments, that is increasing $k$, further improves accuracy since the asymptotic relation (33) is more accurate for larger $k$. This is demonstrated by choosing $k=5$, which gives $\varepsilon_{s}^{(3)}=0.0200=\delta$ (to the four decimal places displayed) for all three test cases.

In Figure 2, we plot the exact local transition time $t_{s}(x)$ (Definition 1) and local transition time estimate $t_{s}^{(3)}(x)$ (36) for all three test cases and four different choices of the tolerance $\delta$. To calculate $t_{s}(x)$ we solve equation (1) (with the first 50 terms used in the eigenfunction solution expansion for $u(x,t)$) using MATLAB’s fzero function with the default tolerance [19]. Recall that $t_{s}(x)$ provides the finite time required for the solution at position $x$ to transition to within a specified tolerance $\delta$ of the steady state solution at position $x$ (as defined in equation 1) and $t_{s}^{(3)}(x)$ is an estimate of that value. The local transition time profiles are therefore useful for determining the position(s) that take the shortest/longest time to reach steady state. For example, for Case C, the steady state is reached in the longest time at $x=0,0.5,1$ and the shortest time near $x=0.25$ and $x=0.75$, where the solution $u(x,t)=0.5=u_{\infty}(x)$ for all $t>0$. Both $t_{s}^{(1)}(x)$ and $t_{s}^{(2)}(x)$ are also useful in this regard as their shape follows closely that of $t_{s}(x)$: for both Case A and Case C, the value of $x$ that maximises/minimises $t_{s}^{(1)}(x)$ and $t_{s}^{(2)}(x)$ also maximises/minimises the exact local transition time $t_{s}(x)$. As observed previously, increasing $k$ leads to a better match with the exact transition time, with the curves almost indistinguishable for $k=10$. Finally, one must take care when using $t_{s}^{(3)}(x)$ for moderately large values of the tolerance $\delta$. As previously remarked in section 4.2, $t_{s}^{(3)}(x)$ is non-physical if $\alpha_{k}(x)<\delta$ and this behaviour is observed in Figure 2 for $\delta=10^{-1}$ near $x=0$ in Case A and near $x=1$ in Case B. This anomaly is tied to the inaccuracy of taking only the leading term in the transient solution $u(x,t)$ (which is effectively being approximated in equation (32)) for moderately small values of time $t$.

In Table 2, the new global transition time estimate $\widehat{t}_{s}^{(3)}$ is compared to the exact global transition time $\widehat{t}_{s}$ for several values of the tolerance $\delta$ and index $k$. Using only the moments and without explicit calculation of the transient solution $u(x,t)$, it is evident that $\widehat{t}_{s}^{(3)}$ is able to very accurately estimate $\widehat{t}_{s}$. Rough rule-of-thumb values are obtained for $k=1$, while for $k=2$, the relative errors are all between $10^{-1}$ and $10^{-3}$ indicating an accuracy of at least one and at most three significant digits. Increasing $k$ increases the accuracy (due to the accuracy of the asymptotic relation (33) improving) with at least 3 and and at most 6 significant digits obtained for $k=5$, which is probably sufficient for most applications. For $k=10$, one obtains the exact value to all four decimal places displayed with the relative errors indicating an accuracy of between 6 and 11 significant digits (inclusive).

To conclude this section, we provide some observations regarding the following generalised version of Case A, where the constant diffusivity $D$ is arbitrary and the boundary conditions are specified as:

[TABLE]

For this problem, computing the three local transition time estimates and evaluating them at $x=L$, where the maximum occurs, yields the following global transition time estimates:

[TABLE]

where $\gamma_{k}$ and $\theta_{k}$ are constants that depend on the chosen moment index $k$.

Note that all three estimates are proportional to the diffusive timescale $L^{2}/D$, differing only by a multiplicative factor, which for $\widehat{t}_{s}^{(3)}$ depends on the tolerance $\delta$. As a result, for two different diffusion processes, the question of which takes longer can be answered by choosing any of the global transition time estimates (39)–(41) and comparing its value for both processes.

Each of the global transition time estimates (39)–(41) were calculated in Maple222with the environment variable Digits set to 50 [17]., where exact fractional expressions for the constants $\gamma_{k}$ and $\theta_{k}$ appearing in equation (41) can be obtained. In Table 3, we give the corresponding values in decimal form and rounded to four decimal places to improve readability. For increasing values of $k$, observe in Table 3 that $\gamma_{k}$ and $\theta_{k}$ are approaching the values $4/\pi^{2}$ ($\approx 0.4053$) and $4/\pi$ ($\approx 1.2732$). Hence, $\widehat{t}^{(3)}$ is approaching:

[TABLE]

which is precisely the global transition time derived by using the leading term of the transient solution, namely

[TABLE]

in the definition of the local transition time (Definition 1), solving equation (1) for $t$ and evaluating the result at $x=L$, where the maximum occurs.

6 Conclusions

In summary, we have derived a simple formula, denoted by $t_{s}^{(3)}(x)$ and given in equation (36), for calculating a finite measure of the time required for a diffusion process to reach steady state. This formula estimates the local transition time defined as the time required for the transient solution to transition to within a prescribed tolerance $0<\delta\ll 1$ of its steady state, at position $x$. A finite measure of the time required for the entire diffusion process to reach steady state is then obtained by evaluating $t_{s}^{(3)}(x)$ at the value of $x$ that produces a maximum. This novel formula is attractive as it (i) avoids explicit calculation of the transient solution (ii) depends only on the prescribed tolerance $\delta$ and the $(k-1)$th and $k$th moments and (iii) can be used to calculate the exact transition time to effectively any number of significant digits by increasing $k$. Our results confirm that even for $k=2$, for which $t_{s}^{(3)}(x)$ utilises only the first and second moments (or, equivalently, the mean and variance of action time), the accuracy is quite remarkable. In all cases, the new approach comprehensively outperforms existing strategies based on the mean action time and variance of action time, which while useful at characterising the associated time-scale, significantly underestimate transition times for diffusion processes.

This paper has presented a proof-of-concept for the simple problem of homogeneous diffusion. Future work will focus on extending the new method to heterogeneous and higher-dimensional problems as well as other transport processes (e.g. advection-diffusion).

Acknowledgments

This research was funded by the Australian Research Council (ARC) via the Discovery Early Career Researcher Award DE150101137.