Sharp L2 Norm Convergence of Variable-Step BDF2 Implicit Scheme for the Extended Fisher – Kolmogorov Equation

,


Introduction
The Fisher-Kolmogorov-type models have been widely used to study physical, material, and biological systems [1][2][3], such as the mesoscopic model of a phase transition in a binary system [1], the propagation of domain walls in liquid crystals [2], and the growth progression of certain types of primary brain tumors [3].Consider the following free energy (Lyapunov) functional where the spatial domain is Ω = ð0, LÞ 2 ⊂ ℝ 2 , the parameter γ > 0 is a positive constant, and FðΦÞ = 1/4ð1 − Φ 2 Þ 2 is the bistable type admitting two local minima.Mathematically, the governing system of the extended Fisher-Kolmogorov (EFK) model could be derived via an L 2 gradient flow associated with the energy functional E½Φ, that is, subjected to periodic boundary conditions with an initial data Φðx, 0Þ ≔ Φ 0 ðxÞ, and the nonlinear bulk f ðvÞ = ðv3Þ 3 − v.One intrinsic property of the EFK model is the energy dissipation law.
in which ð f , gÞ ≔ Ð Ω f gdx, and the associated ðΩÞ.The existing works [4][5][6] with respect to the numerical methods of the EFK model (2) are mainly Crank-Nicolson (CN) type schemes with a uniform time-step.However, it is well known that the evolutions of the initial solutions of the gradient flow models exhibit numerous tortuous diffuse interfaces with a complex topology that coalesces to form gradually less complex configurations and then simplifies toward an eventual steady state, for instance in the spinodal decomposition problem [7][8][9], the epitaxy thin film growth [10,11], and the growth of a polycrystal in a supercooled liquid [12,13].The temporal evolutions of the gradient flow models involve multiple time scales.Naturally, small time steps should be used to capture the fast time dynamics, but the computation is very costly, especially for a later time when the phase structure evolves more slowly and smoothly; while the quick transition may be overlooked if large time steps are used.Consequently, it is essential to adopt some adaptive timestepping strategy to efficiently revolve widely varying time scales.More precisely, small time steps are employed when the solution has a quick transition, and large times are allowed when the solution changes slowly.Several adaptive time step algorithms [7][8][9][10]12] have been suggested for simulating phase field models.Ample numerical results confirm that adaptive time-stepping techniques can greatly reduce the computation time while ensuring that sufficient time accuracy is achieved.
The variable-step BDF2 method is more favorable than the CN-type methods for solving stiff differential problems due to its strong stability property [14][15][16][17][18].However, the stability and convergence theory of adaptive BDF2 timestepping for parabolic equations is rather difficult, especially for general step-size variations, and remains incomplete so far.For the analysis of the variable-step BDF2 method applied to ordinary differential equations, Grigorieff et al. [19,20] proved that it is zero-stable only if the adjacent time-step ratios Y k < 1 + ffiffi ffi 2 p .Becker ([21], Lemma 4.1) proved the variable-step BDF2 method of a linear parabolic problem is stable if Y k < 1:868.Emmrich ( [15], Theorem 3) then improved Becker's bound-up to Y k < 1:910 for the semilinear problems.As for the nonlinear parabolic problems, Emmrich ([22], Theorem 10) showed that the ratios restriction Y k < 1:702 should be imposed to ensure the stability of the variable-step BDF2 method.Moreover, the stability and error estimates reported in [15,21,22] often involve an undesirable prefactor exp(CΓ n ), where Γ n may be unbounded as the time step sizes vanish or require the time meshes appropriately close to the uniform one.Recently, Chen et al. [23] applied a new discrete Grönwall inequality to remove the undesirable prefactor exp(CΓ n ) under the step-ratio condition Y k < 1:534; while the somewhat severe restriction might limit the advantage of using variable time steps.
Very recently, a technique of discrete orthogonal convolution (DOC) kernels was developed in [24] to establish the L 2 norm stability and convergence theories of adaptive BDF2 method for the linear reaction-diffusion problem under an updated time-step ratios restriction S1.The adjacent time-step ratios 0 In what follows, we extend the DOC technique to establish the convergence analysis of full implicit adaptive BDF2 methods for the phase field models, including the molecular beam epitaxial growth model without slope selection [11] and the phase field crystal model [13].However, some nonphysical step-size constraints (which may be unnecessary in numerical simulations) have to be imposed to ensure the unique solvability and energy dissipation property of the BDF2 method.
In this paper, we continue to develop the DOC technique and analyze the variable-step BDF2 time-stepping method for solving the EFK equation.We point out the significant difference between the current work and the earlier ones [11,13]: the nonphysical time step restrictions are removed by introducing a stabilizing term [25,26] as we show below.Consider the nonuniform time levels For v n = vðt n Þ, the BDF2 formula with unequal time steps reads where the discrete convolution kernels b ðnÞ n−k are defined by b ð1Þ 0 ≔ 1/τ 1 for n = 1, and when the time- Without losing the generality, we use the BDF1 formula to compute the first-level solution by putting r 1 ≡ 0 in (4).
To present the fully discrete scheme, we describe briefly the spatial approximation.For the domain Ω = ð0, LÞ 2 , let the grid length h x = h y = h ≔ L/M with an integer M. Let the discrete domains Ω h ≔ fx h = ðih, jhÞj1 ≤ i, j ≤ Mg and The operators Δ y w ij and δ 2 y w ij can be defined similarly so that the discrete gradient vector Ω h .We focus on the adaptability of the numerical method with respect to the variations of time steps by considering a stabilized implicit BDF2 approach for solving the EFK Equation (2).That is, to find the numerical solution subjected to the initial value φ 0 h = Φ 0 ðx h Þ and periodic boundary conditions, where S 1 ≥ 0 is a stabilized parameter.Note that, the stabilized technique is originally proposed to guarantee the unconditional energy stability of some linearized numerical methods [25,26].Currently, we cannot establish the stability and convergence theory for the 2 Journal of Function Spaces second-order semi-implicit stabilized BDF2 method with unequal time steps.The stabilized strategy is introduced in the variable-step BDF2 method to ensure the unconditionally unique solvability, see Theorem 1, and the energy dissipation law in Theorem 3. Then Theorem 4 establishes the boundedness of the solution in the maximum norm.On the other hand, extra errors could be incurred by the stabilized term S 1 τ n ∇ τ ∅ n h and the numerical results reported in Example 12 indicate the stabilization parameter should be selected carefully to maintain numerical accuracy.In Section 3, with the help of a new class of DOC kernels, an optimal L 2 norm error estimate of the numerical scheme ( 6) is established for the nonlinear fourth-order parabolic problem.It is worth mentioning that, our stability (Theorem 4) and error estimates (Theorems 8-10) are almost independent of the step ratios r n , so the variable-step BDF2 scheme is surprisingly robust with respect to the variations of time steps.To the best of our knowledge, this is the first time such stability and convergence theory is built for the variable-step BDF2 scheme to the extended Fisher-Kolmogorov equation.Numerical experiments and comparisons are presented in Section 4 to show the effectiveness of our method.An adaptive time-stepping strategy on the basis of solution accuracy is then used to improve computational efficiency.
Throughout this paper, any subscripted C, such as C u and C Φ , denotes a generic positive constant, not necessarily the same at different occurrences; while, any subscripted c, such as c Ω , c 0 , c 1 and c 2 , denotes a fixed constant.Always, the appeared constants dependent on the given data and the solution, but independent of the spatial length h and time steps τ n .

