Efficient and accurate exponential SAV algorithms with relaxation for dissipative system

TL;DR

This paper introduces two linear, unconditionally energy-stable exponential SAV schemes with relaxation for dissipative systems, including Navier-Stokes equations, ensuring positivity and high-order accuracy with demonstrated effectiveness.

Contribution

The paper develops novel R-ESAV schemes that improve positivity preservation and stability over existing methods, and extends high-order BDF schemes for complex dissipative systems.

Findings

Schemes are linear and unconditionally energy stable.

Numerical examples confirm high accuracy and effectiveness.

Approaches preserve positivity without additional assumptions.

Abstract

In this paper, we construct two kinds of exponential SAV approach with relaxation (R-ESAV) for dissipative system. The constructed schemes are linear and unconditionally energy stable. They can guarantee the positive property of SAV without any assumption compared with R-SAV and R-GSAV approaches, preserve all the advantages of the ESAV approach and satiesfy dissipation law with respect to a modified energy which is directly related to the original free energy. Moreover the second version of R-ESAV approach is easy to construct high-order BDF$k$ schemes. Especially for Navier-Stokes equations, we construct wo kinds of novel schemes based on the R-ESAV method. Finally, ample numerical examples are presented to exhibit that the proposed approaches are accurate and effective.

Tables

Table 1. Table 1: Example 1 (ii). A comparison of L 2 superscript 𝐿 2 L^{2} -error obtained by four approaches based on BDF 2 2 2 scheme for Allen-Cahn equation at T = 200 𝑇 200 T=200 with various time steps.

	ESAV-1	R-ESAV-1	ESAV-2	R-ESAV-2
1E-1	9.49E-04	2.76E-04	0.98	1.01E-04
5E-2	8.43E-04	7.22E-05	1.07E-03	2.74E-05
1E-2	1.18E-04	2.94E-06	4.55E-05	1.01E-06
5E-3	3.53E-05	7.36E-07	2.29E-05	2.49E-07
1E-3	1.63E-06	2.96E-08	1.47E-06	9.88E-09

Table 2. Table 2: Example 5 (i). Convergence test for Navier-Stokes equation using Scheme I/BDF 1 1 1 with and without relaxation.

		ESAV-2				R-ESAV-2
$\delta t$	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate
2.50E-2	4.29E-02	–	1.40	–	1.08E-02	–	2.62E-01	–
1.25E-2	1.96E-02	1.13	6.46E-01	1.11	4.85E-03	1.15	1.15E-01	1.19
6.25E-3	9.40E-03	1.06	3.10E-01	1.06	2.32E-03	1.06	5.35E-02	1.10
3.13E-3	4.61E-03	1.03	1.52E-01	1.03	1.14E-03	1.03	2.58E-02	1.05
1.56E-3	2.28E-03	1.01	7.52E-02	1.01	5.63E-04	1.01	1.26E-02	1.03

Table 3. Table 3: Example 5 (i). Convergence test for Navier-Stokes equation using Scheme I/BDF 2 2 2 with and without relaxation.

		ESAV-2				R-ESAV-2
$\delta t$	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate
2.50E-2	1.94E-03	–	4.47E-02	–	1.54E-03	–	9.25E-03	–
1.25E-2	4.64E-04	2.06	9.96E-03	2.17	3.88E-04	1.99	2.77E-03	1.74
6.25E-3	1.14E-04	2.02	2.40E-03	2.05	9.73E-05	1.99	8.79E-04	1.66
3.13E-3	2.84E-05	2.00	6.19E-04	1.96	2.46E-05	1.99	3.03E-04	1.54
1.56E-3	7.28E-06	1.97	1.70E-04	1.86	6.37E-06	1.95	1.07E-04	1.50

Table 4. Table 4: Example 5 (i). Convergence test for Navier-Stokes equation using Scheme II/BDF 1 1 1 with and without relaxation.

		ESAV-2				R-ESAV-2
$\delta t$	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate
2.50E-2	4.86E-01	–	9.34E-01	–	1.99E-01	–	2.89E-01	–
1.25E-2	3.12E-01	0.64	6.03E-01	0.63	9.51E-02	1.07	1.34E-01	1.11
6.25E-3	1.86E-01	0.75	3.58E-01	0.75	4.60E-02	1.05	6.37E-02	1.07
3.13E-3	1.03E-01	0.85	1.99E-01	0.85	2.26E-02	1.03	3.09E-02	1.04
1.56E-3	5.47E-02	0.91	1.05E-01	0.92	1.15E-02	1.02	1.52E-02	1.02

Table 5. Table 5: Example 5 (i). Convergence test for Navier-Stokes equation using Scheme II/BDF 2 2 2 with and without relaxation.

		ESAV-2				R-ESAV-2
$\delta t$	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate	${\Vert e_{𝐮}\Vert }_{L^2}$	Rate	${\Vert e_p\Vert }_{L^2}$	Rate
2.50E-2	3.83E-02	–	6.23E-02	–	8.25E-03	–	1.59E-02	–
1.25E-2	7.60E-03	2.33	1.30E-02	2.26	2.05E-03	2.01	3.96E-03	2.00
6.25E-3	1.73E-03	2.14	2.98E-03	2.12	5.15E-04	2.00	9.93E-04	1.99
3.13E-3	4.13E-04	2.06	7.17E-04	2.06	1.29E-04	2.00	2.49E-04	2.00
1.56E-3	1.01E-04	2.03	1.76E-04	2.03	2.23E-05	2.00	6.24E-05	2.00

Equations

$$E(\varphi )=\frac{1}{2}(\mathcal{L}\varphi ,\varphi )+{\int }_{\Omega }F(\varphi )\mathrm{d}\mathbf{x},$$

