# Two novel numerical methods for gradient flows: generalizations of the
Invariant Energy Quadratization method

Yukun Yue
Department of Mathematical Sciences
Carnegie Mellon University
5000 Forbes Avenue, Pittsburgh, PA 15213, USA<EMAIL_ADDRESS>
###### Abstract.

In this paper, we conduct an in-depth investigation of the structural
intricacies inherent to the Invariant Energy Quadratization (IEQ) method as
applied to gradient flows, and we dissect the mechanisms that enable this
method to uphold linearity and the conservation of energy simultaneously.
Building upon this foundation, we propose two methods: Invariant Energy
Convexification and Invariant Energy Functionalization. These approaches can
be perceived as natural extensions of the IEQ method. Employing our novel
approaches, we reformulate the system connected to gradient flow, construct a
semi-discretized numerical scheme, and obtain a commensurate modified energy
dissipation law for both proposed methods. Finally, to underscore their
practical utility, we provide numerical evidence demonstrating these methods’
accuracy, stability, and effectiveness when applied to both Allen-Cahn and
Cahn-Hilliard equations.

Y.Y.’s work was supported in part by NSF award DMS-1912854

##  1\. Introduction

Gradient flows are a class of partial differential equations (PDEs) that arise
in various scientific and engineering fields, such as fluid dynamics,
materials science, and optimization, as demonstrated by [4, 10, 14, 18, 19,
30, 31, 45]. These equations model the evolution of a given quantity under the
influence of a driving force, which is derived from an energy functional. The
study of gradient flows has attracted considerable attention in recent years
due to their wide applicability and inherent mathematical structures. In this
paper, we present two innovative numerical methods for gradient flows that
extend the Invariant Energy Quadratization (IEQ) method, which is recently
introduced. These methods offer a generalization of the pre-existing method,
leveraging the mathematical structures possessed by the method while
preserving the favorable properties, such as the energy-stable property and
the ability to construct efficient linear schemes to solve the problems.

To illustrate, let’s consider a free energy functional
$E(\phi)=\int_{\Omega}\rho\left(\phi(x)\right)\,dx$, where
$\rho\left(\phi(x)\right)=\left[\frac{1}{2}\lvert\nabla\phi\rvert^{2}+F(\phi)\right]$
serves as the energy density function of $E(\phi)$. The corresponding gradient
flow can then be formulated as

(1.1a) $\frac{\partial\phi}{\partial t}=\mathcal{G}\mu,$ (1.1b)
$\mu=\frac{\delta E}{\delta\phi}=-\Delta\phi+\frac{\delta
F}{\delta\phi}:=-\Delta\phi+f(\phi),$

with $\phi(0)=\phi_{0}$ to be its initial condition. In this context,
$\mathcal{G}=I$ signifies the identity operator when we are considering a
gradient flow in $L^{2}$, and $\mathcal{G}=\Delta$, which denotes the
Laplacian, when we are investigating a gradient flow in $H^{-1}$ [3, 36].

Gradient flows have a rich mathematical structure and are related to several
important concepts in mathematics, such as the Wasserstein distance, optimal
transport, and convex optimization. There has been a surge of interest in
unearthing its mathematical properties, with notable contributions made by
[13, 15, 35, 37]. Concurrently, researchers have been striving to establish
various numerical methods to resolve gradient flow problems. These efforts
have given rise to techniques such as the convex-splitting method [6, 23, 38,
43], the Invariant Energy Quadratization (IEQ) [7, 25, 26, 46, 48, 52], and
the Scalar Auxiliary Variable (SAV) methods [39, 40, 41, 42]. These methods
allow for the development of efficient schemes for solving gradient flows, and
hold the energy-stable property, assuring the numerical solutions maintain
certain physical and mathematical properties of the continuous problem.

The IEQ method, particularly, has gained popularity as a numerical method in
recent years. It introduces an auxiliary variable whose quadratization equals
the original energy density function when the energy density function has a
lower bound. It is based on the idea of conserving certain invariants of the
continuous problem in the discrete setting. As an illustration, one could
define

$r(\phi)=\sqrt{F(\phi)+A_{1}},$

with $A_{1}$ as a constant ensuring $r$ to be well-defined. Provided that $F$
is bounded from below, an appropriate constant $A_{1}$ can always be found.
Consequently, system (1.1) can be reformulated as

(1.2a) $\frac{\partial\phi}{\partial t}=\mathcal{G}\mu,$ (1.2b)
$\mu=-\Delta\phi+2rP,$ (1.2c) $r_{t}=P\phi_{t},$ (1.2d) $P=\frac{\frac{\delta
F}{\delta\phi}}{2\sqrt{F(\phi)+A_{1}}}$

The chain rule can be conveniently applied to the auxiliary variable $r$ to
validate this reformulation. This approach has successfully facilitated the
construction of linear numerical schemes for various gradient-flow type
problems, including the Cahn-Hilliard equation [47, 49], the Ericksen–Leslie
model [12] and Beris-Edwards model for liquid crystals [24, 50, 53], and the
sine-Gordon equation [21, 27]. For instance, a linear numerical scheme can be
constructed by treating $r$ implicitly and $P$ explicitly. Specifically, we
have:

