Entropy-dissipating finite-difference schemes for nonlinear fourth-order parabolic equations

The design of numerical schemes that preserve the structure of the associated partial differential equations is an important task in numerical mathematics. In this paper, we develop new finite-difference approximations conserving the mass, preserving the positivity, and dissipating the entropy of nonlinear fourth-order parabolic equations of the type

where ${\mathbb{T}}={\mathbb{R}}/{\mathbb{Z}}$ is the one-dimensional torus, $a$, $b\in{\mathbb{R}}$, and $\beta\geq 0$, together with the initial condition $u(0)=u^{0}$ in ${\mathbb{T}}$. We also discuss briefly the discretization of multi-dimensional equations; see Section 4.

A special case of (1) (with $a=b=0$) is the thin-film equation

which models the flow of a thin liquid along a solid surface with film height $u(x,t)$ or the thin neck of a Hele-Shaw flow in the lubrication approximation [1]. The multi-dimensional version is given by $\partial_{t}u=-\operatorname{div}(u^{\beta}\nabla\Delta u)$ and its discretization will be discussed in Section 4. Another example (with $a=-2$, $b=1$, $\beta=0$) is the Derrida–Lebowitz–Speer–Spohn (DLSS) equation

which arises as a scaling limit of interface fluctuations in spin systems [7] and describes the evolution of the electron density $u(x,t)$ in a quantum semiconductor under simplifying assumptions.

Equations (2) and (3) conserve the mass $\int_{\mathbb{T}}u(x,t)dx$ and preserve the positivity of the solution, even for the corresponding multi-dimensional versions [6, 14]. They also dissipate the entropy¹¹1Strictly speaking, equation (1) produces entropy and dissipates energy, but the notion entropy dissipation seems to be common in numerical schemes. in the sense that there exists a Lyapunov functional $S(t)=\int_{\mathbb{T}}s(u(x,t))dx$ with entropy density $s(u)$ such that the entropy production $-dS/dt$ provides gradient estimates. For instance, (2) and (3) dissipate the entropy density $s(u)=s_{\alpha}(u)$, where

if $\frac{3}{2}\leq\alpha+\beta\leq 3$ (thin-film equation) and $0\leq\alpha\leq\frac{3}{2}$ (DLSS equation); see [15, Section 4]. These bounds are optimal. We call $s_{1}(u)$ the Shannon entropy (density), $s_{0}(u)$ the logarithmic entropy, and $s_{\alpha}(u)$ for $\alpha\neq 0,1$ a Rényi entropy. It holds that $s_{\alpha}(u)\to s_{0}(u)$ pointwise for $\alpha\to 0$ and $s_{\alpha}(u)\to u\log u$ pointwise for $\alpha\to 1$. We prefer to define $s_{1}(u)$ as in (4) to avoid the additional term in $d(u\log u)/du=\log u+1$ which would complicate the subsequent computations.

The proof of the dissipation of the entropy is based on suitable integrations by parts and the chain rule [15]. On the discrete level, we face the problem that the chain rule is generally not available. We are aware of two general strategies.

The first strategy is to exploit the gradient-flow structure of the parabolic equation (if it exists). It involves only one integration by parts, and the discrete chain rule can be formulated by means of suitable mean functions. This idea was elaborated as the Discrete Variational Derivative Method for finite-difference approximations [11]. The gradient-flow formulation (with respect to the $L^{2}$-Wasserstein metric) yields a natural semi-discretization in time of the evolution using the minimizing movement scheme in finite-dimensional spaces from finite-volume or finite-difference approximations. These techniques were used in [8, 13, 19, 20] for the multi-dimensional DLSS equation. It allows for the proof of entropy dissipation of the Shannon entropy and the Fisher information $\int_{\mathbb{T}}(\sqrt{u})_{x}^{2}dx$, but not for general Rényi entropies, since no Wasserstein gradient-flow formulation seems to be available for these functionals. The thin-film equation with $0<\beta<1$ is shown in [18] to be a gradient flow with respect to a weighted Wasserstein metric. In the work [22], a finite-difference scheme that dissipates the discrete $H^{1}$ norm of the solution to the one-dimensional thin-film equation was analyzed.

The minimizing movement scheme is based on the implicit Euler method. We mention that higher-order time discretizations were investigated too, in the framework of semi-discrete problems; see [4] using the two-step BDF method and [17] using one-leg multistep generalizations. A generic framework for Galerkin methods in space and discontinuous Galerkin methods in time was presented in [9].

