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

Authors: Yang Li; Qihang Sun; Naidan Feng; Jianjun Liu
Journal: Journal of Function Spaces
 (2023)

---

## 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 ratiosrn≔Τn/Τn−1<3+17/2≈3.561. With the aid of discrete orthogonal convolution kernels, concise L2 norm error estimates are proved, for the first time, under the mild step ratios constraint 0<rn<3.561. Our error estimates are almost independent of the step ratios rn 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.


---

## Body

## 1. Introduction

The Fisher-Kolmogorov-type models have been widely used to study physical, material, and biological systems [1–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
(1)EΦ=∫ΩY2ΔΦ2+12∇Φ2+FΦdx,where the spatial domain is Ω=0,L2⊂ℝ2, the parameter γ>0 is a positive constant, and FΦ=1/41−Φ22 is the bistable type admitting two local minima. Mathematically, the governing system of the extended Fisher-Kolmogorov (EFK) model could be derived via an L2 gradient flow associated with the energy functional EΦ, that is,
(2)∂tΦ=−δEδΦ=ΔΦ−γΔ2Φ−fΦforx∈Ωand0<t≤T,subjected to periodic boundary conditions with an initial data Φx,0≔Φ0x, and the nonlinear bulk fv=v33−v. One intrinsic property of the EFK model is the energy dissipation law.
(3)dEdt=δEδΦ,∂tΦ=−∂tΦL22≤0.in which f,g≔∫Ωfgdx, and the associated L2 norm fL2=f,f for all f,g∈L2Ω.The existing works [4–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–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 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 [7–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–18]. 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 Yk<1+2. Becker ([21], Lemma 4.1) proved the variable-step BDF2 method of a linear parabolic problem is stable if Yk<1.868. Emmrich ([15], Theorem 3) then improved Becker’s bound-up to Yk<1.910 for the semilinear problems. As for the nonlinear parabolic problems, Emmrich ([22], Theorem 10) showed that the ratios restriction Yk<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 Yk<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 L2 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<Yk<rmax≔3+17/2≈3.561 for 2≤k≤N. 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 0=t0≤t1≤⋯≤tN=T with the time-step sizes τk≔tk−tk−1for1≤k≤N. Denote the maximum step size τ≔max1≤k≤Ntk and the local time-step ratio rk≔τk/τk−1 for k≥2. Given a grid function vkk=0N, put ∇τvn≔vn−vn−1 and ∂τvn≔∇τvn/τn for k≥1. For vn=vtn, the BDF2 formula with unequal time steps reads
(4)D2vn=∑k=1nbn−kn∇τvkforn≥1,where the discrete convolution kernels bn−kn are defined by b01≔1/τ1 for n=1, and when the time-level index n≥2,
(5)b0n≔1+2γnτn1+γn,b1n≔−γn2τn1+γnandbjn≔0forj≥2.Without losing the generality, we use the BDF1 formula to compute the first-level solution by puttingr1≡0 in (4).To present the fully discrete scheme, we describe briefly the spatial approximation. For the domainΩ=0,L2, let the grid length hx=hy=h≔L/M with an integer M. Let the discrete domainsΩh≔xh=ih,jh1≤i,j≤M and Ω¯h≔xh=ih,jh0≤i,j≤M. For the function wh=wxh, let Δxwij≔wi+1,j−wi,j/h and δx2wij=wi+1,j−2wij+wi−1,j/h2. The operators Δywij and δy2wij can be defined similarly so that the discrete gradient vector Δhwij≔Δxwij,ΔywijT and the discrete Laplacian Δhuij≔δx2+δy2uij. Denote the space of L-periodic grid functions Vh≔vhvhisL−periodicforxh∈Ω¯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 ∅hn∈Vh such that
(6)D2∅hn+γΔh2∅hn+f∅hn+S1τn∇τ∅hn=Δh∅hnforxh∈Ωhand1≤n≤N,subjected to the initial value φh0=Φ0xh and periodic boundary conditions, where S1≥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 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 S1τn∇τ∅hn 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 L2 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 rn, 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 subscriptedC, such as Cu and CΦ, denotes a generic positive constant, not necessarily the same at different occurrences; while, any subscripted c, such as cΩ, c0, c1 and c2, 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.


## 2. Solvability and Energy Dissipation Law

For any grid functionsv,w∈Vh, define the discrete L2 inner product <v,w>≔h2∑xh∈Ωhvhwh and the associated L2 norm v≔<v,w>. The discrete seminorms ∇hv and Δhv can be defined, respectively, by
(7)∇hv≔h2∑xh∈Ωh∇hvh2andΔhv≔h2∑xh∈ΩhΔhvh2forv∈Vh.The discreteL∞ norm v∞≔maxxh∈Ωh and we have the Sobolev embedding inequality
(8)v∞≤cΩv+Δhvforv∈Vh,where cΩ>0 is a constant dependent on the spatial domain Ωh. For any grid functions v,w∈Vh, the discrete Green’s formula with periodic boundary conditions yields.
(9)<Δh2v,w>=<Δhv,Δhw>and<−Δh2v,w>=<∇hv,∇hw>.


### 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 parameterS1≥1/4, 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 vh,
(10)Gz≔12b0nz−∅n−12+b1n<∇τ∅n−1,z>+γ2Δhz2+12∇hz2+14<1−z22,1>+12S1τnz−∅n−12.
For any time-level indexn≥1, the condition implies b0n+S1τn≥2b0nS1τn≥1. Then, G is strictly a convex function, for any λ∈ℝ and any ψh∈vh,
(11)d2dλ2Gz+λψλ=0=b0n+S1τn−1ψ2+γΔhψ2+γ∇hψ2+3zψ2≥γΔhψ2+∇hψ2+3zψ2.
Thus, the functionalG has a unique minimizer, denoted by ∅hn, if and only if it solves the equation
(12)0=ddλ2Gz+λψλ=0=b0n<z−∅n−1,ψ>+γ<Δhz,Δhψ>+<∇hz,∇hψ>+<fz,ψ>+S1τn<z−∅n−1,ψ>=<b0nz−∅n−1+b1n∇τ∅n−1+λΔh2z−Δh∙Δhz+fz+S1τnz−∅n−1,ψ>.
This equation holds for anyψh∈vhh if and only if the unique minimizer φjn∈vh solves
(13)b0n∇τ∅hn+b1n∇τ∅hn−1+γΔh2∅hn+f∅hn+S1τn∇τ∅hn=Δh∅hn.
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(14)Rz,s≔2+4z−z21+z−s1+sfor0≤z,s<γmax=3+172.Note that,γmax=3+17/2 solves the algebraic equation 2+3γmax–γmax2=0. Also, one has ∂/∂zRz,s=1+z−2z+1+33−1−z so that Rz,s is increasing in 0,3−1 and decreasing in 3−1,max. Rz,s is decreasing with respect to s. Thus
(15)Rz,s>Rz,γmax=2+4z−z21+z−17−14>0for0≤z<γmax.Then, the following result, also see ([24], Lemma 2), shows the BDF2 convolution kernels bn−kn are positive definite if the adjacent time-step ratios γk satisfy S1.Lemma 2.
For any real sequencewkk=1n with n entries, it holds that
(16)2wk∑j=1kbk−jkwj≥γk+11+γk+1wk2τk−γk1+γkwk−12τk−1+Rγk,γk+1wk2τkfork≥2.IfS1 holds, the discrete convolution kernels bn−kn are positive and definite in the sense that
(17)∑k=1nwk∑j=1kbk−jkwj≥12∑k=1nRγk,γk+1wk2τkforn≥1.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),
(18)E∅n≔γ2Δh∅n2+12∇h∅n2+14∅n22forn≥0.Define the following modified discrete energyε∅k by ε∅0≔Eε∅0 and
(19)ε∅k≔E∅k+γk+121+γk+1τk∇τ∅k2fork≥1.Theorem 3.
IfS1 holds and the stabilized parameter S1 fulfills
(20)S1≥1minn≥14,8Rγn,γn+1.
The BDF2 time-stepping scheme (6) preserves an energy dissipation law.
(21)ε∅n≤ε∅n−1≤ε∅0=E∅0forn≥1.Proof.
Taking the inner product of (6) by ∇τ∅n, one has
(22)<D2∅n,∇τ∅n>+γ<Δh∅n,∇τΔh∅n>+<Δh∅n,∇τΔh∅n>+<∅n3−∅n,∇τ∅n>+S1τn∇τ∅n2=0forn≥1,
where the discrete Green’s formula with periodic boundary conditions has been used. With the help of the discrete Green’s formula and2aa−b=a2–b2+a−b2, we have
(23)γ<Δh∅n,∇τΔh∅n>=γ2Δh∅n2−γ2Δh∅n−12+γ2∇τΔh∅n2.(24)<∇h∅n,∇τ∇h∅n>=12∇h∅n2−γ2∇h∅n−12+γ2∇τ∇h∅n2.
It is easy to get the following relationship(25)4a3−aa−b=a2−12−b2−12−21−a2a−b2+a2−b22≥a2−12−b2−12−2a−b2.
Then, we can obtain(26)<∅n3−∅n,∇τ∅n>≥14∅n2−12−14∅n−12−12−12∇τ∅n2.
By inserting (23)-(26) into (22), we have
(27)<D2∅n,∇τ∅n>+S1τn∇τ∅n2−12∇τ∅n2+E∅n−E∅n−1≤0forn≥1.
Applying Lemma2 with the restriction (20) of stabilized parameter S1, one has
(28)<D2∅n,∇τ∅n>≥γn+121+γn+1τn∇τ∅n2−γn21+γnτn−1∇τ∅n−12+116S1τn∇τ∅n2.
Thanks to the fact1/16S1τn+S1τ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 S1≥1/4, so that
(29)<D2∅1,∇τ∅1>=1τ1∇τ∅12≥γ221+γ2τ1∇τ∅12+18∇τ∅12.
Then, we derive from (27) that ε∅1≤E∅0=ε∅0 and complete the proof.In adaptive computations, we can choose a proper step-ratioγk+1 (especially when the current step-ratio γk or time-step τ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
(30)Rγk,γk+1≥Rγk,2≥R2,2=43,if0<γk,γk+1≤2,Rγk,3≥R3,3=12,if0<γk,γk+1≤3,Rγmax,γk+1≥1417−3,if3<γk<γmax,0<γk+1≤1.In this case, the stabilized parameterS1=0.5 is sufficient for energy stability since the Condition (20) only requires S1≥1/217−3≈0.445.Theorem 4.
IfS1 holds and the stabilized parameter S1 fulfills (20), the numerical solution of the BDF2 time-stepping scheme (6) is bounded (stable) in the L∞ norm
(31)∅n∞≤c0≔cΩ4γ−1E∅0+2+γΩhforn≥1.
Similarly, the continuous energy dissipation law givesΦtn∞≤ΦtnL∞≤c1. Note that the fixed constants c0 and c1 are independent of the spatial lengths, the time steps τn, the current time tn, and the step ratios γn.Proof.
Sincea4≥21+γa2−1+γ2 due to the simple fact a2−1+γ2≥0, we apply Theorem 3 to obtain
(32)4E∅0≥4E∅n=2γΔh∅n2+2γ∇h∅n2+∅n4−2∅n2+Ωh≥2γΔh∅n2+2γ∅n2−γ2+γΩh≥γΔh∅n+∅n2−γ2+γΩh.
So the Sobolev embedding inequality leads to(33)∅n∞2≤cΩ2Δh∅n+∅n2≤cΩ24γ−1E∅0+2+γΩh.
It implies the maximum norm estimate and completes the proof.


## 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 parameterS1≥1/4, 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 vh,
(10)Gz≔12b0nz−∅n−12+b1n<∇τ∅n−1,z>+γ2Δhz2+12∇hz2+14<1−z22,1>+12S1τnz−∅n−12.
For any time-level indexn≥1, the condition implies b0n+S1τn≥2b0nS1τn≥1. Then, G is strictly a convex function, for any λ∈ℝ and any ψh∈vh,
(11)d2dλ2Gz+λψλ=0=b0n+S1τn−1ψ2+γΔhψ2+γ∇hψ2+3zψ2≥γΔhψ2+∇hψ2+3zψ2.
Thus, the functionalG has a unique minimizer, denoted by ∅hn, if and only if it solves the equation
(12)0=ddλ2Gz+λψλ=0=b0n<z−∅n−1,ψ>+γ<Δhz,Δhψ>+<∇hz,∇hψ>+<fz,ψ>+S1τn<z−∅n−1,ψ>=<b0nz−∅n−1+b1n∇τ∅n−1+λΔh2z−Δh∙Δhz+fz+S1τnz−∅n−1,ψ>.
This equation holds for anyψh∈vhh if and only if the unique minimizer φjn∈vh solves
(13)b0n∇τ∅hn+b1n∇τ∅hn−1+γΔh2∅hn+f∅hn+S1τn∇τ∅hn=Δh∅hn.
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(14)Rz,s≔2+4z−z21+z−s1+sfor0≤z,s<γmax=3+172.Note that,γmax=3+17/2 solves the algebraic equation 2+3γmax–γmax2=0. Also, one has ∂/∂zRz,s=1+z−2z+1+33−1−z so that Rz,s is increasing in 0,3−1 and decreasing in 3−1,max. Rz,s is decreasing with respect to s. Thus
(15)Rz,s>Rz,γmax=2+4z−z21+z−17−14>0for0≤z<γmax.Then, the following result, also see ([24], Lemma 2), shows the BDF2 convolution kernels bn−kn are positive definite if the adjacent time-step ratios γk satisfy S1.Lemma 2.
For any real sequencewkk=1n with n entries, it holds that
(16)2wk∑j=1kbk−jkwj≥γk+11+γk+1wk2τk−γk1+γkwk−12τk−1+Rγk,γk+1wk2τkfork≥2.IfS1 holds, the discrete convolution kernels bn−kn are positive and definite in the sense that
(17)∑k=1nwk∑j=1kbk−jkwj≥12∑k=1nRγk,γk+1wk2τkforn≥1.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),
(18)E∅n≔γ2Δh∅n2+12∇h∅n2+14∅n22forn≥0.Define the following modified discrete energyε∅k by ε∅0≔Eε∅0 and
(19)ε∅k≔E∅k+γk+121+γk+1τk∇τ∅k2fork≥1.Theorem 3.
IfS1 holds and the stabilized parameter S1 fulfills
(20)S1≥1minn≥14,8Rγn,γn+1.
The BDF2 time-stepping scheme (6) preserves an energy dissipation law.
(21)ε∅n≤ε∅n−1≤ε∅0=E∅0forn≥1.Proof.
Taking the inner product of (6) by ∇τ∅n, one has
(22)<D2∅n,∇τ∅n>+γ<Δh∅n,∇τΔh∅n>+<Δh∅n,∇τΔh∅n>+<∅n3−∅n,∇τ∅n>+S1τn∇τ∅n2=0forn≥1,
where the discrete Green’s formula with periodic boundary conditions has been used. With the help of the discrete Green’s formula and2aa−b=a2–b2+a−b2, we have
(23)γ<Δh∅n,∇τΔh∅n>=γ2Δh∅n2−γ2Δh∅n−12+γ2∇τΔh∅n2.(24)<∇h∅n,∇τ∇h∅n>=12∇h∅n2−γ2∇h∅n−12+γ2∇τ∇h∅n2.
It is easy to get the following relationship(25)4a3−aa−b=a2−12−b2−12−21−a2a−b2+a2−b22≥a2−12−b2−12−2a−b2.
Then, we can obtain(26)<∅n3−∅n,∇τ∅n>≥14∅n2−12−14∅n−12−12−12∇τ∅n2.
By inserting (23)-(26) into (22), we have
(27)<D2∅n,∇τ∅n>+S1τn∇τ∅n2−12∇τ∅n2+E∅n−E∅n−1≤0forn≥1.
Applying Lemma2 with the restriction (20) of stabilized parameter S1, one has
(28)<D2∅n,∇τ∅n>≥γn+121+γn+1τn∇τ∅n2−γn21+γnτn−1∇τ∅n−12+116S1τn∇τ∅n2.
Thanks to the fact1/16S1τn+S1τ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 S1≥1/4, so that
(29)<D2∅1,∇τ∅1>=1τ1∇τ∅12≥γ221+γ2τ1∇τ∅12+18∇τ∅12.
Then, we derive from (27) that ε∅1≤E∅0=ε∅0 and complete the proof.In adaptive computations, we can choose a proper step-ratioγk+1 (especially when the current step-ratio γk or time-step τ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
(30)Rγk,γk+1≥Rγk,2≥R2,2=43,if0<γk,γk+1≤2,Rγk,3≥R3,3=12,if0<γk,γk+1≤3,Rγmax,γk+1≥1417−3,if3<γk<γmax,0<γk+1≤1.In this case, the stabilized parameterS1=0.5 is sufficient for energy stability since the Condition (20) only requires S1≥1/217−3≈0.445.Theorem 4.
IfS1 holds and the stabilized parameter S1 fulfills (20), the numerical solution of the BDF2 time-stepping scheme (6) is bounded (stable) in the L∞ norm
(31)∅n∞≤c0≔cΩ4γ−1E∅0+2+γΩhforn≥1.
Similarly, the continuous energy dissipation law givesΦtn∞≤ΦtnL∞≤c1. Note that the fixed constants c0 and c1 are independent of the spatial lengths, the time steps τn, the current time tn, and the step ratios γn.Proof.
Sincea4≥21+γa2−1+γ2 due to the simple fact a2−1+γ2≥0, we apply Theorem 3 to obtain
(32)4E∅0≥4E∅n=2γΔh∅n2+2γ∇h∅n2+∅n4−2∅n2+Ωh≥2γΔh∅n2+2γ∅n2−γ2+γΩh≥γΔh∅n+∅n2−γ2+γΩh.
So the Sobolev embedding inequality leads to(33)∅n∞2≤cΩ2Δh∅n+∅n2≤cΩ24γ−1E∅0+2+γΩh.
It implies the maximum norm estimate and completes the proof.


## 3.L2 Norm Convergence Analysis

To facilitate the subsequent numerical analysis, we introduce a modified BDF2 formula(34)D~2vn≔D2vn+S1τn∇τvn=∑k=1ndn−kn∇τvkforn≥1,where the modified BDF2 kernels dn−kn are defined by
(35)d0n≔b0n+S1τnanddjn≔bjnfor1≤j≤n−1,Thus the stabilized BDF2 implicit scheme (6) reads
(36)D~2∅hn+γΔh2∅hn+f∅hn=Δh∅hnforxh∈Ωhand1≤n≤N.


### 3.1. A New Class of DOC Kernels

We define a new class of DOC kernelsvn−knk=1n.
(37)v0n≔1d0nforn≥3andvn−kn≔−1d0k∑j=k+1nvn−jndj−kjfor3≤k≤n−1.This type of DOC kernel was first introduced in [24] with respect to the original BDF2 convolution kernels bn−kn . It is easy to check that the following discrete orthogonal identity
(38)∑j=knvn−jndj−kj≡δnkfor∀1≤k≤n.This identity directly leads to the following relationship(39)∑j=1nvk−jkD~2vj=∑j=1kvk−jk∑l=1jdj−lj∇τvl=∑l=1k∇τvl∑j=lkvk−jkdj−lj=∇τvk,where the order of summation has been exchanged in the second equality. This equality (39) plays an essential role in the subsequent L2 norm analysis.Applying Lemma2 and the definition (35) of modified BDF2 kernels dn−kn, it is not difficult to check the following result.Lemma 5.
IfS1 holds and the parameter S1≥0, the modified BDF2 kernels dk=1n in (35) are positive definite in the sense that, for any nonzero sequence wkk=1n,
(40)∑k=1nwk∑j=1kdk−jkwj>0.Lemma 6.
If the discrete convolution kernelsdn−kn defined in (35) are positive definite, then the DOC kernels vn−jn defined in (37) satisfy the following properties:
(I)
The discrete kernelsvn−jn are positive definite, and then the symmetric matrix
Kn×n≔Θ+ΘT is positive definite
where the elements ofΘ are defined by Θij for i>j and Θij≔1/2v≤ji−ji for i≤j.
The discrete kernelsvn−jn are positive and vn−jn=1/d0j∏i=j+1nβi for 1≤k≤n, where
(41)βi≔−d1id0i=γi21+2γi+S1τi1+γi≤γi21+2γifor2≤i≤n.(II)
∑j=1nvn−jn≤τn such that ∑k=1n∑j=1kvk−jk≤tn for n≥1Proof.
Following the proof of ([11], Lemma 6), one can prove the positive definiteness of vn−jn by using Lemma 5, and the orthogonal relationship between the two classes of discrete convolution kernels vn−jn and dn−jn. To prove (II), define
(42)v^n−kn≔vn−knd0kfor1≤k≤n.
The definition (37) yields v^0n=1 for n≥1. For 1≤k≤n−1, we use the definition (37) to find a recursive procedure
(43)v^n−kn=−d1k+1d0k+1v^n−k−1n=βk+1v^n−k−1n
Sinceβi>0, we obtain v^n−kn=∏i=k+1nβi>0 for 1≤k≤n, and confirm (II) immediately. Lastly, taking vn=tn in equality (3.6), we have
(44)∑j=1nvn−knD~tj∇τtn=τnforn≥1.
Thus the result of (III) follows from the definition (34) since D~2tj=1+S1τj2≥1.


### 3.2.L2 Norm Error Estimate

Now, we are to establish theL2 norm error estimate of our BDF2 scheme (6) via the equivalent form (36). Let the local consistency error at the time =tj,
(45)εj≔D~2Φtj−∂tΦtj.Consider a convolutional consistency errorΞk defined by
(46)Ξk≔∑j=1kvn−knD2Φtj−∂tΦtj+S1τj∇τΦtjfork≥1.By adding the error termS1τj∇τΦtj, stemmed from the stabilized term S1τj∇τ∅hn with a proper parameter S1=O1, follow to the proof of ([13], Lemma 3.4), we arrive at the following result.Lemma 7.
IfS1 holds, the convolutional consistency error ΞΦk in (46) satisfies
(47)Ξk≤vk−1k∫0t1Φ′′sds+∑j=1kvk−jkτj∫tj−1tjS1Φ′s+3Φ′′′sds,k≥1.
Such that(48)∑k=1nΞk≤C∅τ12∑k=1n∏i=2kγi21+2γi+C∅τ2forn≥1Theorem 8.
Assume that the EFK problem (2) has a solution Φ∈Cx,t6,3Ω×0,T. If S1 holds, the stabilized parameter S1 fulfills (20), and the maximum step size τ is sufficiently small such that τ≤1/4c2, the stabilized BDF2 implicit time-stepping scheme (6) with unequal time steps is convergent in the L2 norm,
(49)Φn−∅n≤C∅exp4c2tn−1τ12∑k=1n∏i=2kγi21+2γi+tnτ2+h2,1≤n≤N.
Here,c2≔1+c02+c0c1+c12 is dependent on the domain Ω and the initial values Φ0 and ∅0, but independent of the spatial length h, the time tn, the step sizes τn, and the step ratios rn.Proof.
Let the error functionehn≔Φhn–∅hn with eh0≔0 for xh∈Ω¯h. We have the following error equation
(50)D~2ehj+γΔh2ehj−Δhehj+fΦhj−fϕhj=εhj+ηhj,xh∈Ω¯h,1≤j≤N,where εhj and ηhj denote the local consistency error in time and space, respectively. Under the regularity setting of the solution, Lemma 6 (III) gives
(51)∑k=1n∏k≤C∅tnh2for1≤n≤Nwhere∏hk≔∑j=1kvk−jkηhjfor xh∈Ωh and 1≤k≤N, where Ξhk and ∏hkare defined by (46) and (51), respectively. Making the inner product of the above equality with 2ek and summing up the superscript from k=1 to n, we apply the discrete Green’s formula to derive that
(52)en2−e02+2γ∑k=1n∑j=1kvk−jk<Δhej,Δhek>+2∑k=1n∑j=1kvk−jk<∇hej,∇hek>≤2∑k=1n∑j=1kvk−jk<f∅j−fΦj,ek>+2∑k=1n<Ξk+Πk,ek>,for 1≤n≤N. Thanks to the maximum norm estimates in Theorem 3, one gets
(53)fϕhj−fΦhj=ϕhj2+ϕhjΦhj+Φhj2−1ehj≤c2ehj,xh∈Ωh,so that the Cauchy-Schwarz inequality leads to
(54)2∑k=1n∑j=1kvk−jk<f∅j−fΦj,ek>≤2c2∑k=jn∑j=1kvk−jkejek.
Then, it follows from (52) that
(55)en2−e02+2c2∑k=1nek∑j=1kvk−jk+2∑k=1nekΞk+Πk.
Choose some integern00≤n0≤n such that en0=max0≤k≤nek. Let n=n0 in the above inequality (55), such that
(56)en≤en0≤e0+2c2∑k=1n0τkek+2∑k=1nΞk+Πk.
Under the settingτ≤1/4c2, we have
(57)en≤2e0+4c2∑k=1n−1τkek+4∑k=1nΞk+Πk,1≤n≤N.
Obviously, applying the discrete Gronwall inequality ([24], Lemma 5), one can obtain
(58)en≤2exp4c2tn−1e0+2∑k=1nΞk+2∑k=1nΠk,1≤n≤N.
The proof is completed from Lemma7 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 S1 because τ1∑k=1n∏i=2kγi2/1+2γi≤tn. The second-order convergence in time can be recovered if S2. The time-step ratios rk are contained in S1, but almost all of them less than 1+2, or R=N0≪N, where R is an index set R≔k1+2≤γk<3+17/2. Actually, very large steps always result in a loss of numerical accuracy, although the condition S1 permits us to use a series of increasing time steps with amplification factors up to 3.561. The condition S2 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 Cr such that ∑k=1n∏i=2kγi2/1+2γi≤Cr. Thus, we have the following corollary.Corollary 9.
Assume that the EFK problem (2) has a unique smooth solution. If S2 holds, the stabilized parameter S1 fulfills (20) and the maximum step size τ≤1/4c2, the stabilized BDF2 scheme (6) is second-order convergent in the L2 norm,
(59)Φn−ϕn≤C∅exp4c2tn−1Cγτ12+tnτ2+h2for1≤n≤N.At the same time, under the mild step-ratio conditionS1, 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 S1 holds with the maximum step size τ≤1/4c2, the stabilized BDF2 implicit time-stepping scheme (6) without the first-level BDF1 scheme is convergent in the L2 norm,
(60)Φn−ϕn≤2exp4c2tn−1Φ1−ϕ1+2tnΦ1−ϕ1+C∅tnτ2−h2,2≤n≤N.


## 3.1. A New Class of DOC Kernels

We define a new class of DOC kernelsvn−knk=1n.
(37)v0n≔1d0nforn≥3andvn−kn≔−1d0k∑j=k+1nvn−jndj−kjfor3≤k≤n−1.This type of DOC kernel was first introduced in [24] with respect to the original BDF2 convolution kernels bn−kn . It is easy to check that the following discrete orthogonal identity
(38)∑j=knvn−jndj−kj≡δnkfor∀1≤k≤n.This identity directly leads to the following relationship(39)∑j=1nvk−jkD~2vj=∑j=1kvk−jk∑l=1jdj−lj∇τvl=∑l=1k∇τvl∑j=lkvk−jkdj−lj=∇τvk,where the order of summation has been exchanged in the second equality. This equality (39) plays an essential role in the subsequent L2 norm analysis.Applying Lemma2 and the definition (35) of modified BDF2 kernels dn−kn, it is not difficult to check the following result.Lemma 5.
IfS1 holds and the parameter S1≥0, the modified BDF2 kernels dk=1n in (35) are positive definite in the sense that, for any nonzero sequence wkk=1n,
(40)∑k=1nwk∑j=1kdk−jkwj>0.Lemma 6.
If the discrete convolution kernelsdn−kn defined in (35) are positive definite, then the DOC kernels vn−jn defined in (37) satisfy the following properties:
(I)
The discrete kernelsvn−jn are positive definite, and then the symmetric matrix
Kn×n≔Θ+ΘT is positive definite
where the elements ofΘ are defined by Θij for i>j and Θij≔1/2v≤ji−ji for i≤j.
The discrete kernelsvn−jn are positive and vn−jn=1/d0j∏i=j+1nβi for 1≤k≤n, where
(41)βi≔−d1id0i=γi21+2γi+S1τi1+γi≤γi21+2γifor2≤i≤n.(II)
∑j=1nvn−jn≤τn such that ∑k=1n∑j=1kvk−jk≤tn for n≥1Proof.
Following the proof of ([11], Lemma 6), one can prove the positive definiteness of vn−jn by using Lemma 5, and the orthogonal relationship between the two classes of discrete convolution kernels vn−jn and dn−jn. To prove (II), define
(42)v^n−kn≔vn−knd0kfor1≤k≤n.
The definition (37) yields v^0n=1 for n≥1. For 1≤k≤n−1, we use the definition (37) to find a recursive procedure
(43)v^n−kn=−d1k+1d0k+1v^n−k−1n=βk+1v^n−k−1n
Sinceβi>0, we obtain v^n−kn=∏i=k+1nβi>0 for 1≤k≤n, and confirm (II) immediately. Lastly, taking vn=tn in equality (3.6), we have
(44)∑j=1nvn−knD~tj∇τtn=τnforn≥1.
Thus the result of (III) follows from the definition (34) since D~2tj=1+S1τj2≥1.


## 3.2.L2 Norm Error Estimate

Now, we are to establish theL2 norm error estimate of our BDF2 scheme (6) via the equivalent form (36). Let the local consistency error at the time =tj,
(45)εj≔D~2Φtj−∂tΦtj.Consider a convolutional consistency errorΞk defined by
(46)Ξk≔∑j=1kvn−knD2Φtj−∂tΦtj+S1τj∇τΦtjfork≥1.By adding the error termS1τj∇τΦtj, stemmed from the stabilized term S1τj∇τ∅hn with a proper parameter S1=O1, follow to the proof of ([13], Lemma 3.4), we arrive at the following result.Lemma 7.
IfS1 holds, the convolutional consistency error ΞΦk in (46) satisfies
(47)Ξk≤vk−1k∫0t1Φ′′sds+∑j=1kvk−jkτj∫tj−1tjS1Φ′s+3Φ′′′sds,k≥1.
Such that(48)∑k=1nΞk≤C∅τ12∑k=1n∏i=2kγi21+2γi+C∅τ2forn≥1Theorem 8.
Assume that the EFK problem (2) has a solution Φ∈Cx,t6,3Ω×0,T. If S1 holds, the stabilized parameter S1 fulfills (20), and the maximum step size τ is sufficiently small such that τ≤1/4c2, the stabilized BDF2 implicit time-stepping scheme (6) with unequal time steps is convergent in the L2 norm,
(49)Φn−∅n≤C∅exp4c2tn−1τ12∑k=1n∏i=2kγi21+2γi+tnτ2+h2,1≤n≤N.
Here,c2≔1+c02+c0c1+c12 is dependent on the domain Ω and the initial values Φ0 and ∅0, but independent of the spatial length h, the time tn, the step sizes τn, and the step ratios rn.Proof.
Let the error functionehn≔Φhn–∅hn with eh0≔0 for xh∈Ω¯h. We have the following error equation
(50)D~2ehj+γΔh2ehj−Δhehj+fΦhj−fϕhj=εhj+ηhj,xh∈Ω¯h,1≤j≤N,where εhj and ηhj denote the local consistency error in time and space, respectively. Under the regularity setting of the solution, Lemma 6 (III) gives
(51)∑k=1n∏k≤C∅tnh2for1≤n≤Nwhere∏hk≔∑j=1kvk−jkηhjfor xh∈Ωh and 1≤k≤N, where Ξhk and ∏hkare defined by (46) and (51), respectively. Making the inner product of the above equality with 2ek and summing up the superscript from k=1 to n, we apply the discrete Green’s formula to derive that
(52)en2−e02+2γ∑k=1n∑j=1kvk−jk<Δhej,Δhek>+2∑k=1n∑j=1kvk−jk<∇hej,∇hek>≤2∑k=1n∑j=1kvk−jk<f∅j−fΦj,ek>+2∑k=1n<Ξk+Πk,ek>,for 1≤n≤N. Thanks to the maximum norm estimates in Theorem 3, one gets
(53)fϕhj−fΦhj=ϕhj2+ϕhjΦhj+Φhj2−1ehj≤c2ehj,xh∈Ωh,so that the Cauchy-Schwarz inequality leads to
(54)2∑k=1n∑j=1kvk−jk<f∅j−fΦj,ek>≤2c2∑k=jn∑j=1kvk−jkejek.
Then, it follows from (52) that
(55)en2−e02+2c2∑k=1nek∑j=1kvk−jk+2∑k=1nekΞk+Πk.
Choose some integern00≤n0≤n such that en0=max0≤k≤nek. Let n=n0 in the above inequality (55), such that
(56)en≤en0≤e0+2c2∑k=1n0τkek+2∑k=1nΞk+Πk.
Under the settingτ≤1/4c2, we have
(57)en≤2e0+4c2∑k=1n−1τkek+4∑k=1nΞk+Πk,1≤n≤N.
Obviously, applying the discrete Gronwall inequality ([24], Lemma 5), one can obtain
(58)en≤2exp4c2tn−1e0+2∑k=1nΞk+2∑k=1nΠk,1≤n≤N.
The proof is completed from Lemma7 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 S1 because τ1∑k=1n∏i=2kγi2/1+2γi≤tn. The second-order convergence in time can be recovered if S2. The time-step ratios rk are contained in S1, but almost all of them less than 1+2, or R=N0≪N, where R is an index set R≔k1+2≤γk<3+17/2. Actually, very large steps always result in a loss of numerical accuracy, although the condition S1 permits us to use a series of increasing time steps with amplification factors up to 3.561. The condition S2 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 Cr such that ∑k=1n∏i=2kγi2/1+2γi≤Cr. Thus, we have the following corollary.Corollary 9.
Assume that the EFK problem (2) has a unique smooth solution. If S2 holds, the stabilized parameter S1 fulfills (20) and the maximum step size τ≤1/4c2, the stabilized BDF2 scheme (6) is second-order convergent in the L2 norm,
(59)Φn−ϕn≤C∅exp4c2tn−1Cγτ12+tnτ2+h2for1≤n≤N.At the same time, under the mild step-ratio conditionS1, 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 S1 holds with the maximum step size τ≤1/4c2, the stabilized BDF2 implicit time-stepping scheme (6) without the first-level BDF1 scheme is convergent in the L2 norm,
(60)Φn−ϕn≤2exp4c2tn−1Φ1−ϕ1+2tnΦ1−ϕ1+C∅tnτ2−h2,2≤n≤N.


## 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 10−12.


### 4.1. Tests on Random Time Meshes

Example 11.
Consider the exterior-forced EFK equation∂tΦ+YΔ2Φ−ΔΦ+fΦ=gX,t in the domain Ω=0,π2 with γ=0.02 such that it has a solution Φ=costsin2xsin2y.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=∑k=1Nσk. The discrete L2 norm error eN≔ΦT−∅N is recorded in each run and the experimental order of convergence is computed by Order≈logeN/e2N/logτN/τ2N, in which τ (N) denotes the maximum time-step size for total N subintervals.In our computations, we take the stabilization parameterS1=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 N1 in Table 1) of time levels with the step ratio γk≥3+17/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.Table 1
Accuracy of BDF2 scheme (6) on random time mesh.
NτeNOrdermaxγkN1209.57e-023.53e-04—19.223405.09e-025.48e-052.9533.387802.80e-022.09e-061.6267.54151601.37e-025.04e-061.9854.10213206.03e-031.03e-061.9486.0248To exploit the numerical performance of the stabilized termS1τn∇t∅hn in the BDF2 scheme (6), we also run the scheme (6) by varying the value S1. Note that the parameter values used in our examinations are the same as before but the stabilization parameter S1. The error curves obtained using different stabilization parameters S1=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 S1 gives the less accurate solution. Thus the stabilization parameter S1 should be selected as small as possible to maintain the numerical accuracy.Figure 1
The semilog plot of numerical error e(N) of stabilized BDF2 scheme (6) with different stabilization parameters S1.


### 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 Algorithm1 refer to the backward Euler scheme and BDF2 scheme (6), respectively. The adaptive time step τada is given by τadae,τcur=tol/eρτ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 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.Algorithm 1: Adaptive time-stepping strategy.
Require: Set τ1≔τmin. Given ∅n and time step τn1: Compute∅1n+1 by using BDF1 scheme with time step τn.2: Compute∅2n+1 by using BDF2 scheme with time step τn.3: Calculateen+1=∅2n+1−∅1n+1/∅2n+1.4: ifen<tol then5:  Update time-step sizeτn+1⟵minmaxτmin,τada,2.414τn,τmax.6: else7:  Recalculate with time-step sizeτn⟵maxτmin,τada.8:  Goto 19: end if


### 4.3. Numerical Application

Example 12.
Consider the EFK Equation (2) with the following initial data
(61)ΦX,0=0.1sin3xsin2y+sin5xsin5y,in the domain Ω=0,2π2 with parameter γ=10−4 . We use 128×128 uniform mesh to discrete the spatial domain Ω and take the stabilization parameter S1=0.5 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τ=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.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.
(a)(b)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.Figure 3
Solution snapshots of the BDF2 scheme (6) for the EFK equation using adaptive time strategy at t=0,0.1,0.9,1,1.1,15, respectively.Figure 4
Evolutions of energy (a) and adaptive time steps (b) of the BDF2 scheme (6) for the EFK equation using adaptive time strategy.
(a)(b)


## 4.1. Tests on Random Time Meshes

Example 11.
Consider the exterior-forced EFK equation∂tΦ+YΔ2Φ−ΔΦ+fΦ=gX,t in the domain Ω=0,π2 with γ=0.02 such that it has a solution Φ=costsin2xsin2y.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=∑k=1Nσk. The discrete L2 norm error eN≔ΦT−∅N is recorded in each run and the experimental order of convergence is computed by Order≈logeN/e2N/logτN/τ2N, in which τ (N) denotes the maximum time-step size for total N subintervals.In our computations, we take the stabilization parameterS1=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 N1 in Table 1) of time levels with the step ratio γk≥3+17/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.Table 1
Accuracy of BDF2 scheme (6) on random time mesh.
NτeNOrdermaxγkN1209.57e-023.53e-04—19.223405.09e-025.48e-052.9533.387802.80e-022.09e-061.6267.54151601.37e-025.04e-061.9854.10213206.03e-031.03e-061.9486.0248To exploit the numerical performance of the stabilized termS1τn∇t∅hn in the BDF2 scheme (6), we also run the scheme (6) by varying the value S1. Note that the parameter values used in our examinations are the same as before but the stabilization parameter S1. The error curves obtained using different stabilization parameters S1=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 S1 gives the less accurate solution. Thus the stabilization parameter S1 should be selected as small as possible to maintain the numerical accuracy.Figure 1
The semilog plot of numerical error e(N) of stabilized BDF2 scheme (6) with different stabilization parameters S1.


