Abstract

A variable-step BDF2 time-stepping method is investigated for simulating the extended Fisher-Kolmogorov equation. The time-stepping scheme is shown to preserve a discrete energy dissipation law if the adjacent time-step ratios . With the aid of discrete orthogonal convolution kernels, concise norm error estimates are proved, for the first time, under the mild step ratios constraint . Our error estimates are almost independent of the step ratios so that the proposed numerical scheme is robust with respect to the variations of time steps. An adaptive time-stepping strategy based on solution accuracy is then applied to update the computational efficiency. Numerical examples are included to illustrate our theoretical results.

1. Introduction

The Fisher-Kolmogorov-type models have been widely used to study physical, material, and biological systems [13], 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 , the parameter is a positive constant, and is the bistable type admitting two local minima. Mathematically, the governing system of the extended Fisher-Kolmogorov (EFK) model could be derived via an gradient flow associated with the energy functional , that is, subjected to periodic boundary conditions with an initial data , and the nonlinear bulk . One intrinsic property of the EFK model is the energy dissipation law. in which , and the associated norm for all .

The existing works [46] 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 [79], 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 time-stepping 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 [710, 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 [1418]. However, the stability and convergence theory of adaptive BDF2 time-stepping 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 . Becker ([21], Lemma 4.1) proved the variable-step BDF2 method of a linear parabolic problem is stable if . Emmrich ([15], Theorem 3) then improved Becker’s bound-up to for the semilinear problems. As for the nonlinear parabolic problems, Emmrich ([22], Theorem 10) showed that the ratios restriction 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(), where 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() under the step-ratio condition ; 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 norm stability and convergence theories of adaptive BDF2 method for the linear reaction-diffusion problem under an updated time-step ratios restriction . The adjacent time-step ratios for . 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 with the time-step sizes . Denote the maximum step size and the local time-step ratio for . Given a grid function , put and for . For , the BDF2 formula with unequal time steps reads where the discrete convolution kernels are defined by for and when the time-level index ,

Without losing the generality, we use the BDF1 formula to compute the first-level solution by putting in (4).

To present the fully discrete scheme, we describe briefly the spatial approximation. For the domain , let the grid length with an integer M. Let the discrete domains and . For the function , let and . The operators and can be defined similarly so that the discrete gradient vector and the discrete Laplacian . Denote the space of L-periodic grid functions .

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 such that subjected to the initial value and periodic boundary conditions, where 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 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 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 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 , 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 and , denotes a generic positive constant, not necessarily the same at different occurrences; while, any subscripted , such as , , and , 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 .

2. Solvability and Energy Dissipation Law

For any grid functions , define the discrete inner product and the associated norm . The discrete seminorms and can be defined, respectively, by

The discrete norm and we have the Sobolev embedding inequality where is a constant dependent on the spatial domain . For any grid functions , 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.

Theorem 1. If the stabilized parameter , scheme (6) is uniquely solvable.

Proof. We prove the solvability of the BDF2 implicit scheme (6) via a discrete energy functional G on the space , For any time-level index , the condition implies . Then, G is strictly a convex function, for any and any , Thus, the functional G has a unique minimizer, denoted by , if and only if it solves the equation This equation holds for any if and only if the unique minimizer solves Which is just our BDF2 time-stepping scheme (6). It completes the proof.

2.2. Energy Dissipation Law

For the convenience of subsequent analysis, we denote

Note that, solves the algebraic equation . Also, one has so that is increasing in and decreasing in . is decreasing with respect to . Thus

Then, the following result, also see ([24], Lemma 2), shows the BDF2 convolution kernels are positive definite if the adjacent time-step ratios satisfy .

Lemma 2. For any real sequence with entries, it holds that

If holds, the discrete convolution kernels 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, be the discrete version of energy (Lyapunov) functional (1),

Define the following modified discrete energy by and

Theorem 3. If holds and the stabilized parameter fulfills The BDF2 time-stepping scheme (6) preserves an energy dissipation law.

Proof. Taking the inner product of (6) by , 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 , we have It is easy to get the following relationship Then, we can obtain By inserting (23)-(26) into (22), we have Applying Lemma 2 with the restriction (20) of stabilized parameter , one has Thanks to the fact , one combines the above inequality with (27) to obtain the energy dissipation plaw for . For the case , the Condition (20) with gives , so that Then, we derive from (27) that and complete the proof.

In adaptive computations, we can choose a proper step-ratio (especially when the current step-ratio or time-step τk is large) from the sufficient condition S1 to make the function R(, ) have a positive lower bound. It is easy to check the following estimates

In this case, the stabilized parameter is sufficient for energy stability since the Condition (20) only requires .

Theorem 4. If holds and the stabilized parameter fulfills (20), the numerical solution of the BDF2 time-stepping scheme (6) is bounded (stable) in the norm Similarly, the continuous energy dissipation law gives . Note that the fixed constants and are independent of the spatial lengths, the time steps , the current time , and the step ratios .

Proof. Since due to the simple fact , we apply Theorem 3 to obtain So the Sobolev embedding inequality leads to It implies the maximum norm estimate and completes the proof.

3. Norm Convergence Analysis

To facilitate the subsequent numerical analysis, we introduce a modified BDF2 formula where the modified BDF2 kernels are defined by

Thus the stabilized BDF2 implicit scheme (6) reads

3.1. A New Class of DOC Kernels

We define a new class of DOC kernels .

This type of DOC kernel was first introduced in [24] with respect to the original BDF2 convolution kernels . 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 norm analysis.

Applying Lemma 2 and the definition (35) of modified BDF2 kernels , it is not difficult to check the following result.

Lemma 5. If holds and the parameter , the modified BDF2 kernels in (35) are positive definite in the sense that, for any nonzero sequence ,

Lemma 6. If the discrete convolution kernels defined in (35) are positive definite, then the DOC kernels defined in (37) satisfy the following properties: (I)The discrete kernels are positive definite, and then the symmetric matrix is positive definite
where the elements of are defined by for and for .
The discrete kernels are positive and for , where (II) such that for

Proof. Following the proof of ([11], Lemma 6), one can prove the positive definiteness of by using Lemma 5, and the orthogonal relationship between the two classes of discrete convolution kernels and . To prove (II), define The definition (37) yields for . For , we use the definition (37) to find a recursive procedure Since , we obtain for , and confirm (II) immediately. Lastly, taking in equality (3.6), we have Thus the result of (III) follows from the definition (34) since .

3.2. Norm Error Estimate

Now, we are to establish the norm error estimate of our BDF2 scheme (6) via the equivalent form (36). Let the local consistency error at the time ,

Consider a convolutional consistency error defined by

By adding the error term , stemmed from the stabilized term with a proper parameter , follow to the proof of ([13], Lemma 3.4), we arrive at the following result.

Lemma 7. If holds, the convolutional consistency error in (46) satisfies Such that

Theorem 8. Assume that the EFK problem (2) has a solution . If S1 holds, the stabilized parameter fulfills (20), and the maximum step size is sufficiently small such that , the stabilized BDF2 implicit time-stepping scheme (6) with unequal time steps is convergent in the norm, Here, is dependent on the domain Ω and the initial values and , but independent of the spatial length h, the time , the step sizes , and the step ratios .

Proof. Let the error function with for . We have the following error equation where and denote the local consistency error in time and space, respectively. Under the regularity setting of the solution, Lemma 6 (III) gives for and , where and defined by (46) and (51), respectively. Making the inner product of the above equality with and summing up the superscript from to , we apply the discrete Green’s formula to derive that for . 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 such that . Let in the above inequality (55), such that Under the setting , 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 because . The second-order convergence in time can be recovered if . The time-step ratios rk are contained in , but almost all of them less than , or , where is an index set . Actually, very large steps always result in a loss of numerical accuracy, although the condition permits us to use a series of increasing time steps with amplification factors up to 3.561. The condition 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 such that . Thus, we have the following corollary.

Corollary 9. Assume that the EFK problem (2) has a unique smooth solution. If holds, the stabilized parameter fulfills (20) and the maximum step size , the stabilized BDF2 scheme (6) is second-order convergent in the norm,

At the same time, under the mild step-ratio condition , 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 holds with the maximum step size , the stabilized BDF2 implicit time-stepping scheme (6) without the first-level BDF1 scheme is convergent in the norm,

4. 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 .

4.1. Tests on Random Time Meshes

Example 11. Consider the exterior-forced EFK equation in the domain with such that it has a solution .

The time accuracy of the stabilized BDF2 scheme (6) is examined via the random time mesh. Let for , where is the uniformly distributed random number and . The discrete norm error is recorded in each run and the experimental order of convergence is computed by , in which τ () denotes the maximum time-step size for total subintervals.

In our computations, we take the stabilization parameter , and the spatial grid point in each direction such that the temporal error dominates the spatial error in each run. The problem is solved until time . 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 in Table 1) of time levels with the step ratio . 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 in the BDF2 scheme (6), we also run the scheme (6) by varying the value . Note that the parameter values used in our examinations are the same as before but the stabilization parameter . The error curves obtained using different stabilization parameters =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 gives the less accurate solution. Thus the stabilization parameter should be selected as small as possible to maintain the numerical accuracy.

