Optimal finite-differences discretization for the diffusion equation from the perspective of large-deviation theory

When applying the finite-differences method to numerically solve the one-dimensional diffusion equation, one must choose discretization steps Δx, Δt in space and time, respectively. By applying large-deviation theory on the discretized dynamics, we analyze the numerical errors due to the discretization, and find that the (relative) errors are especially large in regions of space where the concentration of particles is very small. We find that the choice Δt=Δx2/(6D) , where D is the diffusion coefficient, gives optimal accuracy compared to any other choice (including, in particular, the limit Δt→0 ), thus reproducing the known result that may be obtained using truncation error analysis. In addition, we give quantitative estimates for the dynamical lengthscale that describes the size of the spatial region in which the numerical solution is accurate, and study its dependence on the discretization parameters. We then turn to study the advection–diffusion equation, and obtain explicit expressions for the optimal Δt and other parameters of the finite-differences scheme, in terms of Δx, D and the advection velocity. We apply these results to study large deviations of the area swept by a diffusing particle in one dimension, trapped by an external potential ∼|x| . We extend our analysis to higher dimensions by combining our results from the one dimensional case with the locally one-dimension method.


I. INTRODUCTION
The diffusion equation and the advection-diffusion equation (in d spatial dimensions) are of central importance in many fields, including hydrodynamics, statistical physics, biology, engineering and many more.In these equations, u (x, t) is the concentration of some substance, D > 0 is the diffusion coefficient and in Eq. ( 2), F is the velocity due to advection.In simple cases these equations can be solved analytically, but in general they can be difficult to solve, and as a result many numerical methods have been developed [1][2][3][4][5][6][7][8][9][10][11][12][13][14].One of the simplest methods involves a discretization of time and space, followed by an approximation of the temporal and spatial derivatives by using the method of finite differences.One is then faced with the problem of choosing the discretization steps ∆x, ∆t in space and time, respectively.It is well known that if this choice is not done carefully, the numerical solution may be unstable even if ∆x and ∆t are both very small (see below).
For simplicity, let us begin by discussing the diffusion equation in d = 1, We choose a rectangular grid in the xt plane, defined by the points ) t j = j∆t, , j = 0, 1, 2, . . ., (5) and denote the concentration at position x i and time t j by u j i .Approximating the derivatives by using the method of finite differences, one obtains which can be rewritten as where p = D∆t/∆x 2 .In principle, one must also deal with the boundary conditions.However, for now they will not play an important role in our analysis so we will assume for simplicity that the system size is infinite.What is the best way to choose the space-time discretization?Generally speaking, as ∆x and ∆t tend to zero, the accuracy of the solution improves but the computational cost increases.However, it is well known that, for the numerical solution to be stable, one must choose p ≤ 1/2, see e.g.[4].Phrasing the question slightly differently, for a given spatial discretization ∆x, how should one choose p (or equivalently ∆t) to obtain good efficiency and accuracy?One may naively expect that, as p decreases, the accuracy of the solution improves but the computational cost increases.However, it turns out that (rather surprisingly and counterintuitively) the best accuracy is obtained for p = 1/6 and, in particular, if p is decreased below this value then both the efficiency and the accuracy of the numerical solution worsen.This result is well-known and is easy to obtain using error-truncation analysis, see e.g.Ref. [2] and references therein.For completeness we reproduce this analysis in Appendix A.
In this paper, we revisit this problem from the point of view of large-deviation theory.We first note that the relative numerical errors are especially large in regions of space in which the concentration of particles u is very small.We characterize the region of space in which the relative numerical errors are small through a dynamical lengthscale ℓ(t).We find that, for generic choices of 0 < p ≤ 1/2, this lengthscale behaves (at t ≫ ∆t) as but for the specific choice p = 1/6, an improvement that becomes especially significant at long times.The remainder of the paper is organized as follows.In section II we analyze the diffusion equation in d = 1 in detail by applying tools from large-deviation theory.In section III we perform an analogous analysis for the one-dimensional advection-diffusion equation, which we then apply in section IV to numerically study large deviations of the area swept by a Brownian particle trapped in an external potential ∼ |x|.In section V we extend our analysis to d > 1.Finally, we briefly summarize and discuss our results in section VI.

