We consider a family of variable time-stepping Dahlquist-Liniger-Nevanlinna (DLN) schemes, which is unconditionally non-linear stable and second order accurate, for the Allen-Cahn equation. The finite element methods are used for the spatial discretization. For the non-linear term, we combine the DLN scheme with two efficient temporal algorithms: partially implicit modified algorithm and scalar auxiliary variable algorithm. For both approaches, we prove the unconditional, long-term stability of the model energy under any arbitrary time step sequence. Moreover, we provide rigorous error analysis for the partially implicit modified algorithm with variable time-stepping. Efficient time-adaptive algorithms based on these schemes are also proposed. Several one- and two-dimensional numerical tests are presented to verify the properties of the proposed time-adaptive DLN methods.
]]>The Allen-Cahn equation was first introduced by Allen and Cahn in [2] to describe the motion of anti-phase boundaries in crystalline solids. Given the domain Ω ⊂ R d (d = 1, 2, 3), u(x, t) models the difference between the concentrations of two mixtures' components at time t, and is governed by
where is the interfacial parameter. The Helmholtz free-energy density F(u) and its derivative takes the form
The model energy of the Allen-Cahn equation ( 1) is
which satisfies the following energy dissipation law:
Phase field models, including the Allen-Cahn equation, have been widely applied to complex moving interface problems in materials science and fluid dynamics. Over the past forty years, numerical simulations of phase field models have been extensively studied, leading to a wide array of spatial discretization techniques including finite element methods, alongside other approaches such as finite difference, finite volume, and Fourier-spectral methods. Various splitting techniques and stabilizers, coupled with classical time integrators, have been proposed to handle the potential functions in phase field models, resulting in numerous efficient and stable algorithms, as summarized in the recent survey paper [12].
Convex splitting technique is widely employed to ensure the discrete energy dissipation law [15,18,38,43,47]. The main idea is to split the potential function into the convex and concave parts, treated implicitly and explicitly, respectively. Linearly-implicit stabilizers are another popular way for gradient flows, leading to linear solvers at each time step and ensuring unconditional energy stability [31,45,49,50]. Exponential integrators, transferring phase field equations into equivalent integral formulations, also achieve energy stability [3,14,25]. Additionally, Lagrange multiplier approach is an alternative way to construct stable numerical schemes [19], and two new classes of time-stepping schemes, invariant energy quadratization (IEQ) [34,51,52] and scalar auxiliary variable (SAV) schemes [1,4], have been recently proposed. The variable step backward differentiation formulation (BDF2) method has been thoroughly investigated to show energy stability for the Allen-Cahn equation, under the ratio of successive time steps satisfying τ < (3 + √ 17)/2 [33]. Fully discrete interior penalty discontinuous Galerkin methods for phase field models were developed in [16,17] to derive error estimates dependent on the polynomial order of the reciprocal model parameter (1/ ) rather than the exponential order. Reduced-order finite element methods, based on the proper orthogonal decomposition method, can improve efficiency by reducing the dimension of the solution space while preserving the discrete energy dissipation law [32]. Explicit hybrid finite difference and operator splitting methods have demonstrated the ability to obtain pointwise boundedness of the numerical solutions [24].
In this paper, we consider a family of two-step time-stepping schemes proposed by Dahlquist, Liniger and Nevanlinna (hence the DLN method) which is unconditionally Gstable (non-linear stable) and second order accurate under non-uniform time grids [8][9][10][11]. To the best of our knowledge, the DLN method is one of the few multi-step schemes possessing these two properties and thus provides a great potential in the simulations of differential equations. Recently the DLN method has been applied to various fluid models and its fine properties have been confirmed by some benchmark test problems [27,[35][36][37]44]. Moreover, the DLN implementation can be simplified by adding a few lines of code to the commonlyused backward Euler scheme [28]. To utilize the unconditional G-stability to a large extent, some adaptive DLN algorithms have been proposed to balance the time accuracy and computational costs [29].
Given the initial value problem:
with y: [0, T ] → R d , g: [0, T ]×R d → R d and y 0 ∈ R d , the two-step, variable time-stepping DLN method, with parameter θ ∈ [0, 1] and starting points y 0 , y 1 , takes the form Here we denote
as the time grid on the interval [0, T ], with k n = t n+1t n as the time step, k n = α 2 k n -α 0 k n-1 as the average time step, and y n as the discrete DLN solution approximating y(t n ). The coefficients in (DLN) take the form
The step variability ε n = (k nk n-1 )/(k n + k n-1 ) represents the variation between two consecutive step sizes. On uniform time grids, step variability ε n reduces to 0, and the coefficients [β
. Various efficient ways to address the non-linear term of the Allen-Cahn model (1) in pursuit of unconditional energy stability have been designed. In the paper, we consider the implicit modified Crank-Nicolson method first, followed by the SAV method, which have been widely studied in the last decade. The implicit modified Crank-Nicolson scheme has been shown to provide unconditional energy stability in the Allen-Cahn model [7,13,42,48]. The SAV scheme for gradient flows [39][40][41] achieves unconditional stability of the model energy by introducing an auxiliary scalar function, which eliminates the need to solve non-linear algebraic equations. In this context, we investigate the partially implicit modified DLN algorithm and the DLN-SAV algorithm, combined with the finite element spatial discretization, for robust simulations of the Allen-Cahn model (1). The main contributions of this work are to:
• Construct the variable step modified DLN and DLN-SAV finite element algorithms based on the DLN refactorization process (simplifying the DLN implementation by adding filters on the backward Euler scheme); • Provide unconditional stability of model energy under arbitrary time grids for both modified DLN and DLN-SAV schemes; • Present a rigorous proof that the fully discrete modified DLN algorithm is 1st order accurate in time under arbitrary time steps and 2nd order accurate in time under uniform time grids;
• Carry out a rigorous error estimate of the fully discrete DLN-SAV algorithm via a novel way to analyze the nonlinear term; • Design time adaptive algorithms for both modified DLN and DLN-SAV schemes by the local truncation error (LTE) criterion: utilizing certain explicit time-stepping methods to estimate the LTE and adjusting time step size based on the estimator.
The rest of the paper is organized as follows. We provide necessary preliminaries and notations in Section 2. In Section 3, we study the partially implicit modified DLN scheme. We prove that the discrete model energy is unconditionally stable under variable time steps. Error estimate of the resulting fully discrete method is also provided. In Section 4, we propose the variable time-stepping DLN-SAV scheme. We prove that the corresponding model energy satisfies the discrete energy dissipation law. Moreover, we discuss how the DLN-SAV algorithm can be implemented by a simplified refactorization process: adding a few lines of code on the backward Euler-SAV (BE-SAV) scheme. Optimal error estimate of the proposed DLN-SAV algorithm in the H 1 norm is also carried out. In Section 5, the corresponding time adaptive algorithms (based on LTE criterion) are presented to improve computational efficiency. We estimate the LTE by the numerical solutions of the explicit AB2-like method (a revised two-step Adams-Bashforth method) with little or almost no extra computational costs. Several one-and two-dimensional numerical tests are presented in Section 6 to confirm the main conclusions of the paper. Conclusion remarks are provided in Section 7. Proofs of some lemmas in Section 2 are given in the appendices.
We start by introducing some necessary notations. Let H r (Ω) be the usual standard Sobolev space W r ,2 (Ω) with the norm • r . When r = 0, this reduces to the inner product space L 2 (Ω) with the L 2 norm • and the standard L 2 inner product (•, •). W r ,∞ (Ω) is the function spaces containing all functions with essentially bounded weak derivatives up to order r . When r = 0, the space is reduced to L ∞ (Ω) space with norm • ∞ . For spatial discretization, we denote X h s ⊂ H 1 (Ω) as a standard finite element space on Ω with the highest polynomial degree s and a mesh diameter h > 0. Throughout the paper, C denotes a generic positive constant independent of h and time step size k n . For any given function u(x, t) on time grids {t n } N n=0 , u n = u n (x) represents the function u(•, t n ) at time t n . For an arbitrary sequence {z n } ∞ n=0 , we define
with α ℓ and β (n) ℓ defined in (4). The following lemmas on stability and consistency of the DLN method (DLN) would play a key role in the numerical analysis.
and the coefficients γ (n) ℓ are given by
The derivation of ( 7) follows from algebraic calculation. We refer to [11,27] for details, and omit the derivation here to save space. Following this, the variable time-stepping DLN method in (DLN) is unconditionally G-stable. (We refer to [8,9] for the definition of Gstability for multi-step, time-stepping methods.) Lemma 2 Let u(•, t) be the mapping from [0, T ] to L 2 (Ω), and u n be the function u(•, t n ) in L 2 (Ω). Assuming the mapping u(•, t) is second-order differentiable in time, then for any θ ∈ [0, 1), we have
When θ = 1, the DLN method reduces to the midpoint rule, and we have
Under uniform time grids with constant time step k, (11) becomes
Lemma 3 Let u(•, t) be the mapping from [0, T ] to H r (Ω), and u n be the function u(•, t n ) in H r (Ω). Assuming the mapping u(•, t) is third-order differentiable in time, then for any θ ∈ [0, 1), we have
When θ = 1, the DLN method reduces to the midpoint rule, and we have
The proof of these two Lemmas can be found in Appendix A and B.
Let u h n ∈ X h s be the numerical approximation of u(x, t) in (1) at time t n . The modified DLN finite element algorithm for the Allen-Cahn equation ( 1) is formulated as the following: given
holds for all v h ∈ X h s , where u h n,α , u h n,β and u h n,θ are defined in (5) and
Remark 1 The function f (u) is approximated by f numerically, which has the following properties.
Note that the solution u(x, t) of the Allen-Cahn model (1) satisfies the maximum principle: if the initial value and boundary conditions are bounded by some constant Γ , then the solution at a later time is also bounded by Γ , i.e.,
Following the maximum principle, there exists
) is a first order approximation to f (u(t n,β )) in time under variable time steps and a second order approximation in time under constant time steps. This can be observed from (11), (14) and
which follows from Taylor's expansion, where ξ n is between u n and u n,β . For θ = 1, f (u n+1,θ , u n,θ ) is always a second order approximation to f (u(t n,β )) in time, which is implied by (13) and (20).
where
The DLN-CSS scheme is obtained by replacing f (u h n+1,θ , u h n,θ ) in (19) with f in (21). We test the DLN-CSS scheme on the 1D traveling wave example in Section 6.1. The numerical results are similar to those of the modified DLN scheme, and are omitted to save space.
We define the corresponding discrete model energy for the modified DLN algorithm (19) at time t n as
The following theorem presents the discrete energy dissipation law satisfied by the numerical solution.
Theorem 1 The model energy of the variable time-stepping modified DLN algorithm (19) satisfies the discrete energy dissipation law:
thus the modified DLN algorithm (19) is unconditionally stable in model energy.
following the G-stability identity (7). Then we use the fact
which implies (23).
This subsection is devoted to the error estimate of the fully discrete modified DLN algorithm (19).
We remind that {β
ℓ=0 are uniformly bounded for any fixed θ ∈ [0, 1], i.e. |β (n) ℓ | < C β (θ ) for all ℓ, n, and define the following discrete Bochner function space for error analysis ℓ 2 θ (0, N ; H r (Ω))
Theorem 2 Let u h n ∈ X h s denote numerical solutions of the modified DLN algorithm in (19) and u n denote exact solutions of the Allen-Cahn equation in (1) at time t n . We assume
there exists L > 0 such that f (u) ∞ < L for all t ∈ [0, T ], and time step sizes and time step ratios satisfy
for some C ℓ , C u > 0. Moreover we assume that the initial two solutions u h 0 and u h 1 satisfy
Then we have: for θ ∈ [0, 1)
Under uniform time grids with constant time step k, (28) becomes
Proof We start by defining the following Ritz projection operator Π:
The following approximation theorem [6,46] holds for the Ritz projection: for any v ∈ H r (Ω) and its Ritz projection Π(v) ∈ X h s , we have
We define the error of u at time t n to be
The exact solution u at t n,β satisfies
Subtracting ( 19) from (33) yields
The definition of Π leads to ∇η n,β , ∇v h = 0, and (34) can be rewritten as
We set v h = φ h n,β in (35) and use the G-stability identity in (7) to obtain
Next, we provide a bound for each term on the right-hand side of (36). Utilizing ( 16), Young's inequality and the fact |β
where δ > 0 will be decided later and C(θ, δ, L) is a constant only depending on θ, δ, L. By Cauchy-Schwarz inequality, Young's inequality and the projection error (32)
By HR older's inequality, we have
therefore, (38) becomes
By Cauchy-Schwarz inequality and ( 15)
For the non-linear term, utilizing Cauchy Schwarz inequality, triangle inequality and Young's inequality leads to
Considering the first term on the right-hand side, we have
By the assumption f (u) ∞ < L, f (•, •) are Lipschitz-continuous for both components [42]. Thus
where the last inequality follows from HR older's inequality as in (39). Similarly, we have
We use ( 10)-( 11) in Lemma 2 and ( 20) to obtain
Combining ( 42) -( 46), we have
Now, we are ready to derive the error estimate. We combine (36), ( 37), ( 40), ( 41) and ( 47), multiply both sides by k n and use the time ratio condition in (27) to obtain
We sum (48) over
2 and use the time ratio condition in ( 27) to obtain
To estimate φ h N 2 , we need
to relax the upper bound for k max and obtain the time step restriction in (26). By the discrete GrR onwall inequality [23] and the fact that the maximum time step k max < 1, we have
which leads to
By (51), triangle inequality and (32), we have (28). Lastly, we comment on the case of uniform time grids with constant time step k. Under such case, utilizing (10) and ( 14) in Lemma 2, Eq. ( 46) becomes
Following similar arguments, we have (30).
We have proven that the modified DLN algorithm is first-order accurate in time under arbitrary time step sequences. However, all numerical tests in Section 6 demonstrate that numerical solutions on both uniform and non-uniform time grids converge at secondorder in time.
In this section, we consider the energy-stable SAV approach, and present the DLN-SAV algorithm. For the Allen-Cahn model ( 1), we introduce the scalar auxiliary variable of the form
where λ > 0 is a small positive number and C 0 > 0 is picked to ensure r (t) > C r > 0 for some constant C r . Then the Allen-Cahn model in (1) can be equivalently written as
where g(u) = f (u) -λu.
Let u h n , r h n be fully-discrete numerical solutions of u(x, t n ), r (t n ) respectively. The fullydiscrete weak formulation for the DLN-SAV algorithm with the finite element spatial discretization is: given
123
where u h n, * defined in ( 6) is the explicit second-order approximation to u(t n,β ) in time. The variable time-stepping DLN scheme can be simplified by the refactorization process: adding a few lines of codes to the backward Euler (BE) scheme to obtain the DLN solutions (see [28] for details). The DLN-SAV algorithm can be equivalently rewritten as the following process:
Step 1. Pre-process
, where the coefficients for the pre-process are
Step 2. BE solver to solve u h,temp n+1 ∈ X h s and r h,temp n+1 such that for all
Step 3. Post-process
, where the coefficients for the post-process are
In practice, the scalar auxiliary variable r need not be refactorized in Step 3, since its effect on accuracy can be negligible [53]. The above refactorization process can be applied to the modified DLN algorithm in (19) as well. However, the refactorization process of the modified DLN algorithm does not simplify the implementation much since the partially implicit, nonlinear term needs to be refactorized as well, and was omitted in Section 3.
Proof Since the DLN-SAV algorithm in ( 53) is equivalent to the above refactorization process, it suffices to show that the BE solver in (54) has a unique solution. We denote
and eliminate r h,temp n+1 in (54) to obtain
Let {ψ ℓ } N d ℓ=1 be the basis for X h s . We have u h,temp n+1
and vectors
Then the algebraic structure for ( 55) is
Since the mass matrix M > 0, the stiffness matrix K ≥ 0 and B (n) (B (n) ) T ≥ 0, the linear system in (56) has a unique solution. In addition, r h,temp n+1 is uniquely obtained from the second equation of (54).
The left matrix in (56) is far away from a sparse matrix. To efficiently implement the DLN-SAV scheme, we solve (55) in the following simpler way. We denote A = I -2 k BE n Δ > 0 and apply integration by parts to (55) to rewrite it as
To solve u h,temp n+1 in (57) efficiently, we first solve for
s to be the unique finite element solution to the following second-order system
Thus A -1 G ∈ X h s satisfies the following weak form
We use the definition (58) and integration by parts to obtain
We use (59) and set v h = A -1 ϕ h n in (57) to solve
To obtain u h,temp n+1 from (57), it suffices to solve
or equivalently,
We summarize the simplified implementation in the following steps.
(i) Solving for 60). (iv) Computing u h,temp n+1 from (61) and r h,temp n+1 from the second equation of (54). ) obtained by the above process is the unique solution to the algorithm in (54), and the details are omitted here.
Remark 6 Solving two second-order equations in (62) and ( 63) is much more efficient than solving (56) directly since the left matrix of the two linear systems in (62) and (63) becomes a sparse matrix M + 2 k BE n K .
We define the discrete model energy for the DLN-SAV algorithm (53) at time t n as
and have the following theorem on its unconditional stability with respect to E SAV n .
thus the DLN-SAV algorithm is unconditionally stable with respect to this model energy.
Proof We set v h = u h n,α in the first equation of the DLN-SAV algorithm (53)
We multiply the second equation of the DLN-SAV algorithm (53) by 2 k n r n,β and obtain
We combine (66) and (67) and use the G-stability identity (7) to obtain
which results in the unconditional stability property (65).
The inequality in (65) is the discrete version of the energy dissipation law (2) and the discrete energy E SAV n in ( 64) is an approximation of the model energy E.
In this subsection, we investigate the error estimate of the fully discrete DLN-SAV algorithm (53). Let u n , r n be exact solutions u(x, t), r (t) evaluated at t n respectively. The errors of numerical solutions of the DLN-SAV algorithm in (53) at time t n are denoted as
Introduce the following discrete Bochner function space for error analysis
Since t n,β is a point between t n-1 and t n+1 , the above discrete norm can be understood as the Riemann sum of √ 2 u L 2 (0,T ;L 2 (Ω)) . To carry out the error analysis, one key idea is to use the mathematical induction to prove the following statement
for any M between 0 and N , under certain time step assumptions. To start, we assume that (70) holds for n = 0, 1, • • • , M -1. Since the exact solution u satisfies the maximum principle and is bounded, it can be derived that { u h n ∞ } M-1 n=0 uniformly bounded following the triangle inequality. Next, we use the Ritz projection operator to decompose the error e u n as
Applying inverse inequality, triangle inequality and approximation theorem in (32) leads to
Thus for n = 0, 1,
and the diameter h is small enough.
Then we combine all these facts to achieve
which will be summarized in Theorem 4 below. Therefore, we can verify (70) for n = M under the time step restriction
The main error analysis of the DLN-SAV method is summarized as follows.
Theorem 4 We assume
the approximations of the DLN-SAV algorithm satisfy
the maximum time step has bound k max ≤ min{1, h 1/2 } and the time step ratios satisfy
for some C ℓ , C u > 0. Moreover we assume that the initial two DLN-SAV solutions u h 0 and u h 1 satisfy
Then we have
Proof The exact solutions of (52) satisfy
We subtract the first equation of (76) evaluated at t n,β from the first equation of DLN-SAV algorithm (53) and apply integration by parts to obtain
where
Following the error decomposition e u n = η n + φ h n , Eq. (77) becomes
where the fact 2 (∇η n,β , ∇φ h n,α ) = 0 is utilized in (79), due to the definition of Ritz projection.
We subtract the second equation of (76) evaluated at t n,β from the second equation of ( 53) and use the notation in (78) to derive
Setting v h = φ h n,α in (79), multiplying (80) by e r n,β , adding these two equations together and utilizing the G-stability identity (7) yield
Now we deal with each term on the right side of (81). By similar argument to ( 38)-( 40)
By the approximation theorem in (32), triangle inequality and Lemma 3
By Cauchy-Schwarz inequality, Young's inequality and Lemma 3
By algebraic calculation,
where
By maximum principle of exact solutions, we have that r (t) is uniformly bounded for all t,
We apply approximation theorem (32), triangle inequality, Lemma 3 to (88)
where u n, * is the second-order approximation to u(t n,β ) in time is used in the last inequality. Similarly,
We combine (89) and (90) to obtain
By the uniform boundedness of r (t), u(t) ∞ for all t, Cauchy-Schwarz inequality and Lemma 3
Similar to (89)-( 92)
By the assumption in (72), { u h n, * ∞ } M-1 n=1 is uniformly bounded and
where δ > 0 is a certain constant to be decided later. We use the assumption in (72) again and the approximation theorem (32) to obtain
By triangle inequality, Young's inequality and the fact that u n, * is the second-order approximation to u(t n,β ) in time
We combine (95), ( 96), (97) to have
By the maximum bound principle for the exact solution and Hölder's inequality
By the maximum bound principle for the exact solution and Lemma 3
We combine (81) -( 86), ( 91) -( 94), ( 98) -(101) and drop the numerical dissipation terms
-1 u t 2 s+1 dt + Ch 2s+2 k 4 max t n+1 t n-1 u tt 2 s+1 dt +(k n +k n-1 ) u(t n,β ) 2 s+1 + Ck 4 max t n+1 t n-1 u ttt 2 dt + Ck 4 max t n+1 t n-1 Δu tt 2 dt + C(δ)k 4 max t n+1 t n-1
We sum (102) over n from 1 to N -1 and use time step ratio bounds in (73)
We set δ = (1 + θ)/(12C 2 β (θ )) and apply the discrete Grönwall inequality [23] with the time step restriction k max < 1 to (103)
We use the triangle inequality, approximation theorem in (32), (104) and initial conditions in (74) to have (75).
Remark 8 From (104), we can conclude that e u M ≤ O(k 2 max , h s+1 ) achieves the optimal error estimate in the L 2 norm, if the initial two step solutions match the exact solutions (leading to zero φ h 0 and φ h 1 ). This is also verified by Table 4 in the one-dimensional traveling wave test problem in Subsection 6.1.
In this section, we discuss adaptive DLN algorithms for solving the Allen-Cahn equation, to better take advantage of this variable time-stepping method. Here we consider adaptive DLN algorithms by using the LTE criterion, which involve two essential parts:
• the estimator for LTE, • the time step controller.
An explicit, variable step AB2-like scheme 1 The deviation of this scheme is similar to twostep Adam-Bashforth scheme, thus we name it an AB2-like scheme. is adopted to estimate the LTE of the DLN method. The AB2-like time-stepping scheme, applying to the initial value problem in (3), takes the form
where g DLN (t n-1,β , y n-1,β ) and g DLN (t n-2,β , y n-2,β ) are calculated by the DLN scheme in (DLN)
Thus the AB2-like solution y AB2-like n+1 in ( 105) is just the linear combination of the previous four DLN solutions {y n , y n-1 , y n-2 , y n-3 }, and thus obtained with low computational costs.
We denote time step ratio τ n = k n /k n-1 and utilize the AB2-like scheme (105) to estimate the LTE of the DLN method via
where G (n) , R (n) are coefficients before the LTEs of the DLN scheme and AB2-like scheme respectively:
We refer to [29] for the derivation of the explicit AB2-like scheme in (105) and the estimator of LTE in (106).
For the time step controller, there are many choices and we consider the following one proposed in [20]
where κ ∈ (0, 1] is the safety factor and Tol is the required tolerance for the LTE. We summarize the adaptive DLN algorithm using the estimator of LTE in (106) and the step controller in (108) in the following pseudocode.
Algorithm 1: Adaptive DLN method (LTE estimator by AB2-like scheme) Input: tolerance Tol, four previous DLN solutions u h n , u h n-1 , u h n-2 , u h n-3 , current time step k n , three previous time step k n-1 , k n-2 , k n-3 , safety factor κ ; Compute the DLN solution u h,DLN n+1 by (19) or (53) ; Compute the AB2-like solution u h,AB2-like n+1 by (105) ; Use k n
AB2-like n+1 u h,DLN n+1 ; // relative estimator if T n+1 < Tol then u h n+1 ⇐ u h,DLN n+1 ; // accept the result k n+1 ⇐ k n • min 1.5, max 0.2, κ Tol T n+1 1/3 ; // adjust step by (108) else // adjust current step to recompute k n ⇐ k n • min 1.5, max 0.2, κ Tol T n+1 1/3 ;
In this section, we test the performance of the proposed modified DLN algorithm and DLN-SAV algorithm on three numerical tests:
• Accuracy test on the 1D traveling wave solution;
• Accuracy test on the 2D test with known solutions;
• Simulation of the 2D test with random initial conditions.
We implement these algorithms with three different values of θ : 2/3, 2/ √ 5, 1. Among these values, θ = 2/3 is proposed in [11] to balance the magnitude of LTE and keep fine stability properties. The value θ = 2/ √ 5 is suggested in [26] to have the best stability at infinity. The algorithm with value θ = 1 reduces to the classical one-step midpoint rule. In the implementation, we use MATLAB for the 1D test problem and the software FreeFem++ [22] for two 2D test problems. For the modified DLN algorithm (19), we use fixed point iteration to solve the non-linear system at each time step.
We use the 1D traveling wave problem [5] with the known solution to test the accuracy of modified DLN and DLN-SAV algorithms. One traveling wave solution of Allen-Cahn equation (1) in 1D takes the form
where the traveling speed is s = 3 / √ 2. We set the final time as T = 2, the model parameter as = 0.01, and use the inhomogeneous Dirichlet boundary condition. We use P 2 finite element space in the spatial discretization.
Table 1 L -norm of error and rate for 1D modified DLN scheme in time
.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99 θ = 2 √ 5 0.04 2.17e-3 -2.41e-3 -4.80e-1 -0.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99 θ = 1 0 .04 2.17e-3 -2.41e-3 -4.80e-1 -0.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99
First, we consider the constant time step k and set h = k 2 to check the accuracy of both algorithms in time. From Tables 1 and 3, we observe that both modified DLN and DLN-SAV algorithms have second-order accuracy in time for the 1D traveling wave problem, confirming the expected time-stepping accuracy of these algorithms. Next, we set k = h 2 to check the accuracy of both algorithms in space. From Tables 2 and 4, we can see that both algorithms demonstrate third-order spatial convergence for ℓ ∞ (L 2 )-norm and second-order spatial convergence for ℓ 2 (H 1 )-norm.
Table 3 L 2 -norm of error and for 1D DLN-SAV scheme in time
.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99 θ = 2 √ 5 0.04 2.17e-3 -2.41e-3 -4.80e-1 -0.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99 θ = 1 0 .04 2.17e-3 -2.41e-3 -4.80e-1 -0.02 2.91e-4 2.90 3.97e-4 2.60 1.32e-1 1.86 0.01 3.69e-5 2.98 5.19e-5 2.94 3.39e-2 1.97 0.005 4.65e-6 2.99 6.57e-6 2.98 8.53e-3 1.99
We utilize constant time step k as a reference time step and set h = k 2 to test the accuracy of the variable time-stepping modified DLN algorithm (19) and DLN-SAV algorithm (53) over time. Two different time step scenarios are considered:
1. Randomized Time Steps: For each time step, we set k n = k + k • rand, where rand is a random number drawn from the uniform distribution in the interval (0,1).
We alternate the time step size as follows:
Table 5 L 2 -norm of error and rate for 1D modified scheme in time (k n
1 0.1919 2.20e-4 -2.06e-4 -3.45e-2 -0.05 0.0991 5.98e-5 1.97 5.33e-5 2.04 3.45e-3 3.48 0.04 0.0771 3.36e-5 2.31 3.02e-5 2.28 1.80e-3 2.61 0.02 0.0399 1.01e-5 1.82 8.58e-6 1.91 4.50e-4 2.10 θ = 2 √ 5 0.1 0.1995 1.85e-4 -1.48e-4 -3.34e-2 -0.05 0.0992 5.61e-5 1.71 4.57e-5 1.68 2.71e-3 3.60 0.04 0.0789 3.56e-5 1.99 2.72e-5 2.26 1.32e-3 3.13 0.02 0.0398 8.16e-6 2.15 6.37e-6 2.12 2.37e-4 2.51 θ = 1 0 .1 0.1738 1.53e-4 -1.30e-4 -3.32e-2 -0.05 0.0994 4.81e-5 2.07 3.55e-5 2.32 2.40e-3 4.70 0.04 0.0787 3.01e-5 2.00 2.28e-5 1.91 1.15e-3 3.22 0.02 0.0391 7.86e-6 1.92 6.03e-6 1.90 2.02e-4 2.47
For both test scenarios, the time steps lie within the range [k, 2k]. The accuracy rate is calculated using the formula:
where Error 1 and Error 2 are the errors corresponding to the maximum time steps k max,1 and k max,2 , respectively. From Tables 5678, our observation indicates that both algorithms with all θ values achieve second-order accuracy in time for the 1D traveling wave problem, even under non-uniform time steps. This is better than the first order convergence in time that we can prove in Theorem 2 for the non-uniform time step case.
To further test the order of accuracy of the modified DLN scheme, we apply the modified DLN method on the corresponding ODE dy/dt + f (y) = 0, assuming y = y(t) is only timedependent. We have carried out extensive numerical investigation of this test with different ways to determine the non-uniform time step size and different choices of f (u). All of these tests demonstrate a second-order convergence numerically.
In this subsection, we compare the modified DLN algorithm with the DLN-CSS algorithm in (21) in terms of accuracy and efficiency. Both algorithms are implicit and need certain non-linear solver at each time step. From Tables 9 and 10, both algorithms reach second-order accuracy in time and take about the same CPU time to complete simulation. However, the modified DLN algorithm is slightly more accurate than the DLN-CSS algorithm, possibly due to the fact that the modified DLN algorithm addresses the whole non-linear term implicitly.
Next, we test the accuracy and efficiency of the adaptive DLN algorithm in Algorithm 1, by comparing it with the corresponding constant time step algorithms. For adaptive algorithms,
Table 6 L 2 -norm of error and rate for 1D modified DLN scheme in time (k n
1 0.1885 6.25e-4 -4.85e-4 -4.16e-2 -0.05 0.1000 1.62e-4 2.13 1.31e-4 2.07 6.75e-3 2.87 0.04 0.0781 1.20e-4 1.20 9.35e-5 1.36 4.80e-3 1.37 0.02 0.0399 3.10e-5 2.02 2.45e-5 1.99 1.22e-3 2.04 θ = 2 √ 5 0.1 0.1931 7.12e-4 -5.64e-4 -4.10e-2 -0.05 0.0987 2.20e-4 1.76 1.66e-4 1.82 7.20e-3 2.59 0.04 0.0744 1.28e-4 1.89 9.68e-5 1.91 4.05e-3 2.03 0.02 0.0400 3.74e-5 1.98 2.80e-5 2.00 1.15e-3 2.03 θ = 1 0 .1 0.1996 8.19e-4 -6.41e-4 -4.13e-2 -0.05 0.0998 2.58e-4 1.67 2.04e-4 1.65 7.94e-3 2.38 0.04 0.0770 1.38e-4 2.42 1.05e-4 2.57 4.04e-3 2.62 0.02 0.0396 3.86e-5 1.91 2.96e-5 1.90 1.11e-3 1.93
the following setup is adopted: the tolerance Tol = 1.e -6, the minimum time step k min = 1.e -5, the maximum time step k max = 0.1, two initial time steps k 0 = k 1 = 1.e -3, and the safety factor κ = 0.8. For constant step algorithms, we set the constant time step as k = 1.e -3. From Tables 11 and 12, we observe that both adaptive algorithms outperform the constant step algorithms. Adaptive algorithms take much fewer time steps to obtain the same accuracy for all three θ values, when compared with the corresponding constant step algorithms.
Table 8 L 2 -norm of error and rate for 1D DLN-SAV scheme in time (k n = k, 2k, k, 2k, ...; h = k 2 )
04 1.84e-5 -1.05e-3 -72.18 0.02 4.64e-6 1.98 1.58e-4 2.73 431.92 0.01 1.17e-6 1.99 3.74e-5 2.08 3494.29 θ = 2 √ 5 0.04 1.45e-5 -9.38e-4 -61.00 0.02 3.68e-6 1.98 1.05e-4 3.16 452.98 0.01 9.27e-7 1.99 2.27e-5 2.21 3506.87 θ = 1 0 .04 1.27e-5 -9.19e-4 -57.17 0.02 3.22e-6 1.98 9.45e-5 3.28 443.27 0.01 8.12e-7 1.99 1.97e-5 2.26 3519.10
3 3.69e-5 3.67e-5 2.39e-2 1000 θ = 2 √ 5 3.69e-5 3.67e-5 2.39e-2 1000 θ = 1 3.69e-5 3.67e-5 2.39e-2 1000 -
3 3.68e-5 3.67e-5 2.39e-2 174 65 θ = 2 √ 5 3.67e-5 3.67e-5 2.39e-2 146 36 θ = 1 3.67e-5 3.67e-5 2.39e-2 140 43 Constant time-stepping DLN-SAV schemes θ u -u h ℓ ∞ (L 2 ) u -u h ℓ 2 (L 2 ) u -u h ℓ 2 (H 1 ) # Steps # Rejections θ = 2 3 3.68e-5 3.67e-5 2.39e-2 1000 θ = 2 √ 5 3.68e-5 3.67e-5 2.39e-2 1000 θ = 1 3.68e-5 3.67e-5 2.39e-2 1000 -
We consider a two-dimensional example to verify that the proposed DLN algorithms are accurate and efficient for 2D test problems. As studied in [30], we consider the exact function of the form u(x, y, t) = 0.05e -0.1t sin(x) sin(y), (x, y)
with the homogeneous boundary condition. Since this function does not satisfy the Allen-Cahn equation (1) exactly, we add the extra source term g(x, y, t) = 0.05(2 2 -1.1)e -0.1t sin(x) sin(y) + (0.05e -0.1t sin(x) sin(y)) 3 such that u(x, y, t) in ( 109) is the exact solution to the revised Allen-Cahn equation
with the model parameter = 0.01. P 2 finite element space is again used.
Table 13 L 2 -norm of error and rate for 2D modified DLN scheme in time (N = 500)
4.98e-5 -2.57e-5 -1.48e-3 -40 6.08e-6 3.03 3.18e-6 3.02 3.78e-4 1.97 60 1.51e-6 3.43 8.03e-7 3.39 1.59e-4 2.14 80 5.74e-7 3.36 3.12e-7 3.29 8.81e-5 2.04 100 2.65e-7 3.46 1.46e-7 3.39 5.48e-5 2.13 θ = 2 √ 5 20 4.97e-5 -2.57e-5 -1.48e-3 -40 6.07e-6 3.03 3.17e-6 3.02 3.78e-4 1.97 60 1.51e-6 3.43 8.02e-7 3.39 1.59e-4 2.14 80 5.73e-7 3.37 3.11e-7 3.29 8.81e-5 2.04 100 2.63e-7 3.48 1.45e-7 3.41 5.47e-5 2.13 θ = 1 20 4.96e-5 -2.56e-5 -1.48e-3 -40 6.06e-6 3.03 3.17e-6 3.02 3.78e-4 1.97 60 1.51e-6 3.43 8.01e-7 3.39 1.59e-4 2.14 80 5.72e-7 3.37 3.11e-7 3.29 8.81e-5 2.04 100 2.63e-7 3.49 1.45e-7 3.42 5.47e-5 2.13
To confirm the time convergence rate, we first set the final time as T = 4 and consider a refined triangular mesh with N = 500 nodes on each side of the boundaries. Tables 13 and 15 show that both modified DLN and DLN-SAV algorithms achieved second-order convergence in time for the 2D test problem in (109). Then we set the final time T = 1 and the constant time step k = 0.01 to check the spatial convergence. We observe from Tables 14 and 16 that both algorithms have third-order accuracy in L 2 -norm and second-order accuracy in H 1 -norm.
Table 15 L 2 -norm of error and rate for 2D DLN-SAV scheme in time (N = 500)
4.97e-5 -2.57e-5 -1.48e-3 -40 6.08e-6 3.03 3.18e-6 3.01 3.78e-4 1.97 60 1.52e-6 3.43 8.07e-7 3.38 1.59e-4 2.14 80 5.80e-7 3.34 3.15e-7 3.27 8.81e-5 2.04 100 2.71e-7 3.41 1.49e-7 3.35 5.48e-5 2.13 θ = 2 √ 5 20 4.97e-5 -2.56e-5 -1.48e-3 -40 6.08e-6 3.03 3.18e-6 3.01 3.78e-4 1.97 60 1.52e-6 3.42 8.06e-7 3.38 1.59e-4 2.14 80 5.84e-7 3.32 3.16e-7 3.25 8.81e-5 2.04 100 2.78e-7 3.32 1.53e-7 3.27 5.48e-5 2.13 θ = 1 20 4.96e-5 -2.56e-5 -1.48e-3 -40 6.07e-6 3.03 3.17e-6 3.01 3.78e-4 1.97 60 1.52e-6 3.42 8.06e-7 3.38 1.59e-4 2.14 80 5.86e-7 3.31 3.17e-7 3.24 8.81e-5 2.04 100 2.83e-7 3.26 1.55e-7 3.22 5.47e-5 2.13 6.2.2 Variable Time Step Similar to Section 6.1.2, we test the accuracy of the variable time-stepping modified DLN algorithm (19) and DLN-SAV algorithm (53) over time for the 2D case, utilizing the randomized time step k n = k + k • rand. The final time is T = 4.0 and the number of nodes N = 500 on each side of the boundaries. The accuracy rate calculated using k max again demonstrates second-order temporal convergence rate numerically, as seen in Tables 17 and 18.
Table 17 L 2 -norm of error and rate for 2D modified DLN scheme in time (k n
3 0.4 0.7818 3.37e-3 -3.27e-3 -4.63e-3 -0.3 0.5851 2.17e-3 1.52 1.87e-3 1.94 2.64e-3 1.94 0.2 0.3683 1.10e-3 1.47 9.02e-4 1.57 1.28e-3 1.57 0.1 0.1987 3.10e-4 2.05 2.33e-4 2.20 3.29e-4 2.20 θ = 2 √ 5 0.4 0.7818 4.25e-3 -4.12e-3 -5.83e-3 -0.3 0.5851 2.87e-3 1.35 2.47e-3 1.76 3.50e-3 1.76 0.2 0.3683 1.52e-3 1.38 1.25e-3 1.48 1.76e-3 1.48 0.1 0.1987 4.48e-4 1.98 3.36e-4 2.13 4.75e-4 2.13 θ = 1 0 .4 0.7818 4.65e-3 -4.51e-3 -6.37e-3 -0.3 0.5851 3.21e-3 1.28 2.76e-3 1.69 3.90e-3 1.69 0.2 0.3683 1.72e-3 1.35 1.41e-3 1.45 2.00e-3 1.45 0.1 0.1987 5.16e-4 1.95 3.87e-4 2.10 5.47e-4 2.10
Next, we compare Algorithm 1 and the corresponding constant time-stepping algorithms to show the accuracy and efficiency of time adaptivity. We set final time T = 1, model parameter = 0.01, the number of nodes N = 50 on each side of the boundaries for both adaptive and constant algorithms. For adaptive algorithms, we set the maximum time step k max = 0.1, the minimum time step k min = 1.e -5, two initial time step k 0 = k 1 = 1.e -3, tolerance Tol = 1.e -8 and safety factor κ = 0.8. We use time step k = 1.e -3 for constant time-stepping algorithms. From Table 19, we can see that adaptive modified DLN algorithms, especially with θ = 2/3 or 1, work more efficiently than the corresponding constant time-stepping algorithms, since the adaptive algorithms take much fewer number of time steps and obtain errors comparable to those of the constant step algorithms. For
Table 19 Errors and number of time steps of the adaptive and constant time-stepping 2D modified DLN schemes Adaptive modified DLN schemes θ
-6 1.79e-6 2.37e-4 280 73 θ = 2 √ 5 4.86e-6 2.96e-6 2.38e-4 93 18 θ = 1 2.96e-6 1.60e-6 2.37e-4 48 18 Constant time-stepping modified DLN schemes θ
6 1.62e-6 2.37e-4 45 9 θ = 2 √ 5 3.08e-6 1.65e-6 2.37e-4 44 9 θ = 1 3.11e-6 1.66e-6 2.37e-4 44 15 Constant time-stepping DLN-SAV schemes θ u -u h ℓ ∞ (L 2 ) u -u h ℓ 2 (L 2 ) u -u h ℓ 2 (H 1 ) # Steps # Rejections θ = 2 3 2.97e-6 1.56e-6 2.37e-4 1000 θ = 2 √ 5 2.97e-6 1.56e-6 2.37e-4 1000 θ = 1 2.97e-6 1.56e-6 2.37e-4 1000 -
DLN-SAV algorithms, we observe from Table 20 that adaptive algorithms have the same error magnitude as constant time-stepping algorithms. However, adaptive algorithms with all three θ values finish the simulation only in 45 time steps while the constant step algorithms take 1000 time steps.
In this example, we simulate 2D Allen-Cahn equation with random initial condition [21]: where the function rand(x, y) generates value uniformly distributed on [0, 1] at point (x, y).
We set the model parameter to be 0.1 and number of nodes on each side of the boundaries to be 100 to generate a triangular mesh. We test the adaptive Algorithm 1 with periodic boundary condition. Fig. 2 2D Adaptive DLN-SAV algorithm with θ = 1 converges to steady state with 5896 time steps For both adaptive modified DLN and DLN-SAV algorithms, we set two initial time step k 0 = k 1 = 0.01, the maximum time step k max = 0.1, the minimum step k min = 1.e -5, the tolerance for LTE Tol = 1.e -6, and safety factor κ = 0.8. We use the fixed point iteration with tolerance 1.e -8 for non-linear solver in the adaptive modified DLN algorithm. The results of three θ values (2/3, 2/ √ 5, 1) are very close, thus we only present the numerical results for the case of θ = 1. From Fig. 1 and Fig. 2, we observe that both modified DLN and DLN-SAV algorithms converge to the steady state at the time T = 320 with around 5800 time steps. We also test the corresponding constant time-stepping DLN algorithms with constant step k = 0.01 and observe that the constant time-stepping algorithms converge to the steady state with more than 30000 time steps. From this test problem, we demonstrate that the adaptive DLN algorithms are stable and possess superior time efficiency when compared to the corresponding constant time-stepping algorithms.
Then Eq. (11) follows by the fact t n,β ∈ [t n-1 , t n+1 ] and using Hölder's inequality: u n+1,θ + u n,θ 2 u(t n,β ) 2 ≤ 1 + θ 4 2 Ω t n+1 t n-1 |u t (x, t)|dt 2 dx ≤ C(θ ) Ω t n+1 t n-1 1dt t n+1 t n-1 |u t (x, t)| 2 dt dx.
For the case of θ = 1, the corresponding conclusions for the midpoint rule are easy to verify and the details are omitted here.
Next, we consider the case of uniform time grids with constant time step k and aim to prove (14). Applying Taylor's theorem with integral remainder u(x, t n+1 ) = u(x, t n ) + u t (x, t n )k + t n+1 t n u tt (x, t)(t n+1t)dt, u(x, t n-1 ) = u(x, t n )u t (x, t n )k + t n-1 t n u tt (x, t)(t n-1t)dt, u(x, t n,β ) = u(x, t n ) + u t (x, t n )(t n,βt n ) + t n,β t n u tt (x, t)(t n,βt)dt = u(x, t n ) + u t (x, t n ) β (n) 2 -β (n) 0 k + t n,β t n u tt (x, t)(t n,βt)dt, (110) with (under uniform time grids) β (n) 2 = 1 4 (2 + θ -θ 2 ), β (n) 0 = 1 4 (2 -θ -θ 2 ), gives u n+1,θ + u n,θ 2 u(x, t n,β ) = 1 + θ 4 u(x, t n+1 ) + 1 2 u(x, t n ) + 1 -θ 4 u(x, t n-1 )u(x, t n,β ) = 1 + θ 4 t n+1 t n u tt (x, t)(t n+1t)dt + 1 -θ 4 t n-1 t n u tt (x, t)(t n-1t)dt t n,β t n u tt (x, t)(t n,βt)dt. Therefore, Eq. (14) follows from u n+1,θ + u n,θ 2 u(t n,β ) 2 ≤ C(θ ) Ω t n+1 t n |u tt (x, t)| 2 dt t n+1 t n (t n+1 -t) 2 dt+ t n t n-1 |u tt (x, t)| 2 dt t n t n-1 (t -t n-1 ) 2 dt + t n,β t n |u tt (x, t)| 2 dt t n,β t n (t n,β -t) 2 dt dx ≤ C(θ )k 3 Ω t n+1 t n-1 |u tt (x, t)| 2 dtdx.
As in the proof of Lemma 2, we only prove the case θ ∈ [0, 1) below. The case for θ = 1 can be shown easily following the same procedure. Also, it suffices to consider the case r = 0 below. Combining Taylor's theorem with integral remainder (110) (but with variable time step sizes k n and k n-1 ) with the fact that β (n)
1 leads to u n,β (x)u(x, n,β ) = β (n) 2 t n+1 t n u tt (x, t)(t n+1t)dt + β (n) 0 t n-1 t n u tt (x, t)(t n-1t)dt t n,β t n u tt (x, t)(t n,βt)dt. (111) Utilizing the Hölder's inequality, and bounding β (n) ℓ by C(θ ), we have u n,β (•)u(•, t n,β ) 2 ≤ C(θ ) Ω t n+1 t n |u tt (x, t)|(t n+1 -t)dt + t n t n-1 |u tt (x, t)|(t -t n-1 )dt + t n,β t n |u tt (x, t)||t n,β -t|dt 2 dx ≤ C(θ ) Ω t n+1 t n |u tt (x, t)| 2 dt t n+1 t n (t n+1 -t) 2 dt 1 2 + t n t n-1 |u tt (x, t)| 2 dt t n t n-1 (t -t n-1 ) 2 dt 1 2 + t n,β t n |u tt (x, t)| 2 dt t n,β t n (t n,β -t) 2 dt 1 2 2 dx ≤ C(θ )(k n + k n-1 ) 3 Ω t n+1 t n-1 |u tt (x, t)| 2 dtdx, which implies (15) with r = 0. Next, we consider Eq. (16). We again consider Taylor's theorem with integral remainder u(x, t n+1 ) = u(x, t n )+u t (x, t n )k n +u tt (x, t n ) k 2 n 2 + t n+1 t n u ttt (x, t) (t n+1 -t) 2 2 dt, u(x, t n-1 ) = u(x, t n )-u t (x, t n )k n-1 +u tt (x, t n ) k 2 n-1 2 + t n-1 t n u ttt (x, t) (t n-1 -t) 2 2 dt, u t (x, t n,β ) = u t (x, t n )+u tt (x, t n )(β (n) 2 k n -β (n) 0 k n-1 )+ t n,β t n u ttt (x, t)(t n,β -t)dt. and combine them with the fact that α 2 + α 1 + α 0 = 0 to obtain u n,α (x) k n u t (t n,β ) = α 2 k n 2 +α 0 k n-1 2 2 k n -β (n) 2 k n -β (n) 0 k n-1 u tt (x, t n )-t n,β t n u ttt (x, t)(t n,β -t)dt + α 2 k n t n+1 t n u ttt (x, t)(t n+1 -t) 2 dt + α 0 k n t n-1 t n u ttt (x, t)(t n-1 -t) 2 dt. It's easy to check α 2 k n 2 + α 0 k n-1 2 2 k n -β (n) 2 k n -β (n) 0 k n-1 = 0. (112) Therefore, u n,α (•) k n u t (x, t n,β ) 2 ≤ C(θ ) k 2 n Ω t n+1 t n |u ttt (x, t)|(t n+1t) 2 dt + t n t n-1 |u ttt (x, t)|(t n-1t) 2 dt + k n t n,β t n |u ttt (x, t)||t n,β -t|dt 2 dx ≤ C(θ ) k 2 n Ω t n+1 t n |u ttt (x, t)| 2 dt t n+1 t n (t n+1 -t) 4 dt 1 2 + t n t n-1 |u ttt (x, t)| 2 dt t n t n-1 (t -t n-1 ) 4 dt 1 2 + k n t n,β t n |u ttt (x, t)| 2 dt t n,β t n (t n,β -t) 2 dt 1 2 2 dx ≤ C(θ ) (k n + k n-1 ) 2 + k 2 n k 2 n (k n + k n-1 ) 3 Ω t n+1 t n-1 |u ttt (x, t)| 2 dtdx ≤ C(θ )(k n + k n-1 ) 3 Ω t n+1 t n-1 |u ttt (x, t)| 2 dtdx,
where the last inequality follows from k n = 1+θ 2 k n + 1-θ 2 k n-1 . The case with r ≥ 0 can be proven in the same way.
+
f (u n+1,θ , u h n,θ )f (u n+1,θ , u n,θ ) 2 + 3 f (u n+1,θ , u n,θ )f (u(t n,β )) 2 . (43)