$$\{\begin{array}{rl}& \frac{\partial \varphi }{\partial t}=-\mathcal{G}\mu ,\\ & \mu =\frac{\delta E(\varphi )}{\delta \varphi }=\mathcal{L}\varphi +F^{\prime }(\varphi ),\end{array}$$

$$\frac{\mathrm{d}E(\varphi )}{\mathrm{d}t}=\left(\frac{\delta E(\varphi )}{\delta \varphi },\frac{\partial \varphi }{\partial t}\right)=-(\mathcal{G}\mu ,\mu )\le 0.$$

$$\frac{\mathrm{d}{H}^{-1}(r)}{\mathrm{d}t}={\int }_{\Omega }{F}^{\prime }(\varphi ){\varphi }_t\mathrm{d}\mathbf{x},$$

$$\{\begin{array}{rl}& {\partial }_t\varphi =-\mathcal{G}\mu ,\\ & \mu =\mathcal{L}\varphi +\frac{r}{H\left({\int }_{\Omega }F(\varphi )\mathrm{d}\mathbf{x}\right)}{F}^{\prime }(\varphi ),\\ & \frac{\mathrm{d}}{\mathrm{d}t}{H}^{-1}(r)=\frac{r}{H\left({\int }_{\Omega }F(\varphi )\mathrm{d}\mathbf{x}\right)}\left(F^{\prime }(\varphi ),{\varphi }_t\right).\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left\{\frac{1}{2}(\mathcal{L}\varphi ,\varphi )+H^{-1}(r)\right\}=\frac{\mathrm{d}E(\varphi )}{\mathrm{d}t}=-(\mathcal{G}\mu ,\mu ).$$

$$\frac{{\varphi }^{n+1}-{\varphi }^n}{\delta t}=-\mathcal{G}{\mu }^{n+1},$$

$$E(\varphi )=\frac{1}{2}(\mathcal{L}\varphi ,\varphi )+{\int }_{\Omega }F(\varphi )\mathrm{d}\mathbf{x},$$

$$\{\begin{array}{rl}& \frac{\partial \varphi }{\partial t}=-\mathcal{G}\mu ,\\ & \mu =\frac{\delta E(\varphi )}{\delta \varphi }=\mathcal{L}\varphi +F^{\prime }(\varphi ),\end{array}$$

$$r(t)=\exp (E_1(\varphi ))=\exp ({\int }_{\Omega }F(\varphi )\mathrm{d}\mathbf{x}),$$

$$\frac{\mathrm{d}r}{\mathrm{d}t}=r{\int }_{\Omega }{F}^{\prime }(\varphi ){\varphi }_t\mathrm{d}\mathbf{x}=\frac{r^2}{\exp \left(E_1(\varphi )\right)}{\int }_{\Omega }{F}^{\prime }(\varphi ){\varphi }_t\mathrm{d}\mathbf{x}.$$

$$\frac{\mathrm{d}\ln (r)}{\mathrm{d}t}=\left(U(r,\varphi ),{\varphi }_t\right),$$

$$\{\begin{array}{rl}& \frac{\partial \varphi }{\partial t}=-\mathcal{G}\mu ,\\ & \mu =\mathcal{L}\varphi +U(r,\varphi ),\\ & \frac{\mathrm{d}\ln (r)}{\mathrm{d}t}=\left(U(r,\varphi ),{\varphi }_t\right).\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{2}(\mathcal{L}\varphi ,\varphi )+\ln (r)\right]=-(\mathcal{G}\mu ,\mu )\le 0.$$

$$\frac{{\alpha }_k{\varphi }^{n+1}-A_k\left({\varphi }^n\right)}{\delta t}=-\mathcal{G}{\mu }^{n+1},$$

$${\mu }^{n+1}=\mathcal{L}{\varphi }^{n+1}+U\left(B_k\left(r^n\right),B_k\left({\varphi }^n\right)\right),$$

$$\frac{{\alpha }_k\ln ({\tilde{r}}^{n+1})-A_k\left(\ln (r^n)\right)}{\delta t}=\left(U\left(B_k\left(r^n\right),B_k\left({\varphi }^n\right)\right),\frac{{\alpha }_k{\varphi }^{n+1}-A_k\left({\varphi }^n\right)}{\delta t}\right),$$

$${\alpha }_1=1,A_1\left({\varphi }^n\right)={\varphi }^n,B_1\left({\varphi }^n\right)={\varphi }^n;$$

$${\alpha }_2=\frac{3}{2},A_2\left({\varphi }^n\right)=2{\varphi }^n-\frac{1}{2}{\varphi }^{n-1},B_2\left({\varphi }^n\right)=2{\varphi }^n-{\varphi }^{n-1};$$

$${\alpha }_3=\frac{11}{6},A_3\left({\varphi }^n\right)=3{\varphi }^n-\frac{3}{2}{\varphi }^{n-1}+\frac{1}{3}{\varphi }^{n-2},B_3\left({\varphi }^n\right)=3{\varphi }^n-3{\varphi }^{n-1}+{\varphi }^{n-2};$$

$${\alpha }_4=\frac{25}{12},A_4\left({\varphi }^n\right)=4{\varphi }^n-3{\varphi }^{n-1}+\frac{4}{3}{\varphi }^{n-2}-\frac{1}{4}{\varphi }^{n-3},$$

$$B_4\left({\varphi }^n\right)=4{\varphi }^n-6{\varphi }^{n-1}+4{\varphi }^{n-2}-{\varphi }^{n-3}.$$

$$r^{n+1}={\theta }_0^{n+1}{\tilde{r}}^{n+1}+(1-{\theta }_0^{n+1})\exp \left(E_1({\varphi }^{n+1})\right),{\theta }_0^{n+1}\in \mathcal{W},$$

$$\mathcal{W}=\left\{\theta \in [0,1]s.t.\ln (r^{n+1})-\ln ({\tilde{r}}^{n+1})\le \delta t\gamma \left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)\right\};$$

$$\mathcal{W}=\left\{\theta \in [0,1]s.t.\frac{3}{2}\ln (r^{n+1})-\frac{3}{2}\ln ({\tilde{r}}^{n+1})\le \delta t\gamma \left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)\right\};$$

$$\left({\tilde{r}}^{n+1}-\exp \left(E_1\left({\varphi }^{n+1}\right)\right)\right){\theta }_0^{n+1}\le \exp \left(\delta t\gamma \left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)+\ln ({\tilde{r}}^{n+1})\right)-\exp (E_1\left({\varphi }^{n+1})\right),$$

$$\left({\tilde{r}}^{n+1}-\exp \left(E_1\left({\varphi }^{n+1}\right)\right)\right){\theta }_0^{n+1}\le \exp \left(\frac{2}{3}\delta t\gamma \left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)+\ln ({\tilde{r}}^{n+1})\right)-\exp \left(E_1({\varphi }^{n+1})\right).$$

$$R_{R-ESAV-BDF1}^{n+1}-R_{R-ESAV-BDF1}^n\le -\delta t(1-\gamma )\left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)\le 0,$$

$$R_{R-ESAV-BDF2}^{n+1}-R_{R-ESAV-BDF2}^n\le -\delta t(1-\gamma )\left(\mathcal{G}{\mu }^{n+1},{\mu }^{n+1}\right)\le 0,$$

$$r^{n+1}\le \exp \left(E_1({\varphi }^{n+1})\right),\forall n\ge 0.$$

$${\theta }_0^{n+1}=\frac{S-\exp \left(E_1({\varphi }^{n+1})\right)}{{\tilde{r}}^{n+1}-\exp \left(E_1({\varphi }^{n+1})\right)}\in [0,1]$$

$$E\left(\Phi \right)=\sum _{i,j=1}^m\frac{1}{2}{d}_{ij}({\varphi }_i,\mathcal{L}{\varphi }_j)+E_1\left(\Phi \right),$$

Full text

{AMS}

35Q40; 65M12; 35Q55; 65M70

1 Introduction

Many significant scientific and engineering problems, such as complex fluids, new composite materials, the non-convex function optimization, etc., can be modeled by the dissipative system. From the numerical perspective, it’s critical to maintain the discrete energy dissipation law so as to obviate non-physics numerical solutions. Over the past few decades, a large effort has been devoted to construct efficient energy stable time discretized schemes for the dissipative system. Existing and popular approaches can be classified into several categories: stabilized linearly implicit approach [35, 26], exponential time differencing (ETD) approach [28, 7, 8], convex splitting approach [9, 10, 22, 1], invariant energy quadratization (IEQ) approach [30, 32, 33, 27], scalar auxiliary variable (SAV) approach [24, 23, 25] and so on.

Among these approaches, SAV method has become a very efficient and popular tool to construct energy stable schemes and has been successfully applied to gradient flow [24, 25, 6, 18, 27] and general dissipative systems [19, 16, 17]. SAV method possesses a lot of attractive superiorities, but there are still some deficiencies that need to be improved. For example, (i) it requires to solve two linear systems at each time step; (ii) it requires that nonlinear part of free energy has a lower bound; (iii) it satisfies unconditional energy stability according to modified energy rather than original energy. Recently there have been some corresponding improvements for these shortcomings. A generalized SAV approach which only requires solving one linear equation with constant coefficients has been proposed by Huang et al. [13, 12]. Cheng et al. [2, 5] proposed a novel Lagrange multiplier approach which dissipates original energy and do not require the nonlinear part of the free energy to be bounded from below. The constructed scheme requires solving a nonlinear algebraic equation which brings some additional computation costs and theoretical analysis difficulties, and may not exist a suitable solution when the time step is insufficiently small. Besides, some SAV approaches in more general form have been developed to extend its applicability in [20, 21, 3]. The constructed schemes, especially for the exponential SAV methods, don’t need to require that the nonlinear part of the free energy has a lower bound and can guarantee the positive property of SAV without any assumption. Very recently, the relaxation technique has been adopted to the SAV approach to improve the accuracy of numerical solutions in [14, 34]. The key idea in this approach is to make the modified energy link to the original energy closely by updating the auxiliary variable.

nother improvement is to extend the applicability of the SAV approach by changing the definition of the SAV. In order to see this more clearly, we briefly represent the SAV approach in more general form below.

To fix the idea, we consider a free energy in the form

[TABLE]

where $\mathcal{L}$ is a linear self-adjoint elliptic operator, $F(\phi)$ is a nonlinear energy density function. Then, the gradient flow derived from the above free energy by energy-variational principle can be written as follows

[TABLE]

where $\mathcal{G}$ is a positive definite operator which signifies the dissipative mechanism of the system, e.g. $\mathcal{G}=I$ for the $L^{2}$ gradient flow and $\mathcal{G}=-\Delta$ for the $H^{-1}$ gradient flow. And the system equipped with periodic or homogeneous Neumann boundary condition. Taking the inner product of the first equation in (2.11) with $\mu$ and the second equation in (2.11) with $\frac{\partial\phi}{\partial t}$, we obtain the energy dissipation law

[TABLE]

Introducing a SAV $r(t)=H(\int_{\Omega}F(\phi)\mathrm{d}\boldsymbol{x})$, where function $H$ is invertible. Then we have $H^{-1}(r)=\int_{\Omega}F(\phi)\mathrm{d}\boldsymbol{x}$ and can drive that

[TABLE]

The system (2.11) can be rewritten into the equivalent form as follows

[TABLE]

Similarly, we can derive the following energy dissipation law

[TABLE]

For example, a first-order semi-discrete scheme can be constructed as follows

[TABLE]

It is also pointed out in [3] that there are many choices for the definition of function $H(x)$, and commonly used are

• •

$H(x)=(x+C_{0})^{\frac{1}{\alpha}}$ with $H^{-1}(x)=x^{\alpha}-C_{0}$;

• •

$H(x)=\tanh(x)$ with $H^{-1}(x)=\operatorname{artanh}(x)$;

• •

