The Regularized Mesh Scheme to Solve Quasilinear Parabolic Equation with Time-Fractional Derivative

A quasilinear parabolic problem with a time fractional derivative of the Caputo type and mixed boundary conditions is considered. The coefficients of the elliptic operator depend on the gradient of the solution, and this operator is uniformly monotone and Lipschitz-continuous. For this problem, unconditionally stable linear regularized semi-discrete scheme is constructed based on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L1$$\end{document}-approximation of the fractional time derivative. Stability estimates are obtained by the variational method. Accuracy estimates are given provided that the initial data and the solution to the differential problem are sufficiently smooth. The proved result of stability of the semi-discrete scheme is valid for the mesh scheme obtained from the semi-discrete problem using the finite element method in spatial variables.


INTRODUCTION
Partial differential equations with time-fractional derivatives arise in mathematical modeling of anomalous diffusion and dynamic processes in materials with memory. Numerous articles are devoted to the approximation of boundary value problems for linear problems with time-fractional derivatives. In recent years, the numerical methods has been actively developed for equations with time-fractional derivatives and the elliptic parts containing the non-linearities. Thus, implicit mesh schemes for the equations with Caputo time-fractional derivatives and with right-hand sides depending on the solution are studied in [1][2][3]. In [4] an implicit scheme is investigated for the nonlinear diffusion equation with a coefficient depending on the solution and with Caputo time-fractional derivative. In [5][6][7][8][9][10], analysis of ADI schemes was performed for linear time-fractional equations and for a class of semilinear timefractional equations.
The article [11] proposes a new mathematical model for a viscoelastic-plastic process by formulating a temporary fractional equation containing an elliptic operator, which depends on the gradient of the solution. In [12], the authors construct and study backward Euler and locally one dimensional schemes approximating Dirichlet problem with Caputo time-fractional derivative and a quasilinear elliptic part without mixed derivatives.
In this article, we consider a quasilinear parabolic problem with a time fractional derivative of the Caputo type and mixed boundary conditions. The quasilinear elliptic part of the problem defines a uniformly monotone and Lipschitz operator, which is energy equivalent to the Laplace operator. For this problem, we construct and investigate a linear regularized semidiscrete scheme in which the regularization term contains the Laplace operator and time-fractional derivative is approximated using well-known L1-approximation. Stability estimates for the semidiscrete problem by the variational method are obtained. Accuracy estimates are given provided that the initial data and the solution to the differential problem are sufficiently smooth. For discrete schemes obtained as a result of mesh approximations of the elliptic part of the semidiscrete problem, the stability estimates remain valid. They are also generalized for factorized schemes.

FORMULATION OF THE PROBLEM
Let Ω ⊂ R 2 be a polygon with the boundary Define a Caputo-type time-fractional derivative: with a kernel G(t) satisfying the assumptions: The conditions (2) are satisfied for See the definitions of the corresponding derivatives, e.g., in [13][14][15]). Let us define a second order quasilinear elliptic differential operator L and the co-normal derivative ∂u ∂ν L by the equalities where n is the unit vector of the outward normal to Σ N . Suppose that the nonlinear coefficients g i (t, p) = g i (x, t, p 0 , p 1 , p 2 ) of the equation satisfy the following assumptions for all (x, t) ∈Q and p, q, r ∈ R 3 : Consider a quasilinear parabolic problem with the time-fractional derivative and mixed boundary conditions: We assume that given functions f = f (x, t) and q = q(x, t) are continuous in t in order to use the simplest mesh approximations and do not specify their properties in details, since we do not investigate the existence of a solution to problem (5).
Let us make only a few remarks about well-known results on the existence of weak solutions of parabolic equations with time-fractional derivatives. The existence of a regular weak solution to Dirichlet problem for the case of linear operator L and classical Caputo fractional derivative was proved in [16]. A priori estimate for the solution in the norm of H α (0, T ; L 2 (Ω)) ∩ L 2 (0, T ; H 1 0 (Ω) ∩ H 2 (Ω)) through the L 2 (Q)-norm of f was derived. The existence of a weak solution from Besov space [12] for a quasilinear Dirichlet boundary value problem with classical Caputo derivative and a particular case of differential operator (3). This result can be generalized for the case of boundary value problem with elliptic operator (3) and Caputo time derivative. In a general case the problem of the existence of a unique solution to (5) seems to be open. In the future, we will assume that problem (5) has a smooth solution, and will focus on studying the mesh approximation of this problem.

APPROXIMATION
Multiplying the differential equation (5) by a sufficiently smooth function v(x) : v| Γ D = 0, we obtain a variational equation that can serve as a basis for the constructing its approximations: Due to the first inequality in (6) the norms ||u|| V and ||u|| 1 are equivalent. By virtue of this fact and the assumptions (4) the form is well-defined on V × V for all t and has the following properties: Let ω τ = {t j = jτ, j = 0, 1, . . . M; Mτ = T } be a uniform mesh on the segment [0, T ] and let us use the notation y j = y(t j ) for a continuous function y(t). The L1-approximation of the time derivative ∂y ∂s (s)ds of a continuous function y(t) at a point t k ∈ ω τ is given by the following relations: The properties (2) imply that the coefficients satisfy the inequalities First, we consider the semidiscrete implicit scheme approximating problem (5): Since ∂ tȳ with right-hand side Due to the properties (7) the left side of (10) defines a continuous and uniformly monotone operator from V to its dual V * . Using the inequalities (6), it is easy to see that linear functional . Thus, the existence of a unique solution to problem (10) follows from [17], Chapter 2, section 2.1. Below, instead of the nonlinear problem (9), a linear regularized problem is constructed. To shorten writing, we will use the notations Let us introduce a bilinear bounded form r(u, v) on the space V × V , which satisfies the inequality A particular case of this bilinear form is The regularized scheme approximating (5) is as follows: for a given y 0 = 0 find y k for k 1 from the system of linear equations For a fixed k the mesh scheme (12) can be written as The linear variational problem (13) has a unique solution for any σ 0 by Lax-Milgram theorem. The essential point used in proving stability estimates is the property of the L1-approximation of the time-fractional derivative, which is given by the following lemma: Then the symmetric matrix 1/2(D + D T ) is positive definite.
Using the matrix D, the approximation of time derivative at a point t k = kτ ∈ ω τ can be written as Introduce also the bilinear, symmetric and positive definite form on

