Stability and errors analysis of two iterative schemes of fractional steps type associated to a nonlinear reaction-diffusion equation

We present the error analysis of two time-stepping schemes of fractional steps type, used in the discretization of a nonlinear reaction-diffusion equation with Neumann boundary conditions, relevant in phase transition and interface problems. We start by investigating the solvability of a such boundary value problems in the class \begin{document}$ W^{1,2}_p(Q) $\end{document} . One proves the existence, the regularity and the uniqueness of solutions, in the presence of the cubic nonlinearity type. The convergence and error estimate results (using energy methods) for the iterative schemes of fractional steps type, associated to the nonlinear parabolic equation, are also established. The advantage of such method consists in simplifying the numerical computation. On the basis of this approach, a conceptual algorithm is formulated in the end. Numerical experiments are presented in order to validates the theoretical results (conditions of numerical stability) and to compare the accuracy of the methods.


1.
Introduction. Consider a one dimensional nonlinear reaction-diffusion equation with respect to the unknown function v(t, x): where: • Ω is a bounded domain in IR with smooth boundary ∂Ω = Γ and T > 0 stands for some final time; • v(t, x) is the unknown function; in particular, v(t, x) is the phase function (used to distinguish between the states (phases) of a material which occupies the region Ω at every time t ∈ [0, T ]); • p 1 , p 2 , p 3 are positive values; • F : IR −→ IR is a real function having the structure: F (ϕ) = ϕ − ϕ 3 ∀ϕ ∈ IR.
Key words and phrases. Nonlinear PDE of parabolic type, reaction-diffusion equations, finite difference methods, fractional steps method, stability and convergence of numerical methods, performance of numerical algorithms, thermodynamics, phase-changes.
Here we have used the standard notation for Sobolev spaces, namely, given a positive integer k and 1 ≤ p ≤ ∞, we denote by W k,2k p (Q) the usual Sobolev space on Q: W k,2k p (Q) = y ∈ L p (Q) : ∂ r ∂t r ∂ s ∂x s y ∈ L p (Q), for 2r + s ≤ 2k , i.e., the space of functions whose t-derivatives and x-derivatives up to the order k and 2k, respectively, belong to L p (Q). Also, we have used the Sobolev spaces W l p (Ω) with nonintegral l for the initial and boundary conditions, respectively (see [19] -Chapter 1 and references therein).
(1) it is a particular instance of the Allen-Cahn equation [2], which was introduced to describe the motion of anti-phase boundaries in crystalline solids, and it has been widely applied to many [9] complex moving interface problems, e.g., the mixture of two incompressible fluids, the nucleation of solids, the vesicle membranes.
The error analysis for the implicit backward Euler and finite elements approximation is presented in [32], while a discontinuous-Galerkin in time method is analyzed in [11]. Computations with several different higher-order time-stepping schemes, such as BDF2-AB2 and Crank-Nicolson, are used for the sharp interface limit in [37].
The outline of the paper is as follows. In Section 2 we prove the existence, regularity, stability and uniqueness of the solution to the nonlinear problem (1) in the presence of homogeneous Neumann boundary conditions, while, in Section 3 we are concerned with the convergence of fractional steps type scheme associated to the nonlinear reaction-diffusion equation (1) 1 . In Section 4 we introduce two semi-discrete in time approximations, proving consistency and stability results for error equations by energy estimates arguments. The numerical experiments in Section 5 confirms the theoretical rates of convergence. The concluding remarks are formulated in Section 6.
2. Well-posedness of solutions to the nonlinear equation (1). In the present Section we will investigate the solvability of the first boundary value problems of the form (1) in the class W 1,2 p (Q). One proves the existence, the regularity and the uniqueness of solutions (Theorem 2.1 below) to the nonlinear parabolic problem (1) considering the cubic nonlinearity p 3 (v − v 3 ) which verifies for N = 1 and r = 3 the general assumptions H 0 and H 2 formulated in [30], [31], that is: H 2 : There exist a functionF : IR 2 → IR and a constant b 0 > 0 verifying the relations: The basic tools in the analysis of the problem (1) are the Leray-Schauder degree theory [13], the L p -theory of linear and quasi-linear parabolic equations [19], as well as the Lions and Peetre embedding Theorem [20], p. 24, which ensures the existence of a continuous embedding W 1,2 p (Q) ⊂ L µ (Q), where the number µ is defined as follows The main result of this section establishes the dependence of the solution v(t, x) in the nonlinear parabolic equation (1) on the term f (t, x) in the right-hand side.
where the constant C depends on |Ω|, T , M , p, p 1 , p 2 and p 3 , but is independent of v and f . If v 1 , v 2 are two solutions of (1) corresponding to the data f 1 , f 2 ∈ L p (Q), respectively, for the same initial conditions, such that where the constant C depends on |Ω|, T , M , p, p 1 , p 2 , p 3 and b 0 but is independent of v 1 , v 2 , f 1 and f 2 .
Proof. The proof of Theorem 2.1 was given in [5] and [27] noting that there formulation differs from this by certain physical parameters, which implies different values for the constant C in (3) and (4). Moreover, corresponding to different boundary conditions (including nonlinear and nonhomogeneous boundary conditions), similar results were proved in [8], [21], [23], [26] and [28]. Proof. Let f 1 = f 2 = f in the Theorem 2.1. Then (5) shows that the conclusion of the corollary is true.
3. Approximating schemes. The aim of this Section is to use the fractional steps method in order to approximate the solution of the nonlinear problem (1) 1 , whose uniqueness is guaranteed by Corollary 2.2. To do that, let us associate to the time-interval [0, T ] the equidistant grid of length ε M = T M , for any integer M ≥ 1. Corresponding to it, the following two approximating schemes can be written in order to approximate the solution of nonlinear boundary value problem (1): where z M (ε M , ·) is the solution of Cauchy problem as well as with Z M (ε M , ·) -solution of the Cauchy problem We point out that the sequence of approximating problems (6)- (7) and (8)- (9), supplies a decoupling method for the original problem (1) into a linear parabolic boundary value problem (6) or (8) and, corresponding, a nonlinear evolution equation (7) or (9). The advantage of this approach consists in simplifying the numerical computation of the approximations to (1), due to that the fractional steps method avoids the iterative process in passing from a time level to the next one.
The main question is the convergence of the sequence {v M } ({V M }, respectively) of solutions to the approximate problems (6)-(7) ((8)- (9), respectively) to the unique solution v ∈ W 1,2 p (Q) of problem (1) as M −→ ∞. We will treat the convergence of this numerical schemes on the basis of an abstract approximation result (see [24], [26] for a detailed discussion).
We have Theorem 3.1 Assume that the function y → p 3 (y 3 − y) + ωy is strictly increasing Proof. The proof of Theorem 3.1 can be found in [28]; we omit the details.
We note that the difference between the above schemes is given by the position of the linear term from the nonlinearities F (v) in the original problem (1) which can appears in Cauchy problem (7) or in the linear problem (8), as was done in previous works: [3], [4], [15], [16], [23]- [29] and [32].
. Moreover, the fractional step method (11) can be equivalently written as Error estimates for the numerical schemes (10)-(11) will be presented in the sequel. In order to facilitate the reader's easy tracking the next computations, let us introduce the following notation We begin by defining the point-wise truncation errors with respect to the time discretization, that is: and the local truncation errors, respectively: By subtracting (10)-(11) respectively from (12)-(13), we obtain the following equations in errors where