$H(x)=\sin(x)$ with $H^{-1}(x)=\arcsin(x)$;

• •

$H(x)=\cos(x)$ with $H^{-1}(x)=\arccos(x)$;

• •

$H(x)=\ln(x+C_{0})$ with $H^{-1}(x)=\exp(x)-C_{0}$;

• •

$H(x)=\exp(x)$ with $H^{-1}(x)=\ln(x)$.

In the definition of above SAV, two kind of schemes based on exponential SAV proposed in [20, 21] has some additional advantages. No matter which type of the definition of the SAV, it can be easily shown that the above scheme is unconditionally energy stable with a modified energy $\tilde{E}(\phi^{n})=\frac{1}{2}(\mathcal{L}\phi^{n},\phi^{n})+\left(H^{-1}\left(r\right)\right)^{n+1}-\left(H^{-1}\left(r\right)\right)^{n}$. However, as described in [34], $r^{n+1}$ is not linked to $H\left(\int_{\Omega}F(\phi^{n})\mathrm{d}\boldsymbol{x}\right)$ directly and may take very different values, thus the ratio $\frac{r^{n+1}}{H\left(\int_{\Omega}F(\phi^{n})\mathrm{d}\boldsymbol{x}\right)}$ converge to a value away from $1$. Hence, it may lead to inaccurate solutions if the time step is not sufficiently small.

Our aim in this paper is to propose two kinds of relaxed exponential SAV (R-ESAV) approaches for dissipative system. The constructed schemes directly link the SAV to the free energy by introducing relaxation technique, and have outstanding advantages in the following aspects:

• •

R-ESAV schemes are unconditionally energy stable with regard to a modified energy which is closer and directly linked to the original free energy, and can improve the accuracy of the solution noticeably compared with the original ESAV approach;

• •

Only one linear system with constant coefficients needs to be solved;

• •

The constructed schemes do not need the bounded below limitation of the nonlinear part of free energy;

• •

The positive property of SAV without any assumption can be guaranteed.

In addition, our numerical results show that, R-ESAV schemes can improve the accuracy of SAV $\xi^{n+1}$, and for the plenty numerical simulations that we tested, the modified energy of our R-ESAV schemes all most always equals to the original free energy. Moreover, two kinds of BDF$k$ ($1\leq k\leq 4$) numerical schemes based on two different decoupled approaches are constructed for Navier-Stokes equations. To the author’s knowledge, it is the first time to apply relaxation technique to the general dissipative systems.

The remainder of this paper is organized as follows. In Section 2, we present the first version of relaxed exponential SAV approach for gradient flow. Then we extend the R-ESAV-1 approach to gradient systems with multiple components or multiple nonlinear potentials. In Section 3, we give the second version of relaxed exponential SAV approach for general dissipative system and construct two BDF$k$ schemes for Navier-Stokes equation. In Section 4, we present several numerical experiments by using the new approaches, and provide ample numerical simulations to validates its generality and efficiency. In Section 5, we provide some concluding remarks.

2 The first version of relaxed exponential SAV approach

In this section, we consider the improvement of ESAV approach [20] by introducing a relaxation factor to modified the SAV, which is abbreviated as R-ESAV-1 approach for convenience.

2.1 The R-ESAV-1 approach for gradient flow

Without losing generality, we consider a free energy given by

[TABLE]

where $\mathcal{L}$ is a self-adjoint linear elliptic operator, $F(\phi)$ is a nonlinear energy density function. Then, the gradient flow derived from the above free energy by energy-variational principle can be written as follows

[TABLE]

where $\mathcal{G}$ is a positive definite operator which signifies the dissipative mechanism of the system, e.g. $\mathcal{G}=I$ for the $L^{2}$ gradient flow and $\mathcal{G}=-\Delta$ for the $H^{-1}$ gradient flow.

We introduce an exponential SAV

[TABLE]

and by taking the derivative of (2.12) with respect to $t$, we have

[TABLE]

Denote $U(r,\phi)=\frac{r}{\exp\left(E_{1}(\phi)\right)}F^{\prime}(\phi),$ and notice that equation (2.13) can be rewritten equivalently into

[TABLE]

then we rewrite the system (2.11) as follows

[TABLE]

By taking the inner product of the first two equations in (LABEL:eq:gradient-flow-single-ESAV) with $\mu$ and $\frac{\partial\phi}{\partial t}$, respectively, we have

[TABLE]

Inspired by the idea of relaxation factor in [14, 34], we can construct R-ESAV-1/BDF$k$ ($1\leq k\leq 6$) schemes:

Given $r^{n-k},...,r^{n},\phi^{n-k},...,\phi^{n}$, we determine $r^{n+1},\phi^{n+1}$ via two steps as follows:

Step 1: Compute an intermediate solution $(\tilde{r}^{n+1},\phi^{n+1})$ by using the ESAV approach:

[TABLE]

where $\alpha_{k}$, $A_{k}$ and $B_{k}$ are related parameter and operators of BDF$k$ schemes, which can be obtained by Taylor expansion. For the convenience of readers, we provide the form of $k=1,2,3,4$ as follows:

First-order:

[TABLE]

Second-order:

[TABLE]

Third-order:

[TABLE]

Fourth-order:

[TABLE]

Step 2: Update $r^{n+1}$ via relaxation factor as follows:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

First-order:

[TABLE]

Second-order:

[TABLE]

with $\gamma\in[0,1]$ is a adjustable parameter.

We explain below how to choose $\theta_{0}^{n+1}$. For BDF$1$ scheme, plugging (2.24) into (2.25), we derive that if we choose $\theta_{0}^{n+1}$ such that

[TABLE]

then $\theta_{0}^{n+1}\in\mathcal{W}$. Similar to BDF$2$ scheme, we need to choose $\theta_{0}^{n+1}$ satiesfy following inequlity

[TABLE]

Denote $S=\exp\left(\delta t\gamma\left(\mathcal{G}\mu^{n+1},\mu^{n+1}\right)+\ln(\tilde{r}^{n+1})\right)$ for BDF$1$ scheme and $S=\exp\left(\frac{2}{3}\delta t\gamma\left(\mathcal{G}\mu^{n+1},\mu^{n+1}\right)\right.$ $\left.+\ln(\tilde{r}^{n+1})\right)$ for BDF$2$ scheme, the next theorem summarizes the choice of $\theta_{0}^{n+1}$.

Theorem 2.1**.**

We choose $\theta_{0}^{n+1}$ in (2.24) as follows:

If $\tilde{r}^{n+1}=\exp\left(E_{1}\left(\phi^{n+1}\right)\right)$, we set $\theta_{0}^{n+1}=0$. 2. 2.

If $\tilde{r}^{n+1}>\exp\left(E_{1}\left(\phi^{n+1}\right)\right)$, we set $\theta_{0}^{n+1}=0$. 3. 3.

If $\tilde{r}^{n+1}<\exp\left(E_{1}\left(\phi^{n+1}\right)\right)$ and $S-\exp\left(E_{1}(\phi)\right)\geq 0$, we set $\theta_{0}^{n+1}=0$. 4. 4.

If $\tilde{r}^{n+1}<\exp\left(E_{1}\left(\phi^{n+1}\right)\right)$ and $S-\exp\left(E_{1}(\phi)\right)<0$, we set $\theta_{0}^{n+1}=\frac{S-\exp\left(E_{1}(\phi^{n+1})\right)}{\tilde{r}^{n+1}-\exp\left(E_{1}(\phi^{n+1})\right)}$.

*Then, (2.27) (resp. (2.28)) for BDF$1$ (resp. BDF$2$) scheme is satisfied in all cases above and $\theta_{0}^{n+1}\in\mathcal{W}$. Moreover, we have $r^{n+1}>0$, and the scheme (2.17)-(2.24) with the above choice of $\theta_{0}^{n+1}$ satiesfies unconditionally energy stability in the sense that:

First-order:*

[TABLE]

*where $R_{R-ESAV-BDF1}^{n+1}=\frac{1}{2}\left(\mathcal{L}\phi^{n+1},\phi^{n+1}\right)+\ln(r^{n+1})$;

Second-order:*

[TABLE]

where $R_{R-ESAV-BDF2}^{n+1}=\frac{1}{4}\left(\left(\mathcal{L}\phi^{n+1},\phi^{n+1}\right)+\left(\mathcal{L}\left(2\phi^{n+1}-\phi^{n}\right),2\phi^{n+1}-\phi^{n}\right)\right)+\frac{1}{2}\left(3\ln(r^{n+1})-\ln(r^{n})\right)$. Furthermore, we have

[TABLE]

Proof.

It can be verified easily that the above choice of $\theta_{0}^{n+1}$ satiesfies (2.27) (resp. (2.28)) for BDF$1$ (resp. BDF$2$) scheme in all cases such that $\theta_{0}^{n+1}\in\mathcal{W}$.