II. DIFFUSION EQUATION IN d = 1
Our starting point is the diffusion equation (3).In what follows, it is very useful to interpret this as a Fokker-Planck equation for the time-dependent probability density function u(x, t) of the position of a single, diffusing particle (so that u (x, t) dx is the probability that, at time t, the position of the particle is between x and x + dx) [15,16].One can then interpret the discretized equation (7) as describing the temporal evolution of the probability vector u j i of the position of a particle undergoing a random walk in discrete space and time.Indeed, interpreting u j i as the probability that, at time t j the position of the particle is x i , the dynamics (7) describe a particle which, at every timestep (of duration ∆t), either jumps a distance of ∆x to the left or to the right, each with probability p, or does not move, with the complementary probability 1 − 2p.Incidentally, within this interpretation the requirement p ≤ 1/2 corresponds to the probability of not moving being nonnegative.An equivalent way to represent the dynamics of the random walk would be to write them as where X j is a random variable denoting the position of the particle at time t j , and ξ 1 , ξ 2 , . . .are i.i.d.random variables each taking the values −1, 0, 1 with probabilities p, 1 − 2p, p respectively.How well does the discrete random walk approximate the (continuous) diffusion process?For simplicity, let us assume that the particle starts from position x = 0 at time t = 0, i.e., the initial condition for Eqs. ( 3) and ( 7) are given by Dirac-u (x, t = 0) = δ (x) and Kronecker-u 0 i = δ i,0 delta functions, respectively.Moreover, we will assume for now that the system is infinite in space, i.e., −∞ < x < ∞ in (3) and i ∈ Z in (7).In this case, the solution to (3) is simply given by the Green's function of the diffusion equation, to which we would now like to compare the solution to the discretized equation (7).One can gain insight by using the representation (10) of the random walk, which we analyze here with the initial condition X 0 = 0.In this representation it is clear that X j = ∆x j−1 j ′ =0 ξ j ′ , and since the ξ j ′ are i.i.d., X j converges to a normal distribution at large j.Moreover, since the mean and variance of ξ j 's are given by ξ j = 0 (where the angular brackets denote ensemble averaging) and Var ξ j = 2p respectively, it immediately follows that X j = 0 and Var X j = 2pj∆x 2 = 2Dt j .Thus, at j ≫ 1, typical fluctuations of X j are well approximated by a Gaussian distribution (due to the central limit theorem) in agreement with Eq. (11).
We are now interested in analyzing the error in the approximate result (12).This can be done by considering the higher cumulants of the distribution of X j : For an (exactly) Gaussian distribution, all cumulants higher than the second cumulant vanish.At long times, the nonvanishing higher cumulants of P (X j ) mostly affect the tails of the distribution.Indeed, at large j, the large-fluctuations regime of this distribution follows the scaling behavior [17][18][19][20][21] P X j ∼ e −jI(X j /j∆x) . ( This scaling behavior is known as a large-deviation principle, and it describes an exponential decay in time with a "rate function" I(a).The rate function vanishes at its global minimum which is at a = 0, and its asymptotic behavior at |a| ≪ 1 is I (a) ≃ a 2 /4p, thus providing a smooth matching with the typical-fluctuations Gaussian regime (12).I(a) is closely related to the higher cumulants of the distribution: It is given by the Legendre-Fenchel transform of the scaled cumulant generating function (SCGF) λ (k) = lim j→∞ 1 j ln e kX j /∆x [20].Due to the mirror symmetry of P (X j ), the leading-order correction to the quadratic approximation of I(a) around its minimum is generically of order O(a 4 ).
λ(k) can be easily found explicitly: Writing this as a power series in k, we find The power series expansion of λ(k), taken up to quadratic order, describes the mean and variance of the Gaussian distribution of X j at j ≫ 1, in agreement with the corresponding quadratic asymptotic behavior of I(a) around its minimum (at a = 0).The terms in the series ( 16) that are of order O(k n ) with n > 2 represent corrections to this Gaussian behavior [in particular, the O(k 4 ) term in (16) gives rise to the O(a 4 ) term in the expansion of I(a)], and their importance decreases as n grows.One immediately notices that the O(k 4 ) term in Eq. ( 16) term vanishes at p = 1/6.I.e., for the choice p = 1/6, the fourth cumulant of X j exactly vanishes and the subleading corrections to the leading-order Gaussian behavior come from the sixth cumulant and higher.Thus, from the point of view of approximating the diffusion process by the random walk, the choice p = 1/6 is expected to provide the best accuracy.We tested this prediction by solving Eq. ( 7) with different values of p = 0.4, 1/6, 0.1, 0.01, see Fig. 1.Indeed, p = 1/6 provides a level of accuracy that is far better than all other choices.The differences are especially pronounced in the tails of the distribution, which are overestimated (underestimated) for p < 1/6 (p > 1/6) by orders of magnitude.For example, for the choice p = 1/6, u (x = 500, t = 400) differs by just 4% from its theoretical value G (x = 500, t = 400) ∼ 10 −70 , whereas for the choices p = 0.4, 0.1, 0.01 the ratio u/G| x=500,t=400 takes the values 3.6 × 10 −4 , 6.5 and 71, respectively.The ratio u/G between the numerical solutions and the exact solution.It is clearly seen in the figure that for p = 1/6, the agreement is significantly better than any other p.
The large-deviation principle (13) enables us to estimate the range of X j 's beyond which the numerical solution is expected to develop very large errors with respect to the exact solution, for different choices of the parameter p. Indeed, for p ̸ = 1/6 the O(a 4 ) term in the power-series expansion of I(a) causes the numerical solution to lose accuracy at X j ∼ ∆xj 3/4 ∼ (t j D) 3/4 / (∆x) 1/2 , corresponding to Eq. ( 8) above.In contrast, for p = 1/6 the O(a 4 ) term is absent, and instead the leading-order correction to the parabolic behavior comes from an O(a 6 ) term.This causes the numerical solution to be accurate in a much wider range, losing accuracy only at X j ∼ ∆xj 5/6 ∼ (t j D) 5/6 / (∆x) 2/3 , corresponding to Eq. ( 9) above.The differences between these lengthscales is clearly seen in Fig. 1.

III. ADVECTION-DIFFUSION EQUATION IN d = 1
Let us now consider the case in which a (constant) advection term is added, i.e., the equation of motion reads A standard approach for discretizing time and space would be to approximate the terms ∂ t u and ∂ 2 x u as described above, and to approximate the advection term using the finite-differences method by ∂ x u ≃ u j i+1 − u j i−1 / (2∆x).If this is done, Eq. ( 7) gives way to where Once again one is faced with the problem of choosing the discretization parameters, which we shall phrase in the form of choosing ∆t for a given ∆x.
We will follow the same general approach that we used above for the advection-free case.We interpret Eq. ( 18) as the evolution of a random walker in discrete time and space, that at every timestep ∆t, jumps a distance of ∆x to the right (left) with probability p (q), or stays at the same position with the complementary probability 1 − p − q.Note that, clearly from the definition of q in (19), one must choose ∆x ≤ 2D/F , otherwise q is negative.It is again convenient to use the description (10) of the random walker's dynamics, but now the ξ j 's take the values −1, 0, 1 with probabilities q, 1 − p − q, p respectively.
The exact solution to the advection-diffusion equation ( 17) with initial condition u(x, We now analyze the error of the numerical solution, which again is expected to be most pronounced in the tails of the distribution P (X j ).Let us begin by obtaining the central limit theorem's prediction, which should give P (X j ) in the typical-fluctuations regime.The mean and variance of each of the ξ j 's are respectively, so it immediately follows that the mean and variance of X j are given by and respectively.The mean (22) is in agreement with the exact result (20), but the variance (23) corresponds to an effective diffusion coefficient that is slightly modified compared to the exact one, D eff = D − F 2 ∆t/2.Although the difference is not large (and vanishes as ∆t → 0), its effect on the tails of the distribution P (X j ) can be nonnegligible.We therefore propose, as a first improvement over this naive method, to modify the choices of p and q so that the mean and variance of X j attain their desired (exact) values.From the first parts of Eqs. ( 22) and ( 23), one finds that the requirements are The solution to these equations is which replace the naive choices given by Eq. ( 19).In Fig. 2, we show that, for a particular set of parameters (D, F, ∆x), this choice of p and q indeed improves the accuracy of the solution, although the effect is not very large.We found similar improvements for other choices of parameters (D, F, ∆x) (not shown).
Having chosen p and q such that the mean and variance of X j are correct, one can now analyze the error of the numerical solution as encoded in the higher-order cumulants of P (X j ).As in the advection-free case, these are again related to the tails of the distribution as described by Eq. ( 13), but with a rate function I(a) whose features are slightly modified, due to the advection, which breaks the mirror symmetry of the distribution P (X j ).As a result, I(a) will also not be mirror symmetric around a = a * .So first of all, the minimum of I(a) is at the nonzero value a = a * = F , with an asymptotic behavior that is again parabolic I (a) ∼ (a − a * ) 2 around a ≃ a * .Moreover, there 500 0 500 1000 1500 2000 (a) Numerical solutions u(x, t) of the advection-diffusion equation obtained using the finite-differences method (18), using the naive p and q from Eq. ( 19) (solid line), and using the improved values ( 26) and ( 27) (dot-dashed lines) that are adjusted so that the variance X j equals its theoretical value, together with the exact solution GF from Eq. ( 20).Parameters are D = ∆x = 1, F = 1/2, t = 2000 and ∆t = 0.1.(b) The ratio u/GF between the numerical solutions and the exact solution.
As seen in the figure, the accuracy is somewhat improved when using the variance-adjusted p and q.
will in general be a cubic term in this asymptotic expansion around a = a * , reflecting the nonvanishing third cumulant of the distribution P (X j ) [To remind the reader, the third cumulant for the exact (Gaussian) solution vanishes].
Nevertheless, as we now show, one can choose ∆t such that this third cumulant exactly vanishes, thus significantly reducing the error of the numerical solution.The SCGF is again easy to calculate: and its power-series expansion around k = 0 is From here, we read off the third cumulant Now, recalling the choices ( 26) and ( 27) of p and q (repsectively) in terms of ∆t, we can determine the optimal choice of ∆t by requiring ξ j 3 c to vanish, which (using that F ̸ = 0 so clearly p ̸ = q) leads to the equation whose solution is This particular choice is expected to give the best accuracy, because for this choice the third cumulant of P (X j ) vanishes so the error comes from higher cumulants.As in the case of the diffusion equation (with no advection), our choices for the parameters p, q and ∆t coincide with the ones that may be obtained by truncation error analysis, see e.g.Ref. [1].Besides reproducing the optimal values of these parameters, our analysis provides additional insight by interpreting the finite-differences method as the dynamics of a random walk in discrete space and time.Moreover, we point out that the numerical errors can be enormous (many orders of magnitude) in the large-deviation regimes, and find a dynamical lengthscale that describes the size of the region in space in which the relative numerical error is small.
We tested the quality of our choice of parameters by solving Eq. ( 18) with different values of ∆t, see Fig. 3. Indeed, ∆t = ∆t * provides a level of accuracy that is far better than all others.One can see that the skewness changes sign as ∆t crosses the optimal value ∆t * .Using a similar argument to the one given in the previous section, one finds here that for ∆t ̸ = ∆t * , the solution is only expected to be accurate (in the sense that the relative error is much smaller than 1) within a spatial regime around x = F t of typical size ℓ(t) ∼ t 2/3 (33) whereas for ∆t = ∆t * , the relevant lengthscale is  26) and ( 27) respectively.The exact, analytic solution GF (x, t) from Eq. ( 20) is also plotted.(b) The ratio u/GF between the numerical solutions and the exact solution.It is clearly seen in the figure that for ∆t = ∆t * , the agreement is significantly better than any other ∆t.
However, note that if the spatial discretization is very fine, ∆x ≪ D/F , then the optimal time step (32) becomes ∆t * ≃ ∆x 2 / (6D), coinciding with that of the advection-free case.Thus, at ∆x ≪ D/F , choosing ∆t = ∆t * is optimal, because for this choice the third and the fourth cumulants of P (X j ) both (nearly) vanish, so the lowest nonvanishing cumulant will be the fifth.This leads to a further enhancement of the accuracy, encompassing a spatial regime whose typical size grows as This appears to be the case in Fig. 3 for the moderately small value ∆x = D/ (3F ).This analysis also suggests that, regardless of the force F , the choice ∆t = ∆x 2 / (6D) should yield good performance in terms of efficiency and accuracy, as long as ∆x ≪ D/F (because the difference between this choice and the choice ∆t = ∆t * is very small).One may thus expect this choice ∆t = ∆x 2 / (6D) to perform quite well even if F is not constant in space and/or time.

IV. APPLICATION: AREA SWEPT BY A BROWNIAN PARTICLE IN A POTENTIAL ∼ |x|
As an application, we consider a model of Brownian motion in an external potential V (x) = µ|x|, as described by the Langevin equation [22][23][24] where µ > 0 and −∞ < x < ∞.The particle starts at the origin, x(t = 0) = 0. We study the full distribution P (A; t) of the area swept by the particle up to time t.This problem is studied theoretically in Ref. [25], and our goal here is to study it numerically.Using that dA/dt = x, the Fokker-Planck equation for the joint distribution P (x, A; t) of x and A at time t is given by Eq. ( 38) is to be solved for −∞ < x < ∞ and −∞ < A < ∞ with the initial condition P (x, A; t = 0) = δ (x) δ (A).
Once P (x, A; t) is found, one can marginalize over x to obtain the area distribution: In Ref. [25], a variety of theoretical tools from large-deviation theory are applied to obtain the full distribution P (A; t) (including its tails) approximately, in the long-time limit t → ∞.In the present work, we obtain P (A; t) numerically by solving Eq. ( 38) using the finite-differences method.The numerical solution is computationally heavy, especially due to the fact that we are interested in obtaining P at long times and large A (and, as a result, rather large x).As a result, one has no choice but to use relatively large values of the discretization parameters ∆x, ∆A, ∆t, and it becomes of great importance to retain the best level of accuracy possible.
Let us now describe our numerical method.We rely on the results of Section III, while noticing that in each of the regimes x > 0 and x < 0, Eq. ( 36) describes a diffusing particle with a constant drift term (F = −µ and F = µ respectively).Let us describe our discretization method in terms of its corresponding stochastic dynamics that are discrete in x, t and A. As above, we consider a random walker X (t) that, at every timestep ∆t, jumps a distance ∆x to the right (left) with probability p (q), or, with probability 1 − p − q, stays in the same position.The optimal choices of p, q and ∆t are given by Eqs. ( 26), ( 27) and (32), respectively, where one must use F = µ (F = −µ) if the position of the particle is to the right (left) of the origin2 .In addition, we consider a discrete-time process A(t) that, at every timestep ∆t, jumps a distance of X (t).Note that A(t) takes discrete values, that are all given by integer multples of ∆A ≡ ∆X∆t.The dynamics of the joint distribution of X = i∆X and A = j∆A at time t = k∆t are described by a discrete version of Eq. ( 38) where i, j ∈ Z, k = 0, 1, 2, . . ., and p and q depend on the sign of i as described above.The three terms on the right-hand side of Eq. (40) correspond to the three possible positions from which the random walker may reach the point X = i∆X.
We implemented this solution with D = 1 and µ = 0.693.Choosing ∆x = 1, Eqs. ( 26), ( 27) and (32) yield ∆t ≃ 0.1645, p ≃ 0.114 and q ≃ 0.228 at x > 0 (at x < 0 the values of p and q are exchanged).In order to probe the large-deviation regime, we used fairly large cutoff values, X max = 80∆x and A max = 10 5 ∆A, and we were able to run the solution up to time t = 5 × 10 4 ∆t ≃ 8225 within reasonable computational time.The resulting P (A; t) is plotted in [25], showing good agreement with the theoretical predictions obtained there.
To test the dependence of the numerical accuracy on the chosen parameters, we also performed shorter runs, up to time t = 200, with different choices of the time discretization, ∆t = 0.4, 0.1645, 0.1, 0.01.We compared them with a more accurate numerical solution which we obtained by using a finer spatial discretization, ∆x = 0.5, and the corresponding optimal time discretization ∆t = 0.0415 from (32).In all cases (i.e., for each choice of ∆x and ∆t), we used the values of p and q from Eqs. (26), (27).The resulting P (A; t) are plotted in Fig. 4. Using the solution with ∆x = 0.5, ∆t = 0.0415 as a reference, we find that of all the solutions with ∆x = 1, the choice ∆t = 0.1645 yields the best accuracy, whereas the other choices yield an error that, in the large-deviation regime, is of order a factor of a few (this includes the choices ∆t < ∆t * which are computationally more expensive).

P(A)
x = 1, t = 0.4000 x = 1, t = 0.1645 x = 1, t = 0.1000 x = 1, t = 0.0100 x = 0.5, t = 0.0415 FIG. 4. Distributions P (A; t) of the area (37) swept by a Brownian particle in an external potential µ|x|, obtained from numerical solutions of the joint Fokker-Planck equation (38) for the particle's position x(t) and the area A(t).Parameters are D = 1 and µ = 0.693.The solid lines correspond to a spatial discretization of ∆x = 1 and different choices of the temporal discretization, ∆t = 0.4, ∆t * , 0.1, 0.01 where ∆t * = 0.1645.We find that the choice ∆t = ∆t * yields the best accuracy, by comparing it with the dashed line, which corresponds to a numerical solution with a finer spatial discretization ∆x = 0.5 and the corresponding optimal temporal discretization ∆t = 0.0415.

V. HIGHER DIMENSION
In d > 1, we can follow a similar approach as we did in d = 1.Namely, we interpret the dynamics obtained from the finite-difference approximation as the evolution of the position distribution of a random walker in discrete space and time, and choose the parameters of the numerical solution such that the lowest cumulants coincide with those of Brownian motion (for simplicity, we consider a delta function initial position distribution and an infinite system).
Let us begin by considering the diffusion equation ∂ t u = D∇ 2 u in d = 2.A standard way to apply the finitedifferences method is to approximate the Laplace operator by where ∆ℓ is the resolution of the spatial grid.Using this in the diffusion equation, and approximating the time derivative by ∂ t u ≃ (u| t+∆t − u| t ) /∆t, one obtains, in analogy to Eq. ( 7), an equation that can be interpreted as describing the dynamics of a discrete random walker, that jumps with probability p = D∆t/∆x 2 a distance of ∆ℓ to each one of the four directions ±x, ±ŷ, or does not move with probability 1 − 4p.One may now try to choose p such that the numerical error due to discretization is minimized.As in the case d = 1, one can choose p = 1/6 so that the fourth cumulants of the x− and y− coordinates of the random walker vanish.However, it is not difficult to see that the x and y coordinates of such a particle are statistically dependent.Thus, the mixed correlation functions such as x 2 y 2 − x 2 y 2 will not vanish, leading to nonvanishing mixed cumulants.This problem is easily remedied if one uses instead the so-called locally one-dimension (LOD) method, in which each time step consists of two substeps, one in the x direction and then one in the y direction (these substeps correspond to the two parts ∂ 2 x u and ∂ 2 y u that together form the Laplacian).The interpretation of these dynamics is as those of a random walker that, at each time step of duration ∆t, first jumps to the left or right each with probability p (or does not move, with probability 1 − 2p) and then jumps up or down each with probability p (or does not move, with probability 1 − 2p).Clearly, the x and y coordinates of such a random walker are exactly statistically independent.One then follows the same analysis as for the d = 1 case and chooses p = 1/6, so that the fourth cumulants vanish.
In order to test this method, we apply it to the diffusion equation in d = 2 equation with a delta-function initial condition, whose analytic solution is where r 2 = x 2 + y 2 .Note in particular that the x-and y-coordinates are statistically independent, because u(x, y, t) may be decomposed as the product of the one-dimensional Green's functions (11) acting on each coordinate separately.As shown in Fig. 5, the LOD method with p = 1/6 is indeed far more accurate than all alternatives.
We used the quantitative tools from large-deviation theory to analyze the accuracy of the numerical solution to the diffusion and advection-diffusion equations with the finite-differences method, and fine-tune the parameters of this method to optimize its accuracy.It would be interesting to see if this approach can be extended to optimize other numerical schemes (for instance, implicit methods [2]) and/or to the numerical solution of related equations such as reaction-advection-diffusion equations [5].In particular, one could consider the case in which the forcing F and/or diffusion D are time-or space-dependent, or even depend on the concentration u itself.It appears unlikely that it would be possible to obtain a method that performs as well as in the relatively simple case of constant F and D considered here.However, in the particular case where D is constant, and the spatial discretization ∆x is chosen to be sufficiently small, one might still expect that choosing ∆t = ∆x 2 / (6D) would result in good performance in terms of efficiency and accuracy, as explained shortly after Eq. ( 35).This analysis could prove useful for studying, e.g., large-deviation statistics in the Ornstein-Uhlenbeck process [26][27][28], in which there are theoretical predictions which still await numerical verification.
Finally, numerical investigation of large-deviation statistics is a notoriously difficult task, due to the extremely low probabilities of the events involved [28][29][30][31][32][33][34][35][36][37][38][39][40][41][42].We pointed out that an optimal choice of parameters for the finitedifference discretization of the advection-diffusion equation can have a huge effect (many orders of magnitude) in the regions of space and time for which the concentration u is small.These regions are often those that are of crucial importance when numerically analyzing large-deviation statistics.It would be useful to apply our results to numerically study large-deviation properties of diffusing particles with or without drifts (either when solving diffusion equations numerically [25], or in Monte-Carlo simulations of Brownian motion [33]).The example given in the present work (area swept by a trapped Brownian particle) together with Ref. [25] is a first step in this direction.

FIG. 1 .
FIG. 1.(a) Numerical solutions u(x, t) of the diffusion equation obtained using the finite-differences method, with ∆t = p∆x 2 /D.Parameters are D = ∆x = 1, t = 400, and the different curves correspond to different values of p as indicated in the figure label.The exact, analytic solution G (x, t) (11) is also plotted.(b)The ratio u/G between the numerical solutions and the exact solution.It is clearly seen in the figure that for p = 1/6, the agreement is significantly better than any other p.
FIG. 2. (a)Numerical solutions u(x, t) of the advection-diffusion equation obtained using the finite-differences method (18), using the naive p and q from Eq. (19) (solid line), and using the improved values (26) and (27) (dot-dashed lines) that are adjusted so that the variance X j equals its theoretical value, together with the exact solution GF from Eq.(20).Parameters are D = ∆x = 1, F = 1/2, t = 2000 and ∆t = 0.1.(b) The ratio u/GF between the numerical solutions and the exact solution.As seen in the figure, the accuracy is somewhat improved when using the variance-adjusted p and q.

FIG. 3 .
FIG. 3. (a)Numerical solutions u(x, t) of the advection-diffusion equation obtained using the finite-differences method.Parameters are D = ∆x = 1, F = 1/3, t = 300, for which Eq.(32) gives ∆t * = 0.16615 . . . .The different curves correspond to different choices of ∆t, as indicated in the figure label, from which the parameters p and q are determined via Eqs.(26) and (27) respectively.The exact, analytic solution GF (x, t) from Eq. (20) is also plotted.(b) The ratio u/GF between the numerical solutions and the exact solution.It is clearly seen in the figure that for ∆t = ∆t * , the agreement is significantly better than any other ∆t.