(1.3a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (1.3b)
$\mu^{n+1}=-\Delta\phi^{n+1}+2r^{n+1}P^{n},$ (1.3c)
$r^{n+1}-r^{n}=P^{n}:(\phi^{n+1}-\phi^{n}),$ (1.3d) $P^{n}=\frac{\frac{\delta
F}{\delta\phi}(\phi^{n})}{2\sqrt{F(\phi^{n})+A_{1}}}.$

The resulting numerical scheme is characterized by a modified energy
dissipation law (see, for instance, [53]). One of the notable strengths of the
IEQ method is its energy-stable property, which ensures a monotonic decrease
of the discrete energy functional along the trajectory of the numerical
solution, akin to the behavior observed in the continuous problem. This
property plays a vital role in the long-term behavior of the numerical
solution and the preservation of inherent physical and mathematical
structures.

The power of the IEQ method fundamentally lies in its innovative decomposition
of the nonlinear term into a product of two functions, specifically
emphasizing which part should be treated implicitly. This quadratization
formulation aligns seamlessly with a usual form of the energy term encountered
in physical processes, as evidenced by numerous scientific studies. In
addition, quadratization formulation induces a linear term to appear, leading
to a linear scheme. However, the energy term in various natural scientific
systems may not be confined to a quadratic form. This observation prompts us
to question the possibility of other formulations capable of preserving the
advantageous properties of the IEQ method. In other words, we are interested
in discovering other indigenous decomposition of the nonlinear term to
construct energy-stable, linear numerical schemes.

In this paper, we propose two novel strategies for extending the Invariant
Energy Quadratization (IEQ) method, namely the Invariant Energy
Convexification (IEC) and the Invariant Energy Functionalization (IEF)
approaches. Both methods are based on introducing an auxiliary variable
originating from using different functions to replace the original energy
density function. Moreover, they both induce a modified energy dissipation law
within the reformulation of the original gradient-flow type system.
Importantly, these two methods retain the advantageous ability to induce
efficient linear numerical schemes that rival the efficacy of the IEQ method.
A noteworthy contribution of these proposed approaches is their potential to
accommodate various functions as auxiliary variables during the transformation
process. This flexibility suggests the possibility of identifying an optimal
form for implementing the method tailored to the specifics of a given physical
system.

The organization of the remainder of this paper is as follows: In Section 2,
we introduce the Invariant Energy Convexification (IEC) method alongside a new
reformulation of the system as described in equation (1.1). We shall proceed
to derive a semi-discrete numerical scheme and verify the stability of these
numerical schemes under suitably modified energy. Likewise, in Section 3, we
outline the Invariant Energy Functionalization (IEF) method, present a
corresponding numerical scheme, and confirm its energy stability. In Section
4, we substantiate the accuracy and efficacy of our proposed methods through a
series of numerical experiments, providing compelling evidence of their
utility.

##  2\. The invariant energy Convexification

In this section, we introduce our first numerical method that expands upon the
core principles of the IEQ method. As it has been stated in the introduction,
the IEQ method has established itself as a potent and effective technique for
addressing a broad class of problems. At the heart of the IEQ method is the
introduction of an auxiliary variable representing the square root of the
energy density function. The incorporation of this auxiliary variable enables
the method to more effectively manage the non-linear terms that emerge in
various problems. Specifically, the IEQ method handles a portion of the
resulting non-linear term implicitly while maintaining the remainder in an
explicit form. This blend of implicit and explicit treatment is a
distinguishing feature of the IEQ method and contributes to its efficacy in
constructing linear numerical schemes. A critical question arising from this
approach is whether the decomposition of the non-linear term in the IEQ method
is unique or if it can be further generalized. This will be discussed in the
following.

###  2.1. Reformulation with L-smooth convex function

To find a more general class of functions that share the advantages of the IEQ
method, we need to figure out the nature of its effectiveness first. We notice
that the energy-conserving property mainly relies on the convexity of the
auxiliary function used to replace the original energy density function.
Meanwhile, the linearity relies on the fact that the derivative of a quadratic
function is a linear function. Therefore, it enlightens us to find a class of
convex functions to introduce the auxiliary variable and look for a possible
linear approximation of its derivative to construct the numerical scheme. To
investigate this question, we propose a new formulation designed to imitate
the IEQ formulation’s structure and broaden its applicability while preserving
its essential advantages. As a result, we discover that quadratization is not
the sole viable option for variable transformation; a particular class of
convex functions can also pave the way to an energy-stable linear scheme.

Specifically, let us presume the energy density function $F$ is bounded from
below. We will take $c:\mathbb{R}\to\mathbb{R}$ as a smooth convex function
that is monotonically increasing on a connected set $K$, and $\mathbb{R}^{+}$
is contained $c(K)$. We further assume that $c$ is $L-$smooth [34], which
means that for any $x,y$ within domain of function $c$, a constant
$0<L<\infty$ exists such that

$|\nabla c(x)-\nabla c(y)|\leq L\lvert x-y\rvert.$

This condition implies that

(2.1) $c(y)\leq c(x)+c^{\prime}(x)(y-x)+\frac{L}{2}\lvert y-x\rvert^{2}.$

Now, consider $r:[0,T]\times\Omega\to\mathbb{R}$ to be a function solving
equation $c\left(r(t,x)\right)=F\left(\phi(t,x)\right)+A_{1}$ for every
$(t,x)\in[0,T]\times\Omega$ where $A_{1}$ is the constant to ensure the non-
negativity of $F$. Given the invertibility of $c$ on $\mathbb{R}^{+}$, this
definition is indeed valid. Differentiating both sides of the equation results
in

$c^{\prime}(r)\,\frac{\delta r}{\delta\phi}=\frac{\delta
F(\phi)}{\delta\phi}=f(\phi).$

Therefore, systems (1.1) will be rewritten as

(2.2a) $\phi_{t}=\mathcal{G}\mu,$ (2.2b)
$\mu=-\Delta\phi+c^{\prime}(r)\frac{\delta r}{\delta\phi},$ (2.2c)
$r_{t}=\frac{\delta r}{\delta\phi}\phi_{t}.$

Taking inner product above with $\mu,\phi_{t}$ and $c^{\prime}(r)$
respectively, one can easily obtain a modified energy dissipation law as:

(2.3)
$\frac{d}{dt}\left[\frac{1}{2}\|\nabla\phi\|^{2}+\int_{\Omega}c(r)\right]=(G\mu,\mu)\leq
0,$

for $\mathcal{G}=-I$ and $\mathcal{G}=\Delta$ where $\|\cdot\|$ is the $L^{2}$
norm. If $c^{\prime}(r)$ is a linear function with respect to $r$, then we can
construct a linear numerical scheme simply by treating $c^{\prime}(r)$
implicitly and $\frac{\delta r}{\delta\phi}$, which is how IEQ works. However,
when $c^{\prime}(r)$ is not linear, treating $c^{\prime}(r)$ implicitly will
result in a nonlinear equation to be solved at each step which generally loses
a very strong advantage of this approach. In addition, we want to point out
that the convexity is not crucial in the continuous case as reformulation
(2.2) and the corresponding modified energy dissipation law actually can work
by simply choosing $c$ as a smooth function, with the convexity condition
being dropped. However, convexity will play an important role in the way of
constructing linear energy-stable schemes in the discrete case. It will be the
main problem that we want to solve in the next subsection.

###  2.2. Numerical scheme

Our newfound approach provides the cornerstone for devising a linear numerical
scheme to resolve (2.2). We merely need a linear approximation for
$c^{\prime}(r)$ in the discretized scenario. With this in mind, we propose the
following first-order semi-discrete scheme:

(2.4a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (2.4b)
$\mu^{n+1}=-\Delta\phi^{n+1}+\left[c^{\prime}(r^{n})+\alpha
L(r^{n+1}-r^{n})\right]P^{n},$ (2.4c)
$r^{n+1}-r^{n}=P^{n}(\phi^{n+1}-\phi^{n}),$

where $P^{n}=\frac{\delta r}{\delta Q}(\phi^{n})$ and $\alpha\geq\frac{1}{2}$
is a chosen parameter. The elegance of this scheme lies in its efficient
solvability at each time step. Specifically, given that $(\phi^{n},r^{n})$ are
known, substituting (2.2c) into (2.2b) enables the replacement of $r^{n+1}$ by
the formula of $\phi^{n+1}$, hence allowing the computation of
$(\phi^{n+1},r^{n+1})$ by first finding $\phi^{n+1}$ from

$(I+\Delta t\,\mathcal{G}\Delta)\phi^{n+1}-\alpha L\Delta
t\,\mathcal{G}\left(P^{n}\phi^{n+1}P^{n}\right)=\phi^{n}+\Delta
t\mathcal{G}\,\left[c^{\prime}(r^{n})-\alpha LP^{n}\phi^{n}\right]P^{n},$

and subsequently updating $r^{n+1}$ by (2.4c). Alternatively, the scheme can
be executed by amalgamating $\mu^{n+1}$ in the solving process. This approach
determines $(\phi^{n+1},\mu^{n+1},r^{n+1})$ by resolving the linear system

(2.5) $\begin{pmatrix}\frac{1}{\Delta t}I&-\mathcal{G}&0\\\ \Delta&I&-\alpha
L\tilde{P}^{n}\\\ -\tilde{P}^{n}&0&I\end{pmatrix}\begin{pmatrix}\phi^{n+1}\\\
\mu^{n+1}\\\ r^{n+1}\end{pmatrix}=\begin{pmatrix}\frac{1}{\Delta t}\phi^{n}\\\
c^{\prime}(r^{n})P^{n}-\alpha Lr^{n}P^{n}\\\
r^{n}-P^{n}\phi^{n}\end{pmatrix},$

at each time step, given that $(\phi^{n},r^{n})$ have been ascertained. Both
implementations offer a linear pathway to updating the numerical results,
thereby providing us a streamlined alternative to tackling nonlinear gradient-
flow type problems, as opposed to implicitly addressing the nonlinear term and
having to resolve a nonlinear equation at each juncture.

######  Remark 2.1.

In the second implementation approach above, we use $\tilde{P}^{n}$ to denote
a functional operator, a multiplication of $P^{n}$ to $r^{n+1}$ point-wisely.
In the fully-discrete case, assume the space has been discretized by a
$N_{x}\times N_{y}$ grid, we will treat $\phi^{n+1},\mu^{n+1},r^{n+1}$ as a
$1\times(N_{x}*N_{y})$ vector, respectively. And $\tilde{P}^{n}$ will be taken
as a $(N_{x}*N_{y})\times(N_{x}*N_{y})$ diagonal matrix by expanding the
values of $P^{n}$ at the $N_{x}*N_{y}$ grid points onto its diagonal elements.
The same notation $\tilde{\cdot}$ will be repeatedly used in the following to
denote a point-wise multiplication operator.

Besides sustaining the linearity property, we can also immediately obtain the
following modified discretized energy dissipation law from scheme (2.4) .

######  Theorem 2.2.

The numerical scheme (2.4) is energy-stable. Specifically, define

(2.6) $E^{n}=\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\int_{\Omega}c(r^{n}),$

then

$E^{n+1}-E^{n}\leq 0.$

###### Proof.

Taking inner product of (2.2a) with $\mu^{n+1}\Delta t$, (2.2b) with
$\phi^{n+1}-\phi^{n}$, (2.2c) with
$\left[c^{\prime}(r^{n})+\frac{L}{2}(r^{n+1}-r^{n})\right]$, and using
property of $L$-smooth function (2.1), we obtain

$\displaystyle(\mathcal{G}\mu^{n+1},\mu^{n+1})\Delta t$
$\displaystyle=(\phi^{n+1}-\phi^{n},\mu^{n+1})$
$\displaystyle=\left(\phi^{n+1}-\phi^{n},-\Delta\phi^{n+1}+\left[c^{\prime}(r^{n})+\alpha
L(r^{n+1}-r^{n})\right]P^{n}\right)$
$\displaystyle=\frac{1}{2}\|\nabla\phi^{n+1}\|^{2}-\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\frac{1}{2}\|\nabla\phi^{n+1}-\nabla\phi^{n}\|^{2}$
$\displaystyle\quad\,\,+\int_{\Omega}\left[c^{\prime}(r^{n})(r^{n+1}-r^{n})+\alpha
L(r^{n+1}-r^{n})^{2}\right]\,dx$
$\displaystyle\geq\frac{1}{2}\|\nabla\phi^{n+1}\|^{2}-\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\frac{1}{2}\|\nabla\phi^{n+1}-\nabla\phi^{n}\|^{2}+\int_{\Omega}\left[c(r^{n+1})-c(r^{n})\right]\,dx$
$\displaystyle\geq E^{n+1}-E^{n}.$

Using the fact that $(\mathcal{G}\mu^{n+1},\mu^{n+1})\leq 0$ independent of
choice of $\mu^{n+1}$, we have finished the proof. ∎

######  Remark 2.3.

From the deduction, we can see that L-smoothness of $c(x)$ plays a key role in
keeping the energy dissipated rule. Specifically, if we modify $c(x)$ to be a
concave function, then simply taking $L=0$ is enough to ensure the
inequalities appeared in the proof to hold and it will also result in a linear
energy-stable numerical scheme.

Following from this theorem, we can immediately obtain the following estimate
holds for $\mu^{n+1}$:

######  Corollary 2.4.

Assume $E^{0}$ is bounded. For fixed $N>0$, we have

$0\leq-\sum_{n=0}^{N}(\mathcal{G}\mu^{n},\mu^{n})\Delta t\leq-E^{N}+E^{0}\leq
E^{0},$

and so it is bounded. Specifically, if $\mathcal{G}=\Delta$, we have

$0\leq\sum_{n=0}^{N}\|\nabla\mu^{n}\|^{2}\Delta t\leq E^{0},$

if $\mathcal{G}=-I$, we have

$0\leq\sum_{n=0}^{N}\|\mu^{n}\|^{2}\Delta t\leq E^{0}.$

###  2.3. An example of the IEC scheme

Having laid the groundwork for the standard formulation of the IEC approach
and its corresponding numerical scheme’s general form, we now face the
necessity of an exact form of the chosen L-smooth convex function for concrete
computational applications. As we conclude this section, we propose a viable
candidate for such a function, demonstrating that it is indeed possible to
construct an IEC scheme using a function that satisfies these criteria. Our
goal is to identify an option that possesses the properties necessary to
induce a valid IEC approach and is also conveniently calculable, thus
facilitating our numerical implementation.

Before delving into that, it is worth mentioning that the IEQ approach can be
considered a specific instance of the IEC formulation. Indeed, the quadratic
function is an $L$-smooth convex function with $L=2$. Thus, by selecting
$\alpha=1$ in (2.4), the resulting numerical scheme aligns perfectly with the
scheme induced by the IEQ formulation. This observation serves to justify our
perspective that the IEC method is a natural generalization of the IEQ method.

As an alternative to the quadratic function, we turn our attention to the
Softplus function [44]. Defining $c(r)=\ln(1+e^{r})$, we quickly find that the
first-order derivative, $c^{\prime}(r)=\frac{e^{r}}{1+e^{r}}$, and the second-
order derivative, $c^{\prime\prime}(r)=-\frac{e^{r}}{(1+e^{r})^{2}}$, are
uniformly bounded. Consequently, $c(r)=\ln(1+e^{r})$ establishes itself as a
non-negative, monotonically increasing, L-smooth convex function with
$L=\frac{1}{4}$. Taking $c(r)=F(\phi)+A_{1}$, this leads to

(2.7) $r(\phi)=\ln(e^{F(\phi)+A_{1}}-1),$

and subsequently to $P^{n}=P(\phi^{n})$ where

(2.8) $P(\phi)=\frac{\delta
r}{\delta\phi}(\phi)=\frac{e^{F(\phi)+A_{1}}}{e^{F(\phi)+A_{1}}-1}f(\phi).$

These equations provide an illustrative example of (2.4), demonstrating the
applicability of the Softplus function in this context. We can conclude this
procedure in the following Algorithm 1.

Algorithm 1 IEC scheme induced by Softplus function

1:procedure

2: Set the step size $\Delta t$

3: Set the total number of iterations $N$

4: Set the value for $\alpha$

5: Set the initial value $\phi^{0}=\phi_{0}$

6: Compute $r^{0}=r(\phi^{0})$ using equation (2.7)

7: Compute $P^{0}=P(\phi^{0})$ using equation (2.8)

8: for $n=0,1,2,...,N-1$ do

9: Construct the matrix $A=\begin{pmatrix}\frac{1}{\Delta
t}I&-\mathcal{G}&0\\\ \Delta&I&-\frac{\alpha}{4}\tilde{P}^{n}\\\
-\tilde{P}^{n}&0&I\end{pmatrix}$ using $P^{n}$

10: Construct the right-hand side $b=\begin{pmatrix}\frac{1}{\Delta
t}\phi^{n}\\\ \frac{e^{r^{n}}}{1+e^{r^{n}}}P^{n}-\frac{\alpha}{4}r^{n}P^{n}\\\
r^{n}-P^{n}\phi^{n}\end{pmatrix}$

11: Update $\phi^{n+1}$, $r^{n+1}$ by solving the linear system
$A\begin{pmatrix}\phi^{n+1}\\\ \mu^{n+1}\\\ r^{n+1}\end{pmatrix}=b$

12: Update $P^{n+1}=P(\phi^{n+1})$ using equation (2.8)

13: end for

14:end procedure

The flexibility in selecting the function $c$ is relatively wide as long as it
conforms to the previously defined prerequisites that ensure a valid numerical
scheme. In practice, this adjustment necessitates an alteration of the
formulas for $r$ and $P$ and an adjustment of the expression for
$c^{\prime}(r)$ used in computing the right-hand side vector $b$ in the
preceding algorithm. An important observation is that the second derivative of
the Softplus function is uniformly bounded, thereby ensuring it is L-smooth
uniformly across the full range of real numbers. Additionally, we may consider
alternative convex functions like the exponential function $e^{r}$ and
polynomials $r^{k}$ where $k\geq 3$, as potential replacements for the
quadratic function, given that some other specific assumptions are met.
Specifically, the exponential and polynomial functions both possess the local
$L-$smoothness. Hence, if the range of $r$ is bounded, a sufficiently large
$L$ can be identified to make these functions locally $L-$smooth and
facilitate the implementation of corresponding numerical schemes. In Section
4, we will display some numerical evidence supporting these formulations’
applicability.

##  3\. Invariant Energy Functionalization

As previously established, the two main benefits of the IEQ method are its
linearity and energy-dissipation preservation. The conservation of energy,
particularly in a discrete case, is primarily ensured through the convexity of
the auxiliary function, which replaces the original energy density function.
The IEC formulation can be seen as an evolution of this concept, placing the
preservation of convexity at its core. In contrast, the linear approximation
employed in (2.4) serves as a mechanism to design a linear scheme following
the assurance of convexity. It reveals a fundamental tenet of the IEC scheme:
the priority of maintaining convexity over linearity.

Nevertheless, it is the linearity that brings efficiency to the IEQ method. If
we shift our focus towards preserving linearity, it becomes evident that this
attribute is primarily derived from the auxiliary variable’s linearity with
respect to the auxiliary function’s derivative. This understanding guides us
to consider decomposing the auxiliary function $c(r)$ into a product of the
auxiliary variable $r$ and $\frac{c(r)}{r}$. Here we can drop the convexity
assumption of $c(r)$ and simply take it as a smooth function. If we consider
$\frac{c(r)}{r}$ as an independent auxiliary variable, denoted as $g(r)$, then
the derivative of $rg(r)$ ends up being linear concerning the auxiliary
variable $r$. This idea provides an alternative means of replacing the
original energy density function to achieve a linear scheme but emphasizes
preserving linearity more than convexity. We introduce this method in the
following section, referring to it as the ”Invariant Energy Functionalization”
method. This name stems from integrating the auxiliary variable $r$ with a
specific function $g$ that fulfills some conditions.

###  3.1. Model reformulation

Let us consider $g$ to be a smooth function on $\mathbb{R}$. Analogous to the
hypothesis in the IEC formulation, we assume that $s(r):=rg(r)$ is invertible
over a connected set $K$ with $\mathbb{R}^{+}\subset s(K)$. We define
$r:[0,T]\times\Omega\to K$ to be a function that resolves the equation
$r(t,x)g\left(r(t,x)\right)=F\left(\phi(t,x)\right)+A_{1}$ for all
$(t,x)\in[0,T]\times\Omega$, where $A$ is a constant employed to ascertain
positivity of $F$. Utilizing the chain rule, we can then derive the following:

$\left[rg^{\prime}(r)+g(r)\right]\frac{\delta r}{\delta\phi}=f(\phi).$

Thus, we can rewrite gradient-flow system (1.1) as

(3.1a) $\phi_{t}=\mathcal{G}\mu,$ (3.1b)
$\mu=-\Delta\phi+\left[rg^{\prime}(r)+g(r)\right]\frac{\delta r}{\delta\phi},$
(3.1c) $r_{t}=\frac{\delta r}{\delta\phi}\phi_{t}.$ (3.1d)
$g_{t}=g^{\prime}(r)\,r_{t}$

Given the equations above, we can deduce a modified energy dissipation law by
computing the inner product of (3.1) with $\mu$, $\phi_{t}$,
$\left[rg^{\prime}(r)+g(r)\right]$, and $r$ respectively, and then summing the
results. This leads us to

(3.2)
$\frac{d}{dt}\left[\frac{1}{2}\|\nabla\phi\|^{2}+\int_{\Omega}rg(r)\right]=(G\mu,\mu)\leq
0,$

which holds for both $\mathcal{G}=-I$ and $\mathcal{G}=\Delta$. It’s worth
pointing out that one of the significant advantages of this formulation is
that the auxiliary variable $r$ appears in a linear form in the reformulation
of $f(\phi)$. This linear structure facilitates the development and
implementation of a linear scheme to solve the equation, a topic that we will
discuss more in the upcoming subsection.

###  3.2. Numerical scheme

Motivated by the formulation proposed above, a natural next step is
introducing an additional variable to represent $g$ within the corresponding
numerical scheme. In the continuous formulation, $g=g(r)$ is a dependent
function of $r$. Thus, once $r$ is established, the value of $g$ is
subsequently determined. Nonetheless, to preserve the linearity of the
numerical scheme, we propose the introduction of a discrete variable for $g$,
and update its value through an explicit discretization. This idea leads to
the following Invariant Energy Functionalization (IEF) numerical scheme:

(3.3a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (3.3b)
$\mu^{n+1}=-\Delta\phi^{n+1}+\left[r^{n+1}g^{\prime}(r^{n})+g^{n+1}\right]P^{n},$
(3.3c) $r^{n+1}-r^{n}=P^{n}(\phi^{n+1}-\phi^{n}),$ (3.3d)
$g^{n+1}-g^{n}=g^{\prime}(r^{n})(r^{n+1}-r^{n}),$

with $P^{n}=\frac{\delta r}{\delta Q}(\phi^{n})$. In addition to the
regularity and the invertibility conditions that should be satisfied by $g$,
we need to assume the derivative of $g$ to be non-negative, namely,

(3.4) $g^{\prime}(r)\geq 0,$

for any $r\in\mathbb{R}$. This is necessary to ensure the modified energy
stability in discrete cases to hold, as stated below in Theorem 3.1. Before
that, we will discuss the procedure to implement this method first.

In a manner analogous to that of the IEC scheme, (3.3c), (3.3d) provide the
path representing $r^{n+1},g^{n+1}$ in terms of $\phi^{n+1}$. Consequently,
one can replace $r^{n+1},g^{n+1}$ in (3.3b), resulting in a linear scheme
pertaining to $\phi^{n+1}$ when considered in conjunction with (3.1a).
Alternatively, an update can be performed on the tuple
$(\phi^{n+1},\mu^{n+1},r^{n+1},g^{n+1})$ collectively by developing a linear
scheme to resolve (3.3) as a unified system. More precisely, at each
iteration, the following linear system could be solved:

(3.5) $\begin{pmatrix}\frac{1}{\Delta t}I&-\mathcal{G}&0&0\\\
\Delta&I&-\widetilde{{g}^{\prime}(r^{n}){P}^{n}}&-\tilde{P}^{n}\\\
-\tilde{P}^{n}&0&I&0\\\
0&0&-\widetilde{{g}^{\prime}(r^{n})}&I\end{pmatrix}\begin{pmatrix}\phi^{n+1}\\\
\mu^{n+1}\\\ r^{n+1}\\\ g^{n+1}\end{pmatrix}=\begin{pmatrix}\frac{1}{\Delta
t}\phi^{n}\\\ 0\\\ r^{n}-P^{n}\phi^{n}\\\
g^{n}-g^{\prime}(r^{n})r^{n}\end{pmatrix},$

given that $(\phi^{n},\mu^{n},r^{n},g^{n})$ is known. Here
$\widetilde{{g}^{\prime}(r^{n})}$ should also be interpreted as a
multiplication operator similar as $\tilde{P}^{n}$.

Now we turn to prove a modified energy dissipation law for the IEF scheme
which is a discretized version of (3.2).

######  Theorem 3.1.

The numerical scheme (3.3) is energy-stable. Specifically, define

(3.6) $\hat{E}^{n}=\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\int_{\Omega}g^{n}r^{n},$

then

$\hat{E}^{n+1}-\hat{E}^{n}\leq 0.$

###### Proof.

Taking inner product of (3.1a) with $\mu^{n+1}\Delta t$, (3.1b) with
$\phi^{n+1}-\phi^{n}$, (3.1c) with $g^{n+1}$, (3.1d), and summing them up with
(3.1d) together, we obtain,

$\displaystyle(\mathcal{G}\mu^{n+1},\mu^{n+1})\Delta t$
$\displaystyle=(\phi^{n+1}-\phi^{n},\mu^{n+1})$
$\displaystyle=\left(\phi^{n+1}-\phi^{n},-\Delta\phi^{n+1}+\left[r^{n+1}g^{\prime}(r^{n})+g^{n+1}\right]P^{n}\right)$
$\displaystyle=\frac{1}{2}\|\nabla\phi^{n+1}\|^{2}-\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\frac{1}{2}\|\nabla\phi^{n+1}-\nabla\phi^{n}\|^{2}$
$\displaystyle\quad\,\,+\int_{\Omega}(g^{n+1}-g^{n})\,r^{n+1}\,dx+\int_{\Omega}g^{n+1}(r^{n+1}-r^{n})\,dx$
$\displaystyle=\frac{1}{2}\|\nabla\phi^{n+1}\|^{2}-\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\frac{1}{2}\|\nabla\phi^{n+1}-\nabla\phi^{n}\|^{2}$
$\displaystyle\quad\,\,+\int_{\Omega}g^{n+1}r^{n+1}\,dx-\int_{\Omega}g^{n}r^{n}\,dx+\int_{\Omega}(g^{n+1}-g^{n})\,(r^{n+1}-r^{n})\,dx$
$\displaystyle=\frac{1}{2}\|\nabla\phi^{n+1}\|^{2}-\frac{1}{2}\|\nabla\phi^{n}\|^{2}+\frac{1}{2}\|\nabla\phi^{n+1}-\nabla\phi^{n}\|^{2}$
$\displaystyle\quad\,\,+\int_{\Omega}g^{n+1}r^{n+1}\,dx-\int_{\Omega}g^{n}r^{n}\,dx+\int_{\Omega}g^{\prime}(r^{n})(r^{n+1}-r^{n})^{2}\,dx$
$\displaystyle\geq\hat{E}^{n+1}-\hat{E}^{n},$

where we have used the assumption (3.4) to see that $g^{\prime}(r^{n})\geq 0$.
Using the fact that $(\mathcal{G}\mu^{n+1},\mu^{n+1})\leq 0$ independent of
choice of $\mu^{n+1}$, we have finished the proof. ∎

Similar as Corollary 2.4, we can also obtain an estimate for $\mu^{n+1}$ for
the IEF formulation.

######  Corollary 3.2.

Assume $\hat{E}^{0}$ is bounded. For fixed $N>0$, we have

$0\leq-\sum_{n=0}^{N}(\mathcal{G}\mu^{n},\mu^{n})\Delta
t\leq-\hat{E}^{N}+\hat{E}^{0}\leq\hat{E}^{0},$

and so it is bounded. Specifically, if $\mathcal{G}=\Delta$, we have

$0\leq\sum_{n=0}^{N}\|\nabla\mu^{n}\|^{2}\Delta t\leq\hat{E}^{0},$

if $\mathcal{G}=-I$, we have

$0\leq\sum_{n=0}^{N}\|\mu^{n}\|^{2}\Delta t\leq\hat{E}^{0}.$

###  3.3. An example of the IEF formulation

Similar as what we have done after establishing the IEC scheme, in this part,
we will provide a detailed example to illustrate the practicality of the IEF
method. A suitable choice is to take $g(r)=r^{2k+1}$, for $k=0,1,2,\cdots$. It
is easy to justify that $rg(r)=r^{2k+2}$ is invertible over $[0,\infty)$ and
assumption (3.4) is satisfied since $g^{\prime}(r)=(2k+1)r^{2k}\geq 0$. This
would result in

(3.7) $r(\phi)=\left(F(\phi)+A_{1}\right)^{\frac{1}{2k+2}}.$

Therefore, the numerical scheme (3.3) can be implemented as

(3.8a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (3.8b)
$\mu^{n+1}=-\Delta\phi^{n+1}+\left[(2k+1)(r^{n})^{2k}r^{n+1}+g^{n+1}\right]P^{n},$
(3.8c) $r^{n+1}-r^{n}=P^{n}(\phi^{n+1}-\phi^{n}),$ (3.8d)
$g^{n+1}-g^{n}=(2k+1)(r^{n})^{2k}(r^{n+1}-r^{n}),$

where $P^{n}=P(\phi^{n})$ and

(3.9) $P(\phi)=\frac{\delta
r}{\delta\phi}(\phi)=\frac{f(\phi)}{(2k+2)(F(\phi)+A_{1})^{\frac{2k+1}{2k+2}}}.$

This procedure can be summarized as the following Algorithm 2.

Algorithm 2 IEF scheme induced by monomial functions

1:procedure

2: Set the step size $\Delta t$

3: Set the total number of iterations $N$

4: Set the initial value $\phi^{0}=\phi_{0}$

5: Compute $r^{0}=r(\phi^{0})$ using equation (3.7)

6: Compute $g^{0}=(r^{0})^{2k+1}$

7: Compute $P^{0}=P(\phi^{0})$ using equation (3.9)

8: for $n=0,1,2,...,N-1$ do

9: Construct the matrix $A=\begin{pmatrix}\frac{1}{\Delta
t}I&-\mathcal{G}&0&0\\\
\Delta&I&-(2k+1)\widetilde{(r^{n})^{2k}P^{n}}&-P^{n}\\\
-\tilde{P}^{n}&0&I&0\\\ 0&0&-(2k+1)\widetilde{(r^{n})^{2k}}&I\end{pmatrix}$
using $P^{n}$

10: Construct the right-hand side $b=\begin{pmatrix}\frac{1}{\Delta
t}\phi^{n}\\\ 0\\\ r^{n}-P^{n}\phi^{n}\\\
g^{n}-(2k+1)(r^{n})^{2k}r^{n}\end{pmatrix}$

11: Update $\phi^{n+1}$, $r^{n+1},g^{n+1}$ by solving the linear system
$A\begin{pmatrix}\phi^{n+1}\\\ \mu^{n+1}\\\ r^{n+1}\\\ g^{n+1}\end{pmatrix}=b$

12: Update $P^{n+1}=P(\phi^{n+1})$ using equation (3.9)

13: end for

14:end procedure

It merits attention that by simply assigning a value of $k=0$, which
consequently yields $g(r)=r$, we can reconstruct the IEQ method. It is
because, in this scenario, the computation $g(r)r=r^{2}$ holds, and the update
rule given by (3.8d) is equivalent to (3.8c). This enforces $g^{n+1}=r^{n+1}$
for each iterative step $n=0,1,2,\ldots,N-1.$ Such a formulation endorses our
claim that the IEF is genuinely an extension of the IEQ method.

When choosing a suitable function $g$ for the IEF scheme, our selection should
not be confined to monomial functions, provided they adhere to the
abovementioned conditions. However, finding functions other than the monomials
that result in an explicit formula for $r$ in terms of $F(\phi)$ is
challenging. In such instances, one might need to employ root-finding
techniques to determine the value of $P^{n}$ once $\phi^{n}$ is known. This
inherent challenge positions our method as an Invariant Energy Monomial method
in the examples provided here and in subsequent numerical instances. Even so,
should an efficient technique be found for computing the root of the equation
$g(r)r=F(\phi)+A_{1}$ for a specific form of $g$, constructing an IEF scheme
using our approach would not pose a significant challenge. This potential
advancement forms an intriguing aspect for future exploration of this method.

##  4\. Numerical Experiments

In this section, we conduct a series of numerical experiments, focusing on the
Allen-Cahn [2] and Cahn-Hilliard equations [9, 10] under the application of
periodic boundary conditions in a two-dimensional domain. The equations above
are commonplace test cases for gradient flow-targeted numerical algorithms, as
illustrated in [8, 17, 22, 40, 41, 51]. Moreover, these equations are
particular manifestations of general gradient flows (1.1) and can be
formulated using the free energy:

$E(\phi)=\int_{\Omega}\left(\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+F(\phi)\right)\,dx,$

where

(4.1) $F(\phi)=\frac{1}{4}\left(\phi^{2}-1\right)^{2}.$

Here $\epsilon$ symbolizes a parameter that introduces stiffness issues into
the PDE system when $\varepsilon\ll 1$ [47].

Employing the variational approach to $E(\phi)$ in $L^{2}$, the Allen-Cahn
equation materializes as:

(4.2a) $\phi_{t}=-M\mu,$ (4.2b) $\mu=-\varepsilon^{2}\Delta\phi+f(\phi).$

Here $M$ is the mobility constant, and $f(\phi)=\frac{\delta
F}{\delta\phi}=\phi(\phi^{2}-1)$ [11].

In a similar manner, applying the variational approach to $E(\phi)$ within the
$H^{-1}$ space yields the Cahn-Hilliard equation:

(4.3a) $\phi_{t}=M\Delta\mu,$ (4.3b) $\mu=-\varepsilon^{2}\Delta\phi+f(\phi).$

We choose a computational domain defined by a square
$\Omega=[0,2\pi]\times[0,2\pi]$ for the subsequent experiments. We will
implement a standard finite difference method for this problem with periodic
boundary conditions to facilitate full discretization. Unless otherwise
stated, the domain will be discretized utilizing a $40\times 40$ grid, and we
will set the parameters $M=0.6$ and $\varepsilon=0.4$.

###  4.1. Accuracy test

In this section, we start our numerical experiments by assessing the
convergence rates of the Invariant Energy Convexification (IEC) scheme in its
application to the Allen–Cahn equation, as given by (4.2). For these
experiments, we adopt

(4.4) $\phi(x,y,t)=\sin(x)\cos(y)\cos(t)$

as the exact solution and include an appropriate right-hand side force field
term to ensure that the designated solution complies with the system (4.2). We
draw attention to the fact that the function $\phi$ employed here aligns with
the choice made for the accuracy test in [47].

In the context of the IEC scheme (2.4), we examine three distinct choices for
the function $c(r)$, namely, ${\ln(1+e^{r}),\ln(r)^{2},r^{2}}$. Each of these
functions are L-smooth and convex, and an elementary analysis would suggest
that $L=2$ serves as an appropriate choice for all three. Maintaining
consistency, we set $\alpha=0.5$ for our tests.

We then measure the $L^{2}$ errors for the variable $\phi$, contrasting our
numerical solution with the exact solution at $T_{end}=1$ under varying time
step sizes. As visualized in Figure 1a, each of the three choices of convex
functions demonstrates first-order accuracy. The error is linearly decreasing
when the time step size is decreasing as well. More specific error data for
varying time step sizes can be found in Table 1.

A noteworthy observation emerges when keeping the parameters $\alpha$ and $L$
constant in the numerical scheme: errors vary depending on the choice of the
corresponding convex function. Specifically, $c(r)=\ln(r)^{2}$ appears to
yield the best performance, while $c(r)=r^{2}$ lags behind. This finding
underscores the likelihood that, for a particular problem, a certain function
may appear to be more suitable to the construction of the numerical scheme in
contrast to a quadratic function to get better precision.

a

b

Figure 1. The $L^{2}$ numerical errors at $T_{end}=1$ for $\phi$ under different formulation. (a) Illustrates the numerical errors of the IEC scheme using different choices of convex functions. (b) Demonstrates the numerical errors of the IEF scheme, where varying monomials are utilized as the multiplying functions. Table 1. The $L^{2}$ numerical errors at $T_{end}=1$ for $\phi$ with three different convex functions at different time step size Function Time step | $\delta t=1.00\times 10^{-1}$ | $\delta t=5.00\times 10^{-2}$ | $\delta t=2.50\times 10^{-2}$ | $\delta t=1.25\times 10^{-2}$ | $\delta t=6.25\times 10^{-3}$ | $\delta t=3.12\times 10^{-3}$
---|---|---|---|---|---|---
$c(r)=\ln(1+e^{r})$ | $0.094\,215$ | $0.047\,625$ | $0.023\,729$ | $0.011\,624$ | $0.005\,532$ | $0.002\,479$
$c(r)=\ln(r)^{2}$ | $0.063\,632$ | $0.032\,332$ | $0.016\,092$ | $0.007\,814$ | $0.003\,635$ | $0.001\,541$
$c(r)=r^{2}$ | $0.114\,006$ | $0.057\,330$ | $0.028\,529$ | $0.014\,009$ | $0.006\,720$ | $0.003\,069$

We further seek to explore the impact of alterations in the values of $\alpha$
on the accuracy of the IEC scheme. With this objective in mind, we apply the
IEC schemes with $c(r)=\ln(1+e^{r})$ and $c(r)=\ln(r)^{2}$ using a range of
values for $\alpha$. We continue to assess the $L^{2}$ error between the
numerical solution and the exact solution at $T_{end}=1$. The outcomes of this
exploration are depicted in Figure 2. For this experiment, we select values
for $\alpha$ that range from its theoretical lower bound of $0.5$—a value
necessary to ensure the energy stability of the numerical scheme—up to $16$.
When the time step size is large, the performance of numerical schemes with
different values of $\alpha$ has obvious differences. In comparison, such a
difference becomes negligible when the time step size is small. Another
intriguing fact is that when the time step size is relatively large, an
increase in the value of $\alpha$ initially enhances the accuracy of the
corresponding numerical scheme but then reduces precision. When constructing
the numerical scheme using $c(r)=\ln(1+e^{r})$, we find that optimal accuracy
is achieved for values of $\alpha$ within the interval $[4,8]$. Contrarily,
for $c(r)=\ln(r)^{2}$, the scheme exhibits optimal performance when $\alpha$
lies within the range $[2,4]$. These findings suggest that different convex
functions may produce distinct optimal values for $\alpha$. As such, when we
use a relatively large time step, careful consideration of $\alpha$’s value is
necessary when implementing the IEC scheme.

a

b

Figure 2. The $L^{2}$ numerical errors at $T_{end}=1$ for $\phi$ with
different values of $\alpha$ under the IEC formulation. (a) Scheme constructed
with $c(r)=\ln{(1+e^{r})}$. (b) Scheme constructed with $c(r)=\ln(r)^{2}$.

We now discuss the experimental findings of the Invariant Energy
Functionalization (IEF) scheme. We will employ the scheme constructed using a
monomial as the auxiliary function $g$, a concept introduced in Section 3.3.
Our evaluation metric remains the $L^{2}$ error for $\phi$, calculated between
the numerical and exact solutions at $T_{end}=1$. The errors are illustrated
in Figure 1b, and the precise data is provided in Table 2. Insight into these
results reveals that the highest level of algorithmic precision is achieved
when $g(r)=r$, corresponding to the Invariant Energy Quadratization (IEQ)
method. Nevertheless, it is worth noting that the discrepancies between
different values of $k$ are tiny. Furthermore, regardless of the choice of
$k$, all options yield first-order accuracy. It implies that the norm of
errors diminishes linearly to the reduction in time step size. This outcome
further validates the robustness and utility of the IEF scheme in this
context.

Table 2. The $L^{2}$ numerical errors at $T_{end}=1$ for $g(r)=r^{k}$ with different values of $k$ $g(r)=r^{k}$ Time step | $\delta t=1.00\times 10^{-1}$ | $\delta t=5.00\times 10^{-2}$ | $\delta t=2.50\times 10^{-2}$ | $\delta t=1.25\times 10^{-2}$ | $\delta t=6.25\times 10^{-3}$ | $\delta t=3.12\times 10^{-3}$
---|---|---|---|---|---|---
$k=0$ | $0.115\,178\,529\,752\,356$ | $0.057\,731\,360\,592\,976\,0$ | $0.028\,558\,396\,507\,218\,5$ | $0.013\,857\,120\,629\,241\,7$ | $0.006\,478\,969\,492\,079\,82$ | $0.002\,787\,056\,598\,031\,50$
$k=1$ | $0.113\,698\,137\,982\,773$ | $0.057\,009\,266\,299\,733\,2$ | $0.028\,202\,337\,568\,806\,5$ | $0.013\,680\,565\,396\,667\,2$ | $0.006\,391\,301\,738\,899\,35$ | $0.002\,743\,746\,442\,533\,42$
$k=3$ | $0.114\,145\,797\,493\,452$ | $0.057\,226\,254\,754\,756\,1$ | $0.028\,309\,093\,017\,014\,1$ | $0.013\,733\,514\,409\,536\,0$ | $0.006\,417\,674\,266\,329\,56$ | $0.002\,756\,894\,527\,407\,49$
$k=5$ | $0.114\,469\,994\,858\,300$ | $0.057\,383\,723\,011\,471\,7$ | $0.028\,386\,621\,117\,126\,5$ | $0.013\,771\,962\,705\,311\,9$ | $0.006\,436\,803\,801\,024\,96$ | $0.002\,766\,401\,714\,214\,71$
$k=7$ | $0.114\,664\,869\,002\,037$ | $0.057\,478\,392\,735\,728\,9$ | $0.028\,433\,232\,743\,144\,9$ | $0.013\,795\,077\,585\,324\,3$ | $0.006\,448\,302\,461\,876\,41$ | $0.002\,772\,113\,944\,248\,55$

###  4.2. Energy stability

In our discussion thus far, we have highlighted a significant advantage shared
by the IEC and the IEF schemes, which they inherit from the Invariant Energy
Quadratization (IEQ) method: the preservation of the energy dissipation law.
In this section, we set out to empirically validate this property through a
series of numerical experiments. We will examine an example with the initial
condition given by:

(4.5) $\phi(x,y,0)=\sin(x)\cos(y),$

and scrutinize the evolution of energy over the time interval $[0,5]$ when
employing both the IEC and the IEF schemes to solve both the Allen-Cahn and
Cahn-Hilliard equations.

We commence with the IEC scheme. We maintain a consistent choice of $c(r)$ as
the Softplus function for these experiments, thereby employing Algorithm 1 for
our implementations. As depicted in Figure 3, the modified energy exhibits a
monotonic decrease independent of the time step size for both the Allen-Cahn
and Cahn-Hilliard equations. The definition of the modified energy, as
referenced here, aligns with equation (2.6). This outcome confirms preserving
the energy dissipation property within the IEC scheme.

a

b

Figure 3.  Evolution of the modified energy over time under the IEC scheme
with $c(r)=\ln(1+e^{r})$: (a) Allen-Cahn equation; (b) Cahn-Hilliard equation.

Another research question is the difference between the modified and original
energies. Specifically, the original energy function $F(\phi)$, defined by
equation (4.1), is a function of the variable $\phi$, whereas its modified
counterpart, $c(r)$, is a function of $r$ in the IEC formulation. An
intriguing question is whether $c(r)$ would converge to $F(\phi)$ in the limit
as the time step size approaches zero. In more technical terms, we aim to
compute the quantity $\lvert\int_{\Omega}(c(r)-F(\phi))|$ at $T_{end}=5$ for
various time step sizes and observe if this quantity approaches zero as the
time step size diminishes.

To investigate this, we employ the same example with the initial condition
specified by equation (4.5), persisting with the use of the Softplus function
for constructing our IEC scheme. As demonstrated in Figure 4a, the difference
between the original energy and the modified energy indeed appears to approach
zero as the time step size shrinks, when tested on the Allen-Cahn equation,
with a first-order convergence rate. Furthermore, figure 4b validates the same
numerical results for the Cahn-Hilliard equation. This result offers
additional evidence in favor of implementing the change of variable procedure,
underscoring that the auxiliary variable computed by the numerical scheme will
converge to the original function it aims to represent when the time step size
is sufficiently small.

a

b

Figure 4.  Graphs depicting $\lvert\int_{\Omega}(c(r)-F(\phi))|$ at
$T_{end}=5$ for varying time step sizes under the IEC scheme, constructed
using the Softplus function: (a) Allen-Cahn equation; (b) Cahn-Hilliard
equation.

We next turn our attention to comparable experiments performed using the IEF
scheme. We shall confine our discussion to applying the IEF scheme to the
Cahn-Hilliard equation for brevity. Employing the same initial condition given
by equation (4.5), we choose $g(r)=r^{7}$ to facilitate our investigations.
The modified energy in this context is defined as per equation (3.6). Figure
5a delineates the evolution of this modified energy over time for varying time
step sizes. The modified energy exhibits a monotonous decrease, in alignment
with our theoretical predictions outlined in Theorem 3.1.

For our second experiment using the IEF scheme, we diverge slightly from the
approach used in the IEC scheme, where we computed the discrepancy between the
modified and original energies. Instead, we focus on the disparity between the
auxiliary variable and its corresponding value as a function of $\phi$ of
differing time step sizes. As illustrated in Figure 5b, the $L^{2}$ norm of
the difference between $\lvert g-g\left(r(\phi)\right)\rvert$ and $\lvert
r-r(\phi)\rvert$ at $T_{end}=5$ shows a linear dependency on the time step
size, again demonstrating first-order accuracy. The results of these two
experiments attest to the robustness and accuracy of the IEF scheme.

a Evolution of the modified energy over time under the IEF scheme for Cahn-
Hilliard equation, constructed using $g(r)=r^{7}$

b Graph depicting $L^{2}$ norm of $\lvert g-g\left(r(\phi)\right)\rvert$ and
$\lvert r-r(\phi)\rvert$ at $T_{end}=5$ for varying time step sizes under the
IEF scheme for Cahn-Hilliard equation, constructed using $g(r)=r^{7}$

Figure 5.  Numerical results of the IEF scheme applied to the Cahn-Hilliard
equation

###  4.3. Coarsening effect

The coarsening effect [16, 29] refers to the phenomenon where small-scale
structures or features in the system tend to merge and form larger-scale
structures over time. This effect is related to the system’s dynamics of phase
separation or pattern formation. Here we will present some experimental
results by applying the IEC and the IEF scheme to the Cahn-Hilliard equation
to simulate such an effect. We firstly consider a benchmark problem for the
corasening effect, for example, see [11, 33, 47]. The initial condition is set
to be

(4.6)
$\phi(x,y,0)=-\sum_{i=1}^{2}\tanh\left(\frac{\sqrt{(x-x_{i})^{2}+(y-y_{i})^{2}}-r_{i}}{1.2\varepsilon}\right)+1,$

where we choose $(x_{1},y_{1},r_{1})=(\pi-0.7,\pi-0.6,)$ and
$(x_{2},y_{2},r_{2})=(\pi+1.65,\pi+1.6,0.8)$. We apply the IEC scheme
formulated by the Softplus function, which is Algorithm 1, to solve this
problem with times step size set to be $0.001$. As we can see from the
snapshots presented in Figure 6, we observe the coarsening effect as the small
circle is generally absorbed by the large circle over the time span $[0,3]$.
The fully absorption happens at around $t=2.00s$.

$t=0.25s$

$t=0.50s$

$t=0.75s$

$t=1.00s$

$t=1.25s$

$t=1.50s$

$t=1.75s$

$t=2.00s$

$t=2.25s$

$t=2.50s$

$t=2.75s$

$t=3.00s$

Figure 6. Snapshots of $\phi$ taken within the time span $[0,3]$, solving the
Cahn-Hilliard equation with initial condition (4.6) by the IEC scheme
formulated by the Softplus function with time step size set to be $0.001$.

Our last exploration entails a new example, which we introduce with the
following initial condition:

(4.7) $\phi(x,y,0)=0.25+0.4\text{rand}(x,y).$

This numerical experiment reproduces a similar one conducted in [32]. We
employ the IEF scheme with the function $g(r)$ set as $r^{7}$ for our
investigations. The time step size adopted is $0.0001$. The dynamic patterns
of the mixing of the phase field, starting from a random initial value, can be
discerned in Figure 7. The observed results, embodying the intriguing
interplay of phase transitions, comply nicely with empirical findings.

a t = 0s

b t = 0.1s

c t = 0.2s

d t = 0.3s

e t = 0.4s

f t = 0.5s

g t = 0.6s

h t = 0.7s

i t = 0.8s

j t = 0.9s

k t = 1.0s

l t = 1.1s

Figure 7. Snapshots of $\phi$ taken within the time span $[0,1.1]$, solving
the Cahn-Hilliard equation with initial condition (4.7) by the IEF scheme
formulated by the $g(r)=r^{7}$ with time step size set to be $0.0001$.

##  5\. Extension to the SAV scheme

In the forthcoming discussion, we present the extensions of the frameworks of
the IEC and the IEF scheme onto the Scalar Auxiliary Variable method (SAV).
These schemes will share similar properties as our generalization of the IEQ
method in the sense that they will also lead to energy-stable linear numerical
schemes.

###  5.1. A brief review of the SAV method

The Scalar Auxiliary Variable (SAV) method, initially presented in [40, 41]
and subsequently refined [1, 5, 20, 28, 39, 42], is a robust tool analogous to
the IEQ method. It is widely employed to construct efficient and accurate time
discretization schemes for a broad spectrum of gradient flows. In terms of
formulation, it takes a notable resemblance to the IEQ method. Both introduce
an auxiliary variable to replace the original energy and then construct an
energy-stable linear scheme based on this reformulation. A distinct
differential point emerges in the scope of the auxiliary variables utilized by
the two methodologies. While the IEQ approach employs a variable over both
time and space (represented as $r=\sqrt{F(\phi)+A_{1}}$ where $F$ is the
energy density function outlined in (1.1)), the SAV method utilizes a variable
defined solely over time ($r=\sqrt{E_{1}(\phi)+A_{2}}$), represented as the
square root of $E_{1}(\phi)$, which is defined as
$E_{1}(\phi)=\int_{\Omega}F(\phi)\,dx$. $A_{2}$ here is a constant to ensure
the auxiliary variable is well-defined given that $E_{1}$ is bounded from
below, which is necessary for the free energy term to be physically sound.
Specifically, as initially outlined in [41], a first-order SAV scheme can be
carried out as;

(5.1a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (5.1b)
$\mu^{n+1}=-\Delta\phi^{n+1}+\frac{r^{n+1}}{\sqrt{E_{1}(\phi^{n})}+A_{2}}F^{n},$
(5.1c)
$r^{n+1}-r^{n}=\frac{1}{c^{\prime}\left(r(\phi)\right)}\int_{\Omega}F^{n}(\phi^{n+1}-\phi^{n})\,dx.$

We take $F^{n}=F(\phi^{n})$ here. Various numerical experiments have shown the
efficiency and practicability of this method. In addition, as it has been
stated in [40, 41], the SAV methodology inherits all the advantages of the IEQ
method and additionally enjoys the benefits of a simpler implementation and
the requirement for only $E_{1}(\phi)$ to be bounded from below (instead of
the need for $F(\phi)$ also to be bounded from below). Given the commonality
in the core construction principles of these methods, it is feasible to
transpose our ideas onto the SAV method, enabling further generalizations. We
will present it in the following.

###  5.2. The C-SAV method

In the following, we will propose a variant of the standard SAV method, the
Convex Scalar Auxiliary Variable method, abbreviated as C-SAV. Guided by the
parallel adaptations we made for the IEC method, a natural proposition is to
supplant $r^{2}$ with an alternative auxiliary function. In this context, we
will adhere to the same assumptions we have put forth for the IEC scheme with
respect to the auxiliary function $c$; namely, $c:\mathbb{R}\to\mathbb{R}$ is
taken as a smooth, L-smooth, and convex function exhibiting a monotonic
increase, with $\mathbb{R}^{+}$ included within the function’s range. Then by
introducing a scalar auxiliary variable $r(t)$ such that

$c\left(r(t)\right)=E_{1}(\phi)+A_{2},$

system (1.1) can be reformulated as

(5.2a) $\frac{\partial\phi}{\partial t}=\mathcal{G}\mu,$ (5.2b)
$\mu=-\Delta\phi+\frac{c^{\prime}\left(r(t)\right)}{c^{\prime}\left(c^{-1}\left(E_{1}(\phi)+A_{2}\right)\right)}\frac{\delta
F}{\delta\phi},$ (5.2c)
$r_{t}=\frac{1}{c^{\prime}\left(c^{-1}\left(E_{1}(\phi)+A_{2}\right)\right)}\int_{\Omega}\frac{\delta
F}{\delta\phi}\phi_{t},$

From this, we propose the following variant of the IEC scheme, through the
application of the same conceptual framework to the SAV method:

(5.3a) $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}=\mathcal{G}\mu^{n+1},$ (5.3b)
$\mu^{n+1}=-\Delta\phi^{n+1}+\frac{c^{\prime}(r^{n})+\alpha
L(r^{n+1}-r^{n})}{c^{\prime}\left(r(\phi^{n})\right)}F^{n},$ (5.3c)
$r^{n+1}-r^{n}=\frac{1}{c^{\prime}\left(r\left(\phi^{n}\right)\right)}\int_{\Omega}F^{n}(\phi^{n+1}-\phi^{n})\,dx.$

By defining modified energy as
$E_{CS}^{n}=\frac{1}{2}|\nabla\phi^{n}|^{2}+c(r^{n})$, we can deduce, through
parallel arguments to Theorem 2.2, that this discrete energy demonstrates a
monotonic decrease. We can similarly extend the IEF formulation to the SAV
methodology, thereby constructing linear, energy-stable numerical schemes. We
omit the details for the sake of brevity, but we encourage interested readers
to explore this subject further.

##  6\. Concluding remarks

In this work, we have introduced two novel numerical methodologies: the
Invariant Energy Convexification (IEC) and the Invariant Energy
Functionalization (IEF) methods. These methods generalize the well-known
Invariant Energy Quadratization (IEQ) method. The inherent ability of these
proposed techniques to construct linear energy-stable numerical schemes, a
crucial feature inherited from the IEQ method, underscores their significance.

The numerical experiments on the Allen-Cahn and Cahn-Hilliard equations
substantiated the theoretical properties deduced for both methods in practical
implementations. The results conclusively demonstrated that both methods are
robust and provide energy stability. Moreover, if an appropriate function is
meticulously selected and carefully chooses the parameters in the method, the
IEC method exhibits the potential to outperform the standard IEQ method.

As we conclude this manuscript, we wish to underscore several notable points
and delineate potential lanes for future research related to this work:

  1. (1)

In the extant version of the IEC method, the selection of uniformly L-smooth
functions, such as the Softplus function, is imperative for designing the
corresponding numerical scheme in theory. However, in our practical
implementations, we observed that numerous functions, which are only locally
L-smooth, can also facilitate the development of numerical schemes with a
commendable performance. It could be attributed to the energy bound of the
auxiliary variable, which may consequently result in a bound of the function’s
second derivative. The necessity for theoretical justification of this aspect
needs to be considered.

  2. (2)

The accuracy of the numerical solutions is influenced by choice of the
auxiliary function $c(x)$ and the parameter $\alpha$, especially when the time
discretization is coarse. It is evident from our discussion in Section 4.
Therefore, deriving a theoretical justification to find the optimal
formulation for specific problems is challenging and immensely useful.

  3. (3)

Our numerical results reveal that both methods lead to convergence to the
exact solution with a desirable convergence rate. However, contrary to the IEQ
method, the convergence analysis of the IEC and the IEF methods may present
formidable challenges due to the non-quadratic structure of their modified
energy. It necessitates using specific analytical tools for the auxiliary
functions used in their formulation—this facet of our work beckons further
investigation.

  4. (4)

It is worth noting that the challenges and problems mentioned above are also
relevant when considering extending these methods to the Scalar Auxiliary
Variable (SAV) method.

## References

  * [1] G. Akrivis, B. Li, and D. Li.  Energy-decaying extrapolated rk–sav methods for the allen–cahn and cahn–hilliard equations.  SIAM Journal on Scientific Computing, 41(6):A3703–A3727, 2019.
  * [2] S. M. Allen and J. W. Cahn.  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening.  Acta metallurgica, 27(6):1085–1095, 1979.
  * [3] L. Ambrosio, N. Gigli, and G. Savaré.  Gradient flows: in metric spaces and in the space of probability measures.  Springer Science & Business Media, 2005.
  * [4] D. M. Anderson, G. B. McFadden, and A. A. Wheeler.  Diffuse-interface methods in fluid mechanics.  Annual Review of Fluid Mechanics, 30(1):139–165, 1998.
  * [5] X. Antoine, J. Shen, and Q. Tang.  Scalar auxiliary variable/lagrange multiplier based pseudospectral schemes for the dynamics of nonlinear schrödinger/gross-pitaevskii equations.  Journal of Computational Physics, 437:110328, 2021.
  * [6] A. Baskaran, J. S. Lowengrub, C. Wang, and S. M. Wise.  Convergence analysis of a second order convex splitting scheme for the modified phase field crystal equation.  SIAM Journal on Numerical Analysis, 51(5):2851–2873, 2013.
  * [7] S. Bilbao, M. Ducceschi, and F. Zama.  Explicit exactly energy-conserving methods for hamiltonian systems.  Journal of Computational Physics, 472:111697, 2023.
  * [8] S. C. Brenner, A. E. Diegel, and L.-Y. Sung.  A robust solver for a mixed finite element method for the cahn–hilliard equation.  Journal of Scientific Computing, 77:1234–1249, 2018.
  * [9] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. iii. nucleation in a two-component incompressible fluid.  The Journal of chemical physics, 31(3):688–699, 1959.
  * [10] J. W. Cahn and J. E. Hilliard.  Free Energy of a Nonuniform System. I. Interfacial Free Energy.  The Journal of Chemical Physics, 28(2):258–267, 08 2004.
  * [11] C. Chen and X. Yang.  Fast, provably unconditionally energy stable, and second-order accurate algorithms for the anisotropic cahn–hilliard model.  Computer Methods in Applied Mechanics and Engineering, 351:35–59, 2019.
  * [12] R. Chen, X. Yang, and H. Zhang.  Second order, linear, and unconditionally energy stable schemes for a hydrodynamic model of smectic-a liquid crystals.  SIAM Journal on Scientific Computing, 39(6):A2808–A2833, 2017.
  * [13] L. Chizat and F. Bach.  On the global convergence of gradient descent for over-parameterized models using optimal transport.  Advances in neural information processing systems, 31, 2018.
  * [14] A. Christlieb, J. Jones, K. Promislow, B. Wetton, and M. Willoughby.  High accuracy solutions to energy gradient flows from material science models.  Journal of computational physics, 257:193–215, 2014.
  * [15] J. Cortés.  Finite-time convergent gradient flows with applications to network consensus.  Automatica, 42(11):1993–2000, 2006.
  * [16] S. Dai and Q. Du.  Computational studies of coarsening rates for the cahn–hilliard equation with phase-dependent diffusion mobility.  Journal of Computational Physics, 310:85–108, 2016.
  * [17] A. E. Diegel, C. Wang, and S. M. Wise.  Stability and convergence of a second-order mixed finite element method for the cahn–hilliard equation.  IMA Journal of Numerical Analysis, 36(4):1867–1897, 2016.
  * [18] Z. Ding, S. Chen, Q. Li, and S. J. Wright.  Overparameterization of deep resnet: zero loss and mean-field analysis.  Journal of machine learning research, 2022.
  * [19] J. Ericksen.  Equilibrium theory of liquid crystals.  In Advances in liquid crystals, volume 2, pages 233–298. Elsevier, 1976.
  * [20] G. Fu and D. Han.  A linear second-order in time unconditionally energy stable finite element scheme for a cahn–hilliard phase-field model for two-phase incompressible flow of variable densities.  Computer Methods in Applied Mechanics and Engineering, 387:114186, 2021.
  * [21] Y. Fu, W. Cai, and Y. Wang.  A linearly implicit structure-preserving scheme for the fractional sine-gordon equation based on the ieq approach.  Applied Numerical Mathematics, 160:368–385, 2021.
  * [22] Y. Gao, D. Han, X. He, and U. Rüde.  Unconditionally stable numerical methods for cahn-hilliard-navier-stokes-darcy system with different densities and viscosities.  Journal of Computational Physics, 454:110968, 2022.
  * [23] K. Glasner and S. Orizaga.  Improving the accuracy of convexity splitting methods for gradient flow equations.  Journal of Computational Physics, 315:52–64, 2016.
  * [24] V. M. Gudibanda, F. Weber, and Y. Yue.  Convergence analysis of a fully discrete energy-stable numerical scheme for the Q-tensor flow of liquid crystals.  SIAM Journal on Numerical Analysis, 60(4):2150–2181, 2022.
  * [25] F. Guillén-González, M. Rodríguez-Bellido, and D. A. Rueda-Gómez.  Unconditionally energy stable fully discrete schemes for a chemo-repulsion model.  Mathematics of Computation, 88(319):2069–2099, 2019.
  * [26] D. Han and N. Jiang.  A second order, linear, unconditionally stable, crank–nicolson–leapfrog scheme for phase field models of two-phase incompressible flows.  Applied Mathematics Letters, 108:106521, 2020.
  * [27] C. Jiang, W. Cai, and Y. Wang.  A linearly implicit and local energy-preserving scheme for the sine-gordon equation based on the invariant energy quadratization approach.  Journal of Scientific Computing, 80:1629–1655, 2019.
  * [28] M. Jiang, Z. Zhang, and J. Zhao.  Improving the accuracy and consistency of the scalar auxiliary variable (sav) method with relaxation.  Journal of Computational Physics, 456:110954, 2022.
  * [29] L. Ju, J. Zhang, and Q. Du.  Fast and accurate algorithms for simulating coarsening dynamics of cahn–hilliard equations.  Computational Materials Science, 108:272–282, 2015.
  * [30] B. Li and J.-g. Liu.  Epitaxial growth without slope selection: Energetics, coarsening, and dynamic scaling.  Journal of Nonlinear Science, 14:429–451, 01 2004.
  * [31] C. Liu and N. J. Walkington.  Approximation of liquid crystal flows.  SIAM Journal on Numerical Analysis, 37(3):725–741, 2000.
  * [32] Z. Liu and X. Li.  Efficient modified techniques of invariant energy quadratization approach for gradient flows.  Applied Mathematics Letters, 98:206–214, 2019.
  * [33] Z. Liu and X. Li.  Step-by-step solving schemes based on scalar auxiliary variable and invariant energy quadratization approaches for gradient flows.  Numerical Algorithms, pages 1–22, 2022.
  * [34] H. Lu, R. M. Freund, and Y. Nesterov.  Relatively smooth convex optimization by first-order methods, and applications.  SIAM Journal on Optimization, 28(1):333–354, 2018.
  * [35] J. Maas.  Gradient flows of the entropy for finite markov chains.  Journal of Functional Analysis, 261(8):2250–2292, 2011.
  * [36] A. Mielke and U. Stefanelli.  Weighted energy-dissipation functionals for gradient flows.  ESAIM: Control, Optimisation and Calculus of Variations, 17(1):52–85, 2011.
  * [37] F. Santambrogio.  Euclidean, metric, and Wasserstein gradient flows: an overview.  Bulletin of Mathematical Sciences, 7:87–154, 2017.
  * [38] J. Shen, C. Wang, X. Wang, and S. M. Wise.  Second-order convex splitting schemes for gradient flows with Ehrlich-Schwoebel type energy: application to thin film epitaxy.  SIAM Journal on Numerical Analysis, 50(1):105–125, 2012.
  * [39] J. Shen and J. Xu.  Convergence and error analysis for the scalar auxiliary variable (sav) schemes to gradient flows.  SIAM Journal on Numerical Analysis, 56(5):2895–2912, 2018.
  * [40] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  Journal of Computational Physics, 353:407–416, 2018.
  * [41] J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  SIAM Review, 61(3):474–506, 2019.
  * [42] J. Shen and X. Yang.  The IEQ and SAV approaches and their extensions for a class of highly nonlinear gradient flow systems.  Contemp. Math., 754:217–245, 2020.
  * [43] J. Shin, H. G. Lee, and J.-Y. Lee.  Unconditionally stable methods for gradient flow using convex splitting Runge-Kutta scheme.  Journal of Computational Physics, 347:367–381, 2017.
  * [44] T. Szandała.  Review and comparison of commonly used activation functions for deep neural networks.  Bio-inspired neurocomputing, pages 203–224, 2021.
  * [45] I. Vardoulakis and E. C. Aifantis.  A gradient flow theory of plasticity for granular materials.  Acta mechanica, 87(3-4):197–217, 1991.
  * [46] X. Yang and L. Ju.  Efficient linear schemes with unconditional energy stability for the phase field elastic bending energy model.  Computer Methods in Applied Mechanics and Engineering, 315:691–712, 2017.
  * [47] X. Yang and G.-D. Zhang.  Convergence analysis for the invariant energy quadratization (IEQ) schemes for solving the Cahn-Hilliard and Allen-Cahn equations with general nonlinear potential.  Journal of scientific computing, 82:1–28, 2020.
  * [48] X. Yang, J. Zhao, and Q. Wang.  Numerical approximations for the molecular beam epitaxial growth model based on the invariant energy quadratization method.  Journal of Computational Physics, 333:104–127, 2017.
  * [49] X. Yang, J. Zhao, Q. Wang, and J. Shen.  Numerical approximations for a three-component Cahn-Hilliard phase-field model based on the invariant energy quadratization method.  Mathematical Models and Methods in Applied Sciences, 27(11):1993–2030, 2017.
  * [50] Y. Yue.  On the convergence of an unconditionally stable numerical scheme for the q-tensor flow based on the invariant quadratization method.  Applied Mathematics Letters, 138:108522, 2023.
  * [51] J. Zhao and D. Han.  Second-order decoupled energy-stable schemes for cahn-hilliard-navier-stokes equations.  Journal of Computational Physics, 443:110536, 2021.
  * [52] J. Zhao, Q. Wang, and X. Yang.  Numerical approximations for a phase field dendritic crystal growth model based on the invariant energy quadratization approach.  International Journal for Numerical Methods in Engineering, 110(3):279–300, 2017.
  * [53] J. Zhao, X. Yang, Y. Gong, and Q. Wang.  A novel linear second order unconditionally energy stable scheme for a hydrodynamic q-tensor model of liquid crystals.  Computer Methods in Applied Mechanics and Engineering, 318:803–825, 2017.