Scalar auxiliary variable approach
for conservative/dissipative partial differential equations
with unbounded energy

In this paper, we consider geometric numerical integration of Partial Differential Equations (PDEs) of the form

where $D$ is a skew-adjoint or negative semidefinite operator defined over an appropriate Hilbert space $V$, $\nabla E$ is the Fréchet derivative of a functional $\ E$, and $\dot{u}$ denotes the temporal derivative. (The details of the setting will be described later.) A class of PDEs of the form (1) involves many physically important equations, e.g., the Korteveg–de Vries equation, the Camassa–Holm equation [7], the Cahn–Hilliard equation [5] and the Swift–Hohenberg equation [35]. This paper is devoted to novel investigation of the SAV approach for these equations, which is a framework to construct linearly implicit geometric numerical integrators recently proposed in [33].

A prominent property of the above PDEs is a conservation/dissipation law:

where $\langle{\cdot},{\cdot}\rangle$ denotes the inner-product of $V$. Therefore, for these PDEs, countless numerical schemes inheriting the conservation/dissipation law have been considered in the literature.

In particular, there are some unified approaches to construct such schemes. The discrete gradient method [18, 30, 31] was originally devised for conservative/dissipative ODEs, which is now also used for PDEs in the form (1) [8]. On the other hand, Furihata [13] independently proposed the Discrete Variational Derivative Method (DVDM) (see also [15]) for the variational PDEs, which is now recognized as a combination of the discrete gradient method and an appropriate spatial discretization. The schemes based on the discrete gradient method (or DVDM) are often superior to general-purpose methods, in particular for a long run or numerically tough problems. Thus, they have been applied to many PDEs (see [14] and references therein).

However, at the price of their superiority, these schemes are unavoidably fully implicit and thus computationally expensive. Therefore, several remedies have been discussed in the literature. Some researchers have constructed structure-preserving and linearly implicit schemes for each specific PDE, for example, Besse [3] and Zhang, Pérez-García, and Vázquez [39] devised such schemes for the nonlinear Schrödinger equation. Moreover, for polynomial energy function, Matsuo and Furihata [28] proposed a multistep linearly implicit version of the DVDM (see also Dahlby and Owren [11]).

Most recently, for energy functions bounded below, Yang and Han [37] proposed an Invariant Energy Quadratization (IEQ) approach (see also [20, 36, 41]). A typical target of the IEQ approach is the Allen–Cahn equation [2]

where $F(u)=\frac{1}{4}\left\lparen u^{2}-1\right\rparen^{2}$ is a potential function. Since $F$ is non-negative, we can introduce an auxiliary function $q=\sqrt{F(u)+a}$ ($a$: a positive constant), and rewrite the Allen–Cahn equation (3). The reformulated one has a dissipation law of the modified energy $\int_{\Omega}\left\lparen\frac{1}{2}\left|\nabla u\right|^{2}+q^{2}\right\rparen\mathrm{d}x$. Since the modified energy is quadratic, we can easily construct dissipative schemes (see, e.g., [16, 40]). For example, we can use the implicit midpoint rule and its extensions, Gauss methods (see, e.g., [19, Chapter IV]). Moreover, it is possible to construct linearly implicit schemes (see, e.g., [17, 20, 36, 37, 41]).

Though the IEQ approach is successful, it is still a bit expensive due to the presence of the auxiliary function. To overcome it, Shen, Xu, and Yang [33] proposed a Scalar Auxiliary Variable (SAV) approach (see also [34]). There, instead of the auxiliary function in the IEQ approach, a scalar auxiliary variable $r=\sqrt{E_{1}(u)+a}$ is introduced, where $E_{1}(u)=\int_{\Omega}F(u)\mathrm{d}x$. Then the Allen–Cahn equation (3) can be rewritten in the form

Again, since the modified energy $\int_{\Omega}\frac{1}{2}\left|\nabla u\right|^{2}\mathrm{d}x+r^{2}$ is quadratic, we can construct dissipative schemes in several ways. The SAV approach often provides us with quite efficient numerical schemes: in addition to being the smaller system than the IEQ approach, in each step, the principal part is to solve two linear equations having the same $d\times d$ constant coefficient matrix.

Due to its attractive properties, the SAV approach has been intensively studied in these years. It has been applied to many PDEs, for example, the two-dimensional sine-Gordon equation [6], the fractional nonlinear Schrödinger equation [12], the Camassa–Holm equation [21], and the imaginary time gradient flow [42]. Shen and Xu [32] and Li, Shen and Rui [22] conducted a convergence analysis of SAV schemes. Akrivis, Li, and Li [1] devised and analyzed linear and high order SAV schemes. Moreover, the SAV approach is extended in several ways, for example, multiple SAV [10] and generalized Positive Auxiliary Variable (gPAV) [38].

Nevertheless, at least, there remain two issues on the SAV approach. First, the essential assumption “the energy functional is bounded below” is restrictive. In particular, when we deal with conservative PDEs, we often encounter an unbounded energy functional: for example, the KdV equation

Second, the scalar auxiliary variable $r$ in the modified equation (4) seems to be introduced in an unnatural way. As a consequence of this unclear derivation, the construction of resulting structure-preserving schemes is ad hoc.

There are some attempts to resolve the first problem [23, 24, 25, 27, 26, 4, 9]. The strategy of these studies is to change the definition of the scalar auxiliary variable. For example, in some of them, the authors replace the square-root function of the definition of $r$, such as the exponential function. Then, the auxiliary variable is always well-defined and thus the assumption for the energy functional can be removed. Although these studies successfully overcome the first problem, the second issue, namely the natural derivation of the SAV schemes, is still remained.