Solvability and Energy Dissipation Law
For any grid functions v, w ∈ V h , define the discrete L 2 inner product <v, w > ≔h 2 ∑ x h ∈Ω h v h w h and the associated L 2 norm kvk ≔ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi <v, w > p .The discrete seminorms k∇ h vk and kΔ h vk can be defined, respectively, by The discrete L ∞ norm kvk ∞ ≔ max x h ∈Ω h and we have the Sobolev embedding inequality where c Ω > 0 is a constant dependent on the spatial domain Ω h .For any grid functions v, w ∈ V h , the discrete Green's formula with periodic boundary conditions yields.
2.1.Unique Solvability.The following result shows that the stabilized BDF2 scheme ( 6) is unconditionally unique and solvable thanks to the stabilized technique.
Proof.We prove the solvability of the BDF2 implicit scheme (6) via a discrete energy functional G on the space v h , For any time-level index n ≥ 1, the condition implies Then, G is strictly a convex function, for any λ ∈ ℝ and any Thus, the functional G has a unique minimizer, denoted by ∅ n h , if and only if it solves the equation This equation holds for any ψ h ∈ vh h if and only if the unique minimizer Which is just our BDF2 time-stepping scheme (6).It completes the proof.

Energy Dissipation Law.
For the convenience of subsequent analysis, we denote Then, the following result, also see ( [24], Lemma 2), shows the BDF2 convolution kernels b ðnÞ n−k are positive definite if the adjacent time-step ratios γ k satisfy S 1 .
Lemma 2. For any real sequence fwkg n k=1 with n entries, it holds that If S 1 holds, the discrete convolution kernels b ðnÞ n−k are positive and definite in the sense that Now, we consider the energy dissipation law of BDF2 scheme (6) to simulate the continuous property (3).Let, E½∅ n be the discrete version of energy (Lyapunov) functional (1), for n ≥ 0: ð18Þ Define the following modified discrete energy ε½∅ k by ε½∅ 0 ≔ E½ε½∅ 0 and Theorem 3. If S 1 holds and the stabilized parameter S 1 fulfills The BDF2 time-stepping scheme (6) preserves an energy dissipation law.
Proof.Taking the inner product of ( 6) by ∇ τ ∅ n , one has where the discrete Green's formula with periodic boundary conditions has been used.With the help of the discrete Green's formula and 2aða − bÞ = a 2b 2 + ða − bÞ 2 , we have ð23Þ ð24Þ It is easy to get the following relationship Then, we can obtain ð26Þ By inserting ( 23)-( 26) into (22), we have Applying Lemma 2 with the restriction (20) of stabilized parameter S 1 , one has Thanks to the fact ð1/16S 1 τ n Þ + S 1 τ n ≥ 1/2, one combines the above inequality with (27) to obtain the energy dissipation plaw ε½∅ n ≤ ε½∅ n−1 for n ≥ 2. For the case n = 1, the Condition (20) with γ 1 ≡ 0 gives S 1 ≥ 1/4, so that Journal of Function Spaces Then, we derive from (27) that ε½∅ 1 ≤ E½∅ 0 = ε½∅ 0 and complete the proof.
In adaptive computations, we can choose a proper stepratio γ k+1 (especially when the current step-ratio γ k or timestep τk is large) from the sufficient condition S1 to make the function R(γ k , γ k+1 ) have a positive lower bound.It is easy to check the following estimates In this case, the stabilized parameter S 1 = 0:5 is sufficient for energy stability since the Condition (20) Similarly, the continuous energy dissipation law gives kΦðt n Þk ∞ ≤ kΦðt n Þk L ∞ ≤ c 1 .Note that the fixed constants c 0 and c 1 are independent of the spatial lengths, the time steps τ n , the current time t n , and the step ratios γ n .
Proof.Since a 4 ≥ 2ð1 + γÞa 2 − ð1 + γÞ 2 due to the simple fact ða 2 − ð1 + γÞÞ 2 ≥ 0, we apply Theorem 3 to obtain So the Sobolev embedding inequality leads to It implies the maximum norm estimate and completes the proof.