The second strategy uses time-continuous Markov chains on finite state spaces. Birth-death processes that define the Markov chain can be interpreted as a finite-volume discretization of a one-dimensional Fokker–Planck equation, and the dissipation of the discrete Shannon entropy can be proved. The nonlinear integrations by parts are reduced to a discrete Bochner-type inequality [5, 10, 16], which is obtained by identifying the Radon–Nikodym derivative of a measure involving the jump rates of the Markov chain [3, Section 2]. It seems that this idea is restricted to linear diffusion equations.

In this paper, we suggest a third strategy. The idea is to write the flux $J$ as a combination of derivatives of the function $s^{\prime}(u)$. This allows for integrations by parts that can be extended to the discrete level and it avoids the application of the chain rule. More precisely, we determine two functions $A$ and $B$ depending on $v:=s^{\prime}(u)$ and its derivatives such that $J=A_{x}-vB_{x}$. The function $s^{\prime}(u)$ is known in thermodynamics as the chemical potential, and the formulation of the flux in terms of the chemical potential seems to be natural from a thermodynamic viewpoint. We apply this idea to fourth-order parabolic equations for the first time. It turns out that for $s=s_{\alpha}$ with $\alpha\neq 1$, we can write

where $\xi_{1}=v_{x}/v$, $\xi_{2}=v_{xx}/v$, and the coefficients $\lambda_{i}$ depend on $a$, $b$, $\alpha$, and $\beta$; see (10) and (11). Integrating by parts twice and using equation (1) gives for $S_{\alpha}=\int_{\mathbb{T}}s_{\alpha}(u)dx$:

The task is to show that the integrand, written as a polynomial in $(\xi_{1},\xi_{2})$, is nonnegative for all values of $(\xi_{1},\xi_{2})\in{\mathbb{R}}^{2}$. It follows from the product rule $(vv_{x})_{x}=vv_{xx}+v_{x}^{2}$ that, for $\alpha\neq 1$,

Under certain conditions on the parameters, we expect that the integrand is bounded from below by $u^{\alpha+\beta}(\xi_{2}^{2}+\xi_{1}^{4})$, up to a factor, which yields some gradient estimates.

On the discrete level, we imitate this idea: The flux $J=A_{x}-vB_{x}$ and the variables $\xi_{1}$, $\xi_{2}$ of the polynomials $A$ and $B$ is discretized by central finite differences. For this, let ${\mathbb{T}}_{h}={\mathbb{T}}/(h{\mathbb{Z}})$ be a discrete torus with grid size $h>0$ and define the scheme

with the initial condition $u_{i}(0)=u^{0}(i)$ for $i\in{\mathbb{T}}_{h}$, where $u_{i}=u(i)$, $v_{i}=s^{\prime}(u_{i})$, and $A_{i}$ and $B_{i}$ are the polynomials $A$ and $B$, evaluated at $i\in{\mathbb{T}}_{h}$, respectively; see (19). We show below that with the discrete entropy $S_{\alpha}^{h}=h\sum_{i\in{\mathbb{T}}_{h}}s_{\alpha}(u)$, the discrete analog of (5) becomes

where $\xi_{1,i}$, $\xi_{2,i}$ approximate $v_{x}(i)/v(i)$, $v_{xx}(i)/v(i)$, respectively. This yields exactly the polynomial of the continuous case. Thus, we obtain the same conditions on the parameters as for the continuous equation.

We still need a discrete analog of the product rule $(vv_{x})_{x}=vv_{xx}+v_{x}^{2}$ to conclude. This is done by carefully choosing $\xi_{1,i}$ and $\xi_{2,i}$. Definition (20) ensures that $v_{i}^{2}(\xi_{2,i}+\xi_{1,i}^{2})=(v_{i+1}^{2}-2v_{i}+v_{i-1}^{2})/(2h^{2})$ which approximates $\frac{1}{2}(v^{2})_{xx}=vv_{xx}+v_{x}^{2}$. This choice is used in the central scheme (6) for $J_{i+1/2}$. Noncentral schemes require different definitions of $\xi_{1,i}$ and $\xi_{2,i}$; see Remark 7.

A drawback of our technique is that scheme (6) depends on the entropy to be dissipated. The scheme does not dissipate all admissible entropies. In applications, however, one usually wants to dissipate only that entropy which is physically relevant.

Our main results can be sketched as follows:

• •

  Lyapunov functional: Let $\alpha\geq 0$, $\alpha\neq 1$ and assume that

  (7)

  $$K(\alpha,\beta):=-2\alpha^{2}+(3a-4\beta+9)\alpha-2\beta^{2}+(3a+9)\beta-9(a+b+1)\geq 0.$$

  Then the continuous entropy $S_{\alpha}$ and the discrete entropy $S_{\alpha}^{h}$ are dissipated in the sense that $dS_{\alpha}/dt\leq 0$ and $dS_{\alpha}^{h}/dt\leq 0$, i.e., $S_{\alpha}$ and $S_{\alpha}^{h}$ are Lyapunov functionals for $u(t)$ and $u_{i}(t)$, respectively; see Theorems 2 and 5. Condition (7) is optimal for the thin-film and DLSS equations.

• •

  Entropy dissipation: Let $\alpha\geq 0$, $\alpha\neq 1$ and assume that $K(\alpha,\beta)>0$. Then there exists a constant $c(\alpha,\beta)>0$ such that for all $t>0$,

  (8)		$\displaystyle\frac{dS_{\alpha}}{dt}+c(\alpha,\beta)\int_{\mathbb{T}}u^{\alpha+\beta}(\xi_{2}^{2}+\xi_{1}^{4})dx$	$\displaystyle\leq 0,$
  (9)		$\displaystyle\frac{dS_{\alpha}^{h}}{dt}+c(\alpha,\beta)h\sum_{i\in{\mathbb{T}}_{h}}\bar{u}_{i}^{\alpha+\beta}(\xi_{2,i}^{2}+\xi_{1,i}^{4})$	$\displaystyle\leq 0,$

  where $\xi_{1}=v_{x}/v$, $\xi_{2}=v_{xx}/v$, $v=s_{\alpha}(u)$, $\bar{u}_{i}$ is an arbitrary average of $(u_{i})$, and $\xi_{1,i}$, $\xi_{2,i}$ are defined in (20). This result is proved in Theorems 2 and 5.

• •

  Case $\alpha=1$: For the case $\alpha=1$, we need the formulation $J=A_{x}-wB_{x}$ instead of $J=A_{x}-vB_{x}$, where $w=s_{0}^{\prime}(u)=-u^{-1}$ and $v=\log u$, since $J$ generally does not depend on the logarithm. For details, see Proposition 8.

• •

  Case $\alpha=0$: We show that scheme (6) with $\alpha=0$ possesses global positive solutions. This result is a consequence of the discrete entropy inequality and mass conservation, which imply that $h\sum_{i\in{\mathbb{T}}_{h}}(u_{i}(t)-\log u_{i}(t))$ is bounded for all $t>0$. Consequently, $u_{i}(t)-\log u_{i}(t)$ is bounded for all $i\in{\mathbb{T}}_{h}$ and $t>0$, and since the function $s\mapsto s-\log s$ tends to infinity if either $s\to 0$ or $s\to\infty$, this proves that $u_{i}(t)$ is bounded from below and above. We refer to Proposition 6 for details.

• •

  Multi-dimensional case: In principle, the multi-dimensional case can be treated using functions $A$ and $B$ with many variables. Practically, however, the computations are becoming too involved and it may be unclear how to discretize mixed derivatives. One idea to overcome this issue is to use scalar variables only, like $\xi_{1}=|\nabla v|/v$ and $\xi_{2}=\Delta v/v$. This allows us to treat the multi-dimensional thin-film equation; see Proposition 10.

The paper is organized as follows. We prove the continuous entropy inequality (8) in Section 2 and the discrete entropy inequality (9) in Section 3. A scheme for the multi-dimensional thin-film equation is proposed and analyzed in Section 4. Numerical simulations for the thin-film and DLSS equations in one space dimension and for a thin-film equation in two space dimensions are presented in Section 5.

To prepare the discretization, we need to analyze the entropy dissipation properties of the continuous equation (1). We show first that $J$ can be written as $A_{x}-vB_{x}$ with $v=s_{\alpha}^{\prime}(u)$ and functions $A$ and $B$ which depend on $v$, $v_{x}$, and $v_{xx}$.

The lemma shows in particular that the formulation $J=A_{x}-vB_{x}$ is not unique. This fact is used to optimize later the range of admissible parameters $\alpha$, $\beta$, $a$, and $b$.

For the following theorem, we recall definition (4) of $s_{\alpha}$ and set $S_{\alpha}(u)=\int_{\mathbb{T}}s_{\alpha}(u)dx$.