Consistency result.
In what follows we will assume that there exists a maximum principle result for solutions of the fractional steps methods (10)- (11) and v 0 (x) ∈ [−1, 1] for a.e. x ∈ Ω. We have the following result regarding consistency.
Lemma 4.1. Assuming that the weak solution of (1) also satisfies v ∈ W 1,2 p (Q). Then the local truncation errors satisfy: which substituting into (12) and using the original equation (1) evaluated at t n+1 , i.e., we obtain Using the maximum principle for the exact solution we have and therefore using Hölder's inequality we obtain concluding the proof of (18). Further, we will focus our attention on the inequality (19) expressing the consistency of the fractional steps method (11). Thus, making use of the following computation (see the original equation (1) then the local truncation error (13) writes Let's observe that the square brackets from the last term in the relationship above is smaller or equal to 1. Consequently, in accordance with our work assumption regarding the existence of a maximum principle result for the exact solution in the fractional steps method (11), we have Now, using (a + b) m ≤ 2 m−1 (a m + b m ) as well as the Hölder's inequality, we obtain

COSTICȂ MOROŞANU
Finally, summing the previous inequality for n = 0, ..., M − 1 and then multiplying the result with ε M , we get (19) which conclude the proof of the Lemma 4.1.

Stability result.
For future reference, we recall here the following form of the discrete Grönwall lemma. Assume w n , α n , q n ≥ 0, β ∈ [0, 1) satisfy w n + q n ≤ α n + β n κ=1 w κ ∀κ ≥ 0, where {α n } is non-decreasing. Then To obtain convergence results for the numerical methods (10) and (11), we will prove a stability result using energy estimates. In this respect, we begin by testing (14) and (15) with e n+1 |e n+1 | p−2 .
We get for any p ≥ 2. In the following, we will analyze each term in (21)-(22) separately.
T4 & T6. Using the Young inequality with a = |E n+1 |, b = |e n+1 | p−1 , m = p and n = p p−1 , the local truncation error term gives, where E n+1 = E n+1 old or E n+1 = E n+1 new . As we have proposed previously, at this point we have the estimations (23)- (27), corresponding to the terms T1 -T6. Now we will continue with the evaluation of relations (21)- (22). We first substitute in (21) the relations (23)- (24) and (27) corresponding to E n+1 = E n+1 old to obtain In the fractional step case (10) we obtain from (25) and (29) e M p L p (Ω) + p(p − 1) Therefore, assuming that the time step is small enough we obtain from Grönwall's inequality (20) the following stability estimate e M p L p (Ω) + p(p − 1) In the fractional steps method (11) we obtain from (23)- (24) and (27) corre- Therefore, assuming that the time-step is small enough using Grönwall's inequality (20), we obtain from (31) the following stability estimate e M p L p (Ω) + p(p − 1) Corresponding to Ω, already chosen in one dimension, the boundary ∂Ω is reduced to the set {0, b}. To approximate V x (0) (V x (b)) we will use a backward (forward) finite differences; this leads to where V i 0 and V i N +1 are dummy variables. Finally we refer to the term p 3 V (t i , x j ) in (10). To approximate this quantity (the reaction term), we will involve an implicit formula, i.e.: We are now ready to build the approximation scheme, mentioned at the beginning.
If we set than the system (43), coupled with (40), can be rewritten in matrix form as where Therefore, the general design of the algorithm to calculate the approximate solution of nonlinear system (1), via fractional steps method and 1-IMBDF, is the following one Compute V i solving the linear system (44); End-for; End.
As it is well known, most initial value problems reduce to solving large sparse linear systems of the form (44). For later use (e.g., numerical implementation of conceptual algorithms), we will proof the following then the matrix coefficients in linear system (44) can be factored into the product of a lower-triangular matrix and an upper-triangular matrix (LU -factorization).
Proof. Let denote by a mn , m, n = 1, 2, · · · , N , the elements of matrix coefficients in linear system (44). Analyzing the main diagonal elements of matrix A, we finding that c 1 + c 2 = 0 reflect the assumptions expressed in (45), as well as that c 2 = 0. Thus, we find easily that a mm = 0 ∀m = 1, 2, · · · , N and so the Gaussian elimination can be performed on the system (44) without interchanges; consequently A has an LU factorization.