## 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 Algorithm1 refer to the backward Euler scheme and BDF2 scheme (6), respectively. The adaptive time step τada is given by τadae,τcur=tol/eρτ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 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.Algorithm 1: Adaptive time-stepping strategy.
Require: Set τ1≔τmin. Given ∅n and time step τn1: Compute∅1n+1 by using BDF1 scheme with time step τn.2: Compute∅2n+1 by using BDF2 scheme with time step τn.3: Calculateen+1=∅2n+1−∅1n+1/∅2n+1.4: ifen<tol then5:  Update time-step sizeτn+1⟵minmaxτmin,τada,2.414τn,τmax.6: else7:  Recalculate with time-step sizeτn⟵maxτmin,τada.8:  Goto 19: end if


## 4.3. Numerical Application

Example 12.
Consider the EFK Equation (2) with the following initial data
(61)ΦX,0=0.1sin3xsin2y+sin5xsin5y,in the domain Ω=0,2π2 with parameter γ=10−4 . We use 128×128 uniform mesh to discrete the spatial domain Ω and take the stabilization parameter S1=0.5 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τ=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.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.
(a)(b)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.Figure 3
Solution snapshots of the BDF2 scheme (6) for the EFK equation using adaptive time strategy at t=0,0.1,0.9,1,1.1,15, respectively.Figure 4
Evolutions of energy (a) and adaptive time steps (b) of the BDF2 scheme (6) for the EFK equation using adaptive time strategy.
(a)(b)



---

*Source: 1869660-2023-03-21.xml*