L 2 Norm Convergence Analysis
To facilitate the subsequent numerical analysis, we introduce a modified BDF2 formula where the modified BDF2 kernels d ðnÞ n−k are defined by Thus the stabilized BDF2 implicit scheme (6) reads .
This type of DOC kernel was first introduced in [24] with respect to the original BDF2 convolution kernels b ðnÞ n−k .It is easy to check that the following discrete orthogonal identity This identity directly leads to the following relationship where the order of summation has been exchanged in the second equality.This equality (39) plays an essential role in the subsequent L 2 norm analysis.
Applying Lemma 2 and the definition (35) of modified BDF2 kernels d ðnÞ n−k , it is not difficult to check the following result. where Proof.Following the proof of ( [11], Lemma 6), one can prove the positive definiteness of v ðnÞ n−j by using Lemma 5, and the orthogonal relationship between the two classes of discrete convolution kernels v ðnÞ n−j and d ðnÞ n−j .To prove (II), define The definition (37) yields Since β i > 0, we obtain and confirm (II) immediately.Lastly, taking v n = t n in equality (3.6), we have Thus the result of (III) follows from the definition (34) since D2 t j = 1 + S 1 τ 2 j ≥ 1.

L 2
Norm Error Estimate.Now, we are to establish the L 2 norm error estimate of our BDF2 scheme (6) via the equivalent form (36). Let the local consistency error at the time = t j , Consider a convolutional consistency error Ξ k defined by By adding the error term S 1 τ j ∇ τ Φðt j Þ, stemmed from the stabilized term S 1 τ j ∇ τ ∅ n h with a proper parameter S 1 = Oð1Þ, follow to the proof of ( [13], Lemma 3.4), we arrive at the following result.