Stability conditions.
To establish conditions of stability for the linear difference equation (44), introduced in the previous subsection, we will use in our analysis the Lax-Richtmyer definition of stability, expressed in terms of norm · ∞ (see [29] and references therein). The equation (44) may be rewritten in a more convenient form as (the existence of A −1 will be proved in the Proposition 3.1 below). In addition, the matrix A can be written in the form where D = diag(c 1 +c 2 , c 2 , · · · , c 2 , c 1 +c 2 ) and G = A−D. Thus, noting a 1 = c 1 +c 2 , we have and a simple analysis of all lines in matrix D -1 G allows us to deduce that we only have two distinct lines. The sum of each such line is written in vector w below Let's denote by v min = min{|a 1 |, |c 2 |}. Now we are able to prove the following result with respect to the stability in matrix equation (46).
then the equation (46) is stable. Otherwise, it is unstable.
Proof. The proof is reduced to checking the condition of stability which, based on the Lax-Richtmyer definition mentioned above and taking into account the relation (46), it reduces to check the inequality We begin our analysis by determining an estimate for D −1 G ∞ . As we have already noted (see relation (48)), this is equivalent with the following equality: D −1 G ∞ = max |w|, wherefrom we easily derive the estimate The estimate (49) allows us now to prove the existence of A -1 . Indeed, since by hypothesis we have assumed that 2|c 1 | < v min than D −1 G ∞ < 1 which guarantees that there exist (I +D -1 G) −1 . Moreover, there exist A -1 and A -1 = (I +D -1 G) -1 D -1 . Using the well known inequality: (I + D -1 G) -1 ∞ ≤ 1 1− D −1 G ∞ and making use of (47), it follows that vmin > 0, we easily deduce from this that Since D −1 ∞ ≤ 1 vmin and involving the above estimate, from (50) we finally obtain Now we turn our attention to matrix B. Analyzing the matrix B lines, it follows that Summing up and making use of (51)-(52) we derive the following estimate which leads to the estimate A −1 B ∞ < 1 as we claimed at beginning. We will compare the two numerical solutions obtained by (10) and (11) with the following exact solution to (1) v e (t, with the forcing term In numerical tests we will consider a particular case of the nonlinear reactiondiffusion equation (1), namely, the Allen-Cahn equation ([1]), which means p 1 = α * ξ, p 2 = ξ and p 3 = 1 2ξ . The initial values v 0 (x j ), j = 1, 2, ..., N , were computed via Matlab function csapi(v0) -cubic spline interpolant to the given data: For beginning, taking T = 1, ω = 0.5, b = 1, α = 1.0e + 2, ξ = .5, N = 31, ε M = 0.1, M = T /ε M , the shape of the graphs ploted in Figures 1 shows the stability and accuracy of the numerical results obtained by implementing the conceptual algorithm Alg 1-IMBDF (see (45)), while, the errors v e − V N j ∞ produced by three methods analyzed in [32] (the Newton method, the linearized method and the old fractional steps method -(10)), as well as the new fractional steps method (11), are shown in Figure 2.

Conclusions.
• As a novelty of this work we refer firstly to the new scheme of fractional steps type, introduced by (11) in order to approximate the solution to the nonlinear reaction-diffusion problem (1). Corresponding, we have considered an IMEX scheme and we have proved stability result for the error equation,    (11) which shows us that the new fractional steps method depend linearly on the small parameter, like as in the methods studied in [32].
• Secondly, in the numerical experiments we focus our attention on a particular case of (1) -the Allen-Cahn equation, which serves as a mathematical model for many complex moving interface problems and in which the challenge in terms of numerical analysis is due to the thickness of the interface separating different phases. • Not least, let's remark from the graphical representation of the errors (see Figures 2-4), produced by those four methods analyzed, that conditions of consistency and stability are sustained by both theory and numerical experiment and that are significantly influenced by the parameters of time and space (see [29] for detailed discussions regarding the dependence on the physical parameters). In addition, let's remark that conditions of stability are sustained by both theory and numerical experiment (see Figure 1) and that are significantly dependent on all parameters (see Proposition 5.1).