Some Aspects of Numerical Analysis for a Model Nonlinear Fractional Variable Order Equation

: The article proposes a nonlocal explicit ﬁnite-difference scheme for the numerical solution of a nonlinear, ordinary differential equation with a derivative of a fractional variable order of the Gerasimov–Caputo type. The questions of approximation, convergence, and stability of this scheme are studied. It is shown that the nonlocal ﬁnite-difference scheme is conditionally stable and converges to the ﬁrst order. Using the fractional Riccati equation as an example, the computational accuracy of the numerical method is analyzed. It is shown that with an increase in the nodes of the computational grid, the order of computational accuracy tends to unity, i.e., to the theoretical value of the order of accuracy.

Fractional calculus studies fractional derivatives and integrals, as well as their applications. It is thanks to the applications that fractional calculus is intensively used in physics [7,8], mechanics [9][10][11], biology [12,13], economics [14,15] and other sciences. We will not dwell here on the interpretation of the meaning of the fractional derivative. We refer the reader to other works, such as [16], for discussions on this problem.
In the 1990s, within the framework of fractional calculus, there appeared a direction devoted to fractional operators of variable order (FO-VO) [17]. This direction received a powerful impetus for development, which is associated with the emergence of various complex problems in physics, mechanics, biology, economics, and other sciences, which are well described by FO-VO. For example, it was shown in [12] that the kinetics of the reaction of proteins demonstrates relaxation mechanisms that are properly described by the temperature-dependent fractional-order operator. The survey article [18] described the development of fractional calculus as fractional operators of constant and variable orders, and considered analytical and numerical methods for their study.
Of great practical importance is the application of numerical methods of analysis to the solution of nonlinear equations with derivatives of fractional orders [19][20][21][22][23][24][25][26][27]. This is due to the fact that numerical methods provide a quantitative description of the parameters of the process under study. However, as a rule, questions arise here related to the stability and convergence of numerical methods. In this paper, we investigate these issues, using an explicit finite-difference scheme in order to find a numerical solution to a nonlinear equation with a fractional variable order derivative. This equation can be a model and describe various processes, for example, processes with saturation.
In the paper, the necessary concepts and definitions are first given; then, the formulation of the problem is revealed, the method of solution is given, the issues of approximation, stability and convergence of the scheme are considered, and specific examples are given.

Preliminary Information
Definition 1. The Euler gamma function is determined according to the following integral: (1)

Remark 1.
It should be noted that the Euler gamma function is monotonically decreasing over the interval 0 < x < 1.

Definition 2.
A continuous function f (u, t) ∈ C[0, T] satisfies the Lipschitz condition with a constant L in the variable u(t) ∈ C[0, T]: Definition 3. The operator of the fractional variable order 0 < α(t) < 1 acts on a function u(t) ∈ C[0, T]: where Γ(·) is the Euler gamma function (1) andu(t) = du dt will be called the derivative of a fractional variable order of the Gerasimov-Caputo type.
The properties of operator (3) can be found in [28].

Definition 4.
Equations containing derivatives of a fractional variable order of the Gerasimov-Caputo type (3) will be called fractional equations.

Remark 2.
In the case when operator (3) has a constant order, i.e., α does not depend on t, we obtain the Gerasimov-Caputo operator [29,30].

Problem Statement and Solution Method
Consider the following Cauchy problem for a nonlinear fractional equation: where u(t) ∈ C[0, T] is the solution function, t ∈ [0, T] is the time, T is the model time, u 0 is a given constant, 0 < b(t) is a continuous function, f (u, t) is a nonlinear function satisfying the Lipschitz condition (2) with a constant L in the variable u(t), and the operator of a fractional variable order has view (1).

Remark 3.
The Cauchy problem (1) describes a wide class of dynamic processes with variable memory [11].
Due to the nonlinearity of the Cauchy problem (4), we will seek its solution using the numerical method of finite-difference schemes [31][32][33]. Consider a uniform mesh.
To do this, we divide the segment [0, T] into N equal parts-grid nodes with a step τ = T/N. We assume that the function u(t) satisfies the necessary smoothness for constructing a finite-difference scheme, i.e., u(t) ∈ C 2 [0, T]. Then, the solution function u(t) goes over to the grid solution function u(t k ) or u k , where t k = kτ, k = 0, . . . , N − 1.
An approximation of the derivative of a fractional variable order of the Gerasimov-Caputo type in Equation (4) can be written as follows:
Substituting (3) into (4), we obtain a discrete analogue of the Cauchy problem: where is a known constant The following lemma is true for the Cauchy problem (6).

Lemma 1.
The coefficients A k and w k j of the discrete Cauchy problem (6) for any fixed k have the following properties: Proof. The first property follows from the property of the Euler gamma function for any fixed k: 1 Γ(2 − α k ) > 0, as well as τ > 0, then A k > 0. The second property of the weighting coefficients w k j follows from the disclosure of the following sum: We prove the third property of the weight coefficients w k j as follows. Consider the following function: Therefore, the function η(z) is monotonically decreasing and the weight coefficients w k j have a property of 3. Let us now investigate the order of approximation of the fractional operator ∂ 0t u(t). Then, the following lemma is true.

Lemma 2. An approximation∂
0t u(t) satisfies the following estimate: Proof. The proof follows from Definition 3, approximation (3), and the second property of Lemma 1. Note that the approximationū of the first-order derivative has the form The estimate (7) follows.

