A positivity‐preserving and energy stable scheme for a quantum diffusion equation

We propose a new fully‐discretized finite difference scheme for a quantum diffusion equation, in both one and two dimensions. This is the first fully‐discretized scheme with proven positivity‐preserving and energy stable properties using only standard finite difference discretization. The difficulty in proving the positivity‐preserving property lies in the lack of a maximum principle for fourth order partial differential equations. To overcome this difficulty, we reformulate the scheme as an optimization problem based on a variational structure and use the singular nature of the energy functional near the boundary values to exclude the possibility of non‐positive solutions. The scheme is also shown to be mass conservative and consistent.

for higher-order diffusion equations using a standard finite difference discretization, yet with proven positivity of numerical solutions.
Here we focus on the quantum diffusion equation on the domain Ω ⊂ R . To make the presentation simple, we consider Ω = T ( ≥ 1 be any positive integer) and periodic boundary conditions: where the energy functional F ≔ F(u) is given by Note that this functional is also called Fisher information in the literature. A rigorous justification of the gradient flow structure of this equation was given in Reference [13] using the Wasserstein distance W 2 2 (u 0 , u 1 ) = inf ∈Π(u 0 ,u 1 ) ∫ T ×T |x − y| 2 d (x, y), (1.4) where (x, y) is the joint probability measure of x and y, Π(u 0 , u 1 ) is the set of all probability measures on T × T with first marginal u 0 and second marginal u 1 , and the symbol | ⋅ | denotes the usual Euclidean norm in T . Positivity of numerical solutions is highly desirable for any numerical method to solve Equation (1.1). One way to obtain positivity is to introduce auxiliary variables and enforce the solution positivity accordingly. For example, one can introduce v = log u [6,19] or v = √ u [4] to ensure a positive solution. Another way to obtain positivity is to use the Wasserstein gradient flow structure. In the one dimensional case, the Wasserstein distance (1.4) corresponds to the L 2 distance between Lagrangian maps. Positivity of u then follows from the monotonicity of the Lagrangian map [11,28,29]. One can also utilize the Eulerian formulation of the Wasserstein gradient flow structure. More precisely, by the Benamou-Brenier formulation of the Wasserstein distance [2], which is defined by the scheme proposed in Reference [22] for (1.1) takes the form Here Δt is the time step and u n is the solution at the nth time step. Positivity of u n+1 was proved since the objective function in the above optimization problem becomes infinite when u touches zero [22]. Some other numerical approximations such as finite volume methods were also developed to solve this equation, see, for example, Reference [26]. However, the existing positivity-preserving numerical methods for Equation (1.1) do not seem to follow simply from a direct finite difference discretization. In Reference [28], the authors even argued that there is no reason to expect the positivity-preserving property from a standard discretization approach.