Now the equation (13) can be rewritten as
For generality, we derive an a priori estimate for the case of a nonzero initial value. For i = 1, 2 we denote by y i the solutions of the equation (13) with the right-hand sides f i and the initial values y 0 i . Let also y = y 1 − y 2 , y 0 = y 0 1 − y 0 2 and ψ = f 1 − f 2 .

Lemma 2. Let the conditions (7) be satisfied and the regularization parameter satisfies the condition
Then there is a constant μ = μ(δ) such that for all m = 1, 2, . . . M the following a priori estimate for the regularized mesh scheme (13) holds: Proof. From the equations (24) written for y k 1 and y k 2 with test function v = 2(y k 1 − y k 2 ) = 2y k we obtain: In what follows, we use the following estimates: with a constant δ 0 ∈ (0, 2c 0 ) in the second estimate. Using the derived estimates and the inequality (14), after summing the equation (18) over k from 1 to m we obtain Now, due to the definition of the bilinear form b(., .) and the inequalities (11) and (16) we have the inequality With an appropriate choice of the constant δ 0 = δ 0 (δ) ∈ (0, 2c 0 ), we obtain 2 1 , then the last inequality implies (17) P The a priori estimate (17) can be used to obtain various stability estimates for the regularized mesh scheme (13). We give one such estimate in the case of zero initial condition. Let and let the following norms be defined for the functions from ω τ to V or to its subspace W : Theorem 1. Let the conditions (7) be satisfied, y i be the solutions of the equation (13) with the right-hand sides f i , homogeneous initial conditions y 0 i = 0 and the right-hand sides in Neuman boundary conditions q i . If the regularization parameter satisfies the condition (16): then the following stability estimate holds: Proof. The inequality (17) implies with For all k and v ∈ V and due to (6) Now the estimate (20) immediately follows from (21) and the last two inequalities. P Using the proven stability estimate, we can obtain an accuracy estimate under the assumption that the input data and the solution to the differential problem are sufficiently smooth. Let y k , k = 1, 2, . . . , M, be the solution of a semidiscrete problem (12) and u k = u(x, t k ) be the exact solution u(x, t) of the differential problem (5). Then the estimate (18) applied to the difference y − u contains in the right-hand side the function ψ k which is defined by the equality Let us consider the case of classical Caputo derivative. The following estimate is well-known ( [19,20]): Thus, the most significant term in the approximation error is associated with the introduction of the regularizing term, and the accuracy estimate is Note that for the implicit scheme (backward Euler scheme) and locally one-dimensional scheme (see [12]), the terms in the accuracy estimate related to the approximation in time have the following orders respectively: Both of these schemes are nonlinear and require the use of iterative methods when implemented at each time level. Remark 1. A fully discrete approximation of problem (5) can be constructed using an approximation with respect to the spatial variables of the semidiscrete problem (12). All the results on the existence of a unique solution and the proved stability estimates remain valid for mesh problems constructed using the finite element method with finite elements of the first order, as well as for finite element schemes with quadrature formulas. Remark 2. The explicit (forward Euler) mesh scheme is a particular case of fully discrete regularized scheme with σ = 0. A sufficient condition of its stability follows from the estimate (19): The last inequality is equivalent to the following estimate for the maximum eigenvalue λ max of the mesh Laplace operator: As is known, , where h characterizes the mesh step in the spatial variable. Since, for example, for the L1-approximation of the Caputo derivative, d 1 − d 2 = O(τ −α ), then (22) leads to a very restrictive condition for the smallness of the time step: τ < ch 2/α with a constant c, which depends on the problem data.

Remark 3.
Let r 1 (u, v) be a bounded bilinear form in V × V and r 1 (u, u) 0. Then the following scheme satisfies all the proven results on unique solvability and stability: This fact allows to apply these results to the factorized mesh schemes. More precisely, let for a fully discrete scheme constructed using finite element or finite difference scheme r(u, v) = (Ru, v) and the mesh Laplace operator R = −Δ h with corresponding boundary conditions can be splitted into the sum R = R 1 + R 2 such that R 1 R 2 is a non-negative mesh operator: (R 1 R 2 u, u) 0. In this case we set r 1 (u, v) = (R 1 R 2 u, v) and take the regularization term in the form y k , v). The corresponding mesh scheme becomes the mesh scheme with a factorized operator on each time level: The implementation of (24) is reduced to sequential inversion of the mesh operators I + σd −1 1 R 1 and I + σd −1 1 R 2 and is very easy if these operators are locally one-dimensional (that corresponds to alternating direction method) or triangle (that corresponds to alternating triangular method). The term σ 2 d −1 1 r 1 (u k , v) is added to the approximation error. In case of classical Caputo derivative it is asymptotically O(τ α ). So, the accuracy estimate is ||y − u|| L 2 (ωτ ;H 1 ) = O(τ α ).

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.