The aim of this paper is twofold. First we propose a novel approach to deal with unbounded energy functions by a decomposition of the energy functional. Second, we present a new interpretation of the SAV approach using gradient system expression. Furthermore, combining the above two results, we propose several numerical schemes in the framework of Crank-Nicolson method and the Runge-Kutta methods.

In order to address general energy functionals, we assume that the energy $E$ is expressed as

for some linear operator $L$ and lower-bounded functions $E_{\mathrm{L}}$ and $E_{\mathrm{U}}$ (see Assumption 1). We then introduce two auxiliary variables $r_{\mathrm{X}}$ ($\mathrm{X}=\mathrm{L}$ or $\mathrm{U}$) by

for some real numbers $a_{\mathrm{X}}$. In fact, the above decomposition of $E$ is always possible (Remark 1), and thus we can address general (possibly unbounded) energy functionals by this approach.

We then address the second problem: what is a natural way to construct SAV schemes? To answer this question, we begin with the modified energy

which is a quadratic function. The gradient of $\tilde{E}$ is then $\nabla\tilde{E}=[Lu,2r_{\mathrm{L}},-2r_{\mathrm{U}}]^{T}$. With these notations, we will show that the equation of the SAV approach (such as (4)) is expressed as the gradient system for the modified functional $\tilde{E}$. Further, this modified system inherits the conservation/dissipation property (i.e., variational structure) of the original gradient system (1) (see Lemma 1). From this point of view, existing SAV schemes are interpreted as a kind of discrete gradient method.

It is advantageous to keep the modified energy quadratic for the conservative case. Indeed, it is known that a certain class of Runge-Kutta methods preserves quadratic invariants [19]. Therefore, it is easy to construct higher-order time marching methods in the framework of Runge-Kutta methods and the SAV aproach, with the aid of our gradient system expression.

The remainder of this paper is organized as follows. In Section 2, we present the novel gradient system expression of the SAV approach, and show that the novel gradient system inherits the variational structure of the original system. Then, we propose second and fourth order schemes in the framework of our expression in Section 3. Finally, in Section 4, we present some numerical examples of our scheme to confirm the efficiency and we conclude this paper in Section 5.

In this section, we present a novel interpretation of the SAV approach, namely, gradient system expression. We state our result in an abstract setting.

Let $V$ be a Hilbert space over $\mathbb{K}$ equipped with an inner-product $\langle{\cdot},{\cdot}\rangle$, where $\mathbb{K}=\mathbb{R}$ or $\mathbb{K}=\mathbb{C}$, and $V^{*}$ be its dual. Let $E$ be a $\mathbb{K}$-valued $C^{1}$ function defined on some open subset of $V$, which is possibly unbounded, and $\nabla E(u)\colon V\to\mathbb{K}$ be the Fréchet derivative of $E$ at $u\in V$. We further let $D$ be a skew-adjoint or negative semidefinite linear operator on $V$ (not necessarily bounded). Then, we consider the gradient system of the form

where $u\colon(0,T)\to V$ is an unknown function and $T\in(0,\infty]$. We assume that all of the terms in (9) is well-defined.

We present SAV approach for the general gradient flow (9). To begin with, we make an assumption on the energy functional $E$.

We now introduce two scalar auxiliary variables $r_{\mathrm{L}}$ and $r_{\mathrm{U}}$ by

where $\mathrm{X}=\mathrm{L}$ or $\mathrm{U}$ and $a_{\mathrm{X}}>-\inf_{u}E_{\mathrm{X}}(u)$ are real constants. Then we define a modified energy functional by

which is quadratic, and thus $\nabla\tilde{E}=[Lu,2r_{\mathrm{L}},-2r_{\mathrm{U}}]^{T}$. Differentiating the scalar auxiliary variables, we have

where

Although $\phi_{\mathrm{X}}$ depends on $u$, we abbreviate the dependency to simplify the notation. With this notation, the time derivative of the modified energy functional becomes

According to the above observation, we consider the following equation

We show that the system (20) can be expressed as a gradient system. For each time $t$, define a linear operator $\Phi_{\mathrm{X}}(t)\in V^{*}$ by $\Phi_{\mathrm{X}}(t)(v)=\left\langle{\phi_{\mathrm{X}}(t)},{v}\right\rangle$ $(v\in V)$. The adjoint operator of $\Phi_{\mathrm{X}}$ is a multiplication operator $\Phi^{*}_{\mathrm{X}}(t)(r)=r\phi_{\mathrm{X}}(t)$. Then, (20) becomes

where $I$ is the identity operator on $V$. Since it is easy to see that

we obtain

Hence, letting $Z=V\times\mathbb{R}\times\mathbb{R}$, $z(t)=[u(t),r_{\mathrm{L}}(t),r_{\mathrm{U}}(t)]^{T}$, and

we can write (23) as

which is nothing but a gradient system in $Z$. In fact, the system (25) inherits the structure of the original gradient system (9).

The new expression (25) gives insights to construct structure-preserving schemes for the gradient system (9). In this section, we propose a unified approach to construct numerical schemes for gradient systems based on (25). Let ${\Delta t}>0$ be time increment, $t_{n}=n{\Delta t}$, $z^{n}\approx z(t_{n})$, and $z^{n}=[u^{n},r^{n}_{\mathrm{L}},r^{n}_{\mathrm{U}}]^{T}$. We use the same notation as in the previous section.