Our contributions
The main objectives of this paper are to present a novel numerical scheme to approximate (1.1) using a standard finite difference discretization and prove how positivity of numerical solutions can be ensured by such a discretization. We begin with the semi-discretized scheme u n+1 − u n Δt = ∇ ⋅ (u n ∇H n+1 ), H n+1 = u F(u n+1 ), (1.6) where H n+1 is obtained by the variational derivative of energy functional F = F(u). This scheme is formally equivalent to the optimization problem Here with  u n ≔ −∇ ⋅ (u n ∇ ) we define The above optimization problem can be rewritten as The optimization problem (1.8) is similar in structure to the scheme (1.5). Note that the optimal conditions for this optimization problem lead to precisely the scheme (1.6), while the optimal conditions for (1.5) link only to a first order approximation of (1.6) or (2.3), see Reference [22]. The singular nature of F(u) at u = 0 may prevent the minimizer of (1.8) from reaching zero, so it is expected that solution positivity can be deduced at least after some spatial discretization. In this work we will prove solution positivity for the fully-discretized numerical scheme of form Hereû is defined to be an average value of u on the adjacent grid points given by (3.3) and (5.3), and h , D h are finite difference operators, given by (3.2) and (5.1) and (5.2), both for one and two dimensions, respectively. The function F h is a discretized version of the functional (1.3), given by (3.9) and (5.8) for the one and two dimensions, respectively. Our approach in proving the positivity-preserving property using optimization formulations was motivated by the works [7,10], where the authors used the finite difference discretization to study the Cahn-Hillard equation with the energy functional E given by, for example in Reference [7]. Our problem here differs from the Cahn-Hillard equation in that the functional includes a higher order term which develops singularity at zero, and the possibility of ∇u and u being zero simultaneously also poses additional challenges.
The main results of this paper include a new fully discretized finite difference scheme (1.9) to solve (1.1), and rigorous proofs of scheme properties such as positivity-preserving, mass conservation, energy stability, and scheme consistency. One most remarkable contribution of this work is to show that standard finite difference discretization of the quantum diffusion Equation (1.1) can meet the requirement of both solution positivity and energy dissipation simultaneously.
The paper is organized as follows. In the next section we present four semi-discrete schemes with only time discretization for Equation (1.1) and briefly discuss their properties. Section 3 is devoted to the fully discretized scheme in the one dimensional case, where the scheme is shown to be consistent. The optimization formulation is given in Section 4, where the scheme (1.9) is shown to be positivity-preserving, energy dissipative and mass conservative. In Section 5, we present and analyze the scheme in higher dimensions, taking the two dimensional case as an example. Finally, some numerical examples are presented in Section 6.
Notations. We use L ∞ to denote the space of bounded sequences. We use bold symbol f = (f x , f y ) to denote vector in R 2 . We take h to be the mesh size and Δt to be the time step in discretization.

TIME DISCRETIZATION
In this section we present four semi-discrete schemes with only time discretization, and comment on both pros and cons of each scheme.

Explicit scheme
The simplest scheme is the explicit scheme: .
This scheme is easy to implement. However, the sign of the right hand side is not certain and the positivity of solutions is difficult to show. The scheme can also be written in the following form ) . (

2.2)
Assume u n > 0 and u n+1 > 0, we calculate the relative energy defined by The above relation shows that the second term on the right hand side will be negative while the first term be positive. The right hand side is not guaranteed to have a definite sign unless Δt is very small. So for the explicit scheme, the positivity-preserving and energy dissipation properties are not guaranteed.

Fully implicit scheme
For stability reasons, we prefer to use an implicit scheme instead of an explicit one. The implicit scheme for Equation (1.1) reads as Such a scheme was studied in Reference [4], where the existence of weak solutions and the dissipation of Fisher information were shown. In other words, assuming u n > 0 and u n+1 > 0, the scheme dissipates the energy. For the convenience of reference, we present this result and a proof as follows.
Proof. We calculate the relative energy F(u n |u n+1 ) by We learned from the gradient structure that Equation (2.3) can be written as due to the positivity of u n+1 . Substituting it into (2.5) gives Notice that the positivity property is crucial in establishing the energy stability, but it seems difficult to directly prove such a property. Given u n , the equation for u n+1 is a fourth order nonlinear elliptic equation, the positivity of u n+1 does not seem to be derivable from a maximum principle.

A positivity-preserving scheme
Drawing ideas from References [23][24][25] in the design of unconditionally positive schemes for second-order Fokker-Planck equations, we set

Equation (1.1) can be rewritten in the form
which can be approximated by This scheme is linear in u n+1 and hence easy to code. Also it is unconditionally positivity-preserving. Proof. Set G n+1 = u n+1 ∕M n so that Note that both existence and regularity of the solution G n+1 are ensured by the classical elliptic theory. Assume G n+1 achieves its minimum at x 0 , then ∇G n+1 (x 0 ) = 0, and ΔG n+1 (x 0 ) ≥ 0. From the equation when evaluated at x 0 it follows that Hence G n+1 > 0 for all x ∈ T , so is u n+1 . ▪ It seems less obvious to verify the energy dissipation property.