We can obtain that $\tilde{r}^{n+1}>0$ from (2.19), and thanks to $\exp\left(E_{1}(\phi^{n+1})\right)>0$, we have $r^{n+1}>0$.

Taking the inner product of (2.17)-(2.18) with $\mu^{n+1}$ and $\frac{\phi^{n+1}-\phi^{n}}{\delta t}$ (resp. $\frac{3\phi^{n+1}-4\phi^{n}+\phi^{n-1}}{2\delta t}$), respectively, combining them with the equation (2.19) and (2.25) (resp. (2.26)), we can derive the discrete energy dissipation law (2.29) (resp. (2.30)) for BDF$1$ (resp. BDF$2$) scheme.

For Cases 1-3, we have $\theta_{0}^{n+1}=0$ so $r^{n+1}=E_{1}(\phi^{n+1})$. For Case 4, since

[TABLE]

and $\tilde{r}^{n+1}<\exp\left(E_{1}(\phi^{n+1})\right)$, we can obtain that $r^{n+1}\leq\exp\left(E_{1}(\phi^{n+1})\right)$ from (2.24). The proof is complete.

Remark 2.1**.**

*To prevent the solution “blowing up” due to the rapid growth of exponential function, we can add a positive constant $C$ to redefine the exponential SAV $r(t)=\exp\left(\frac{E_{1}(\phi)}{C}\right)$, which is similar to [20]. *

2.2 The application of R-ESAV-1 approach for gradient flows of multiple functions

We provide below the R-ESAV-1 approach for gradient flows of multiple functions by considering the following energy functional

[TABLE]

where $\Phi=\left[\phi_{1},\phi_{2},\cdots,\phi_{m}\right]^{T}$, $\mathcal{L}$ is a self-adjoint linear positive definite operator, and the constant matrix $(d_{ij}),i,j=1,\cdots,m$ is symmetric positive definite. We set $W_{i}\left(\Phi\right)=\frac{\delta E_{1}\left(\Phi\right)}{\delta\phi_{i}}$, then the associated gradient flow is given by

[TABLE]

where $\mathcal{G}$ is a nonnegative operator. Taking the inner products of (2.33) with $\mu_{i}$ and $\frac{\partial\phi_{i}}{\partial t}$ respectively, summing over $i$, and thaks to the self-adjoint of $\mathcal{L}$ and $d_{ij}=d_{ji}$, we have the energy dissipation law as follows

[TABLE]

We construct numerical schemes for gradient flow of multiple functions by using R-ESAV-1 approach. Introducing an exponential SAV

[TABLE]

then (2.33) can be rewritten as

[TABLE]

Denote

[TABLE]

then the system (2.36) can be simplified into

[TABLE]

Taking the inner product of the first two equations in (2.37) with $\mu_{i}$ and $\frac{\partial\phi_{i}}{\partial t}$, respectively, combining them with the third equation in (2.37), and summing over $i$, we obtain the equivalent energy dissipation law

[TABLE]

Next we can construct numerical scheme based on first version of R-ESAV approach and Crank-Nicolson formula (R-ESAV-1/CN) as follows.

Given $r^{n-1},r^{n},\Phi^{n-1},\Phi^{n}$, we determine $r^{n+1},\Phi^{n+1}$ via two steps as follows:

Step 1: Compute an intermediate solution $(\tilde{r}^{n+1},\Phi^{n+1})$ by using the ESAV-1 approach:

[TABLE]

where $g^{n+\frac{1}{2}}=\frac{g^{n+1}+g^{n}}{2}$ and $g^{*,n+\frac{1}{2}}$ can be any second-order explicit approximation of $g(t^{n+1/2})$, such as $g^{*,n+\frac{1}{2}}=\frac{3}{2}g^{n}-\frac{1}{2}g^{n-1}$.

Step 2: Update $r^{n+1}$ via relaxation factor as follows:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

[TABLE]

with $\gamma\in[0,1]$ is a adjustable parameter.

Inserting (2.42) into the inequality of (2.43), we observe that if we choose $\theta_{0}^{n+1}$ satiesfies following condition

[TABLE]

then $\theta_{0}^{n+1}\in\mathcal{W}$. Denote $S=\exp\left(\delta t\gamma\sum_{i=1}^{m}\left(\mathcal{G}\mu_{i}^{n+\frac{1}{2}},\mu_{i}^{n+\frac{1}{2}}\right)+\ln(\tilde{r}^{n+1})\right)$, then we can choose $\theta_{0}^{n+1}$ according to Theorem 2.1 and the scheme (2.39)-(2.42) satisfies the unconditional energy stability in the sense that

[TABLE]

where $R_{R-ESAV-CN}^{n+1}=\frac{1}{2}\sum_{i,j=1}^{m}d_{ij}\left(\phi_{i}^{n+1},\mathcal{L}\phi_{j}^{n+1}\right)+\ln(r^{n+1})$.

2.3 Extension to the multiple ESAV-1 approach

In this subsection, we present how to construct relaxed multiple ESAV-1 (R-MESAV-1) schemes for gradient flow, where the model may include disparate terms such that original schemes with only one SAV has limitation on describing the different disparate evolution processes and may require overly small time steps to obtain accurate numerical solution [4].

Without losing generality, we study gradient flow with two disparate nonlinear terms as follows and it can be easily extended to more than two disparate nonlinear terms

[TABLE]

where $\mathcal{L}$ is a self-adjoint linear elliptic operator, $F_{1}(\phi),F_{2}(\phi)$ are nonlinear potential function, $\mathcal{G}$ is a linear positive definite operator. The system (2.46) satisfies the following energy dissipation law

[TABLE]

where the total energy is

[TABLE]

We first consider relaxed MESAV approach based on first ESAV approach. Setting $E_{1}(\phi)=\int_{\Omega}F_{1}(\phi)\mathrm{d}\boldsymbol{x},\,E_{2}(\phi)=\int_{\Omega}F_{2}(\phi)\mathrm{d}\boldsymbol{x}$ and introducing two SAVs $r_{1}(t)=\exp\left(E_{1}(\phi)\right),\,r_{2}(t)=\exp\left(E_{2}(\phi)\right)$, we can rewrite the system (2.46) as

[TABLE]

Denote $U_{1}(r_{1},\phi)=\frac{r_{1}}{\exp\left(E_{1}(\phi)\right)}F_{1}^{\prime}(\phi)$ and $U_{2}(r_{2},\phi)=\frac{r_{2}}{\exp\left(E_{2}(\phi)\right)}F_{2}^{\prime}(\phi)$, then (2.49) can be transformed as

[TABLE]

Then we can construct R-MESAV-1/CN schemes inspired by the idea of relaxation factor.

Given $\phi^{n-1},\phi^{n},r_{1}^{n-1},r_{1}^{n},r_{2}^{n-1},r_{2}^{n}$, we determine $\phi^{n+1},r_{1}^{n+1},r_{2}^{n+1}$ via two steps as follows:

Step 1: Compute an intermediate solution $(\phi^{n+1},\tilde{r}_{1}^{n+1},\tilde{r}_{2}^{n+1})$ by using the MESAV approach:

[TABLE]

Step 2: Update the SAVs $r_{1}^{n+1},r_{2}^{n+1}$ via relaxation factor as follows:

[TABLE]

where $\mathcal{W}$ is a set defined as follows:

[TABLE]

with $\gamma\in[0,1]$ is a adjustable parameter.

We explain below how to choose $\theta_{0}^{n+1}$. Plugging (2.55) into (2.56), we derive that if we choose $\theta_{0}^{n+1}$ such that

[TABLE]

then $\theta_{0}^{n+1}\in\mathcal{W}$. And $\theta_{0}^{n+1}$ can be regarded as a solution of the optimization problem as follows

[TABLE]

where the coefficients are

[TABLE]

Denote

[TABLE]

and the next theorem summarizes the choice of $\theta_{0}^{n+1}$.

Theorem 2.2**.**

We choose $\theta_{0}^{n+1}$ in (2.55) as follows:

If $a_{1}=0$ and $a_{2}=0$, we set $\theta_{0}^{n+1}=0$. 2. 2.

If $a_{1}>0$ and $a_{2}>0$, we set $\theta_{0}^{n+1}=0$. 3. 3.

If $a_{1}<0$ and $a_{2}<0$, we set $\theta_{0}^{n+1}=\max\left\{0,\theta_{1}\right\}$. 4. 4.

If $a_{1}a_{2}<0$ and $b_{1}b_{2}-c\leq 0$, we set $\theta_{0}^{n+1}=0$. 5. 5.

If $a_{1}a_{2}<0$ and $b_{1}b_{2}-c>0$, we set $\theta_{0}^{n+1}=\theta_{1}$. 6. 6.