Lemma 7.
If S 1 holds, the convolutional consistency error Ξ k Φ in (46) satisfies Such that Theorem 8. Assume that the EFK problem (2) has a solution x,t ðΩ × ð0, TÞ.If S1 holds, the stabilized parameter S 1 fulfills (20), and the maximum step size τ is sufficiently small such that τ ≤ 1/ð4c 2 Þ, the stabilized BDF2 implicit timestepping scheme (6) with unequal time steps is convergent in the L 2 norm, Here, c 2 ≔ 1 + c 2 0 + c 0 c 1 + c 2 1 is dependent on the domain Ω and the initial values Φ 0 and ∅ 0 , but independent of the spatial length h, the time t n , the step sizes τ n , and the step ratios r n .

Proof. Let the error function e
We have the following error equation where ε j h and η j h denote the local consistency error in time and space, respectively.Under the regularity setting of the solution, Lemma 6 (III) gives Q k h are defined by ( 46) and (51), respectively.Making the inner product of the above equality with 2e k and summing up the superscript from k = 1 to n, we apply the discrete Green's formula to 6 Journal of Function Spaces derive that for 1 ≤ n ≤ N. Thanks to the maximum norm estimates in Theorem 3, one gets so that the Cauchy-Schwarz inequality leads to Then, it follows from (52) that Choose some integer n 0 ð0 ≤ n 0 ≤ nÞ such that ke n0 k = max 0≤k≤n ke k k.Let n = n 0 in the above inequality (55), such that Under the setting τ ≤ 1/ð4c 2 Þ, we have Obviously, applying the discrete Gronwall inequality ( [24], Lemma 5), one can obtain The proof is completed from Lemma 7 and the error estimate (51).
As noted early in ([13], Remark 2), Theorem 8 confirms the (at least, first-order) convergence of numerical solution under the step-ratio condition The second-order convergence in time can be recovered if S 2 .The time-step ratios rk are contained in S 1 , but almost all of them less than . Actually, very large steps always result in a loss of numerical accuracy, although the condition S 1 permits us to use a series of increasing time steps with amplification factors up to 3.561.The condition S 2 is reasonable in practice because large amplification factors of time-step size rarely appeared continuously in long-time simulations.As shown in ( [24], Lemma 7), there is a step-ratio-dependent constant C r such that Thus, we have the following corollary.
Corollary 9. Assume that the EFK problem (2) has a unique smooth solution.If S 2 holds, the stabilized parameter S 1 fulfills (20) and the maximum step size τ ≤ 1/ð4c 2 Þ, the stabilized BDF2 scheme ( 6) is second-order convergent in the L 2 norm, ð59Þ  At the same time, under the mild step-ratio condition S 1 , the second-order temporal accuracy can also be achieved theoretically by employing some high-order starting scheme (with very smooth initial data) other than the first-order BDF1 scheme.By following the discussions in ( [24], Remark 4 and Theorem 3.3), it is not difficult to verify the following result.
Theorem 10.Assume that the EFK problem (2) has a unique smooth solution.If S 1 holds with the maximum step size τ ≤ 1/ð4c 2 Þ, the stabilized BDF2 implicit time-stepping scheme (6) without the first-level BDF1 scheme is convergent in the L 2 norm,