An explicit-implicit scheme
The fourth scheme is One can view this scheme as an intermediate one between (2.3) and (2.1) or (2.8). The scheme is energy stable assuming u n and u n+1 are both positive.
Proof. The proof is similar to the proof of Lemma 2.2, except that Notice that here the assumption u n+1 > 0 is also needed in getting the result, since the function In the next section, we shall further discretize (2.9) in space and prove the positivity of numerical solutions in the subsequent sections.

Notations
We follow the notations in Reference [30]. First we define the following two grids on the torus T = [0, L] with mesh size h = L∕N, where N is the number of mesh intervals: We treat  and  as periodic. For example, we write x i as the ith element in , then x N = L and x N+1 = x 1 = h. The elements in  then can be written as x i+ 1 2 with i ∈ Z. We define the discrete N-periodic function space as Here we call  per the space of cell centered functions and  per the space of edge centered functions. We also define a subspace of  per by The discrete gradient D h and h are defined to be We define the average of the function values of nearby points bŷ The inner products are defined by The corresponding norms in each space are defined by Suppose f , g ∈  per and ∈  per , the following summation-by-parts formulas hold: We next introduce a weighted norm for g ∈ then the norm of g on • per is defined by Note that  is a symmetric and nonnegative operator, and satisfies where a discrete Poincaré inequality was used, with c = c 0 min > 0 being a constant. This implies by the Lax-Milgram theorem the existence and uniqueness of (3.5). The norm (3.6) can also be induced by the bilinear form on  per : with the following property:

The scheme
We proceed to study the time discretization scheme (2.9). We adopt the following fully discrete scheme as where D h , h are the difference operators defined in (3.2),û is the average operator defined in (3.3), and H n+1 is taken to be where the discretized energy functional The scheme can be written explicitly as Note that there are different ways to solve the above nonlinear scheme, we choose to use the Newton iteration combined with the Krylov subspace method in our numerical implementation.

Consistency of the scheme
Here we show the consistency of the scheme. Suppose u = u(x, t) is a smooth solution to Equation (1.1) and u n j is a finite difference solution to (3.7) and (3.8), we will show that the local truncation error defined by converges to 0 as Δt, h → 0. The local truncation error n j can be computed using Taylor's expansion. A careful calculation gives the following result, in which each term is evaluated at (x j , t n ): We have the following theorem.
where C depends on u.