If $a_{1}a_{2}=0$ and $a_{1}b_{2}+a_{2}b_{1}>0$, we set $\theta_{0}^{n+1}=0$. 7. 7.

If $a_{1}a_{2}=0$ and $a_{1}b_{2}+a_{2}b_{1}<0$, we set $\theta_{0}^{n+1}=\max\left\{0,\frac{c-b_{1}b_{2}}{a_{1}b_{2}+a_{2}b_{1}}\right\}$.

Then, (2.57) for scheme (2.52)-(2.55) is satisfied in all cases above and $\theta_{0}^{n+1}\in\mathcal{W}$. Moreover, we have $r_{1}^{n+1}>0,r_{2}^{n+1}>0$, and the scheme (2.52)-(2.55) with the above choice of $\theta_{0}^{n+1}$ is unconditionally energy stable in the sense that:

[TABLE]

*where $R_{R-MESAV-CN}^{n+1}=\frac{1}{2}\left(\mathcal{L}\phi^{n+1},\phi^{n+1}\right)+\ln(r_{1}^{n+1})+\ln(r_{2}^{n+1})$. *

Proof.

We find optimal relaxation $\theta_{0}^{n+1}$ by discussing the coefficient of (2.58) case by case.

• •

If $a_{1}=0$ and $a_{2}=0$, we have $\theta_{0}^{n+1}=0$ obviously.

• •

If $a_{1}a_{2}\neq 0$, notice that

[TABLE]

we have $1\in\mathcal{W}$ which means $\mathcal{W}\neq\emptyset.$

• –

If $a_{1}a_{2}>0$, and thanks to $f(1)\leq 0$, we have $\theta_{0}^{n+1}=\max\left\{0,\theta_{1}\right\}$, then we can easily derive Case 2, 3.

• –

If $a_{1}a_{2}<0$, and $f(0)=b_{1}b_{2}-c\leq 0$, we have $\theta_{0}^{n+1}=0$.

• –

If $a_{1}a_{2}<0$, and $f(0)=b_{1}b_{2}-c>0,f(1)\leq 0$, we have $\theta_{0}^{n+1}=\max\left\{\theta_{1},\theta_{2}\right\}$. And since $\theta_{1}\geq\theta_{2}$, then $\theta_{0}^{n+1}=\theta_{1}$.

• •

If $a_{1}a_{2}=0$, and $a_{1}b_{2}+a_{2}b_{1}>0$, thanks to $f(1)\leq 0$, we have $\theta_{0}^{n+1}=0$.

• •

If $a_{1}a_{2}=0$, and $a_{1}b_{2}+a_{2}b_{1}<0$, since $f(1)\leq 0$, then we have $\theta_{0}^{n+1}=\max\left\{0,\frac{c-b_{1}b_{2}}{a_{1}b_{2}+a_{2}b_{1}}\right\}$.

We derive from (2.54)-(2.54) that $\tilde{r}_{1}^{n+1}>0,\tilde{r}_{2}^{n+1}>0$, and thanks to $\exp\left(E_{1}(\phi^{n+1})\right)>0,\exp\left(E_{2}(\phi^{n+1})\right)>0$, we have $r_{1}^{n+1}>0,r_{2}^{n+1}>0$.

Taking the inner product of (2.52)-(2.52) with $\mu^{n+1}$ and $\frac{\phi^{n+1}-\phi^{n}}{\delta t}$, respectively, combining them with the equation (2.54)-(2.54) and (2.56), we can obtain the desired result (2.59).

The R-ESAV-1 approach for Navier-Stokes equation We consider the following classic incompressible Navier-Stokes equation

[TABLE]

where $\Omega$ is an open bounded domain in $\mathbb{R}^{d}(d=2,3)$ with a sufficiently smooth boundary $\partial\Omega,J=(0,T]$, $\mathbf{u},p$ represent the unknown velocity and pressure respectively, $\mathbf{f}$ is an external body force, $\nu>0$ is the viscosity coefficient and $\mathbf{n}$ is the unit outward normal of the domain $\Omega$. The system (3.115) satisfies the following law

[TABLE]

where $E(\mathbf{u})=\frac{1}{2}\|\mathbf{u}\|^{2}$ is the total energy.

Introduce an exponential SAV $r(t)=\exp\left(E\left(\mathbf{u}\right)\right)$, then we can rewrite the governing system (3.115) into the equivalent form as follows

[TABLE]

Next we construct first-order and second-order pressure correction schemes for the above system (2.62). Firstly, the first-order semi-implicit version can be written as

Given $\mathbf{u}^{n},p^{n},r^{n}$, we compute $\mathbf{u}^{n+1},p^{n+1},r^{n+1}$ via the following four steps:

Step 1: Determine an intermediate solution $\tilde{\mathbf{u}}^{n+1}$:

[TABLE]

Step 2: Solve solution $\left(\tilde{\mathbf{u}}^{n+1},p^{n+1}\right)$:

[TABLE]

Step 3: Solve an intermediate solution $\tilde{r}^{n+1}$:

[TABLE]

Step 4: Update the SAV $r^{n+1}$ via the following relaxation:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

[TABLE]

with $\gamma\in[0,1]$ is a adjustable parameter.

Then the second-order semi-implicit version can be written as

Given $\mathbf{u}^{n-1},\mathbf{u}^{n},p^{n-1},p^{n},r^{n-1},r^{n}$, we compute $\mathbf{u}^{n+1},p^{n+1},r^{n+1}$ via the following four steps:

Step 1: Determine an intermediate solution $\tilde{\mathbf{u}}^{n+1}$:

[TABLE]

where $g^{*,n+1}$ can be any second-order explicit approximation of $g(t^{n+1})$, such as $g^{*,n+1}=2g^{n}-g^{n-1}$.

Step 2: Solve solution $\left(\tilde{\mathbf{u}}^{n+1},p^{n+1}\right)$:

[TABLE]

Step 3: Solve an intermediate solution $\tilde{r}^{n+1}$:

[TABLE]

Step 4: Update the SAV $r^{n+1}$ via the following relaxation:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

[TABLE]

with $\gamma\in[0,1]$ is a adjustable parameter.

3 The second version of relaxed exponential SAV approach

Inspired by the total energy based on exponential SAV approach in [21], we construct the second version of relaxed exponential SAV approach(R-ESAV-2) for dissipative system in this section.

3.1 The R-ESAV-2 approach for dissipative system

More generally, we consider dissipative system

[TABLE]

where $\mathcal{A}$ is a positive operator and $g(\phi)$ is a semi-linear or quasi-linear operator. Assume it satisfies the following energy dissipation law

[TABLE]

where $E(\phi)>-C_{0}$ is a free energy and $\mathcal{K}(\phi)>0$ for all $\phi$. Introducing a SAV $R(t)=\exp\left(E(\phi)\right)$, we transform the equation (3.79) into the equivalent system as follows

[TABLE]

where $V(\xi)$ is a function related to $\xi$ and $V(\xi)\equiv 1$ at a continuous level.

Then we can construct R-ESAV-2/BDF$k$ ($1\leq k\leq 6$) schemes inspired by the idea of relaxation factor.

Given $R^{n-k},...,R^{n},\phi^{n-k},...,\phi^{n}$, we determine $R^{n+1},\phi^{n+1}$ via two steps as follows:

Step 1: Compute an intermediate solution $(\tilde{R}^{n+1},\phi^{n+1})$ by using the ESAV approach:

[TABLE]

where $\alpha_{k}$, $A_{k}$ and $B_{k}$ are defined as above and $V_{k}(\xi^{n+1})$ for $k$th-order ($1\leq k\leq 4$) scheme are following:

First-order: $V_{1}(\xi^{n+1})=\xi^{n+1}$;

Second-order: $V_{2}(\xi^{n+1})=\xi^{n+1}(2-\xi^{n+1})$;

Third-order: $V_{3}\left(\xi^{n+1}\right)=\xi^{n+1}\left(3-3\xi^{n+1}+\left(\xi^{n+1}\right)^{2}\right)$;

Forth-order: $V_{4}\left(\xi^{n+1}\right)=\xi^{n+1}\left(2-\xi^{n+1}\right)\left(2-2\xi^{n+1}+\left(\xi^{n+1}\right)^{2}\right)$.

Step 2: Update the SAV $R^{n+1}$ via relaxation factor as follows:

[TABLE]

where

[TABLE]

and $\gamma^{n+1}\geq 0$ is to be determined such that the set $\mathcal{W}$ is not empty.

Then we describe how to choose $\theta_{0}^{n+1}$ and $\gamma^{n+1}$. Insertting (3.85) into the equality of the set (3.86), we derive that if we choose $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ such that

[TABLE]

then, $\theta_{0}^{n+1}\in\mathcal{W}$. The summation of the choice of $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ are provided in next theorem.