Numerical Example
We run the variable-step BDF2 scheme (6) for solving the EFK equation ( 2) numerically.The BDF1 formula is used to obtain the first-level solution, and a simple iteration is employed to solve the nonlinear algebra equations at each time level.Always, we use the solution at the previous level as the initial guess for each iteration, and the numerical tolerance is taken as 10 −12 .The time accuracy of the stabilized BDF2 scheme ( 6) is examined via the random time mesh.Let τ k ≔ Tσ k /S for 1 ≤ k ≤ N, where σ k ∈ ð0, 1Þ is the uniformly distributed random number and S = ∑ N k=1 σ k .The discrete L 2 norm error eðNÞ ≔ kΦðTÞ − ∅ N k is recorded in each run and the experimental order of convergence is computed by Order ≈ log ðeðNÞ/eð2 NÞÞ/log ðτ ðNÞ/τ ð2 NÞÞ, in which τ (N) denotes the maximum time-step size for total N subintervals.
In our computations, we take the stabilization parameter S 1 = 0:5, and the spatial grid point M = 1024 in each direction such that the temporal error dominates the spatial error in each run.The problem is solved until time T = 1.The numerical results are reported in Table 1, in which we also record the maximum time-step size, the maximum step ratio and the number (denote by N 1 in Table  the step ratio γ k ≥ ð3 + ffiffiffiffiffi 17 p Þ/2. From these numerical results, one can observe that the stabilized BDF2 scheme ( 6) is robustly stable and has second-order accuracy on nonuniform time meshes.
To exploit the numerical performance of the stabilized term S 1 τ n ∇ t ∅ n h in the BDF2 scheme (6), we also run the scheme (6) by varying the value S 1 .Note that the parameter values used in our examinations are the same as before but the stabilization parameter S 1 .The error curves obtained using different stabilization parameters S 1 =0, 0.5, 1 are plotted in Figure 1.As can be seen, the extra errors are incurred by the stabilized term and the bigger value S 1 gives the less accurate solution.Thus the stabilization parameter S 1 should be selected as small as possible to maintain the numerical accuracy.

Adaptive Time-Stepping Strategy.
The temporal evolution of the solution of the EFK Equation (2) may involve multiple time scales, such as the numerical example discussed in Example 12, the initial solutions change dramatically while later ones change slowly.In the following, we shall take the following variant time adaptive strategy of [7,8] to select the time steps.
Here, the BDF1 and BDF2 schemes used in Algorithm 1 refer to the backward Euler scheme and BDF2 scheme (6), respectively.The adaptive time step τ ada is given by τ ada ðe, τ cur Þ = ffiffiffiffiffiffiffiffiffi tol/e p ρτ cur , where ρ is a default safety coefficient, tol is a reference tolerance, e is the relative error at each time level, and τ cur is the current time step.Also, τmax and τmin are the predetermined maximum and minimum time steps.In the following  9 Journal of Function Spaces computations, if not explicitly specified, we take the safety coefficient ρ = 0:9, the reference tolerance tol = 10 −3 , the maximum time step τ max = 0:1, and the minimum time step τ min = 10 −4 .
We first examine the numerical performance of the uniform and adaptive time strategies.We here begin with the examination using a constant time step τ = 10 −3 until time T = 5 (i.e., N = 5000).Also, we compute the numerical solution by employing the adaptive time strategy reported in Algorithm 1.The numerical results involving the discrete energies and time steps are plotted in Figure 2. From these diagrams, one can observe that the energy curve obtained using an adaptive timing strategy practically coincides with that obtained using a small uniform step τ = 10 −3 .However, the total number of adaptive time steps is 509, which is far less than the constant time steps of 5000.Thus, the adaptive time strategy increases the computational efficiency but does not affect the numerical results.
In what follows, we employ the stabilized BDF2 scheme (6) with Algorithm 1 to solve the EFK problem until time T = 15.The numerical solutions are depicted in Figure 3.One can find that the initial solution evolves on a fast time scale while later ones evolve on a very slow time scale.In addition, Figure 4 shows the discrete energy and adaptive time steps.One can see that the evolution of energy undergoes large variations initially, while it changes very slowly in the remaining time.Correspondingly, some small time steps are selected when the energy decays quickly, and large steps are adopted when the energy dissipates slowly.

3. 1 .
A New Class of DOC Kernels.We define a new class of DOC kernels fv ðnÞ n−k g n k=1

Lemma 5 .Lemma 6 .
If S 1 holds and the parameter S 1 ≥ 0, the modified BDF2 kernels d ðnÞ k=1 in (35) are positive definite in the sense that, for any nonzero sequence fw k g n k=1 , If the discrete convolution kernels d ðnÞ n−k defined in (35) are positive definite, then the DOC kernels v ðnÞ n−j defined in (37) satisfy the following properties: (I) The discrete kernels v ðnÞ n−j are positive definite, and then the symmetric matrix where the elements of Θ are defined by Θ ij for i > j and Θ ij ≔ 1/2v ≤ j ðiÞ i−j for i ≤ j.The discrete kernels v ðnÞ n−j are positive and v ðnÞ n−j = 1/d ðjÞ 0

1 Figure 1 :
Figure 1: The semilog plot of numerical error e(N) of stabilized BDF2 scheme (6) with different stabilization parameters S 1 .

Figure 2 :
Figure 2: Evolutions of energy (a) and comparison of time steps (b) of the stabilized BDF2 scheme (6) for the EFK equation using different time strategies until time T = 5.

Figure 4 :
Figure 4: Evolutions of energy (a) and adaptive time steps (b) of the BDF2 scheme (6) for the EFK equation using adaptive time strategy.