SOLUTION PROPERTIES VIA OPTIMIZATION FORMULATION
In this section, we show that the finite difference scheme can be reformulated as an optimization problem by using the gradient flow structure. Inspired by the formal equivalence of (1.7) with the scheme (1.6), we will prove the equivalence after discretized in space. More precisely, we first establish the following theorem. .7) and (3.8) if and only if it is the unique minimizer to the following optimization problem Before proving this theorem, we first prepare a lemma, which is due to Reference [7] and holds for both one and higher dimensions. Here we present the result for the one dimensional case.
where min is the minimum of over the grid points and C does not depend on h.
Proof. We adapt the proof of lemma 3.2 in Reference [7]. By the definition (3.6), Since f ∈ • per, the use of the discrete Poincaré inequality gives which when inserted into the previous inequality leads to Using an inverse inequality (e.g., see References [7,8]) leads to The proof is divided into three steps: Step 1. We prove the existence and uniqueness of the optimization problem (4.1) for any > 0. According to the definition (3.6), we have Set v =û n D h f , we get The optimization problem (4.1) is equivalent to the following problem Notice that V h, is a region given by linear constraints and thus is convex. The objective function  1 (u, v) is also convex due to the convexity of [ 1 in v and the convexity of F h (u) in u, which can be seen from that each term in the summation in (3.9) is a convex function of u i and u i+1 . Therefore, the problem (4.3) is a convex optimization problem and thus has a unique solution. Since (4.3) and (4.1) are equivalent, there also exists a unique solution to the optimization problem (4.1) for any > 0.
Step 2. We show that u n+1 > 0 is the solution to the numerical scheme (3.7) and (3.8) if and only if it is the minimizer of (4.1) in • h, with ≤ 0 for some 0 > 0. First, suppose u n+1 is the minimizer in • h, of the convex optimization problem (4.1), it follows that Hence, we must have which is the numerical scheme (3.7) and (3.8). On the other hand for numerical solution u n+1 > 0, we take ≤ 0 < min u n+1 , and reverse the above calculation to get (4.4), which when combined with the convexity of  (u) and the assumption that the minimizer does not touch the boundary of  h, , implies that it is indeed the minimizer.
Step 3. We proceed to show that there exists 0 > 0 such that the minimizer does not touch the boundary of  h, for all 0 < ≤ 0 . We use a contradiction argument: suppose there exists a minimizer u * to the optimization problem (4.1) touching the boundary of  h, at some grid points i 0 , … , i m−1 with 1 ≤ m ≤ N − 1, that is We calculate the directional derivative along (4.5) In order to find a direction along which the above derivative is negative, we distinguish two cases characterized by a parameter q to be determined later: i. There exists a grid index First we consider case (i). Suppose u * reaches its maximum at the grid point i m , we take the direction v in Equation (4.5) to be It is easy to check that u * + sv ∈ • h, for small s. We illustrate by taking the case N = 3 as an example. When m = 1, suppose K∕h = 1 and i 0 = 2, the region of  h, is plotted in Figure 1(a). We use a ternary diagram to represent the set {(u 1 , u 2 , u 3 ) ∶ u 1 + u 2 + u 3 = K∕h}. The region inside the outer trianlge is  h,0 and the region inside the inner triangle is  h, . The red point represents the minimizer (u * 1 , u * 2 , u * 3 ). Here u * 1 > u * 3 > u * 2 so i 1 = 1. The vector v = (−1, 1, 0) is chosen to be in the direction as shown by the red arrow in the figure. When m = 2, suppose i 0 = 1, i 1 = 2, the vector v = (1, 1, −2) is chosen to be in the direction as shown by the red solid arrow in Figure 1(b). The two dashed arrows represent the vectors (1, 0, −1) and (0, 1, −1), respectively.
The directional derivative along the vector v is (4.6) Since u * j > for j ≠ i 0 , … , i m−1 , the value of u * at the adjacent grid points i p + 1 and i p − 1 of i p are either greater than or equal to , so Since u * i m is the maximum point and the summation Taking the above inequality and (4.7) into (4.6) and applying Lemma 4.2, we can deduce that 1 h ds which contradicts to the assumption that u * is a minimizer. Therefore, case (i) cannot occur.
Next we consider the case (ii). We first claim that there exists at least one index i m such that Otherwise, we would have u * j+1 − u * j ≤ qu * j for all j ∈ {1, … , N}. Using this we obtain Taking the summation of above equations gives if we take This is contradictory to the relation of ∑ N j=1 u * j = K. We can further assume because if this does not hold, we can move to the next grid point until such an inequality holds (notice that (4.10) holds for all grid points in the set {i 0 , … , i m−1 }).
We suppose u * reaches it maximum at the grid point i m+1 . Notice that i m+1 ≠ i 0 , … , i m−1 and i m+1 ≠ i m because i m cannot be a maximum point due to (4.9). We take the direction v to be It is easy to check u * + sv ∈

•
Ah, for small s. Taking N = 3 as an example. When m = 1, suppose K∕h = 1 and i 0 = 2, i 1 = 3, i 2 = 1, the region of  h, is plotted in Figure 1(c). The vector v = (−2,1,1) is chosen to be in the direction as show by the red solid arrow in the figure where the two dotted arrows represent the directions of (−1,1,0) and (−1,0,1).
The directional derivative of the objective function  (u) along the above direction v is (4.11) From the inequalities (4.9) and (4.10), we have Since u * j > for all j ∈ {1, … , N} and j ∈{i 0 , … , i m−1 } and u * i m+1 is the maximum point so (4.8) holds, we have which contradicts to the assumption that u * is a minimizer. Therefore, the situation (ii) cannot occur, either.
Combining the above estimates, we set so that for ≤ 0 neither (i) nor (ii) can occur. Hence no minimizer u * of the optimization problem (4.1) can touch the boundary of  h, for any ≤ 0 . Thus we finish the proof. ▪ The mass conservation, positivity-preserving, and energy stability of the scheme (3.7) and (3.8) can be obtained directly from Theorem 4.1. We state the results in the following theorem. (1) (Mass conservation) For any n ≥ 1, n ∈ N, (2) (Positivity-preserving) For any i = 1, … , N, u n i > 0 holds for any n ≥ 1, n ∈ N.
Proof. (1) This property is ensured by the conservative form of the numerical scheme, leading to (2) Starting from u 0 > 0, we apply Theorem 4.1 recursively to obtain for some constant n > 0, together we have for any n, This ensures that u n > 0.
(3) Since the solution of the numerical scheme (3.7) and (3.8) is the minimizer of the optimization problem (4.1), we have Therefore inequality (4.12) holds. ▪ Remark 4.1 In the proof of Theorem 4.1, is dependent of h, Δt andû n . Moreover, q is not guaranteed to be finite when h → 0 and Δt → 0, the above analysis cannot be carried over to the continuous equation.

Remark 4.2 A standard central discretization of Equation (1.1) is given by
which corresponds to the derivative of a discretized energy With this discrete energy Theorems 4.1 and 4.3 can be proven in similar manner. Moreover, Theorems 4.1 and 4.3 also hold for schemes with other forms of the discrete energy such as . (4.14) Notice that the expression (4.13) leads to a first order approximation of H, and (4.14) leads to a second order approximation.
See, for example, Reference [6]. Suppose we take the grid to be x i = a + (b − a)∕N ⋅ i, i = 0, … , N and the value of the numerical solution at x i equals u i , we incorporate the above boundary conditions by setting As a result, the scheme still features the structure as given in (3.7)-(3.9). Therefore Theorems 4.1 and 4.3 can also be established for the reflecting boundary conditions. This remark also applies to Theorem 5.1 when rectangular domains are considered.

THE SCHEME IN HIGHER DIMENSIONS
The numerical scheme (3.7) and (3.8) based on the gradient flow structure of Equation (1.1) can be generalized to higher dimensions. Theorem 4.1 and 4.3 can be shown to hold too. We will illustrate this by only considering the two dimensional case.

Notations
We take the domain to be T 2 = [0, L] × [0, L], and the mesh size h = L∕N for both x and y directions. We label the grid points with i, j for i = 1, … , N and j = 1, … , N. Since we are taking a periodic grid, we extend the label to the whole Z × Z by using  ×  with  defined in (3.1) and also label the mid-point grids using  ×  with  defined in (3.1). We define the periodic function spaces: and also We define the difference operators D x , D y and x , y as and the discrete divergence h as We also define the averagesf The inner products are defined by ) , and for vector functions: The corresponding norms in these spaces are defined accordingly. Suppose f , g ∈  per , and is defined on all the edge-center points  ×  ∪  × , then the following summation-by-parts formulas hold: Suppose ∈  ×  ∪  ×  and > 0, we introduce the following operator  on It can also be induced by the bilinear form