Lemma 3.
The discrete Cauchy problem (6) approximates the original differential problem (4) to the first order, that is, has an error as follows: Proof. Let us compose the difference between the Cauchy differential problem (4) and its discrete analogue (6), we obtain the following: Passing to the absolute value and taking into account estimates (2) and (7), we obtain the following: b t j u t j − u j + L u t j − u j ≤ Cτ, Next, we pass to the following norm: whence the estimate (8) follows.
We write the Cauchy problem (6) in the form of a nonlocal explicit finite-difference scheme: Consider the issues of convergence and stability.
We substitute r = k − 1 into relation (15) and, taking into account estimate (12), we obtain the following: e k+1 ≤ s k · s k−1 · . . . · s 1 e 1 + s(s k s k−1 · . . . · s 2 + . . . Therefore, if for any index condition: b k A k ≤ 1 is satisfied, then the numerical solution (6) converges to the exact one to the first order. It is easy to obtain estimate (10) from this condition. The theorem is proved.
Consider the issues of sustainability. Let U k , W k a be two different solutions of the matrix equation (12) with initial conditions U 0 , W 0 . Then, the theorem on the stability of scheme (11) is valid.

Theorem 2.
A nonlocal explicit finite-difference scheme (11) is conditionally stable if condition (10) is satisfied and the estimate |U k − W k | ≤ C|U 0 − W 0 | is valid for any k, where C > 0 is a constant independent of the step τ.
Proof of Theorem 2. Let us introduce the notation e k+1 = U k+1 − W k . Then, Equation (11) can be written in the following form: e k+1 = Me k + F e,k . Then, the estimate is valid: A k e k = s k e k ≤ s k s k−1 e k−1 ≤ s k s k−1 s k−2 e k−2 ≤ . . . ≤ s k s k−1 · . . . · s k−r e k−r .
For these values of the parameters, condition (10) of Theorem 1 is satisfied. Indeed, in this example, A = 9.505 and b/A = 0.07364 < 1. Then, according to the explicit finite-difference scheme (17), we obtain the calculated curve in Figure 1. From Figure 1, it can be seen that the calculated curve has an s-shape, which is characteristic of dynamic processes in saturated media. With a decrease in the values of the parameter a, the s-shape becomes more pronounced.
Let us estimate the computational accuracy of scheme (17). To do this, we calculate the maximum error ξ by the double recalculation method, i.e., ξ = max i (|u i − u 2i |), where u i and u 2i are the calculated values obtained by Formula (17) at step τand τ/2, respectively. The results are shown in Table 1.  Table 1, it is seen that with an increase in the nodes of the computational grid, the maximum error decreases, and the computational accuracy tends to unity, which corresponds to the condition of Lemma 3.
For this example, consider the case when condition (10) is violated. To do this, it is enough to take the values of the parameter by an order of magnitude, for example, b = 18, b/A = 1.8937 > 1, and we leave the values of other parameters unchanged. The calculated curve obtained by Formula (17) is shown in Figure 2. It should be noted that for this example, condition (10) is satisfied. The simulation results using Formula (17) give the calculated curve in Figure 3.  Figure 3 shows that the shape of the calculated curve has become noticeably more complicated in comparison with the calculated curve in Figure 1. This is primarily due to changes in the parameters of the model Equation (16) with time t.
However, despite the more complex form of the calculated curve, we can use it to describe cyclical dynamics, for example, in economics, to describe cycles and crises.
An estimate of the computational accuracy is given in Table 2.  Table 2, we see that with an increase in the nodes of the computational grid, the error decreases, and the computational accuracy tends to unity, which corresponds to the results of Lemma 3.

Discussion
In the case of the stiff Cauchy problem (4), the method may not work correctly. For example, this can be seen in Figure 2, the case when the b parameter has a sufficiently large value (b = 18). Here, the stability condition (10) is violated. Therefore, the question of the condition for the appearance of the stiff Cauchy problem (4) deserves a separate study.
Another very interesting question also deserves attention: how does the variable order of the fractional derivative affect the stability of the explicit finite-difference scheme (6)? It should be noted that the variable order α(t) is included in estimate (10) for the stability of the method. Here, it is worthwhile to study the stability of the proposed scheme (6) depending on the behavior of the function α(t), for example, monotonicity or its periodicity, etc. Additionally, the function α(t) can have a singularity at zero.

Conclusions
An explicit finite-difference scheme for a fractional equation with a derivative of the variable fractional order of the Gerasimov-Caputo type was investigated. The questions of convergence and stability were investigated. It was shown that the scheme is conditionally stable and converges to the first order. The results were confirmed by examples and estimates of computational accuracy.
It can be concluded that an explicit finite-difference scheme can be used to calculate dynamic processes in environments with saturation. Further continuation of the work lies is the practical application of an explicit finite-difference scheme for calculating the curves of the volumetric activity of radon in a storage chamber by analogy with work [36].

Funding:
The work was carried out within the framework of the research project of the Vitus Bering Kamchatka State University Natural disasters of Kamchatka-earthquakes and volcanic eruptions (monitoring forecast, study, psychological support of the population) No. AAAA-A19-119072290002-9 and subject research IKIR FEB RAS "Physical processes in the system of near space and geospheres under solar and lithospheric influences" No. AAAA-A21-121011290003-0.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results