4.2. 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 is given by , where is a default safety coefficient, is a reference tolerance, is the relative error at each time level, and is the current time step. Also, τmax and τmin are the predetermined maximum and minimum time steps. In the following computations, if not explicitly specified, we take the safety coefficient , the reference tolerance , the maximum time step , and the minimum time step .

Require: Set . Given and time step
1: Compute by using BDF1 scheme with time step .
2: Compute by using BDF2 scheme with time step .
3: Calculate .
4: if then
5:  Update time-step size .
6: else
7:  Recalculate with time-step size .
8:  Goto 1
9: end if
4.3. Numerical Application

Example 12. Consider the EFK Equation (2) with the following initial data in the domain with parameter . We use uniform mesh to discrete the spatial domain and take the stabilization parameter in our simulations.

We first examine the numerical performance of the uniform and adaptive time strategies. We here begin with the examination using a constant time step until time (i.e., ). 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 . 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 . 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.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

Authors’ Contributions

Yang Li provided the main idea of the article and wrote the first draft. Qiang Sun performed the validation and formal analysis. Naidan Feng reviewed and edited the manuscript as it was submitted. Jianjun Liu contributed to the original draft and editing of the methods. All authors read and approved the final draft.

Acknowledgments

This work is supported by the Natural Science Foundation of China (61877065), Liaoning Province Natural Science Foundation (No.2019-KF-23-08), and 111 Project (No. D18003).