The scheme
The scheme in abstract form is where H n+1 is given by with .
Its explicit form is .
Again, the numerical scheme (5.6) and (5.7) is mass conservative, positivity-preserving and unconditionally energy stable. We state the main result in the following.
Proof of Theorem 5.1 We prove Theorem 5.1 for the two dimensional case in the same way as for the one dimensional case, that is, by first proving the result of Theorem 4.1 for 2D case. Since the scheme (5.6) and (5.7) follows the same gradient flow structure, we can prove that u n+1 > 0 is a solution to the numerical scheme (5.6) and (5.7) if and only if it is an interior minimizer of the optimization problem over the set with ≤ 0 for some 0 > 0. Here we only need to use ‖ ⋅ ‖  −1 u n and F h defined in (5.5) and (5.8) for the two dimensional case to replace the corresponding one dimensional definitions in Step 1 and 2 of the previous proof of Theorem 4.1.
What is left is to prove that there exists a constant 0 > 0 such that the minimizer of (5.11) does not touch the boundary. We follow the previous proof for the one dimensional case and calculate the directional derivative at the minimizer u * which touches the boundary at m grid points (1 ≤ m ≤ N 2 − 1) .

(5.13)
Likewise we distinguish the cases based on a positive constant q to be determined: i. There exists an index r ∈ {0, … , m − 1} such that u * satisfies u * i r +1,j r − u * i r ,j r > qu * i r ,j r or First we consider the case (i). Suppose u * reaches its maximum at i = i m , j = j m with value u * i m ,j m . We take and get This contradicts to the assumption that u * is a minimizer. Therefore, none of the assumptions (i) and (ii) hold for Therefore, u * cannot touch the boundary of  h, for ≥ 0 and Theorem 4.1 holds for the two dimensional case. Theorem 5.1 then follows the same way as in the one dimensional case. The above proof is easy to extend to the general dimensions and Theorem 5.1 holds for any dimension ≥ 1. ▪

