Arbitrarily High-order Maximum Bound Preserving Schemes with Cut-off Postprocessing for Allen-Cahn Equations
AARBITRARILY HIGH-ORDER MAXIMUM BOUND PRESERVING SCHEMES WITHCUT-OFF POSTPROCESSING FOR ALLEN-CAHN EQUATIONS ∗ JIANG YANG † , ZHAOMING YUAN ‡ , AND
ZHI ZHOU § Abstract.
We develop and analyze a class of maximum bound preserving schemes for approximately solving Allen–Cahn equations. We apply a k th-order single-step scheme in time (where the nonlinear term is linearized by multi-stepextrapolation), and a lumped mass finite element method in space with piecewise r th-order polynomials and Gauss–Lobattoquadrature. At each time level, a cut-off post-processing is proposed to eliminate extra values violating the maximum boundprinciple at the finite element nodal points. As a result, the numerical solution satisfies the maximum bound principle (atall nodal points), and the optimal error bound O ( τ k + h r + ) is theoretically proved for a certain class of schemes. Thesetime stepping schemes under consideration includes algebraically stable collocation-type methods, which could be arbitrar-ily high-order in both space and time. Moreover, combining the cut-off strategy with the scalar auxiliary value (SAV)technique, we develop a class of energy-stable and maximum bound preserving schemes, which is arbitrarily high-order intime. Numerical results are provided to illustrate the accuracy of the proposed method. Keywords:
Allen-Cahn equation, single step methods, lumped mass FEM, cut off, high-order, maximum bound pre-serving, energy-stable.
AMS subject classifications 2010:
1. Introduction.
The aim of this paper is to design and analyze a high-order maximum boundpreserving (MBP) scheme for solving the Allen–Cahn equation:(1.1) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ u t = ∆ u + f ( u ) in Ω × ( , T ) ,u ( x, t = ) = u ( x ) in Ω × { } ,∂ n u = ∂ Ω × ( , T ) where Ω is a smooth domain in R d with the boundary ∂ Ω. f ( u ) = − F ′ ( u ) with a double-well potential F that has two wells at ± α , for some known parameter α >
0. For two popular choices of potentials. Itis well-known that the Allen–Cahn equation (1.1) has the maximum bound principle [7]: ∣ u ( x )∣ ≤ α (cid:212)⇒ ∣ u ( x, t )∣ ≤ α for all ( x, t ) ∈ Ω × ( , T ] . (1.2)As a typical L gradient flow associating with the following free energy E ( u ) = ∫ Ω ∣∇ u ∣ + F ( u ) d x, the nonlinear energy dissipation law holds in the sense(1.3) dd t E ( u ) = − ∫ Ω ∣ u t ∣ d x ≤ . The Allen–Cahn equation was originally introduced by Allen and Cahn in [2] to describe the motionof anti-phase boundaries in crystalline solids. In the context, u represents the concentration of one of the ∗ The work of J. Yang is supported by National Natural Science Foundation of China (NSFC) Grant No. 11871264,Natural Science Foundation of Guangdong Province (2018A0303130123), and NSFC/Hong Kong RRC Joint ResearchScheme (NFSC/RGC 11961160718), and the research of Z. Yuan and Z. Zhou is partially supported by Hong Kong RGCgrant (No. 15304420). † Department of Mathematics & SUSTech International Center for Mathematics, Southern University of Science andTechnology, Shenzhen 518055, China. ( yangj7sustech.edu.cn ) ‡ Department of Applied Mathematics, The Hong Kong Polytechnic University, Kowloon, Hong Kong.( [email protected] ) § Department of Applied Mathematics, The Hong Kong Polytechnic University, Kowloon, Hong Kong.( [email protected], [email protected] ) 1 a r X i v : . [ m a t h . NA ] F e b wo metallic components of the alloy and the parameter ε involved in the nonlinear term represents theinterfacial width, which is small compared to the characteristic length of the laboratory scale. Recentdecades, the Allen–Cahn equation has become one of basic phase-field equations, which has been widelyapplied to many complicated moving interface problems in materials science and fluid dynamics througha phase-field approach coupled with other models [3, 5, 31]. The development and analysis of MBP method have been in-tensively studied in existing references. It was proved in [23, 27] that the stabilized semi-implicit Eulertime-stepping scheme, with central difference method in space, preserves the maximum principle uncon-ditionally if the stabilizer satisfies certain restrictions. In [6], a stabilized exponential time differencingscheme was proposed for solving the (nonlocal) Allen–Cahn equation, and the scheme was proved to beunconditionally MBP. See also [7] for the generalization to a class of semilinear parabolic equations. Thesecond-order backward differentiation formula (with nonuniform meshes) was applied to develop an MBPscheme in [17] under the usual CFL condition τ = O ( h ) .High-order strong stability preserving (SSP) time-stepping methods are widely used in the devel-opment of MBP scheme for both parabolic equations and hyperbolic equations (see e.g., [10–12, 18, 19,21, 30, 32]). Recently, an SSP integrating factor Runge–Kutta method of up to order four was proposedand analyzed in [14] for semilinear hyperbolic and parabolic equations. For semilinear hyperbolic andparabolic equations with strong stability (possibly in the maximum norm), the method can preserve thisproperty and can avoid the standard parabolic CFL condition τ = O ( h ) , only requiring the stepsize τ to be smaller than some constant depending on the nonlinear source term, also referring to [15]. Anonlinear constraint limiter was introduced in [29] for implicit time-stepping schemes without requiringCFL conditions, which can preserve maximum principle at the discrete level with arbitrarily high-ordermethods by solving a nonlinearly implicit system.Very recently, a new class of high-order MBP methods was proposed in [16]. The method consists ofa k th-order multistep exponential integrator in time, and a lumped mass finite element method in spacewith piecewise r th-order polynomials. At every time level, the extra values exceeding the maximumbound are eliminated at the finite element nodal points by a cut-off operation. Then the numericalsolution at all nodal points satisfies the MBP, and an error bound of O ( τ k + h r ) was proved. However,numerical results in [16, Table 4.1] indicates that the error bound is not sharp in space, and how toimprove the estimate it is still open. Besides, the aforementioned scheme requires to evaluate someactions of exponential functions of diffusion operators, which might be relatively expensive comparedwith solving poisson problems, and the generalization to other time stepping schemes is a nontrivial task.Finally, the proposed scheme (with relatively coarse step sizes) might produce a numerical solution withobviously increasing and oscillating energy. These motivate our current project. The first contribution of thecurrent paper is to develop and analyze a class of MBP schemes, which could be arbitrarily high-orderin both space and time, for approximately solving the Allen–Cahn equation (1.1). In time, we apply asingle-step method, which is (strictly) accurate of order k , and apply multistep extrapolation to linearizethe nonlinear term. In space, we apply the lumped mass FEM with piecewise r th-order polynomials andGauss–Lobatto quadrature, as in [16]. At each time level, we apply a cut-off operation to remove theextra value exceeding the maximum bound at the nodal points. We estabilish the error estimate of order O ( τ k + h r + ) , which fills the gap between the numerical results in [16, Theorem 3.2] showing optimalconvergence rate O ( h r + ) and the theoretical result in [16, Table 4.1] providing only a suboptimal errorestimate of order O ( h r ) . The improvement follows from a careful examination of quadrature errors (seeRemark 2.4 and [16, eq. (2.6) and (3.22)]). To the best of our knowledge, this is the first work derivingoptimal error estimates of arbitrarily high-order MBP schemes for the Allen–Cahn equation (1.1).Nevertheless, the optimal estimate of the fully discrete scheme (with the cut-off postprocessing)requires the L-stability of the time stepping scheme, which excludes some popular and practical singestep method, e.g. Gauss–Legendre method belonging to algebraically stable collocation Runge–Kutta ethod. Therefore, we revisit this class of time stepping methods and prove the same error estimate byusing the energy argument without using the L-stability. This is the second contribution of the paper.In case of relative coarse step sizes, the proposed time stepping scheme (with the cut-off operationat each time level) might result in oscillating and increasing energy (see e.g. Figure 2 (middle)), whichviolates the energy dissipation law (1.3) of Allen–Cahn equation (1.1). This motivates us to developa class of energy-stable and MBP schemes, by combining the cut-off strategy with the scalar auxiliaryvalue (SAV) method [26]. The scheme is second order in space but could be arbitrarily high-order intime. As far as we know, this is the first scheme that is unconditionally energy-dissipative, maximumbound preserving, and arbitrarily high-order in time scheme with a provable error bound. In fact, ournumerical results show that the use of SAV regularizes the numerical solution, stabilizes the energy, andsignificantly reduces the cut-off values at each time level (see e.g. Figure 2).The rest of the paper is organized as follows. In section 2, we consider the single step methods (ina general framework) with cut-off postprocessing and multistep extrapolation. Error estimate for bothsemidiscrete and fully discrete scheme are provided, where the optimal error estimate of the fully discretescheme requires the L-stability of the time stepping scheme. In section 3, we analyze the algebraicallystable collocation scheme and show the same error estimate without using the L-stability. In section 4,combining the cut-off strategy with the scalar auxiliary value (SAV) method, we develop a class of energy-stable and maximum bound preserving schemes, which is arbitrarily high-order in time. In section 5, wepresent numerical results to illustrate the accuracy and effectiveness of the method in solving the Allen–Cahn equation. Throughout, the notation C , with or without subscripts, denotes a generic constant,which may differ at different occurrences, but it is always independent of the mesh size h and the timestep size τ .
2. Cut-off single-step methods with multi-step extrapolation.
In this section, we shall de-velop and analyze a class of MBP scheme for the Allen–Cahn equation (1.1). Optimal error estimatewill also be provided, which fill the gap in the preceding work [16]. Besides, the argument presented inthis section also builds the foundation of developing MBP scheme which also satisfies energy dissipationproperty (in section 4).
To begin with, we consider the time discretization for theAllen–Cahn equation (1.1). To this end, we split the interval ( , T ) into N subintervals with the uniformmesh size τ = T / N , and set t n = nτ , n = , , . . . , N . On the time interval [ t n − , t n ] , we approximate thenonlinear term f ( u ( s )) by the extrapolation polynomial k ∑ j = L j ( s ) f ( u n − j ) , with known u n − k , . . . , u n − . where L j ( s ) is the Lagrange basis polynomials of degree k − L j ( t n − i ) = δ ij , i, j = , . . . , k. Thus, on [ t n − , t n ] , the linearization of (1.1) states as˜ u t = ∆˜ u + k ∑ j = L j ( s ) f ( u n − j ) . Following Duhamel’s principle yields˜ u ( t n ) = e τ ∆ u ( t n − ) + ∫ τ e ( τ − s ) ∆ k ∑ j = L j ( t n − + s ) f ( u n − j ) d s. hen a framework of a single step scheme of approximating ˜ u ( t n ) reads:(2.1) ˜ u n = σ (− τ ∆ ) u n − + τ m ∑ i = p i (− τ ∆ )( k ∑ j = L j ( t ni ) f ( u n − j )) , for all n ≥ k, with t ni = t n − + c i τ . Here, σ ( λ ) and { p i ( λ )} mi = are rational functions and c i are distinct real numbers in [ , ] . For simplicity, we assume that the scheme (2.1) satisfies the following assumptions. (P1) ∣ σ ( λ )∣ < ∣ p i ( λ )∣ ≤ c , for all i = , . . . , m , uniformly in τ and λ >
0. Besides, the numerator of p i ( λ ) is of lower degree than its denominator. (P2) The time stepping scheme (2.1) is accurate of order k in sense that σ ( λ ) = e − λ + O ( λ k + ) , as λ → . and, for 0 ≤ j ≤ k m ∑ i = c ji p i ( λ ) − j ! (− λ ) j + ( e − λ − j ∑ (cid:96) = (− λ ) (cid:96) (cid:96) ! ) = O ( λ k − j ) , as λ → . (P3) The time discretization scheme (2.1) is strictly accurate of order q in sense that m ∑ i = c ji p i ( λ ) − j ! (− λ ) j + ( σ ( λ ) − j ∑ (cid:96) = (− λ ) (cid:96) (cid:96) ! ) = , for all 0 ≤ j ≤ q − . Remark 2.1.
In practice, it is convenient to choose p i ( λ ) s that share the same denominator of σ ( λ ) ,for instance: σ ( λ ) = a ( λ ) g ( λ ) , and p i ( λ ) = a i ( λ ) g ( λ ) , for i = , , . . . , m, where a i ( λ ) and g ( λ ) are polynomials. Then the time stepping scheme (2.1) could be written as g (− τ ∆ ) ˜ u n = a (− τ ∆ ) u n − + τ m ∑ i = a i (− τ ∆ )( k ∑ j = L j ( t ni ) f ( u n − j )) , for all n ≥ k. See e.g. [28, pp. 131] for the construction of such rational functions satisfying the Assumptions (P1)-(P3).
Unfortunately, the time stepping scheme (2.1) does not satisfy the maximum bound principle. There-fore, at each time step, we apply the cut-off operation: for n ≥ k , we find u n such thatˆ u n = σ (− τ ∆ ) u n − + τ m ∑ i = p i (− τ ∆ )( k ∑ j = L j ( t ni ) f ( u n − j )) , (2.2) u n = min ( max ( ˆ u n , − α ) , α ) , (2.3)where α is the maximum bound given in (1.2). The accuracy of this cut-off semi-discrete method isguaranteed by the next theorem. Theorem 2.1.
Suppose that the Assumptions (P1) and (P2) are fulfilled, and (P3) holds for q = k .Let u ( t ) be the solution to the Allen–Cahn equation, and u n be the solution to the time stepping scheme (2.2) - (2.3) . Assume that ∣ u ∣ ≤ α and the maximum principle (1.2) holds, and assume that the startingvalues u j , j = , . . . , k − , are given and ∣ u j ∣ ≤ α, for all j = , . . . , k − . hen the semi-discrete solution given by (2.2) - (2.3) satisfies for all n ≥ k ∣ u n ∣ ≤ α, and ∥ u n − u ( t n )∥ ⩽ Cτ k + C k − ∑ j = ∥ u j − u ( t j )∥ , provided that f is locally Lipschitz continuous, ∆ u ∈ C k ([ , T ] ; L ( Ω )) , u ∈ C k + ([ , T ] ; L ( Ω )) and f ( u ) ∈ C k ([ , T ] ; L ( Ω )) .Proof . Due to the cut-off operation (2.3), the discrete maximum bound principle follows immediately.Then it suffices to show the error estimate.Let e n = u n − u ( t n ) and ˆ e n = ˆ u n − u ( t n ) . Since the exact solution satisfies the maximum bound (1.2),we have ∥ e n ∥ L ( Ω ) ≤ ∥ ˆ e n ∥ L ( Ω ) . Then it is easy to note that ˆ e n = σ (− τ ∆ ) e n − + ϕ n , n ≥ k. where ϕ n can be written as ϕ n = − u ( t n ) + σ (− τ ∆ ) u ( t n − ) + τ m ∑ i = p i (− τ ∆ )( k ∑ j = L j ( t ni ) f ( u n − j ))= τ m ∑ i = p i (− τ ∆ )( k ∑ j = L j ( t ni ) f ( u n − j ) − f ( t ni )))+ ( − u ( t n ) + σ (− τ ∆ ) u ( t n − ) + τ m ∑ i = p i (− τ ∆ )( ∂ t u − ∆ u )( t ni ))=∶ I + II.
Then the bound of I follows from the approximation property of Lagrange interpolation, the maximumbound of u n − j and u ( t n − j ) , j = , . . . , k , the locally Lipschitz continuity of f , and the Assumption (P1): ∥ I ∥ L ( Ω ) ≤ τ m ∑ i = ∥ p i (− τ ∆ )∥ L ( Ω )→ L ( Ω ) ∥ k ∑ j = L j ( t ni ) f ( u ( t n − j )) − f ( u ( t n − + c i τ ))∥ L ( Ω ) + τ m ∑ i = ∥ p i (− τ ∆ )∥ L ( Ω )→ L ( Ω ) k ∑ j = ∣ L j ( t ni )∣ ∥ f ( u n − j ) − f ( u ( t n − j ))∥ L ( Ω ) ≤ Cτ k + ∥ f ( u )∥ C k ([ t n − k ,t n ] ; L ( Ω )) + Cτ k ∑ j = ∥ e n − j ∥ L ( Ω ) . Now we term to the second term II , which can be rewritten by Taylor’s expansion at t n − II = − k ∑ j = τ j j ! u ( j ) ( t n − ) + σ (− τ ∆ ) u ( t n − )+ τ m ∑ i = p i (− τ ∆ ) k − ∑ j = ( c i τ ) j j ! ( u ( j + ) − ∆ u ( j ) )( t n − ) + R + R . here the remainders R and R are R = ∫ t n t n − ( t n − s ) k k ! u ( k + ) ( s ) d s and R = τ m ∑ i = p i (− τ ∆ ) ∫ t n − + c i τt n − ( t n − + c i τ − s ) k − ( k − ) ! ( u ( k + ) − ∆ u ( k ) )( s ) d s respectively. Hereafter, we use u ( j ) to denote the j th derivative in time. Then Assumption (P1) implies ∥ R + R ∥ L ( Ω ) ≤ Cτ k + (∥ u ∥ C k + ([ t n − ,t n ] ; L ( Ω )) + ∥ ∆ u ∥ C k ([ t n − ,t n ] ; L ( Ω )) ) . Now we revisit the three leading terms of II . Note that − k ∑ j = τ j j ! u ( j ) ( t n − ) + σ (− τ ∆ ) u ( t n − ) + τ m ∑ i = p i (− τ ∆ ) k − ∑ j = ( c i τ ) j j ! ( u ( j + ) − ∆ u ( j ) )( t n − )=( − I + σ (− τ ∆ ) − τ m ∑ i = p i (− τ ∆ ) ∆ ) u ( t n − )+ k − ∑ j = τ j j ! ( − I + j m ∑ i = c j − i p i (− τ ∆ ) − τ m ∑ i = c ji p i (− τ ∆ ) ∆ ) u ( j ) ( t n − )+ τ k k ! ( − I + k m ∑ i = c k − i p i (− τ ∆ )) u ( k ) ( t n − ) = ∑ (cid:96) = II (cid:96) . Since the time stepping scheme is strictly accurate of order q = k (by Assumption (P3)), we have II = II =
0. Meanwhile, we apply Assumption (P3) again to arrive at for λ > − + k m ∑ i = c k − i p i ( λ ) = λ k ! (− λ ) k + ( σ ( λ ) − k ∑ (cid:96) = (− λ ) (cid:96) (cid:96) ! ) =∶ λγ ( λ ) . Note that ∣ γ ( λ )∣ = O ( ) for λ → ∣ γ ( λ )∣ → λ → +∞ . Hence ∣ γ ( λ )∣ isbounded uniformly in [ , ∞) . Then we arrive at ∥ II ∥ L ( Ω ) ≤ Cτ k + ∥ ∆ u ( k ) ( t n − )∥ ≤ Cτ k + ∥ ∆ u ∥ C k ([ t n − ,t n ] ; L ( Ω )) . In conclusion, we obtain the following estimate ∥ e n ∥ L ( Ω ) ≤ ∥ σ (− τ ∆ ) e n − ∥ L ( Ω ) + Cτ k + + Cτ k ∑ j = ∥ e n − j ∥ L ( Ω ) . Then the assumption (P1) leads to ∥ e n ∥ L ( Ω ) ≤ ∥ e n − h ∥ L ( Ω ) + Cτ k + + Cτ k ∑ j = ∥ e n − j ∥ L ( Ω ) . Finally, the desired assertion follows immediately by using discrete Gronwall’s inequality ∥ e n ∥ L ( Ω ) ≤ Ce cT τ k + Ce cT k − ∑ j = ∥ e j ∥ L ( Ω ) . emark 2.2. Theorem 2.1 implies that the cut-off operation preserves the maximum bound withoutlosing global accuracy. However, the Assumption (P3) is restrictive. It is well-known that a single stepmethod with a given m ∈ Z + could be accurate of order m (Gauss–Legendre method) [8, Section 2.2],but at most strictly accurate of order m + [4, Lemma 5]. In general, a collocation-type method is onlystrictly accurate of order m + .Without the assumption of strict accuracy, one may still show the error estimate, provided that f ( u ) satisfies certain compatibility conditions, e.g., f ( u ) ∈ C (cid:96) ([ , T ] ; Dom ( ∆ k − (cid:96) )) for all (cid:96) = , , . . . , k, that requires ∂ n ∆ q f ( u ) = for (cid:96) = , , . . . , k − . Unfortunately, those compatibility conditions cannot befulfilled in general for semilinear parabolic problems. Remark 2.3.
The same error estimate could be proved by assuming that the scheme satisfies theassumption (P3) with q = k − and some additional conditions (see e.g. [28, Theorem 8.4] and [20]).However, the proof is not directly applicable when we apply the cut-off operation at each time step. Itwarrants further investigation to show the sharp convergence rate O ( τ k ) with weaker assumptions. In this part, we discuss the fully discrete scheme. To illustrate themain idea, we consider the one-dimensional case Ω = [ a, b ] , and the argument could be straightforwardlyextended to multi-dimensional cases, see Remark 2.5. We denote by a = x < x < ⋅ ⋅ ⋅ < x Mr = b a partitionof the domain with a uniform mesh size h = x ir − x ( i − ) r = ( b − a )/ M , and denote by S rh the finite elementspace of degree r ≥
1, i.e., S rh = { v ∈ H ( Ω ) ∶ v ∣ I i ∈ P r , i = , . . . , M } , where I i = [ x ( i − ) r , x ir ] and P r denotes the space of polynomials of degree ≤ r .Let x ( i − ) r + j and ω j , j = , . . . , r , be the quadrature points and weights of the ( r + ) -point Gauss–Lobatto quadrature on the subinterval I i , and denote w ( i − ) r + j = { ω j for 1 ≤ j ≤ r − , ω j for j = , r. Then we consider the piecewise Gauss–Lobatto quadrature approximation of the inner product, i.e., ( f, g ) h ∶= Mr ∑ j = w j f ( x j ) g ( x j ) . This discrete inner product induces a norm ∥ f h ∥ h = √( f h , f h ) h ∀ f h ∈ S rh . Then we have the following lemma for norm equivalence. The proof follows directly from the positivityof Gauss–Lobatto quadrature weights [22, p. 426].
Lemma 2.2.
The discrete norm ∥ ⋅ ∥ h is equivalent to usual L norm ∥ ⋅ ∥ L ( Ω ) in sense that C ∥ v h ∥ L ( Ω ) ≤ ∥ v h ∥ h ≤ C ∥ v h ∥ L ( Ω ) , ∀ v h ∈ S rh . where C and C are independent of h . To develop the fully discrete scheme, we introduce the discrete Laplacian − ∆ h ∶ S rh → S rh such that (− ∆ h v h , w h ) h = (∇ v h , ∇ w h ) for all v h , w h ∈ S rh . (2.4) hen at n -th time level, with given u n − kh , . . . , u n − h ∈ S rh , we find an intermediate solution ˆ u nh ∈ S rh suchthat(2.5) ˆ u nh = σ (− τ ∆ h ) u n − h + τ m ∑ i = p i (− τ ∆ h )( k ∑ j = L j ( t ni ) Π h f ( u n − jh )) where t ni = t n − + c i τ , and Π h ∶ C ( Ω ) → S rh is the Lagrange interpolation operator. In order to imposethe maximum bound, we apply the cut-off postprocessing: find u nh ∈ S rh such that u nh ( x j ) = min ( max ( ˆ u nh ( x j ) , − α ) , α ) , j = , . . . , M r. (2.6)It is equivalent to u nh = Π h min ( max ( ˆ u nh , − α ) , α ) . Essentially, the cut-off operation (2.6) only works on the finite element nodal points.Next, we shall prove the optimal error estimate of the fully discrete scheme (2.5)-(2.6). To this end,we need the following stability estimate of operators σ (− τ ∆ h ) and p i (− τ ∆ h ) . Lemma 2.3.
Let ∆ h be the discrete Laplacian defined in (2.4) , and σ ( λ ) and p i ( λ ) are rationalfunctions satisfying the Assumption (P1). Then there holds that for all v h ∈ S rh (2.7) ∥∇ q σ (− τ ∆ h ) v h ∥ h ≤ ∥∇ q v h ∥ h and ∥∇ q p i (− τ ∆ h ) v h ∥ h ≤ C ∥∇ q v h ∥ h with i = , . . . , m and q = , . Meanwhile, (2.8) τ ∥∇ q ∆ h p i (− τ ∆ h ) v h ∥ h ≤ C ∥∇ q v h ∥ h i = , . . . , m, q = , Proof . Let {( λ j , ϕ hj )} Mr + j = be eigenpairs of − ∆ h , where { ϕ hj } Mr + j = forms an orthogonal basis of S rh insense that ( ϕ hi , ϕ hj ) h = δ i,j . Then by the Assumption (P1), we have for any v h ∈ S rh and q = , ∥∇ q σ (− τ ∆ h ) v h ∥ h = Mr + ∑ j = ( λ hj ) q ∣ σ ( τ λ j )∣ ∣( v h , ϕ hj ) h ∣ ≤ Mr + ∑ j = ( λ hj ) q ∣( v h , ϕ hj ) h ∣ = ∥∇ q v h ∥ h . This shows the first estimate. The estimate for p i follows analogously.Moreover, the numerator of p i ( λ ) is of lower degree than its denominator (by Assumption (P1)), andhence there exists constants C , C > ∣ p i ( λ )∣ ≤ C + C λ , for all λ > . Then we derive that for any v h ∈ S rh and q = , τ ∥∇ q ∆ h p i (− τ ∆ h ) v h ∥ h = τ Mr + ∑ j = ( λ hj ) q + ∣ p i ( τ λ j )∣ ∣( v h , ϕ hj ) h ∣ ≤ Cτ Mr + ∑ j = ( λ hj ) q + ( + Cτ λ hj ) ∣( v h , ϕ hj ) h ∣ ≤ C Mr + ∑ j = ( λ hj ) q ∣( v h , ϕ hj ) h ∣ = C ∥∇ q v h ∥ h , here the constant C only depends on C and C . This proves the assertion (2.8). Lemma 2.4.
Let v ∈ H r + ( Ω ) with the homogeneous Neumann boundary condition and ϕ h ∈ S rh .Then we have the following estimate ( Π h ∆ v − ∆ h Π h v, ϕ h ) h ≤ Ch r + ∥ v ∥ H r + ∥ ϕ h ∥ H ( Ω ) . Proof . Using the homogeneous Neumann boundary condition and (2.4), we obtain(2.9) ( Π h ∆ v − ∆ h Π h v, ϕ h ) h = ( Π h ∆ v, ϕ h ) h − ( ∆ h Π h v, ϕ h ) h = (( ∆ v, ϕ h ) h − ( ∆ v, ϕ h )) + (( ∆ v, ϕ h ) − ( ∆ h Π h v, ϕ h ) h )= (( ∆ v, ϕ h ) h − ( ∆ v, ϕ h )) + (( ∂ x v, ∂ x ϕ h ) − ( ∂ x Π h v, ∂ x ϕ h )) Since the ( r + ) -point Gauss–Lobatto quadrature on each subinterval I i is exact for polynomials of degree2 r − ∣( ∆ v, ϕ h ) h − ( ∆ v, ϕ h )∣ = ∣ M ∑ i = ( r ∑ j = ω j ( ∆ vϕ h )( x ( i − ) r + j ) − ∫ I i ( ∆ v ) ϕ h d x )∣≤ Ch r M ∑ i = ∥ ∆ vϕ h ∥ W r, ( I i ) ≤ Ch r M ∑ i = ∥ v ∥ H r + ( I i ) ∥ ϕ h ∥ H r ( I i ) ≤ Ch r + M ∑ i = ∥ v ∥ H r + ( I i ) ∥ ϕ h ∥ H ( I i ) ≤ Ch r + ∥ v ∥ H r + ( Ω ) ∥ ϕ h ∥ H ( Ω ) . Similar argument also leads to the estimate for the second term in (2.9) for r ≥ ∣( ∂ x ( v − Π h v ) , ∂ x ϕ h )∣ = ∣ M ∑ i = ∫ I i ∂ x ( v − Π h v ) ∂ x ϕ h d x ∣ = ∣ M ∑ i = ∫ I i ( v − Π h v ) ∂ x ϕ h d x ∣= ∣ M ∑ i = ∫ I i v∂ x ϕ h d x − r ∑ j = ω j ( v∂ x ϕ h )( x ( i − ) r + j )∣≤ Ch r M ∑ i = ∥ v∂ x ϕ h ∥ W r, ( I i ) ≤ Ch r M ∑ i = ∥ v ∥ H r + ( I i ) ∥ ϕ h ∥ H r ( I i ) ≤ Ch r + M ∑ i = ∥ v ∥ H r + ( I i ) ∥ ϕ h ∥ H ( I i ) ≤ Ch r + ∥ v ∥ H r + ( Ω ) ∥ ϕ h ∥ H ( Ω ) . Finally, in case that r =
1, it is easy to observe that ( ∂ x ( v − Π h v ) , ∂ x ϕ h ) = M ∑ i = ∫ I i ∂ x ( v − Π h v ) ∂ x ϕ h d x = − M ∑ i = ∫ I i ( v − Π h v ) ∂ x ϕ h d x = . To derive an error estimate for the fully discrete scheme (2.5)-(2.6). We need the following extraassumptions on the rational function σ ( λ ) . (P4) The rational function σ ( λ ) satisfies ∣ σ ( λ )∣ → λ → ∞ . ote that the Assumption (P4) immediately implies [28, eq. (7.37)] ∣ σ ( λ )∣ ≤ + c λ for any λ ≥ , with a generic constant c >
0. This further implies1 − ∣ σ ( λ )∣ − ≤ − c λ for any λ ≥ . Therefore, we have for any v h ∈ S rh ∥ σ (− τ ∆ h ) v h ∥ h = Mr + ∑ j = ∣ σ ( τ λ j )∣ ( v h , ϕ hj ) h = ∥ v h ∥ h + Mr + ∑ j = (∣ σ ( τ λ j )∣ − )( v h , ϕ hj ) h = ∥ v h ∥ h + Mr + ∑ j = ( − ∣ σ ( τ λ j )∣ − )∣ σ ( τ λ j )∣ ( v h , ϕ hj ) h ≤ ∥ v h ∥ h − c τ Mr + ∑ j = λ j ∣ σ ( τ λ j )∣ ( v h , ϕ hj ) h = ∥ v h ∥ h − c τ ∥∇ σ (− τ ∆ h ) v h ∥ . Then we are ready to state following main theorem.
Theorem 2.5.
Suppose that the Assumptions (P1), (P2) and (P4) are fulfilled, and (P3) holds for q = k . Assume that ∣ u ∣ ≤ α and the maximum principle (1.2) holds, and assume that the starting values u lh , l = , . . . , k − , are given and ∣ u lh ( x j )∣ ≤ α, j = , . . . , M r, l = , . . . , k − . Then the fully discrete solution given by (2.5) - (2.6) satisfies ∣ u nh ( x j )∣ ≤ α, j = , . . . , M r, n = k, . . . , N, and for n = k, . . . , N ∥ u ( t n ) − u nh ∥ L ( Ω ) ≤ C ( τ k + h r + ) + C k − ∑ l = ∥ u ( t l ) − u lh ∥ L ( Ω ) , provided that u ∈ C k + ([ , T ] ; C ( ¯Ω )) ∩ C k ([ , T ] ; Dom ( ∆ )) ∩ C ([ , T ] ; H r + ( Ω )) , f is locally Lipschitzcontinuous and f ( u ) ∈ C k ([ , T ] ; L ( Ω )) ∩ C ([ , T ] ; H r + ( Ω )) .Proof . In [ t n − , t n ] , we note that Π h u satisfies ∂ t Π h u ( t ) − ∆ h Π h u ( t ) = Π h f ( u ( t )) + g h ( t ) , t ∈ ( t n − , t n ] , with Π h u ( t n − ) given , and g h ( t ) = ( Π h ∆ − ∆ h Π h ) u ( t ) . Then we define its time stepping approximation w nh satisfying w nh = σ (− τ ∆ h ) Π h u ( t n − ) + τ m ∑ i = p i (− τ ∆ h )( Π h f ( u ) + g h )( t n + c i τ ) . Then the argument in Theorem 2.1 implies that ∥ Π h u ( t n ) − w nh ∥ h ≤ Cτ k + ( sup t n − ≤ t ≤ t n ∥ Π h u ( k + ) ( t )∥ h + sup t n − ≤ t ≤ t n ∥ ∆ h Π h u ( k ) ( t )∥ h ) . he first term of the right hand side is bounded by ∥ u ∥ C k + ([ ,T ] ; C ( ¯Ω )) , while the second one is boundedas ∥ ∆ h Π h u ( k ) ( t )∥ h = sup ϕ h ∈ S rh ( ∆ h Π h u ( k ) ( t ) , ϕ h ) h ∥ ϕ h ∥ h = sup ϕ h ∈ S rh (∇( Π h u ( k ) ( t ) − u ( k ) ( t )) , ∇ ϕ h ) + (∇ u ( k ) ( t ) , ∇ ϕ h )∥ ϕ h ∥ h ≤ Ch − ∥∇( Π h u ( k ) ( t ) − u ( k ) ( t ))∥ L ( Ω ) + ∥ ∆ u ( k ) ( t )∥ L ( Ω ) ≤ C ∥ u ( k ) ∥ H ( Ω ) . Therefore, we conclude that ∥ Π h u ( t n ) − w nh ∥ h ≤ Cτ k + (∥ u ∥ C k + ([ t n − ,t n ] ; C ( ¯Ω )) + ∥ u ∥ C k ([ t n − ,t n ] ; H ( Ω )) ) . Then the simple triangle inequality leads to ∥ ˆ u nh − Π h u ( t n )∥ h ≤ (∥ ˆ u nh − w nh ∥ h + ∥ w nh − Π h u ( t n )∥ h ) ≤ ( + Cτ )∥ ˆ u nh − w nh ∥ h + Cτ k + . (2.10)Let ρ nh = ˆ u nh − w nh and e nh = u nh − Π h u ( t n ) , then ρ nh satisfies ρ nh = σ (− τ ∆ h ) e n − h + I n + I n (2.11)where I n = τ m ∑ i = p i (− τ ∆ h )( k ∑ j = L j ( t n − + c i τ ) Π h f ( u n − jh ) − Π h f ( u ( t n − + c i τ ))) , and I n = − τ m ∑ i = p i (− τ ∆ h ) g h ( t n − + c i τ ) . Now take the discrete inner product between (2.11) and ρ nh ∥ ρ nh ∥ h = ( σ (− τ ∆ h ) e n − h , ρ nh ) h + ( I n , ρ nh ) h + ( I n , ρ nh ) h . Then first term, we apply the Assumption (P4) to obtain that ( σ (− τ ∆ h ) e n − h , ρ nh ) ≤ ∥ σ (− τ ∆ h ) e n − h ∥ h + ∥ ρ nh ∥ h ≤ ∥ e n − h ∥ h − c τ ∥∇ σ (− τ ∆ h ) e n − h ∥ + ∥ ρ nh ∥ h ≤ ∥ e n − h ∥ h − c τ ∥∇( ρ nh − I n − I n )∥ + ∣∣ + ∥ ρ nh ∥ h ≤ ∥ e n − h ∥ h − c τ ∥∇ ρ nh ∥ − c τ ∥∇( I n + I n )∥ + c τ (∇ ρ nh , ∇( I n + I n )) + ∥ ρ nh ∥ Then applying the definition of ∆ h , we arrive at(2.12) 12 ∥ ρ nh ∥ h ≤ ∥ e n − h ∥ h − c τ ∥∇ ρ nh ∥ − c τ ( ρ nh , ∆ h ( I n + I n )) h + ( I n , ρ nh ) h + ( I n , ρ nh ) h . y using the approximation property of interpolation I kτ , Lemma 2.3, and the fact that u n − kh , . . . , u n − h satisfies the maximum bound, we bound the fourth term in (2.12) as ∣( I n , ρ nh ) h ∣ ≤ τ m ∑ i = ∣( k ∑ j = L j ( t n − + c i τ ) Π h f ( u ( t n − j )) − Π h f ( u ( t n − + c i τ ) , p i (− τ ∆ h ) ρ nh ) h ∣+ τ m ∑ i = ∣( k ∑ j = L j ( t n − + c i τ )( Π h f ( u ( t n − j )) − Π h f ( u n − jh )) , p i (− τ ∆ ) ρ nh ) h ∣≤ Cτ m ∑ i = ∥ p i (− τ ∆ h ) ρ nh ∥ h k ∑ j = ∥ Π h f ( u ( t n − j )) − Π h f ( u n − j )∥ h + Cτ k + m ∑ i = ∥ p i (− τ ∆ h ) ρ nh ∥ h ∥ Π h f ( u )∥ C k ([ t n − k ,t n ] ; L ( Ω )) ≤ Cτ k + ∥ Π h f ( u )∥ C k ([ t n − k ,t n ] ; L ( Ω )) + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ∥ ρ nh ∥ h ≤ Cτ k + ∥ f ( u )∥ C k ([ t n − k ,t n ] ; C ( ¯Ω )) + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ∥ ρ nh ∥ h . The fifth term in (2.12) can be bounded by using lemmas 2.3 and 2.4, i.e., ∣( I n , ρ nh ) h ∣ ≤ Cτ m ∑ i = ∣( g h ( t n − + c i τ ) , p i (− τ ∆ h ) ρ nh ) h ∣≤ Cτ m ∑ i = h r + ∥ u ( t n − + c i τ )∥ H r + ( Ω ) ∥ p i (− τ ∆ h ) ρ nh ∥ H ( Ω ) ≤ Cτ h r + η ∥ u ∥ C ([ t n − ,t n ] ; H r + ( Ω )) + Cτ η ∥ ρ nh ∥ H ( Ω ) . (2.13)For the third term in the right hand side of (2.12), we shall apply the preceding argument again, togetherwith the stability estimate (2.8), and obtain that τ ∣( ρ nh , ∆ h ( I n + I n )) h ∣ ≤ Cτ m ∑ i = ∥ ∆ h p i (− τ ∆ h ) ρ nh ∥ h k ∑ j = ∥ Π h f ( u ( t n − j )) − Π h f ( u n − j )∥ h + Cτ k + m ∑ i = ∥ ∆ h p i (− τ ∆ h ) ρ nh ∥ h ∥ Π h f ( u )∥ C k ([ t n − k ,t n ] ; L ( Ω )) + Cτ m ∑ i = h r + ∥ u ( t n − + c i τ )∥ H r + ( Ω ) ∥ ∆ h p i (− τ ∆ h ) ρ nh ∥ H ( Ω ) ≤ Cτ k + ∥ f ( u )∥ C k ([ t n − k ,t n ] ; C ( ¯Ω )) + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ∥ ρ nh ∥ h + Cτ h r + η ∥ u ∥ C ([ t n − ,t n ] ; H r + ( Ω )) + Cτ η ∥ ρ nh ∥ H ( Ω ) . (2.14)Then by choosing η small, we arrive at ( − Cτ )∥ ρ nh ∥ h ≤ ∥ e n − h ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ( τ k + h r + ) . his together with (2.10) and the property of the cut-off operation lead to ∥ e nh ∥ h ≤ ∥ ˆ u nh − Π h u ( t n )∥ h ≤ ( + Cτ )∥ ρ nh ∥ h + cτ k + ≤ ∥ e n − h ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ( τ k + h r + ) , and hence we rearrange terms and obtain ∥ e nh ∥ h − ∥ e n − h ∥ h τ ≤ C ( τ k + h r + ) + C k ∑ j = ∥ e n − jh ∥ h . Then the discrete Gronwall’s inequality implies ∥ e nh ∥ h ≤ Ce cT ( τ k + h r + ) + Ce cT k − ∑ j = ∥ e jh ∥ h , and the desired error estimate follows from the equivalence of different norms by Lemma 2.2. Remark 2.4.
In [16], an error estimate O ( τ k + h r ) , which is suboptimal in space, was derived for themultistep exponential integrator method by using energy argument. The loss of the optimal convergencerate is due to the suboptimal estimate of the term ( ∂ x ( Π h u − u ) , ∂ x v h ) in [16, eq. (2.6) and (3.22)]. Theoptimal rate could be also proved by using Lemma 2.4.The Assumption (P4) , called L-stability , is useful when solving stiff problems. It is also essentialin the proof of Theorem 2.5 to derive the optimal error estimate of the extrapolated cut-off single stepscheme. In particular, Assumption (P4) immediately leads to the estimate ∥ σ (− τ ∆ h ) v h ∥ h ≤ ∥ v h ∥ h − c τ ∥∇ σ (− τ ∆ h ) v h ∥ , where the second term in the right side is used to handle the term involving ∥ ρ nh ∥ H ( Ω ) in (2.13) and (2.14) . Many single step methods, e.g., Lobatto IIIC and Radau IIA methods are L-stable [8, 13]. Forboth classes, arbitrarily high-order methods can be constructed. Nevertheless, it is not clear how to removethe restriction (P4) in general. Remark 2.5.
It is straightforward to extend the argument to higher dimensional problems, e.g., Ω isa multi-dimensional rectangular domain ( a, b ) d ⊂ R d , with d ≥ . Then we can divide Ω in to some smallsub-rectangles, called partition K , and apply the tensor-product Lagrange finite elements on the partition K . As a result, Lemma 2.4 is still valid, which implies the desired error estimate. See more details aboutthe setting for multi-dimensional problems in [16, Section 2.2].
3. Collocation-type methods with the cut-off postprocessing.
Note that the Assumption(P4) excludes some popular methods, e.g., Gauss–Legendre methods. This motivates us to discuss thecollocation-type schemes, which belong to implicit Runge–Kutta methods, and derive error estimatewithout Assumption (P4). This class of time stepping methods is easy to implement, and plays anessential role in the next section to develop an energy-stable scheme. For simplicity, we only present theargument for one-dimensional case, and it can be extended to multi-dimensional cases straightforwardlyas mentioned in Remark 2.5.Now we consider an m -stage Runge–Kutta method, described by the Butcher tableau 1. Here { c i } mi = denotes m distinct quadrature points. Definition 3.1.
We call a Runge–Kutta method is algebraically stable if the method satisfies (P5)(a)
The matrix A = ( a ij ) , with i, j = , . . . , m is invertible; (P5)(b) The coefficients b i satisfy b i > for i = , , . . . , m ; (P5)(c) The symmetric matrix
M ∈ R m × m with entries m ij ∶= b i a ij + b j a ji − b i b j , i, j = , . . . , m is positivesemidefinite. . . . a m c ⋮ ⋮ ⋮ a m . . . a mm c m b . . . b m Table 1
Butcher tableau for Runge–Kutta scheme.
Here we assume that the Runge–Kutta scheme described by Table 1 associates with a collocationmethod, i.e., coefficients a ij , b i , c i satisfy m ∑ i = b i c l − i = l , l = , ⋯ , p, (3.1) m ∑ j = a ij c l − j = c li l , l = , ⋯ , m, (3.2)with some integers p ≥ m . Two popular families of algebraically stable Runge–Kutta methods of collo-cation type satisfying (2.6) of orders p = m and p = m − n , with given u n − kh , . . . , u n − h ∈ S rh , we find an intermediate solution ˆ u nh ∈ S rh suchthat(3.3) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ u nih = ∆ h u nih + ∑ k(cid:96) = L (cid:96) ( t n − + c i τ ) Π h f ( u n − (cid:96)h ) for i = , , . . . , m,u nih = u n − h + τ ∑ mj = a ij ˙ u njh for i = , , . . . , m, ˆ u nh = u n − h + τ ∑ mi = b i ˙ u nih , where k = min ( p, m + ) , and Π h ∶ C ( Ω ) → S rh is the Lagrange interpolation operator. Then we apply thecut-off operation: find u nh ∈ S rh such that(3.4) u nh ( x j ) = min ( max ( ˆ u nh ( x j ) , − α ) , α ) , j = , . . . , M r. Remark 3.1.
Note that the scheme (3.3) is equivalent to (2.5) with ( p ( λ ) , . . . , p m ( λ )) = ( b , . . . , b m )( I + λA ) − , σ ( λ ) = − λ m ∑ j = b j p j ( λ ) . Then the Assumption (P5), and (3.1) - (3.2) imply Assumptions (P1), (P2) with order k = min ( p, m + ) and (P3) with order q = min ( p, m + ) . Hence Theorem 2.5 indicates the temporal error O ( τ min ( p,m + ) ) .This is the reason why we choose k -step extrapolation, where k = min ( p, m + ) , in the time steppingscheme (3.3) . Next, we shall derive an error estimate for the fully discrete scheme (3.3)-(3.4). To begin with, weshall examine the local truncation error. We define the local truncation error η ni and η n + as(3.5) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ u ni ∗ = ∆ u ( t ni ) + ∑ k(cid:96) = L (cid:96) ( t ni ) f ( u ( t n − (cid:96) )) for i = , , . . . , m,u ( t ni ) = u ( t n − ) + τ ∑ mj = a ij ˙ u nj ∗ + η ni for i = , , . . . , m,u ( t n ) = u ( t n − ) + τ ∑ mi = b i ˙ u ni ∗ + η n where t ni = t n − + c i τ and k = min ( p, q + ) . Then the next lemma give an estimate for the local truncationerror η ni and η n . We sketch the proof in Appendix for completeness. Lemma 3.2.
Suppose that the Assumption (P5), and relations (3.1) and (3.2) are valid. Then the ocal truncation error η ni and η n , given by (3.5) , satisfy the estimate ∥ η n ∥ H ( Ω ) + τ m ∑ i = ∥ η ni ∥ H ( Ω ) ≤ Cτ k + . with k = min ( p, q + ) , provided that u ∈ C k + ([ , T ] ; H ( Ω )) and f ( u ) ∈ C k ([ , T ] ; H ( Ω )) . Then we are ready to present the following theorem, which gives the error estimate for the cut-offRunge–Kutta scheme (3.3)-(3.4).
Theorem 3.3.
Suppose that the Runge–Kutta method given by Table 1 satisfies Assumption (P5),and relations (3.1) and (3.2) are valid. Assume that ∣ u ∣ ≤ α and the maximum principle (1.2) holds, andassume that the starting values u nh , l = , . . . , k − , are given and ∣ u lh ( x j )∣ ≤ α, j = , . . . , M r, l = , . . . , k − . Then the fully discrete solution given by (3.3) - (3.4) satisfies ∣ u nh ( x j )∣ ≤ α, j = , . . . , M r, n = k, . . . , N, and for n = k, . . . , N ∥ u ( t n ) − u nh ∥ L ( Ω ) ≤ C ( τ k + h r + ) + C k − ∑ l = ∥ u ( t l ) − u lh ∥ L ( Ω ) , provided that u ∈ C k + ([ , T ] ; H ( Ω ))∩ C ([ , T ] ; H r + ( Ω )) , f is locally Lipschitz continuous and f ( u ) ∈ C k ([ , T ] ; H ( Ω )) ∩ C ([ , T ] ; H r + ( Ω )) .Proof . Due to the cut-off operation (2.3), the discrete maximum bound principle follows immediately.With the notation e nih = Π h u ( t ni ) − u nih , ˙ e nih = Π h ˙ u ni ∗ − ˙ u nih , e nh = Π h u ( t n ) − u nh , ˆ e nh = Π h u ( t n ) − ˆ u nh , we derive the error equations(3.6) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ e nih = ∆ h e nih + ( Π h ∆ − ∆ h Π h ) u ( t ni ) + ∑ k(cid:96) = L (cid:96) ( t ni ) Π h ( f ( u ( t n − (cid:96) )) − f ( u n − (cid:96)h )) for i = , , . . . , m,e nih = e n − h + τ ∑ mj = a ij ˙ e njh + Π h η ni for i = , , . . . , m, ˆ e nh = e n − h + τ ∑ mi = b i ˙ e nih + Π h η n . Take the square of discrete L norm of both sides of the last relation of (3.6), we obtain ∥ ˆ e nh ∥ h = ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h + ( η n , e n − h + τ m ∑ i = b i ˙ e nih ) h + ∥ Π h η n ∥ h . (3.7)For the first term on the right hand side, we expand it and apply the second equation of (3.6) to obtain ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h = ∥ e n − h ∥ h + τ m ∑ i = b i ( ˙ e nih , e nih − η ni ) h − τ m ∑ i,j = m ij ( ˙ e nih , ˙ e njh ) h ⩽ ∥ e n − h ∥ h + τ m ∑ i = b i ( ˙ e nih , e nih − η ni ) h , here in the last inequality we use the positive semi-definiteness of the matrix M in the Assumption(P5). Next, we note that the first relation of (3.6) implies ( ˙ e nih , e nih − η ni ) h = ( ∆ h e nih + k ∑ (cid:96) = L (cid:96) ( t ni )( f ( u ( t n − (cid:96) )) − f ( u n − (cid:96)h )) + ( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h = −∥∇ e nih ∥ L ( Ω ) + (∇ e nih , ∇ Π h η ni ) + ( k ∑ (cid:96) = L (cid:96) ( t ni )( f ( u ( t n − (cid:96) )) − f ( u n − (cid:96)h )) , e nih − η ni ) h + (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h The bound of second term of the right hand side can be derived via Cauchy-Schwarz inequality ∣(∇ e nih , ∇ Π h η ni )∣ ≤ ∥∇ e nih ∥ L ( Ω ) + C ∥ η ni ∥ H ( Ω ) . Meanwhile, using the fact that f is locally Lipschitz and the fully disctete solutions satisfy maximumbound principle at the Gauss–Lobatto points, the third term can be bounded as ( k ∑ (cid:96) = L (cid:96) ( t ni )( f ( u ( t n − (cid:96) )) − f ( u n − (cid:96)h )) , e nih − η ni ) h ≤ C (∥ e nih ∥ h + ∥ η ni ∥ H ( Ω ) + k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h ) The bound of the last term follows from Lemma 2.4 (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h ≤ Ch r + ∥ e nih − Π h η ni ∥ H ( Ω ) ≤ ∥∇ e nih ∥ L ( Ω ) + C (∥ e nih ∥ h + ∥ η ni ∥ H ( Ω ) + h r + ) . Therefore, we arrive at2 ( ˙ e nih , e nih − η ni ) h ≤ −∥∇ e nih ∥ L ( Ω ) + C ( k ∑ j = ∥ e n − jh ∥ h + ∥ e nih ∥ h + ∥ η ni ∥ H ( Ω ) + h r + ) , and hence by Lemma 3.2, we derive ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h ⩽ ∥ e n − h ∥ h − τ m ∑ i = b i ∥∇ e nih ∥ L ( Ω ) + Cτ m ∑ i = ∥ e nih ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ( h r + + τ k ) . In view of the first relation of the error equation (3.6), we have the estimate ( η n , e n − h + τ m ∑ i = b i ˙ e nih ) h ≤ ∥ η n ∥ H ( Ω ) (∥ e n − h ∥ h + Cτ m ∑ i = b i (∥∇ e nih ∥ h + k ∑ j = ∥ e n − jh ∥ h + h r + ))≤ Cτ ( h r + + τ k ) + τ m ∑ i = b i ∥∇ e nih ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h which gives a bound of the second term in (3.7). In conclusion, we obtain that ∥ ˆ e nh ∥ h + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + Cτ m ∑ i = ∥ e nih ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h . (3.8) ext, we shall derive a bound for ∑ mi = ∥ e nih ∥ h on the right-hand side. To this end, we test the secondrelation of (3.6) by e nih . This yields m ∑ i = ∥ e nih ∥ h ≤ C ∥ e n − h ∥ h + Cτ m ∑ i,j = a ij ( ˙ e njh , e nih ) h + C m ∑ i = ∥ Π h η ni ∥ h ≤ C ∥ e n − h ∥ h + Cτ m ∑ i,j = a ij ( ˙ e njh , e nih ) h + Cτ k . Then, we apply the first relation of (3.6) and Lemma 2.4 to derive m ∑ i,j = a ij ( ˙ e njh , e nih ) h = − m ∑ i,j = a ij (∇ e njh , ∇ e nih ) + m ∑ i,j = a ij ( k ∑ (cid:96) = L (cid:96) ( t ni )( f ( u ( t n − (cid:96) )) − f ( u n − (cid:96)h )) , e nih ) h + m ∑ i,j = a ij (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih ) h ≤ C m ∑ i = (∥∇ e nih ∥ L ( Ω ) + ∥ e nih ∥ h ) + Ch r + + C k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h . Therefore, we obtain m ∑ i = ∥ e nih ∥ h ≤ C ( τ h r + + τ k ) + C ∥ e n − h ∥ h + Cτ k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h + Cτ m ∑ i = (∥∇ e nih ∥ L ( Ω ) + ∥ e nih ∥ h ) . Then for sufficiently small τ , Cτ ∑ mi = ∥ e nih ∥ h on the right-hand side can be absorbed by the left-hand side.Then, we obtain m ∑ i = ∥ e nih ∥ h ≤ C ( τ h r + + τ k ) + C ∥ e n − h ∥ h + Cτ k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h + Cτ m ∑ i = ∥∇ e nih ∥ L ( Ω ) . Now substituting the above estimate into (3.8), there holds for sufficiently small τ ∥ ˆ e nh ∥ h ≤ Cτ ( h r + + τ k ) + ∥ e n − h ∥ h + Cτ k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h . Noting that ∥ e nh ∥ h ≤ ∥ ˆ e nh ∥ h and rearranging terms, we obtain ∥ e nh ∥ h − ∥ e n − h ∥ h τ ≤ C ( h r + + τ k ) + C k ∑ (cid:96) = ∥ e n − (cid:96)h ∥ h . Then the discrete Gronwall’s inequality impliesmax k ≤ n ≤ N ∥ e nh ∥ h ≤ C ( h r + + τ k ) + C k − ∑ j = ∥ e jh ∥ h . This completes the proof of the theorem.
Remark 3.2.
In Theroem 3.3, we discuss the algebraically stable collocation-type method with cut-offtechnique. We still prove the optiaml error estimate O ( τ k + h r + ) , without the L-stability, i.e. Assumption(P4). Note that this class of methods includes Gauss–Legendre and Radau IIA methods [13, Theorem12.9], while the first one is not L-stable [13, Table 5.13]. . Fully discrete scheme based on SAV method. In the preceding section, we develop andanalyze a class of maximum bound preserving schemes. Unfortunately, the proposed scheme (with rel-atively large time steps) might produce solutions with increasing and oscillating energy, see Figure 2.This violates another essential property of the Allen–Cahn model, say energy dissipation. The aim forthis section is to develop a high-order time stepping schemes via combining the cut-off strategy and thescalar auxiliary variable (SAV) method.SAV method is a common-used method for gradient flow models. It was firstly developed in [25, 26]and have motived a sequence of interesting work on the development and analysis of high-order energy-decayed time stepping scheme in recent years [1, 9, 24].In particular, assuming that E ( u ( t )) = ∫ Ω F ( u ( x, t )) d x is globally bounded from below, i.e., E ( u ( t )) >− C . we introduce the following scalar auxiliary variable [25](4.1) z ( t ) = √ E ( u ( t )) + C and W ( u ) = f ( u )√ E ( u ) + C Then the Allen–Cahn equation in (1.1) can be reformulated as(4.2) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ u t = ∆ u + z ( t ) W ( u ) in Ω × ( , T ) ,u ( x, t = ) = u ( x ) in Ω × { } ,∂ n u = ∂ Ω × ( , T ) and the scalar auxiliary variable r ( t ) satisfies(4.3) ⎧⎪⎪⎪⎨⎪⎪⎪⎩ z ′ ( t ) = − ( W ( u ( t )) , u t ( t )) , in ( , T ) ,z ( ) = √ E ( u ) + C . One can easily show that the coupled problem (4.2)-(4.3) is equivalent to the original equation (1.1).Meanwhile, simple calculation leads to the SAV energy dissipation:(4.4) dd t ( ∥∇ u ∥ + ∣ z ( t )∣ ) = −∥ u t ( t )∥ ≤ . Inspired by [1], we discretize the coupled problem (4.2)-(4.3) by using the m -stage Runge–Kuttamethod in time (described by Table 1) and lumped mass finite element method with r = m -stage Runge–Kutta method (described by Table 1) satisfies the Assump-tion (P5) and relations (3.1) and (3.2). Then at n -th time level, with known u n − kh , . . . , u n − h ∈ S rh and z n − ∈ R , we find ˆ u nh ∈ S rh and z n ∈ R such that(4.5) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ u nih = ∆ h u nih + z ni W nih for i = , , . . . , m,u nih = u n − h + τ ∑ mj = a ij ˙ u njh for i = , , . . . , m, ˆ u nh = u n − h + τ ∑ mi = b i ˙ u nih , and(4.6) ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩ ˙ z ni = − ( W nih , ˙ u nih ) h for i = , , . . . , m,z ni = z n − + τ ∑ mj = a ij ˙ z nj for i = , , . . . , m,z n = z n − + τ ∑ mi = b i ˙ z ni , here Π h ∶ C ( Ω ) → S rh is the Lagrange interpolation operator, and the linearized term W ni is defined by W nih = k ∑ (cid:96) = L (cid:96) ( t n − + c i τ ) Π h W ( u n − jh ) , with k = min ( p, m + ) . Then we apply the cut-off operation: find u nh ∈ S rh such that u nh ( x j ) = min ( max ( ˆ u nh ( x j ) , − α ) , α ) , j = , . . . , M r. (4.7) Lemma 4.1.
For r = , the cut-off operation (4.7) indicates ∥∇ u nh ∥ L ( Ω ) ≤ ∥∇ ˆ u nh ∥ L ( Ω ) . (4.8) Proof . Since both ˆ u nh and u nh are piecewise linear, it is easy to see that ∥∇ u nh ∥ L ( Ω ) = h M ∑ j = ∣ u nh ( x j ) − u nh ( x j − )∣ , ∥ ˆ u nh ∥ L ( Ω ) = h M ∑ j = ∣ ˆ u nh ( x j ) − ˆ u nh ( x j − )∣ . Obviously, the cut-off operation (4.7) derives ∣ u nh ( x j ) − u nh ( x j − )∣ ≤ ∣ ˆ u nh ( x j ) − ˆ u nh ( x j − )∣ , for j = , ⋯ , M, which completes the proof.The next theorem shows that the cut-off SAV-RK scheme (4.5)-(4.7) satisfies the energy decay prop-erty and discrete maximum bound principle. Theorem 4.2.
Suppose that the Runge–Kutta method in Table 1 satisfies Assumption (P5), and weapply the lumped mass finite element method with r = in space discretization. Then, the time steppingscheme (4.5) - (4.7) satisfies the energy decay property: (4.9) 12 ∥∇ u nh ∥ L ( Ω ) + ∣ z n ∣ ≤ ∥∇ u n − h ∥ L ( Ω ) + ∣ z n − ∣ , for all n ≥ k. Meanwhile, the fully discrete solution (4.5) - (4.7) satisfies the maximum bound principle (4.10) max k ≤ n ≤ N ∣ u nh ( x )∣ ≤ α, for all x ∈ Ω . Proof . Due to the cut-off operation in each time level, we know thatmax k ≤ n ≤ N ∣ u nh ( x j )∣ ≤ α, for all j = , , . . . , M. Since the finite element function is piecewise linear, then for any x ∈ ( x j − , x j )∣ u nh ( x )∣ ≤ max (∣ u nh ( x j − )∣ , ∣ u nh ( x j )∣) ≤ α. Next, we turn to the energy decay property (4.9). According to the third relation of (4.5), we have ∇ ˆ u nh = ∇ u n − h + τ m ∑ i = b i ∇ ˙ u nih . Squaring the discrete L -norms of both sides, yields ∥∇ ˆ u nh ∥ = ∥∇ u n − h ∥ + τ m ∑ i = b i (∇ ˙ u nih , ∇ u n − h ) + τ m ∑ i,j = b i b j (∇ ˙ u nih , ∇ ˙ u njh ) . y the second relation in (4.5), we arrive at ∥∇ ˆ u nh ∥ = ∥∇ u n − h ∥ + τ m ∑ i = b i (∇ ˙ u nih , ∇ u nih − τ m ∑ j = a ij ∇ ˙ u nih ) + τ m ∑ i,j = b i b j (∇ ˙ u nih , ∇ ˙ u njh )= ∥∇ u n − h ∥ + τ m ∑ i = b i (∇ ˙ u nih , ∇ u nih ) − τ m ∑ i,j = m ij (∇ ˙ u nih , ∇ ˙ u njh )≤ ∥∇ u n − h ∥ + τ m ∑ i = b i (∇ ˙ u nih , ∇ u nih ) , where we apply the Assumption (P4) in the last inequality. Then we apply the first relation in (4.5) toderive ∥∇ ˆ u nh ∥ = ∥∇ u n − h ∥ − τ m ∑ i = b i ∥ ˙ u nih ∥ + τ m ∑ i = b i z ni ( ˙ u nih , W nih ) h On the other hand, the similar argument also leads to ∣ z n ∣ ≤ ∣ z n − ∣ − τ m ∑ i = b i z ni ( ˙ u nih , W nih ) h Therefore we conclude that12 ∥∇ ˆ u nh ∥ h + ∣ z n ∣ ≤ ∥∇ u n − h ∥ h + ∣ z n − ∣ − τ m ∑ i = b i ∥ ˙ u nih ∥ h ≤ ∥∇ u n − h ∥ h + ∣ z n − ∣ . which together with (4.8) implies the desired energy decay property immediately. Remark 4.1.
Note that the energy dissipation law holds valid only if r = , since in this case thecut-off operation does not enlarge the H semi-norm, which is present as (4.8) in Lemma 4.1. Thisproperty is not clear for finite element method with high degree polynomials. Hence, how to design anspatially high-order (unconditionally) energy dissipative and maximum bound preserving scheme is stillunclear and warrants further investigation. Next, we shall derive an error estimate for the fully discrete scheme (4.5)-(4.7). To begin with, weshall examine the local truncation error. We define the local truncation error η ni and η n as(4.11) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ u ni ∗ = ∆ u ( t ni ) + z ( t ni ) W ni ∗ for i = , , . . . , m,u ( t ni ) = u ( t n − ) + τ ∑ mj = a ij ˙ u nj ∗ + η ni for i = , , . . . , m,u ( t n ) = u ( t n − ) + τ ∑ mi = b i ˙ u ni ∗ + η n where t ni = t n − + c i τ and W ni ∗ denotes the extrapolation W ni ∗ = m ∑ (cid:96) = L (cid:96) ( t n − + c i τ ) W ( u ( t n − j )) . Similarly, we define d ni and d n as(4.12) ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩ ˙ z ni ∗ = − ( W ni ∗ , ˙ u ni ∗ ) for i = , , . . . , m,z ( t ni ) = z ( t n − ) + τ ∑ mj = a ij ˙ z nj ∗ + d ni for i = , , . . . , m,z ( t n ) = z ( t n − ) + τ ∑ mi = b i ˙ z ni ∗ + d n , Provided the assumption (P5) and relations (3.1) and (3.2), the local truncation errors η ni , η n , d ni , d n satisfy the estimate(4.13) ∥ η n ∥ H ( Ω ) + ∣ d n ∣ + τ m ∑ i = (∥ η ni ∥ H ( Ω ) + ∣ d ni ∣) ≤ Cτ k + . e omit the proof, since it is similar to the one of Lemma 3.2, given in Appendix. See also [1, Lemma3.1]. Theorem 4.3.
Suppose that the Runge–Kutta method satisfies Assumption (P4) and the relations (3.1) and (3.2) . Assume that ∣ u ∣ ≤ α and the maximum principle (1.2) holds, and assume that thestarting values u lh and z l , l = , . . . , k − , are given and ∣ u lh ( x j )∣ ≤ α, j = , . . . , M, l = , . . . , k − . Then the fully discrete solution given by (4.5) - (4.7) satisfies for n = k, . . . , N ∥ u ( t n ) − u nh ∥ L ( Ω ) ≤ C ( τ k + h ) + C k − ∑ l = ∥ u ( t l ) − u lh ∥ L ( Ω ) + C ∣ z ( t k − ) − z k − ∣ , (4.14) provided that u, f and f ( u ) are sufficiently smooth in both time and space variables.Proof . Subtracting (4.5)-(4.6) from (4.11)-(4.12), and with the notation e nih = Π h u ( t ni ) − u nih , ˙ e nih = Π h ˙ u ni ∗ − ˙ u nih , e nh = Π h u ( t n ) − u nh , ˆ e nh = Π h u ( t n ) − ˆ u nh ,ξ ni = z ( t ni ) − z ni , ˙ ξ ni = ˙ z ni ∗ − ˙ z ni , ξ n = z ( t n ) − z n . we have the following error equations(4.15) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ˙ e nih = ∆ h e nih + ( z ( t ni ) Π h W ni ∗ − z ni W nih ) + ( Π h ∆ − ∆ h Π h ) u ( t n − ) for i = , , . . . , m,e nih = e n − h + τ ∑ mj = a ij ˙ e nj + Π h η ni for i = , , . . . , m, ˆ e nh = e n − h + τ ∑ mi = b i ˙ e nih + Π h η n and(4.16) ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎩ ˙ ξ ni = − ( W ni ∗ , ˙ u ni ∗ ) + ( W nih , ˙ u nih ) h for i = , , . . . , m,ξ ni = ξ n − + τ ∑ mj = a ij ˙ ξ nj + d ni for i = , , . . . , m,ξ n = ξ n − + τ ∑ mj = b i ˙ ξ ni + d n , Now, take the square of discrete L norm of both sides of the last relation of equation (4.15), we canget ∥ ˆ e nh ∥ h = ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h + ( η n , e n − h + τ m ∑ i = b i ˙ e nih ) h + ∥ Π h η n ∥ h . (4.17)For the first term on the right hand side, we expand it and apply the second equation of (4.15) to obtain ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h = ∥ e n − h ∥ h + τ m ∑ i = b i ( ˙ e nih , e nih − η ni ) h − τ m ∑ i,j = m ij ( ˙ e nih , ˙ e njh ) h ⩽ ∥ e n − h ∥ h + τ m ∑ i = b i ( ˙ e ni , e nih − η ni ) h , where in the last inequality we use the positive semi-definiteness of the matrix M in Assumption (P4).Next, we note that the relation of (4.15) implies ( ˙ e nih , e nih − η ni ) h = ( ∆ h e nih + ( z ( t ni ) Π h W ni ∗ − z ni W nih ) + ( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h = −∥∇ e nih ∥ L ( Ω ) + (∇ e nih , ∇ Π h η ni ) + ( z ( t ni ) Π h W ni ∗ − z ni W nih , e nih − η ni ) h + (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h he bound of second term of the right hand side can be derived via Cauchy-Schwarz inequality ∣(∇ e nih , ∇ Π h η ni )∣ ≤ ∥∇ e nih ∥ L ( Ω ) + C ∥ η ni ∥ H ( Ω ) . Then the third term can be bounded as ( z ( t ni ) Π h W ni ∗ − z ni W nih , e nih − η ni ) h ≤ z ( t ni )( Π h W ni ∗ − W nih , e nih − η ni ) h + ξ ni ( W nih , e nih − η ni ) h ≤ C ( k ∑ j = ∥ e n − jh ∥ h + ∥ e nih ∥ h + ∥ Π h η ni ∥ L ( Ω ) + ∣ ξ ni ∣ ) . The bound of the last term follows from Lemma 2.4 (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih − η ni ) h ≤ Ch ∥ e nih − η ni ∥ H ( Ω ) ≤ ∥∇ e nih ∥ L ( Ω ) + C (∥ e nih ∥ h + ∥ η ni ∥ H ( Ω ) + h ) Therefore, we arrive at2 ( ˙ e nih , e nih − η ni ) h ≤ −∥∇ e nih ∥ L ( Ω ) + C ( k ∑ j = ∥ e n − jh ∥ h + ∥ e nih ∥ h + ∣ ξ ni ∣ + ∥ η ni ∥ H ( Ω ) + h ) , and hence ∥ e n − h + τ m ∑ i = b i ˙ e nih ∥ h ⩽ ∥ e n − h ∥ h − τ m ∑ i = b i ∥∇ e nih ∥ L ( Ω ) + Cτ m ∑ i = (∣ ξ ni ∣ + ∥ e nih ∥ h )+ Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ ( h + τ k ) . In view of the first relation of the error equation (4.15), we have the estimate ( η n , e n − h + τ m ∑ i = b i ˙ e nih ) h ≤ ∥ η n ∥ h ∥ e n − h ∥ h + Cτ ∥ η n ∥ H ( Ω ) m ∑ i = b i (∥∇ e nih ∥ h + k ∑ j = ∥ e n − jh ∥ h + ∣ ξ ni ∣ + h )≤ Cτ ( h + τ k ) + τ m ∑ i = b i (∥∇ e nih ∥ h + ∣ ξ ni ∣ ) + Cτ k ∑ j = ∥ e n − jh ∥ h which gives a bound of the second term in (4.17). In conclusion, we obtain that ∥ ˆ e nh ∥ h + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + Cτ m ∑ i = (∥ e nih ∥ h + ∣ ξ ni ∣ ) + Cτ k ∑ j = ∥ e n − jh ∥ h . (4.18)Similarly, from (4.16) and (4.13) we can derive ∣ ξ n ∣ ≤ Cτ ( h + τ k ) + ( + cτ )∣ ξ n − ∣ + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) + Cτ m ∑ i = (∥ e nih ∥ h + ∣ ξ ni ∣ ) + Cτ k ∑ j = ∥ e n − jh ∥ h , here we use the estimate that ( W ni ∗ , ˙ u ni ∗ ) − ( W nih , ˙ u nih ) h = ( W ni ∗ , ˙ u ni ∗ ) − ( W ni ∗ , ˙ u ni ∗ ) h + ( W ni ∗ − W nih , ˙ u ni ∗ ) h + ( W nih , ˙ e nih ) h ≤ Ch + C k ∑ j = ∥ e n − jh ∥ h ∥ Π h ˙ u ni ∗ ∥ h + (∇ W nih , ∇ e nih ) h + ( W nih , z ( t ni ) Π h W ni ∗ − z ni W nih ) h + ( W nih , ( Π h ∆ − ∆ h Π h ) u ( t n − )) h ≤ Ch + C k ∑ j = ∥ e n − jh ∥ h + C ∥∇ e nih ∥ + C ∣ ξ ni ∣ , where we use the fact that ∥∇ u nh ∥ ≤ C (by Theorem 4.2) in the last inequality. To sum up, we arrive at ∥ ˆ e nh ∥ h + ∣ ξ n ∣ + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + ( + cτ )∣ ξ n − ∣ + Cτ m ∑ i = (∥ e nih ∥ h + ∣ ξ ni ∣ ) + Cτ k ∑ j = ∥ e n − jh ∥ h . Note that ∣ e nh ( x j )∣ ≤ ∣ ˆ e nh ( x j )∣ for all j = , , . . . , M , which implies ∥ e nh ∥ h + ∣ ξ n ∣ + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + ( + cτ )∣ ξ n − ∣ + Cτ m ∑ i = (∥ e nih ∥ h + ∣ ξ ni ∣ ) + Cτ k ∑ j = ∥ e n − jh ∥ h . (4.19)Next, we shall derive a bound for ∑ mi = (∥ e nih ∥ h + ∣ ξ ni ∣ ) on the right-hand side. To this end, we testthe second relation of (4.15) by e nih . This yields m ∑ i = ∥ e nih ∥ h ≤ C ∥ e n − h ∥ h + Cτ m ∑ i,j = a ij ( ˙ e njh , e nih ) + C m ∑ i = ∥ Π h η ni ∥ h ≤ C ∥ e n − h ∥ h + Cτ m ∑ i,j = a ij ( ˙ e njh , e nih ) h + Cτ k . Then, we apply the first relation of (4.15) and Lemma 2.4 to derive m ∑ i,j = a ij ( ˙ e njh , e nih ) h = − m ∑ i,j = a ij (∇ e njh , ∇ e nih ) + m ∑ i,j = a ij ( z ( t ni ) Π h W ni ∗ − z ni W nih , e nih ) h + m ∑ i,j = a ij (( Π h ∆ − ∆ h Π h ) u ( t n − ) , e nih ) h ≤ C m ∑ i = (∥∇ e nih ∥ L ( Ω ) + ∥ e nih ∥ h + ∣ ξ ni ∣ ) + Ch + C k ∑ j = ∥ e n − jh ∥ h . Therefore, we obtain m ∑ i = ∥ e nih ∥ h ≤ C ( τ h + τ k ) + C ∥ e n − h ∥ h + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ m ∑ i = (∥∇ e nih ∥ L ( Ω ) + ∥ e nih ∥ h + ∣ ξ ni ∣ ) . imilarly, from (4.16) we can derive m ∑ i = ∣ ξ ni ∣ ≤ C ∣ ξ n − ∣ + Cτ m ∑ i,j = a ij ˙ ξ nj ξ ni + C m ∑ i = ∣ d ni ∣ ≤ C ( τ h + τ k ) + C ∣ ξ n − ∣ + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ m ∑ i = (∥∇ e nih ∥ L ( Ω ) + ∥ e nih ∥ h + ∣ ξ ni ∣ ) Sum up these two estimates and note that, for sufficiently small τ , m ∑ i = (∥ e nih ∥ h + ∣ ξ ni ∣ ) ≤ C ( τ h + τ k ) + C ∣ ξ n − ∣ + Cτ k ∑ j = ∥ e n − jh ∥ h + Cτ m ∑ i = ∥∇ e nih ∥ L ( Ω ) . Now substituting the above estimate into (4.19), we have ∥ e nh ∥ h + ∣ ξ n ∣ + τ m ∑ i = ∥∇ e nih ∥ L ( Ω ) ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + ( + Cτ )∣ ξ n − ∣ + Cτ m ∑ i = ∥∇ e nih ∥ L ( Ω ) + Cτ k ∑ j = ∥ e n − jh ∥ h . Then for sufficiently small τ , there holds ∥ e nh ∥ h + ∣ ξ n ∣ ≤ Cτ ( h + τ k ) + ∥ e n − h ∥ h + ( + Cτ )∣ ξ n − ∣ + Cτ k ∑ j = ∥ e n − jh ∥ h . Rearranging terms, we obtain (∥ e nh ∥ h + ∣ ξ n ∣ ) − (∥ e n − h ∥ h + ∣ ξ n − ∣ ) τ ≤ C ( h + τ k ) + C ∣ ξ n − ∣ + C k ∑ j = ∥ e n − jh ∥ h . Then the discrete Gronwall’s inequality impliesmax k ≤ n ≤ N (∥ e nh ∥ h + ∣ ξ n ∣ ) ≤ C ( h + τ k ) + C ∣ ξ k − ∣ + C k − ∑ j = ∥ e jh ∥ h . This completes the proof of the theorem.
5. Numerical Results.
In this section, we present numerical results to illustrate the the theoreticalresults with a one-dimensional example:(5.1) ⎧⎪⎪⎪⎪⎨⎪⎪⎪⎪⎩ ∂ t u = ∂ xx u + f ( u ) , in Ω × ( , T ] ,∂ x u = , on ∂ Ω × ( , T ] u ( x, t = ) = u ( x ) in Ω , where Ω = ( , ) and f ( u ) = ε − ( u − u ) with ε = . u ( x ) = ⎧⎪⎪⎨⎪⎪⎩ , if 0 < x < / , cos ( π ( x + )) , if 1 / ⩽ x < . The smooth initial value is chosen to satisfy the Neumann boundary condition. e solve the problem (5.1) with spatial mesh size h = / N x and temporal mesh size τ = T / N t , with T = ε and 5 ε . Throughout the section, we shall apply the Gauss–Legendre methods with m = , , k = , ,
4. We compute the numerical solution at the first k − k − r = m =
3) withsmall mesh sizes. In particular, the temporal error e τ is computed by fixing the spatial mesh size h = / τ = T / e h is computed to by fixing the temporal step size τ = T / h = / r = , , r =
1. Numerical results show the optimal rate O ( h r + ) ,which fully supports our theoretical results in Theorems 3.3 and 4.3. Temporal errors are presented in 3and 4, both of which show the empirical convergence rate O ( τ m + ) and hence coincidence to Theorems3.3 and 4.3. Table 2 e h of cut-off RK (3.3) - (3.4) and cut-off SAV-RK (4.5) - (4.7) . r / N x T
10 20 40 80 160 rateRK 0 .
01 3.03e-2 7.42e-3 1.84e-3 4.60e-4 1.14e-4 ≈ .
05 1.49e-1 1.03e-2 2.32e-3 5.71e-4 1.43e-4 ≈ .
01 4.37e-3 4.99e-4 5.90e-5 7.27e-6 9.05e-7 ≈ .
05 6.15e-2 1.64e-3 1.73e-4 2.09e-5 2.60e-6 ≈ .
01 5.10e-4 3.19e-5 1.99e-6 1.23e-7 7.74e-9 ≈ .
05 5.89e-3 1.21e-4 8.12e-6 5.03e-7 3.14e-8 ≈ .
01 3.03e-2 7.42e-3 1.84e-2 4.62e-4 1.17e-4 ≈ .
05 1.49e-1 1.03e-2 2.34e-3 5.85e-4 1.56e-4 ≈ Table 3 e τ of cut-off RK scheme (3.3) - (3.4) , with τ = T / N t . m / N t T
10 20 40 80 160 320 rate1 0 .
01 3.76e-4 9.61e-5 2.43e-5 6.10e-5 1.53e-6 3.82e-7 ≈ .
05 8.01e-4 5.36e-5 1.16e-5 2.71e-6 6.56e-7 1.61e-7 ≈ .
01 4.92e-5 6.20e-6 7.74e-7 9.65e-8 1.21e-8 1.51e-9 ≈ .
05 1.73e-2 3.60e-5 1.78e-6 2.08e-7 2.51e-8 3.08e-9 ≈ .
01 1.05e-5 6.83e-7 4.31e-8 2.71e-9 1.69e-10 1.05e-11 ≈ .
05 2.88e-2 3.66e-3 3.82e-7 1.56e-8 9.61e-10 6.06e-11 ≈ ρ n = max ≤ j ≤ Mr + ∣ u nh ( x j ) − ˆ u nh ( x j )∣ and the error of the numerical solution e ( x ) = u Nh ( x ) − u ( x, T ) . Our numerical results show that thecut-off operation is active in the computation. Meanwhile, we observe that a coarse step mesh will resultin a larger cut-off value, without affecting the convergence rate.Finally, we test the numerical results in case of relatively large time steps, and compare the numericalsolutions of extrapolated RK, cut-off RK (3.3)-(3.4), and cut-off SAV-RK schemes (4.5)-(4.7), with r = able 4 e τ of cut-off SAV-RK scheme (4.5) - (4.7) , with τ = T / N t . m / N t T
10 20 40 80 160 320 rate1 0 .
01 8.08e-3 2.23e-3 5.96e-4 1.53e-4 3.79e-5 8.78e-6 ≈ .
05 7.94e-4 1.79e-4 4.80e-5 1.24e-5 3.09e-6 7.17e-7 ≈ .
01 5.56e-9 5.95e-4 8.82e-5 1.11e-5 1.37e-6 1.65e-7 ≈ .
05 1.47e-2 5.17e-5 7.17e-6 1.00e-6 1.31e-7 1.63e-8 ≈ .
01 6.97e-11 2.56e-4 2.47e-5 1.66e-6 1.06e-7 6.60e-9 ≈ .
05 2.45e-2 2.86e-3 7.73e-7 6.16e-8 4.38e-9 2.93e-10 ≈ Fig. 1 . Error at T = . and maximal cut-off value at each time level. see Figure 2. Without the cut-off postprocessing, the numerical solutions of RK scheme significantlyexceed the maximum bound, and present oscillating solution profiles. With the cut-off operation ateach time step, the numerical solutions satisfy the maximum bound, and present reasonable solutionprofiles. However, numerical results show that the cut-off RK scheme might produce a solution with aobviously increasing and oscillating energy curve. This issue could be significantly improved by applyingthe cut-off SAV-RK method, whose solution satisfy the maximum bound and the numerical energy ismore stable. Moreover, the numerical results show that the cut-off SAV-RK scheme will produce a more m = , ε = . , T = , τ = / m = , ε = . , T = , τ = / Fig. 2 . Left: solution profiles of numerical solutions of RK, cut-off RK and cut-off SAV-RK scheme. Middle: solutionenergy of cut-off RK and cut-off SAV-RK scheme. Right: cut-off values of cut-off RK and cut-off SAV-RK scheme. regular numerical solution and smaller cut-off values, compared with the cut-off RK scheme.
Appendix A.
In this part, we sketch a proof for Lemma 3.2.
Proof . We note that the second relation in equation (3.5) implies u ( t ni ) − u ( t n − ) − τ m ∑ j = a ij u njt = τ m ∑ j = a ij ( ˙ u nj ∗ − u t ( t nj )) + η ni for i = , , . . . , m. Then we substitute the first relation of (3.5) and derive that for i = , , . . . , mu ( t ni ) − u ( t n − ) − τ m ∑ j = a ij u njt = τ m ∑ j = a ij ( k ∑ (cid:96) = L (cid:96) ( t n − + c j τ ) f ( u ( t n − (cid:96) )) − f ( t nj )) + η ni . Define ˜ η ni as the left hand side of the above relation. Now we apply Taylor’s expansion at t n − and use(3.2) to derive˜ η ni = m ∑ l = τ l ( l − ) ! ⎛⎝ c li l − m ∑ j = a ij c l − j ⎞⎠ u ( (cid:96) ) ( t n ) + m ! ∫ t ni t n − ( t ni − s ) m u ( m + ) ( s ) d s + τ ( m − ) ! m ∑ j = a ij ∫ t nj t n − ( t nj − s ) m − u ( m + ) ( s ) d s = m ! ∫ t ni t n − ( t n − s ) m u ( m + ) ( s ) d s + τ ( m − ) ! m ∑ j = a ij ∫ t nj t n − ( t nj − s ) m − u ( m + ) ( s ) d s Then we obtain the estimate for ˜ η ni , with i = , , . . . , m , that ∥ ˜ η ni ∥ H ( Ω ) ⩽ Cτ m + ∥ u ( m + ) ∥ C ([ t n − ,t n ] ; H ( Ω )) . This together with the approximation property of Lagrange interpolation lead to ∥ η ni ∥ H ( Ω ) ⩽ C ( τ k + ∥ f ( u )∥ C k ([ t n − k ,t n ] ; H ( Ω )) + τ m + ∥ u ∥ C ( m + ) ([ t n − ,t n ] ; H ( Ω )) ) . or i = , , . . . , m . Similarly, we have u ( t n ) − u ( t n − ) − τ m ∑ i = b i u nit = τ m ∑ i = b i ( k ∑ (cid:96) = L (cid:96) ( t n − + c i τ ) f ( u ( t n − (cid:96) )) − f ( t ni )) + η n . Take the left hand side as ˜ η n . Then Taylor expansion and (3.1) imply˜ η n = p ! ∫ t n t n − ( t n − s ) p u ( p + ) ( s ) d s + τ ( p − ) ! m ∑ i = b i ∫ t ni t n − ( t ni − s ) p − u ( p + ) ( s ) d s. This together with the approximation property of Lagrange interpolation leads to ∥ η ni ∥ H ( Ω ) ⩽ C ( τ k + ∥ f ( u )∥ C k ([ t n − k ,t n ] ; H ( Ω )) + τ p + ∥ u ∥ C p + ([ t n − ,t n ] ; H ( Ω )) ) . Using the choice that k = min ( p, m + ) , we derive the desired result. REFERENCES[1]
G. Akrivis, B. Li, and D. Li , Energy-decaying extrapolated RK-SAV methods for the Allen-Cahn and Cahn-Hilliardequations , SIAM J. Sci. Comput., 41 (2019), pp. A3703–A3727.[2]
S. M. Allen and J. W. Cahn , A microscopic theory for anti-phase boundary motion and its application to anti-phasedomain coarsening , Acta Metall, 27 (1979), pp. 1085–1095.[3]
D. M. Anderson, G. B. McFadden, and A. A. Wheeler , Diffuse-interface methods in fluid mechanics , Annualreview of fluid mechanics, 30 (1998), pp. 139–165.[4]
P. Brenner, M. Crouzeix, and V. Thom´ee , Single-step methods for inhomogeneous linear differential equations inBanach space , RAIRO Anal. Num´er., 16 (1982), pp. 5–26.[5]
L.-Q. Chen , Phase-field models for microstructure evolution , Annual review of materials research, 32 (2002), pp. 113–140.[6]
Q. Du, L. Ju, X. Li, and Z. Qiao , Maximum principle preserving exponential time differencing schemes for thenonlocal Allen-Cahn equation , SIAM J. Numer. Anal., 57 (2019), pp. 875–898.[7] ,
Maximum bound principles for a class of semilinear parabolic equations and exponential time differencingschemes , arXiv preprint: 2005.11465, to appear in SIAM Review, (2020).[8]
B. L. Ehle , On Pad´e approximations to the exponential function and A-stable methods for the numerical solution ofinitial value problems , ProQuest LLC, Ann Arbor, MI, 1969. Thesis (Ph.D.)–University of Waterloo (Canada).[9]
Y. Gong, J. Zhao, and Q. Wang , Arbitrarily high-order unconditionally energy stable schemes for thermodynamicallyconsistent gradient flow models , SIAM J. Sci. Comput., 42 (2020), pp. B135–B156.[10]
S. Gottlieb, D. I. Ketcheson, and C.-W. Shu , Strong stability preserving Runge-Kutta and multistep time dis-cretizations , World Scientific Press, 2011.[11]
S. Gottlieb and C.-W. Shu , Total variation diminishing runge-kutta schemes , Mathematics of computation, 67(1998), pp. 73–85.[12]
S. Gottlieb, C.-W. Shu, and E. Tadmor , Strong stability-preserving high-order time discretization methods , SIAMRev., 43 (2001), pp. 89–112.[13]
E. Hairer and G. Wanner , Solving ordinary differential equations. II , vol. 14 of Springer Series in ComputationalMathematics, Springer-Verlag, Berlin, 2010. Stiff and differential-algebraic problems, Second revised edition,paperback.[14]
L. Isherwood, Z. J. Grant, and S. Gottlieb , Strong stability preserving integrating factor Runge-Kutta methods ,SIAM J. Numer. Anal., 56 (2018), pp. 3276–3307.[15]
L. Ju, X. Li, Z. Qiao, and J. Yang , Maximum bound principle preserving integrating factor runge-kutta methods forsemilinear parabolic equations , arXiv preprint arXiv:2010.12165, (2020).[16]
B. Li, J. Yang, and Z. Zhou , Arbitrarily high-order exponential cut-off methods for preserving maximum principleof parabolic equations , SIAM J. Sci. Comput., 42 (2020), pp. A3957–A3978.[17]
H. L. Liao, T. Tang, and T. Zhou , On energy stable, maximum-principle preserving, second order bdf scheme withvariable steps for the allen-cahn equation , SIAM J. Numer. Math, 58 (2020), pp. 2294–2314.[18]
H. Liu and H. Yu , Maximum-principle-satisfying third order discontinuous Galerkin schemes for Fokker-Planckequations , SIAM J. Sci. Comput., 36 (2014), pp. A2296–A2325.[19]
X.-D. Liu and S. Osher , Nonoscillatory high order accurate self-similar maximum principle satisfying shock capturingschemes i , SIAM J. Numer. Anal., 33 (1996), pp. 760–779.[20]
A. Ostermann and M. Roche , Runge-Kutta methods for partial differential equations and fractional orders ofconvergence , Math. Comp., 59 (1992), pp. 403–420.2821]
J. Qiu and C.-W. Shu , Runge–Kutta discontinuous Galerkin method using WENO limiters , SIAM J. Sci. Comput.,26 (2005), pp. 907–929.[22]
A. Quarteroni, R. Sacco, and F. Saleri , Numerical mathematics , vol. 37 of Texts in Applied Mathematics,Springer-Verlag, New York, 2000.[23]
J. Shen, T. Tang, and J. Yang , On the maximum principle preserving schemes for the generalized Allen-Cahnequation , Commun. Math. Sci., 14 (2016), pp. 1517–1534.[24]
J. Shen and J. Xu , Convergence and error analysis for the scalar auxiliary variable (SAV) schemes to gradient flows ,SIAM J. Numer. Anal., 56 (2018), pp. 2895–2912.[25]
J. Shen, J. Xu, and J. Yang , The scalar auxiliary variable (sav) approach for gradient flows , Journal of Computa-tional Physics, 353 (2018), pp. 407–416.[26]
J. Shen, J. Xu, and J. Yang , A new class of efficient and robust energy stable schemes for gradient flows , SIAMReview, 61 (2019), pp. 474–506.[27]
T. Tang and J. Yang , Implicit-explicit scheme for the Allen-Cahn equation preserves the maximum principle , J.Comput. Math., 34 (2016), pp. 471–481.[28]
V. Thom´ee , Galerkin Finite Element Methods for Parabolic Problems , Springer-Verlag, Berlin, second ed., 2006.[29]
J. J. W. van der Vegt, Y. Xia, and Y. Xu , Positivity preserving limiters for time-implicit higher order accuratediscontinuous Galerkin discretizations , SIAM J. Sci. Comput., 41 (2019), pp. A2037–A2063.[30]
Z. Xu , Parametrized maximum principle preserving flux limiters for high order schemes solving hyperbolic conserva-tion laws: one-dimensional scalar problem , Math. Comp., 83 (2014), pp. 2213–2238.[31]
P. Yue, J. J. Feng, C. Liu, and J. Shen , A diffuse-interface method for simulating two-phase flows of complexfluids , Journal of Fluid Mechanics, 515 (2004), p. 293.[32]