Theorem 3.1**.**

The choice of $\theta_{0}^{n+1}$ in (3.85) and $\gamma^{n+1}$ in (3.86) are shown as follows:

*If $\tilde{R}^{n+1}=\exp\left(E(\phi^{n+1})\right)$, we set $\theta_{0}^{n+1}=0$ and $\gamma^{n+1}=\frac{\mathcal{K}\left(B_{k}(\phi^{n})\right)}{\mathcal{K}({\phi}^{n+1})}$. * 2. 2.

If $\tilde{R}^{n+1}>\exp\left(E(\phi^{n+1})\right)$, we set $\theta_{0}^{n+1}=0$ and

[TABLE] 3. 3.

If $\tilde{R}^{n+1}<\exp\left(E(\phi^{n+1})\right)$ and $\tilde{R}^{n+1}-\exp\left(E(\phi^{n+1})\right)+\delta t\tilde{R}^{n+1}\mathcal{K}\left(B_{k}(\phi^{n})\right)\geq 0$, we set $\theta_{0}^{n+1}=0$ and $\gamma^{n+1}$ the same as (3.88). 4. 4.

If $\tilde{R}^{n+1}<\exp\left(E(\phi^{n+1})\right)$ and $\tilde{R}^{n+1}-\exp\left(E(\phi^{n+1})\right)+\delta t\tilde{R}^{n+1}\mathcal{K}\left(B_{k}(\phi^{n})\right)<0$, we set $\theta_{0}^{n+1}=1-\frac{\delta t\tilde{R}^{n+1}\mathcal{K}\left(B_{k}(\phi^{n})\right)}{\exp\left(E(\phi^{n+1})\right)-\tilde{R}^{n+1}}$ and $\gamma^{n+1}=0$.

Then, (3.87) is satisfied in all cases above and $\theta_{0}^{n+1}\in\mathcal{W}$. Moreover, we have $R^{n+1}>0,\,\xi_{k}^{n+1}>0$, and the scheme (3.82)-(3.85) with the above choice of $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ satiesfies unconditionally energy stability in the sense that

[TABLE]

and more importantly we have

[TABLE]

*Furthermore, we have *

[TABLE]

Proof.

It can be verified easily that the above choice of $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ satiesfies (3.87) in all cases such that $\theta_{0}^{n+1}\in\mathcal{W}$.

Since $R^{0}>0$. It follows from (3.83) that

[TABLE]

Then we derive from (3.84) that $\xi_{k}^{1}>0$, and we derive from (3.85) that $R^{1}>0$. Therefore, it is easy to obtain $\xi_{k}^{n+1}>0$ and $R^{n+1}>0$ by induction method.

Then we obtain (3.89) by combining (3.83) and (3.86).

For Cases 1-3, it can obtain that $\theta_{0}^{n+1}=0$, then we have $R^{n+1}=\exp\left(E(\phi^{n+1})\right)$. For Case 4, thanks to $\theta_{0}^{n+1}=1-\frac{\delta t\tilde{R}^{n+1}\mathcal{K}\left(B_{k}(\phi^{n})\right)}{\exp\left(E(\phi^{n+1})\right)-\tilde{R}^{n+1}}\in(0,1]$ and $\tilde{R}^{n+1}<\exp\left(E(\phi^{n+1})\right)$, we derive that $R^{n+1}\leq\exp\left(E(\phi^{n+1})\right)$ from (3.85) .

Remark 3.1**.**

As a further extension, firstly, we can also construct numerical schemes for gradient flows of multiple functions (2.33) by using the R-ESAV-2 approach, and $\theta_{0}^{n+1},\gamma^{n+1}$ can be chosen similarly according to Theorem 3.1. Secondly, the R-ESAV-2 approach can be also extended to multiple ESAV form based on gradient flow with two disparate nonlinear terms (2.46). Setting $E_{1}(\phi)=\frac{1}{2}(\mathcal{L}\phi,\phi)+\int_{\Omega}F_{1}(\phi)\mathrm{d}\boldsymbol{x},\,E_{2}(\phi)=\int_{\Omega}F_{2}(\phi)\mathrm{d}\boldsymbol{x}$ and introducing two SAVs $R_{1}(t)=\exp\left(E_{1}(\phi)\right),\,R_{2}(t)=\exp\left(E_{2}(\phi)\right)$, we can rewrite the equation (2.46) as

[TABLE]

*The forms of the third and fourth energy equations of (3.104) are different from that in (3.81). Then we can construct BDF$k$ numerical schemes for above system. And we can choose $\theta_{0}^{n+1}$ according to Theorem 2.2 if we set $0\leq\gamma\leq\frac{\left(\mathcal{G}\mu\left(B_{k}(\phi^{n})\right),\mu\left(B_{k}(\phi^{n})\right)\right)}{\left(\mathcal{G}\mu^{n+1},\mu^{n+1}\right)}$. *

The application of R-ESAV-2 approach for gradient flows of multiple functions We can also construct numerical schemes for gradient flows of multiple functions (2.33) by using the R-ESAV-2 approach. Introduce an exponential SAV

[TABLE]

and rewrite (2.33) as

[TABLE]

Then we construct R-ESAV-2/CN schemes as follows.

Given $\Phi^{n-1},\Phi^{n},r^{n-1},r^{n}$, we compute $\Phi^{n+1},r^{n+1}$ via the following two steps:

Step 1: Determine an intermediate solution $(\Phi^{n+1},\tilde{r}^{n+1})$ by using the ESAV-2 method:

[TABLE]

Step 2: Update the SAV $R^{n+1}$ via the following relaxation:

[TABLE]

where

[TABLE]

with $\gamma^{n+1}\geq 0$ to be determined so that $\mathcal{W}$ is not empty. Plugging (3.100) into the equality of (3.101), we find that if we choose $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ such that the following condition is satisfied:

[TABLE]

then, $\theta_{0}^{n+1}\in\mathcal{W}$. Then we can choose $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ according to Theorem 3.1, and the scheme (3.96)-(3.100) satisfies energy dissipation law as follows

[TABLE]

3.2 Extension to the multiple ESAV-2 approach

The R-ESAV-2 approach can be also extended to multiple ESAV form based on gradient flow with two disparate nonlinear terms (2.46). The method is also applicable to general dissipative systems, but the algorithm may need to be adjusted accordingly. Setting $E_{1}(\phi)=\frac{1}{2}(\mathcal{L}\phi,\phi)+\int_{\Omega}F_{1}(\phi)\mathrm{d}\boldsymbol{x},\,E_{2}(\phi)=\int_{\Omega}F_{2}(\phi)\mathrm{d}\boldsymbol{x}$ and introducing two SAVs $R_{1}(t)=\exp\left(E_{1}(\phi)\right),\,R_{2}(t)=\exp\left(E_{2}(\phi)\right)$, we can rewrite the equation (2.46) as

[TABLE]

Then we can construct R-MESAV-2/BDF$k$ ($1\leq k\leq 6$) schemes for (3.104).

Given $\phi^{n-k},...,\phi^{n},R_{1}^{n-k},...,R_{1}^{n},R_{2}^{n-k},...,R_{2}^{n}$, we compute $\phi^{n+1},R_{1}^{n+1},R_{2}^{n+1}$ via the following two steps:

Step 1: Determine an intermediate solution $(\phi^{n+1},\tilde{R}_{1}^{n+1},\tilde{R}_{2}^{n+1})$ by using the MESAV method:

[TABLE]

Step 2: Update the SAV $R_{1}^{n+1},R_{2}^{n+1}$ via the following relaxation:

[TABLE]

where

[TABLE]

with $\gamma\geq 0$ is a adjustable parameter.

We explain below how to choose $\theta_{0}^{n+1}$. Plugging (3.110) into the equality of (3.111), we find that if we choose $\theta_{0}^{n+1}$ such that the following condition is satisfied:

[TABLE]

then $\theta_{0}^{n+1}\in\mathcal{W}$. And $\theta_{0}^{n+1}$ can be chosen as a solution of the following optimization problem,

[TABLE]

where the coefficients are

[TABLE]

Then we can choose $\theta_{0}^{n+1}$ according to Theorem 2.2 if set $0\leq\gamma\leq\frac{\left(\mathcal{G}\mu\left(B_{k}(\phi^{n})\right),\mu\left(B_{k}(\phi^{n})\right)\right)}{\left(\mathcal{G}\mu^{n+1},\mu^{n+1}\right)}$. And scheme (3.105)-(3.110) satisfies unconditional energy stable in the sense that

[TABLE]

3.3 The fully decoupled R-ESAV-2 approach for Navier-Stokes equation

In this subsection, we consider the fully decoupled R-ESAV-2 approach for the following incompressible Navier-Stokes equation, which is a classic dissipative system

[TABLE]