NUMERICAL EXAMPLES
In this section, we present both one and two dimensional numerical examples to verify the mass conservation, positivity-preserving, energy dissipation properties, and the convergence order for scheme (3.7) and (3.8).

One dimensional example
First we consider a numerical example from Reference [3]. We take T = [0, 1] and the initial condition Here is taken to be 0.001 and m = 1 and 8. We take h = 0.01 and Δt = 10 −7 . Plots are made for t = 0.8 × 10 −6 , 3.2 × 10 −5 , 1 × 10 −4 and 7.2 × 10 −4 with solid, dashed, dotted, dash-dot and dash-dash lines, respectively. The solution u h is plotted in Figure 2. Mass variations, energy, and minimum values over time are plotted in Figure 3. These three figures verify the mass conservation, energy stability, and positivity of solutions, respectively. It can be seen that the minimum value remains positive.
We also preform error computations for different mesh sizes. We take Δt = 1 × 10 −6 and take the number of mesh intervals to be 10, 20, 40, 80, 160, with mesh size h = 0.1, 0.05, 0.025, 0.0125, 0.00625, respectively. Here we take Δt = 1×10 −6 , which appears to give accurate numerical solutions. Notice that the scheme was shown to be energy stable for any Δt > 0. However, for a larger Δt, the convergence of the nonlinear solver tends to be slower. We compare the solutions at time t = 7.2×10 −4 . We interpolate u h using the nearest neighbor method in space to find its difference withũ and calculate the L 2 error. The result is plotted in Figure 4. The solid line is the error with respect to the reference solution and the dashed line is the line 1∕h. Thus the numerical scheme is approximately of order 1 in h, that is, p = 1. This initial condition is chosen resembling the corresponding initial condition (6.1) in 1D with m = 8. However, there is one peak inside the domain T 2 at point x = y = 1 2 . So the diffusion of this inner point and the four peaks at the boundary dominate the behavior of the solutions.
We take h = 0.01 and Δt = 10 −6 to and the numerical solutions and their energy functions as well as minimum values are plotted over time in Figure 5. The minimum value decreases first and then increases over time but remains positive. We also perform error computations for different mesh sizes in 2D. The result is plotted in Figure 6. It is similar with the 1D example with same Δt = 10 −6 and h and the numerical convergence order is approximately one in h, that is, p = 1 in (6.2).