where $\Omega$ is an open bounded domain in $\mathbb{R}^{d}(d=2,3)$ with a sufficiently smooth boundary $\partial\Omega,\mathcal{T}=(0,T]$, $\mathbf{u},p$ are the unknown velocity and pressure respectively, $\mathbf{f}$ represents an external body force, $\nu>0$ is the viscosity coefficient and $\mathbf{n}$ is the unit outward normal of the domain $\Omega$. The system (3.115) satisfies the following law

[TABLE]

where $E(\mathbf{u})=\frac{1}{2}\|\mathbf{u}\|^{2}$ is the total energy.

Then we can consider the R-ESAV-2 approach for Navier-Stokes equation (3.115).

Introduce an exponential SAV $R(t)=\exp\left(E\left(\mathbf{u}\right)\right)$, then we can rewrite the governing system (3.115) into the equivalent form as follows

[TABLE]

Next, we construct two BDF$k$ schemes for (3.117). First one is based on pressure correction approach.

**Scheme I: ** Given $\mathbf{u}^{n-k},\cdots,\mathbf{u}^{n},p^{n-k},\cdots,p^{n},r^{n-k},\cdots r^{n}$, we solve $\mathbf{u}^{n+1},p^{n+1},r^{n+1}$ via four steps as follows:

Step 1: Determine solution $\tilde{R}^{n+1}$ and compute $\xi^{n+1}$:

[TABLE]

Step 2: Compute an intermediate solution $\tilde{\mathbf{u}}^{n+1}$:

[TABLE]

Step 3: Solve solution $\left(\mathbf{u}^{n+1},p^{n+1}\right)$:

[TABLE]

where operator $\hat{B}_{1}=B_{1}$ for BDF$1$ scheme and $\hat{B}_{k}=B_{k-1}$ for BDF$k$ ($k\geq 2$) scheme.

Step 4: Update the SAV $R^{n+1}$ via relaxation factor as follows:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

[TABLE]

with $\gamma^{n+1}\geq 0$ to be determined such that $\mathcal{W}$ is not empty.

Inspired by projection method in [29, 11], we can construct the following fully decoupled R-ESAV-2 scheme:

**Scheme II: ** Given $\mathbf{u}^{n-k},\cdots,\mathbf{u}^{n},p^{n-k},\cdots,p^{n},r^{n-k},\cdots r^{n}$, we determine $\mathbf{u}^{n+1},p^{n+1},r^{n+1}$ via two steps as follows:

Step 1: Solve solution $\left(\mathbf{u}^{n+1},p^{n+1},\tilde{R}^{n+1}\right)$:

[TABLE]

where $\mathbf{n}$ is the outward normal of $\partial\Omega$.

Step 2: Update the SAV $R^{n+1}$ via relaxation factor as follows:

[TABLE]

where, $\mathcal{W}$ is a set defined as follows:

[TABLE]

and $\gamma^{n+1}\geq 0$ is to be determined such that the set $\mathcal{W}$ is not empty.

Remark 3.2**.**

In the case of periodic boundary condition, the operators $\nabla,\nabla\cdot$ and $\Delta^{-1}$ can commute with each other by defining them in the Fourier space. In the absence of $\mathbf{f}$, taking the divergence on both sides of first equation of (3.115), we obtain

[TABLE]

Then the first equation of (3.115) can be rewritten as

[TABLE]

where $\mathbf{J}$ is defined by

[TABLE]

Moreover, in addition to (3.116), the system (3.115) satisfies energy dissipation law as follows

[TABLE]

Thus, (3.3)-(3.3) can be replaced by

[TABLE]

Setting $K(\mathbf{u})=\nu\|\nabla\mathbf{u}\|^{2}$ or $K(\mathbf{u})=\nu\|\Delta\mathbf{u}\|^{2}$, we can choose $\theta_{0}^{n+1}$ and $\gamma^{n+1}$ in Scheme I and Scheme II accordding to Theorem 3.1. Similarly, above schemes satiesfy energy dissipation law (3.89) and (3.90).

4 Numerical simulations

In this section, we demonstrate ample numerical results to verify that the constructed R-ESAV-1 and R-ESAV-2 approaches are accurate and efficient. Besides, we also give detailed comparisons between the original ESAV schemes with the constructed R-ESAV schemes. We consider the numerical examples with periodic boundary condition and use the Fourier spectral method for spatial discretization in what follows unless explicitly given, and the dissipation rate parameter $\gamma$ is set to $1$ for the R-ESAV-1 schemes by default.

Example 1. We first consider the Allen-Cahn equation

[TABLE]

(i) We give an exact solution

[TABLE]

and $f$ is the external force satiesfying (4.141). We set computational domain $\Omega=[0,2]^{2}$ and the model parameter $\sigma_{0}=0.01^{2}$. For spatial discretization, we use Fourier mode $N^{2}=64^{2}$, so that compared with the time discretization error, the spatial discretization error is negligible.

The convergence rates of the $L^{2}$ error at $T=1$ obtained by various schemes are presented in Fig. 1, where we can observe that

(a) Numerical results are all consistent with the expected convergence rates;

(b) The errors of of R-ESAV-1 (resp. R-ESAV-2) schemes for BDF$1$ scheme are obviously smaller than that of ESAV-1 (resp. ESAV-2) schemes;

(c) The improvement in the accuracy for the ESAV schemes with relaxation for higher-order schemes is not as notable as for first-order scheme.

We also demonstrate the evolution of relaxation factor $\theta_{0}^{n+1}$ obtained by R-ESAV-1/BDF$2$ and R-ESAV-2/BDF$2$ scheme with time step $\delta t=1e-3$ in Fig. 2. It can be obviously show that $\theta_{0}^{n+1}$ always takes the value zeros except at an initial time interval for R-ESAV-2/BDF$2$ scheme.

(ii) We choose the initial condition as

[TABLE]

where $(\lambda,\rho)$ are the polar coordinates of $(x,y)$. We set computational domain as $\Omega=[0,1]^{2}$, the other parameters are $\sigma_{0}=0.01^{2}$, and Fourier modes are $N^{2}=128^{2}$. The computational solution of the semi-implicit/BDF2 scheme obtained by time step $\delta t=1e-5$ is regarded as the reference solution. It represents the $L^{2}$-norm error of four numerical schemes we constructed above at $T=200$ with various time steps in Table 1. It can be observed that, compared with ESAV-1 (resp. ESAV-2) schemes, R-ESAV-1 (resp. R-ESAV-2) schemes can noticeably reduce the error of the solution. The error of solution is large when the time step of ESAV-2 scheme is not sufficiently small, while R-ESAV-2 scheme can improve accuracy obviously. It also shows a comparison of energy (first), errors of energy (second) and the evolution of errors of energy at different time steps (third) for the constructed schemes in Fig. 3. Moreover the evolution of error of $\xi^{n+1}$ is presented in Fig. 4, which indicates that the R-ESAV-1 (resp. R-ESAV-2) scheme can improve the accuracy compared with ESAV-1 (resp. ESAV-2) scheme and the error of $\xi^{n+1}$ for R-ESAV-2 scheme will reach the machine accuracy after the simulation reaching the steady state.

figure[htbp]

(Example 1B.) Dynamics driven by the Allen-Cahn equation using R-GSAV-BDF2 scheme.

Example 2. The Cahn-Hilliard equation

[TABLE]

(i) We also choose (4.142) as the exact solution, and set model parameter to be $\sigma_{0}=0.04,M=0.005,\epsilon=1$. Fig. 6 shows the convergence rates of different schemes. We can observe similar results as those of the Allen-Cahn equation.

(ii) We consider a rectangular array of $9\times 9$ circles as the initial condition

[TABLE]

where $r_{0}=0.085,x_{m}=0.2\times m,y_{n}=0.2\times n$ for $m,n=1,2,\cdots,9$. We set computational domain as $[0,2]^{2}$, the other parameters are $M=1e-6,\sigma_{0}=1,\epsilon=0.01$ and Fourier modes are $N^{2}=512^{2}$ in the simulations. The evolutions of an array of circles governed by Cahn-Hilliard equation using the R-ESAV-2/BDF$2$ scheme with $\delta t=1e-3$ are shown in Fig. 7(b).

Example 3. In this example, we consider the phase-field crystal model

[TABLE]

which is a $H^{-1}$-gradient flow associated with the total free energy as follows

[TABLE]

where $M>0$ is the mobility coefficient. We set $M=1,\zeta=1$ in the following simulations.

(i) We first simulate the crystal growth in a super-cooled liquid in $2D$. We adopt the following initial condition

[TABLE]

where $x_{j}$ and $y_{j}$ define a local system of Cartesian coordinates, which is oriented with the crystallite lattice. We choose the constant parameters as $\bar{\phi}=0.285,\alpha_{1}=0.446,\alpha_{2}=0.66$. Then, we define five crystallites in five small square patches, each with a length of $40$, located at $(200,200)$, $(150,600)$, $(350,400)$, $(600,300)$ and $(550,700)$ respectively. In order to produce crystallites with different orientations, we use the following affine transformation to generate rotation by five diffferent angles $\rho=-\frac{3\pi}{4},-\frac{\pi}{4},0,\frac{\pi}{4},\frac{3\pi}{4}$ respectively

[TABLE]

We choose Fourier modes $N^{2}=1024^{2}$ to discretize the space and $\delta t=0.02$ to discretize the time. And we set the other parameters as $\epsilon=0.25,T=800$. The relaxation parameter $\theta_{0}^{n+1}$ is also always zero in this example. We present the crystal growth in a super-cooled liquid by using the R-ESAV-1/BDF$2$ scheme for the PFC equation in Fig. 8(b), which indicate that the different alinement of the crystallites leads to defects and dislocations, just similar with those in [31, 18].

(ii) The initial condition is set to be

[TABLE]

where $\delta=rand(x,y)$ is the uniformly distributed random number in $[-1,1]$ with zeros mean. Set the computational domain to $\Omega=[0,128]^{2}$. The model parameter is chosen as $\epsilon=0.025$, time step is $\delta t=0.1$ and $N^{2}=256^{2}$ Fourier modes to discretize the space. It presents the configuration evolution in Fig. 9 and can observe from that uniform phase separation is formed finally. Similar computation results can be found in [20, 21].

(iii) Then we simulate the crystal growth in a super-cooled liquid in $3D$. We choose the computational domain $[0,100]^{3}$ and the initial condition are two crystallites generated by

[TABLE]

The other parameters are chosen as $\epsilon=0.25,\delta t=0.02,T=2000$. We adopt Fourier modes $N^{3}=128^{3}$ to discretize space. It can be observed that the effects of different arrangement of crystallites on the growth of the crystalline phase and the motion of crystal-liquid interfaces in Fig. 10(b), where these results are also consistent with those in [15].

(iv) Finally we study phase transition behaviors in $3D$. The initial data are chosen as

[TABLE]

and set computational domains $\left[0,50\right]^{3}$. Other parameters are chosen as $\epsilon=0.56,\delta t=0.02,T=3000$ and Fourier modes $N^{3}=64^{3}$. We present the steady state microstructure of the phase transition behavior for $\bar{\phi}=0.20,0.35$ and $0.43$, respectively in Fig. 11.

Example 4. In this example, we simulate the phase-field vesicle membrane (PFVM) model [4, 5] to demonstrate the effectiveness of the relaxed MESAV schemes. We consider the following penalized free energy to preserve the area and volume of vesicle membrane,

[TABLE]

where $\sigma_{1},\sigma_{2}$ are two small parameters, and $v_{0},s_{0}$ are the initial volume and initial surface area, respectively. The definition of bending energy $E_{b}(\phi)$, volume $V(\phi)$ and surface area $S(\phi)$ of the vesicle are as follows

[TABLE]

[TABLE]

where

[TABLE]

Then, the dynamic equation based on the above total energy can be described by

[TABLE]

with the periodic boundary condition, and $M$ is the mobility constant. Then, it can easily obtain that the system (4.154) satisfy the energy law as follows

[TABLE]

It can be observed that the system (4.154) contains two nonlinear terms associated with two small parameters $\epsilon$ and $\sigma_{2}$ respectively. Therefore, two SAVs are needed to introduce to deal with the different nonlinear terms. In the following simulations, we set computational domain as $\Omega=(-\pi,\pi)^{3}$, and the model parameters are $\sigma_{1}=\sigma_{2}=0.01,\epsilon=\frac{6\pi}{128},M=1$. We compute the results using the R-MESAV-1/BDF$2$ scheme with time step $\delta t=1e-4$ and $N^{3}=128^{3}$ Fourier modes.

(i) We first consider the interaction of four closeby spheres as the initial condition given by

[TABLE]

where $R_{j}=\frac{\pi}{6},x_{j}=0,\left(y_{1},y_{2},y_{3},y_{4}\right)=\left(-\frac{\pi}{4},\frac{\pi}{4},-\frac{3\pi}{4},\frac{3\pi}{4}\right),$ and $z_{j}=0$ for $j=1,2,3,4$.

Snapshots of iso-surfaces of $\phi=0$ at $t=0,0.02,1$ are presented in Fig. 12, which indicates that four small spheres gradually linked the shape of ‘ice sugar gourd’, and finally merge into a cylinder shape. The results are consistent with those presented in [5]. We also plot the evolution of relaxation factor $\theta_{0}^{n+1}$ in Fig. 13, and observe that, except the first several steps, $\theta_{0}^{n+1}$ always takes the value zeros.

(ii) We simulate the evolution of five close-by spherical vesicles by choosing the initial condition

[TABLE]

where $R_{j}=\frac{\pi}{6},z_{j}=0$ for $j=1,2,\ldots,5$, $\left(x_{1},x_{2},x_{3},x_{4},x_{5}\right)=\left(-\frac{\pi}{3},\frac{\pi}{3},0,-\frac{\pi}{3},\frac{\pi}{3}\right)$, and

$\left(y_{1},y_{2},y_{3},y_{4},y_{5}\right)=\left(-\frac{\pi}{3},-\frac{\pi}{3},0,\frac{\pi}{3},\frac{\pi}{3}\right)$.

We represent the evolution process in Fig. 14(b). It can be observed that five spheres connect within a small time interval, gradually form a doughnut shape which is a final state.

(iii) Then we consider a more complicated initial condition which is nine close-by spherical vesicles given by

[TABLE]

where $R_{j}=\frac{\pi}{6},z_{j}=0$ for $j=1,2,\ldots,9$, $\left(x_{1},x_{2},x_{3},x_{4},x_{5},x_{6},x_{7},x_{8},x_{9}\right)=\left(-\frac{\pi}{2},0,\frac{\pi}{2},-\frac{\pi}{2},0,\frac{\pi}{2},-\frac{\pi}{2},0,\frac{\pi}{2}\right)$, and $\left(y_{1},y_{2},y_{3},y_{4},y_{5},y_{6},y_{7},y_{8},y_{9}\right)=\left(-\frac{\pi}{2},-\frac{\pi}{2},-\frac{\pi}{2},0,0,0,\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2}\right)$.

The evolutions of nine close-by spherical vesicles are demonstrated in Fig. 15(b), which represents that the initially nine spheres gradually connect with each other and finally form a big vesicle.

Example 5. In this numerical simulation, we test the Navier-Stokes equation.

(i) We start with the accuracy test. The right hand side is computed according to the following analytic solution

[TABLE]

the computational domain is $\Omega=(-1,1)^{2}$, other parameters are chosen as $\nu=0.1,T=1$. We use Legendre spectral method to discretize space and $N^{2}=64^{2}$. The numerical results for BDF$1$ and BDF$2$ of Schemes I and II are presented in Tables 2-5 respectively. It can be seen that the errors of the solution will be greatly reduced by using the relaxation factor.

(ii) Next we simulate double shear layer problem. We consider the Navier-Stokes equation with periodic boundary condition, and the initial condition provided by

[TABLE]

where $\sigma$ is the parameter of the shear layer width and $\epsilon$ is the size of the perturbation. We choose $\sigma=30,\epsilon=0.05$ and computational domain $\Omega=(0,2)^{2}$ in the next simulations. We use Fourier modes $N^{2}=128^{2}$ and set $\delta t=6.7e-4$ to test the Navier-Stokes equation with $\nu=5e-5$. The vorticity contours of velocity field $\mathbf{u}$ at $T=1.2$ using BDF$k$ ($k=1,2,3,4$) schemes are presented in Fig. 16. It can be observed that the BDF$3$ and BDF$4$ schemes give correct solutions, while the BDF$1$ scheme leads to a totally wrong result and the result of the BDF$2$ scheme is also inaccurate. The numerical phenomenon demonstrates the superiorities of high-order schemes. The evolution of contours of vorticity with $\nu=1e-4$, and $\delta t=6e-4$ obtained by BDF$2$ scheme is shown in Fig. 17(b). We can observe from the results that the vortex increases gradually.

5 Conclusions

In this paper, we constructed two kind of R-ESAV/BDF$k$ schemes which can improve the accuracy significantly by introducing a relaxation factor to improve the accuracy for the classical SAV method. The constructed schemes are linear and unconditionally energy stable. They can guarantee the positive property of SAV without any assumption compared with R-SAV and R-GSAV approach. Moreover the constructed R-ESAV-2 approach is easy to construct high-order BDF$k$ schemes and can be applied to general dissipative system. Finally we proved that the constructed RESAV scheme with relaxation are more accuracy and closer to original energy by using ample numerical examples.