Solving the Boundary Value Problems for Differential Equations with Fractional Derivatives by the Method of Separation of Variables

: This paper is devoted to solving boundary value problems for differential equations with fractional derivatives by the Fourier method. The necessary information is given (in particular, theorems on the completeness of the eigenfunctions and associated functions, multiplicity of eigenvalues, and questions of the localization of root functions and eigenvalues are discussed) from the spectral theory of non-self-adjoint operators generated by differential equations with fractional derivatives and boundary conditions of the Sturm–Liouville type, obtained by the author during implementation of the method of separation of variables (Fourier). Solutions of boundary value problems for a fractional diffusion equation and wave equation with a fractional derivative are presented with respect to a spatial variable.


Boundary Value Problems for Differential Equations of Fractional Order
The paper is devoted to the method of separation of variables (the Fourier method). This method, which is so widely used in solving boundary value problems for partial differential equations of integer order, until recently remained unsuitable for solving boundary value problems for differential equations with fractional derivatives. The main reason, of course, is that the spectral theory of non-self-adjoint operators generated by the corresponding differential expressions of fractional order and boundary conditions of the Sturm-Liouville type has been supplemented with the necessary information quite recently.
Almost all of the author's papers, to varying degrees, are devoted to the study of the spectral structure of these operators, which constitute the theoretical basis of the method of separation of variables, and the author of this paper, in a sense, summarizes his work in this area.
So, in this paper, a method of separation of variables is presented for solving boundary value problems for differential equations with fractional derivatives of the form First of all, we note that anomalous diffusion or dispersion we can describe using the fractional space derivatives, and some processes with 'memory' effects-using the fractional time derivatives.
It should be noted that depending on the modeled process, the fractional differentiation operators appearing in these equations can be both the Riemann-Liouville fractional differentiation operators and the fractional differentiation operators in the Caputo sense. One of the most important problems in modeling physical processes using differential equations with fractional derivatives is the problem of establishing in what sense the fractional derivative is taken and the identification of the order of this fractional derivative.
Undoubtedly, the most significant, fundamental point in the study of boundary value problems for these equations by the method of separation of variables is the question of completeness of systems of eigenfunctions of boundary value problems for the equations (these equations arise when the variables are separated in Equations (2) and (3)). Therefore, we present basic results from the spectral theory of operators generated by differential equations of the form (2) and boundary conditions of the Sturm-Liouville type.
The relationship between eigenvalues and zeros of a Mittag-Leffler function is shown. The Green's functions of boundary value problems for equations of the form (2) are considered in detail (it should be noted that these Green's functions were first obtained by the author in his post-graduate student paper [2]), the study of which made it possible to approach problems of the distribution of zeros of a function of the Mittag type from completely new positions-Leffler and reveal the deeply hidden properties of these functions, which for many years have not been possible for specialists in the theory of functions. First of all, we note that the asymptotic properties of the Mittag-Leffler function have been sufficiently well studied [1], but the study of the non-asymptotic properties of the zeros of the Mittag-Leffler function or, similarly, the eigenvalues of operators generated by boundary value problems for Equation (2), is conjugate with large analytical difficulties (in particular, M. M. Dzhrbashian wrote in [1] that "the question about the completeness of the eigenfunctions of boundary value problems for Equation (2) or a finer question about whether these systems compose a basis in L 2 (0,1) has a certain interest. However, their solution is apparently associated with significant analytic difficulties".). Therefore, the author gives these properties in sufficient detail.

Boundary Value Problems for the Fractional Order Diffusion Equation
In this section we present the necessary information from the spectral theory of operators generated by differential equations of fractional order and boundary conditions of the Sturm-Liouville type.
2.1. Spectral Analysis of Operators, Generated by Fractional Differential Equations of Order More than 1 but Less than 2 and Boundary Conditions of Sturm-Liouville Type and on One Method for Identifying the Order of the Fractional Derivative We devote this subsection to the spectral analysis of two boundary value problems [3,4] L(u; 1, 1 − α, These problems are the focus of many researchers. First, we note that in [4] (and references therein), the following problem was considered with studying of the spectrum of the operator The operator D (β) arose great interest after F. Mainardi's paper [3]. In this paper, the following equation was considered where ω is a positive constant and 1 < γ < 2, which Mainardi called a fractional oscillatory equation. This paper has been, without exaggeration, very interesting for a lot of researchers. First of all, note that: 1. If λ = 0, then any solution u(x) ∈ S 2 [0, 1] (where S 2 [0, 1] is the class of summable (integrable) on [0, 1] functions u(x) including their derivatives of first and second order) for the equation coincides with the solution for Equation (4); 2. Equations (4) and (7) are equal if Of course, the fractional oscillatory equation, or equation for fractional oscillator (as an equation, which describes an oscillatory physical system), will have at least the main oscillatory properties.
Hereafter, the following integral equations will play the main role: is the Green's function of the problem (4)- (5), which was constructed in [2] and G 1 (x, ζ) is the Green function of the problem which was considered in [5] for the first time (see also references therein). Important note: the operators L(u; 1, 1 − α, 1, 0) and L(u; 1, 1, 1 − α, 0) have the same orders, but the γ 0 , γ 1 , γ 2 , of those orders are different.
It is easy to show [6] that G 0 (x, ζ) is not with a fixed sign, and this fact says that Equation (7) was incorrectly chosen as the oscillatory equation. Physically, it is clear that the order of the operator L(u; γ 0 , γ 1 , γ 2 , q(x)) is close to 2 (or when γ 0 + γ 1 + γ 2 − 1 is close to 2), then the operator L(u; γ 0 , γ 1 , γ 2 , q(x)) has the main oscillatory properties. We have the following result.
then the first eigenvalue of the problem (4)-(5) is positive and simple (the multiplicity of this eigenvalue is equal to 1), and basic (main) tone has no nodes (i.e., the first eigenfunction corresponding to the first eigenvalue does not vanish in (0, 1)).
Next, we consider in detail the function G 1 (x, ζ). As it was shown in [7], this function has many useful properties, in particular G 1 (x, ζ) = G 1 (1 − ζ, 1 − x) and G 1 (x, ζ) > 0, for any x, t ∈ (0, 1) (i.e., this Green's function is a persymmetric function). Part of the results of the theorem below follow from the well-known Perron's theorem.
Theorem 2. The first eigenvalue λ 1 of the problem (8)-(9) is real, simple, and satisfies the condition and the basic (main) tone has no nodes for all α ∈ (0, 1).

Proof.
As it was written above, from Perron's theorem, it follows that the first eigenvalue is real and simple, and the basic (main) tone has no nodes. Let us show that As it was mentioned above, the problem (8)-(9) is equivalent to the integral equation of Fredholm (II kind) and the value λ is an eigenvalue of the problem (8)- (9) if and only if it is a zero of the Mittag-Leffler function E 1 1+α (−λ, 1 +α) [5,7]. Let us rewrite the operator As operators A 0 and A 1 are trace class operators [10], then Since A 0 is Volterra's operator, then sp(A 0 ) = 0, and so It is easy to find the trace of the operator A 1 (A 1 it is one-dimensional operator). Let us consider the equation The Fredholm determinant of this equation is From this we obtain Thus Since the kernel of the operator −A is non-negative, then λ 1 is a positive number, and Theorem 2 is proved.

On Completeness of System of Eigenfunctions and Associated Functions of Operator, Generated by Model Fractional Differential Equation and Boundary Conditions of Sturm-Liouville Type
Let us start from the equation where At first, Equation (10) was studied in [5] as a model equation of the fractional order 1 < σ 2 < 2. In particular, it was established in [5] that the two-point Dirichlet problem for Equation (10) with q(x) = 0 is equivalent to the integral equation We have: Then the system of eigenfunctions and associated functions of the problem (10)- (11) is complete in L 2 (0, 1).
A close result (for a semibounded potential q(x)) was obtained in [11]. It should be noted that the proof of these statements are based on the fact that the operator, generated by the problem (10)-(11), is sectorial [12]. Theorem 4. All eigenvalues of the problem (10)-(11) for q(x) ≡ 0 are in the angle |argz| < π(1−γ) where α = 1 − γ and 0 J α x is the operator of fractional integration of order α: Since the number λ is an eigenvalue of the problem (10)- (11) if λ is a zero of the function E 1/µ (−λ; µ) (µ = 1 + γ) [5], the following proposition is valid.

Methods of the Theory of Perturbations in Fractional Calculus and the Questions of Localization and Multiplicity of Eigenvalues
To prove that the studied operator does not have the associated functions, we present the main points of the method presented in [14].
In L 2 (0, 1) we consider the operator which was for the first time studied in [7]. Here, 0 < ρ < 2, and is the Green function of the following problem S (for λ = 0): In this case [7], if γ 0 = γ 1 = · · · = γ n = 1 then the problem S takes the form The last function was studied very well, and we will use it in the sequel. The operator A ρ was investigated in [5,7,15]. Let us study this operator carefully, because it turns out that the Mainardi equation [3] (fractional oscillatory equation) does not have many basic oscillatory properties. The search for such a differential equation that has these properties led us to study the operator A ρ . Now, let us introduce the most significant properties of this operator established by the author earlier: 1. For ρ > 1 the operator A ρ is completely nonself-adjoint [5,7,15]; 2. For ρ ≤ 1 the operator A ρ is sectorial [6] (and see the references therein); 3. For 0 < ρ < 2 the system of eigenfunctions of the operator A ρ is complete in L 2 (0, 1) [5,16]; Now, let us study integral operators corresponding to boundary value problems for fractional differential equations using methods of the theory of perturbations.
The holomorphic dependence of these operators on the order of fractional differentiation is proved. There are several useful criteria for holomorphy. In accordance with this, various types of holomorphic families are considered. We will use type (A). Type (A) is defined in terms of the boundedness of the perturbation with respect to the unperturbed operator.
Let us formulate a very important criterion, which we will use later [17].
Theorem 6. (Criterion of holomorphy (A)). Let T be a closable operator from X in Y, and let T (n) , n = 1, 2, ..., be operators from X in Y, of which the domains of definition contain D(T) = D. Assume that there exists constants a, b, c ≥ 0, such that Then for |κ| < 1/c the series then the operator T(κ) is closable, and the closures T(κ) form a holomorphic family of type (A) [7].
We shall note that the holomorphic families of this type and, in particular bounded-holomorphic families, were studied since Rellich's papers [18] (and references therein). A wide list of references is presented in papers of M.K. Gavurin and V.B. Loginov [18] (and references therein). Theorem 7. If |ε| < 1, then the operator forms a holomorphic family of type (A), i.e., the unperturbed operator, and Theorem 8. If |ε| < 3/2, then the operator forms a holomorphic family of type (A) where is the unperturbed operator, and Since [15] the Fredholm spectrum of the operators under study coincides with the zeros of the appropriate function of the Mittag-Leffler type, the presented method allows to efficiently study the problem of distribution of zeros for functions of Mittag-Leffler type. To confirm this assertion, we give two examples. Following [8], we introduce the following notation: λ n (α) are the eigenvalues of the problem (4)-(5). In [8], it was written that "... in the limiting case α = 0 the problem (4)-(5) becomes the Sturm-Liouville boundary value problem with the sequence of eigenvalues λ n (α) = (πn) 2 . Is it true that lim α→0+ λ n (α) = (πn) 2 for any fixed n? The answer will be positive." Let us prove a stronger proposition.
Proof. Theorem 8 is a trivial corollary of Theorem 4.2 (see [10], p. 35) and the fact that the operator function B(ε) is strongly continuous for |ε| < 1.
Finally, we consider one more significant question of the multiplicity of eigenvalues of the operator B(ε) (as was mentioned above, this question is related to the question of the multiplicity of zeros of a corresponding function of the Mittag-Leffler type [7]).
It is known ( [1], theorem of Dzhrbashian-Nersesian) that all zeros of a function of the Mittag-Leffler type E ρ (z, µ) (where ρ > 1/2, ρ = 1, Im(µ) = 0) that are sufficiently large in the modulus are simple. Therefore, we mainly pay attention to the multiplicity of the first eigenvalues of the operator B(ε). The following theorem holds [14] Theorem 10. Let |ε| < 32π 2 9 + 2 3 −1 . Then the first eigenvalue λ 1 (ε) of the operator B(ε) is simple [7]. Proof. It is known [7] that if the spectrum of the operator B(0) is divided into two parts by a closed curve Γ, then the spectrum of the operator B(ε) is also divided by the curve Γ for sufficiently small ε.
In this case, the estimate of the smallness of ε is as follows [7]: (where a, b, and c are parameters that enter inequality (12)). As the contour Γ in Formula (13), we take the circumference |ζ − 1 π 2 | = ρ 2 , where ρ is the distance from 1 π 2 to the set of the rest eigenvalues of the operator B(0). The parameters a, b and c are already calculated [7]. Theorem 9 is proved.
We note that it can be shown in the same way that the second eigenvalue of the operator B(ε) is simple too. It is the principal point that this method gives the possibility to include the study of nonselfadjoint operators of the form A  In the last few years, fractional integro-differentiation has been the focus of many researchers of science and engineering [19,20] (and see references therein). We can describe anomalous diffusion or dispersion using the fractional space derivatives, and some processes with 'memory' effects using the fractional time derivatives. In this paragraph, we solve the problem of finding the radon flux density [14] by its concentration at different depths of the earth's surface by the method of approximate solution of the first boundary value problem for the fractional differential equation of advection-diffusion [21].
It is known [21,22] that the problem of finding the radon flux density by its concentration at different depths of the earth's surface is set as follows: to find a solution to the boundary value problem where D α 0+ u(x, t)-the Riemann-Liouville fractional derivative of the order α, with boundary conditions u(0, t) = u(1, t) = 0, Using the method of separation of variables [21], we can write out the solution to this problem here ..., λ −3 , λ −2 , λ −1 , λ 1 , λ 2 , λ 3 , ...
-Zeros of the function E α,α (λ), arranged in the appropriate order according to [21], and δ n = {δ(x), z n (x)} L 2 (0,1) , n = 1, 2, ...  For an approximate solution of this problem, one can use the formula In [23], using this formula, the problem of finding the radon flux density by its concentration at different depths of the earth's surface was solved. It was shown that rather well approximates the exact solution u(x, t) and also there are algorithms for finding the eigenvalues λ n and the Fourier coefficients δ n . Remark 1. Note the paper [24], where the same problem is solved by numerical methods.

Method of Separation of Variables for Time-Space Fractional Vibration Equations-The Basic Theory
In this paragraph we present the necessary information from the spectral theory of operators generated by differential equations of the second order with fractional derivatives in the lowest terms with boundary conditions of the Sturm-Liouville type.
Many problems of mathematical physics [25][26][27] associated with perturbations of normal operators with discrete spectrum lead to the consideration in Hilbert space H of the compact operator A = (I + S)H, called a weak perturbation H (for a compact S) or as the operator of Keldysh type (the information about of such operators and last investigations in this field were published in our brief [28]).
In [28] the basis property of the system of root vectors and localization of root vectors and eigenvalues for the investigated operators were established (we shall also note paper [29] of one of our co-authors E. Larionov, in which the spectral theory of the operators of the Keldysh type is very strongly developed).
In the present paper, we consider the operator of Keldysh type B, generated by the differential expression and the boundary conditions of the Sturm-Liouville type Note, that for 0 < α < 2, the spectral structure of operatorB generated by problem was considered in detail in our paper [12]. In particular, the following theorem was shown: Theorem 11. If |ε| < 10 20 , then all eigenvalues of operatorB are simple and real.
From this theorem, it follows that the operatorB generated no associated functions.
Theorem 12. The number λ is an eigenvalue of problem (19)- (20) if λ is a zero of the function The eigenfunctions of problem (A) take the form where λ i are zero of the function ω(λ). In [12] it was proved that the system of eigenfunctions (22) is complete in L 2 (0, 1). However, this system is not orthogonal. Therefore, in paper [12] they were considered together with problems (19)- (20), the problem conjugate to it. Let us consider operator B, generated by the differential Equation (19) and boundary conditions (20). The following theorem holds [14] Theorem 13. Let 0 < α < 2, then, the system of eigenfunctions of operator B is complete in L 2 (0, 1).

1.
Next, let us denote n(r, B) the exact number of characteristic values of the operator B lying in circle |λ| ≤ r. The problem of allocation of characteristic values of the operator B formulates an investigation of asymptotic properties of n(r, B) for r → ∞. In [25], this problem was solved when the order of fractional derivative D α 0x was less than 1. In [25], the study of the function n(r, B) was reduced to the one of the spectra for the linear beam operator L(λ) = I + M − λN.
Since Obviously, in our case as in [25] we may take the function √ r as ϕ(r). Any linearized mechanical system in which there is energy dissipation is described by a linear operator A, densely defined in H, with values of the form (A f , f ) in the left half-plane: In quantum mechanics, energy dissipation is characterized by the fact that the form of the linear operator describing the physical system lies in the upper half-plane, i.e., For definiteness, when speaking of dissipative operators, we shall have in mind the operators of the latter type; dissipative operators of quantum mechanics.

2.
Since (see [14] and references therein) any linearized mechanical system, which has an energy dissipation described by a linear operatorÃ, is densely defined in a Hilbert space H with values of the form (Ã f , f ) in the left half-plane As the operatorB describes oscillations of the mechanical system, then it should be dissipative (see [30] and references therein).
In this paragraph we show that the operatorB is dissipative. First, we shall note papers of F. Tricomi, Matsaev and Palant [25] (and references therein) (where it was shown, that values of the form (I α f , f ) are lying in the angle | arg λ| ≤ απ 2 , here I α -is fractional integral in the Riemann-Liouville sense of order α) and papers of authors [25] (and references therein), where it was established and Theorem 14. If 0 < α < 1 and ε > 0, then the operator, generated by the problem Proof. This theorem follows from the relation (23) and the fact that the operator Theorem 15. If 1 < α < 2 and ε < 0, then the operator, generated by the problem Proof. The scheme of the proof of Theorem 14 is the same as the one of Theorem 13.

Remark 2.
Let D = {0 < x < 1, 0 < t < 1}, and consider the first boundary value problem for equations of vibration of a string with a fractional derivative of order α with respect to partial variable If C 0 is negative, the numerical methods can work well, and the related numerical theoretical results can also be well established. However, when C 0 is positive, our numerical methods can not compute well, and the convergence and stability of the proposed methods can not be proved. This is a consequence of the fact that the term C 0 D β describes dissipation, as it was established by the Theorems 13 and 14 (in the case when C 0 is negative, the physical sense of the term C 0 D β is incomprehensible). Thus, proved Theorems 13 and 14 allow to correctly formulate the boundary value problem for Equation (25).
3. Let us consider operator A, generated by the problem where 0 < α < 2.
Finally, let us show that operator A is oscillatory (if the operator describes the oscillation motions, then it should have a whole complex of the oscillatory properties).
It is known that [25] (and references therein) if 0 ≤ ε ≤ 1 3 , and 1 < α < 2, then the Green function of the problem (29)-(30) is of fixed sign (we shall note, that Green's function of problem (29)-(30) was firstly constructed by one of the authors in his paper [25] (and references therein)). Unfortunately, this very important property of Green's function is possible to get only for a small enough ε. This is primarily due to the fact that Green's function G 2 (x, τ) [4] of the problem (29)-(30), for 1 < α < 2, has the following complex structure For 0 < α < 1, |ε| < 1/4, Green's function of the problem (29)-(30) was constructed in [25] (and references therein). Let us show how this function was constructed. Since the problem (29)-(30), for 0 < α < 1, is equivalent to the equation We can show that Thus, Since, for |ε| < 1/4 the kernel k(x, t) of the operator K satisfies the condition |k(x, t)| < 2, we have that the Green's function of problem (29)-(30) is of fixed-sign for 0 < α < 1 too.
Note the paper [29], which is very important in our opinion, which contains the proof of the basis property for the eigenfunctions.

Parametric Identification for Time-Space Fractional Vibration Equations
In numerous publications of the last decades, the problem of identifying the parameters of fractional models is mainly solved at a theoretical level, for example, by methods of spectral analysis. In our paper [31] (and see references therein), the model parameters are determined based on several characteristic points obtained in the experiment, by substituting the deformation values in the analytical solutions of the corresponding problem. We will use the same technique in what follows to identify the order of the fractional derivative in model (1).

The Bagley-Torvik Equation and the Laplace Transform
We consider the problem where D α u(x) is a fractional derivative of the order α. When 1 < α < 2, by the Riemann-Liouville definition, this problem is presented as follows [32]: Equation (31) was proposed in papers [33,34] (and see references therein) for modeling the oscillatory properties of various viscoelastic materials (polymers, glass, etc.). We shall note one recent paper [35] (and see references therein) where this scheme is used to model changes of the stress-strain characteristics of polymer concrete when subjected to loadings. In these papers, the polymer concrete samples based on polyester resin (dian-and dihloangidrid-1,1-dichloro-2,2-diethylene) were studied. Polymer concrete is represented as the set of granules of a mineral filler that is located in a viscoelastic medium. In this case, the motion of a granule is described by Equation (31), where λ is the rigidity modulus of resin, α is the viscoelasticity parameter of the medium and c is the viscosity modulus of resin.
Note that physical systems modeled by Equation (31) are very sensitive to changes in the order of fractional damping and it lead us to the very important task of the parametric identification of this value. We shall note that the problem of the parametric identification [36] (see and the references therein) remains poorly understood. The paper [37] (see and the references therein) is devoted to solving this important problem. Let us briefly give the technique presented in this paper.
We integrate Equaiton (31) from 0 to x and transform the resulting expression. We have Obtained expression (32) we integrate from 0 to x again: We solve the latest Equation (33) using the Laplace transform. Let us designate by U(s) the image of the function u(x); i.e., U(s) = u(t) or [38] (which is the same),

It is clear that the function
represents the convolution of the functions u(x) and x 1−α . For the convolution of functions, there exists the simple formula of images where F 1 (s) and F 1 (s) are the images of the functions f 1 (x) and f 2 (x), respectively.
It is clear that After some clearly transformations we obtain From this we obtain the formula for the image Formula (36) makes it possible to express the solution of problem (31) using Laplace's integral

Numerical Construction of the Solution
The obtained Formula (37) allows to numerically construct the graphs of solutions. The calculations are performed using the package Mathcad 14. Figure 1 presents the graphs of solutions for various values of parameter α. As in [35] we take c = 1.8 and λ = 93. Note that used parameter values were obtained in the experiments on the samples of polymer concrete [35]. The presented numerical check proves the validity of the limit behavior of the solution, which transforms into harmonic oscillations for values of α that are close to 2. To proof the correctness of the formulation of the problem of the parametric identification, it is necessary to investigate the solution's stability to the inaccuracy of the parameter α (α has an arbitrary value from (0,1)). For this purpose, in the neighborhood of the point α, consider the relative increment  To prove the correctness of the formulation of the problem of the parametric identification, it is necessary to investigate the solution's stability to the inaccuracy of the parameter α (α has an arbitrary value from (0,1)). For this purpose, in the neighborhood of the point α, consider the relative increment of this parameter by δ (i.e., α = α(1 + δ)) and determine the deviation function ρ(α, δ) with respect to the norm in the space L 1 (the of the summable functions) where u(x, α) is the solution of (31) with the parameter α. Here, the function ε(α, δ) = ∂ρ ∂δ determines the sensitivity of the solution to a possible error of the parameter α. The values of the function ε(α, δ) for various values of parameter α and the levels, 0.1, and 0.15 are found numerically; they are presented graphically in Figure 2.

284
It is known that [50] the solution of the Cauchy's problem (31) is determined by the following formula Comparing the graphs of the solutions numerically obtained by formulas (37) and (39), we can draw a 285 conclusion about their identity (Fig. 3). The obtained values of the function ε(α, δ) show the growth of the sensitivity with increasing of parameter α. The maximum value of the function ε(α, δ) does not surpass 0.2; this allows us to draw a conclusion about the stability of the solutions of problem (31) in relation to a small error of parameter α and the correctness of the problem of identifying this parameter.
It is known that [31] the solution of Cauchy's problem (31) is determined by the following formula Comparing the graphs of the solutions numerically obtained by Formulas (37) and (39), we can draw a conclusion about their identity ( Figure 3).

Parametric Identification of the Model by the Experimental Data
The following technique for the parametric identification of parameter α is based on the experimental data, assuming that the rest of the parameters of the equation are known (with some degree of accuracy). This technique was developed due to the possibility of finding a solution at any point. conclusion about their identity (Fig. 3).   So, assume that we know several experimental points u(x i ) = U i , i = 1, ..., N. To find the parameter α we minimize the deviation of the theoretical curves from the experimental curves. For calculating of the theoretical points we use Expression (33) for u(x i , α). Using the least-squares method we determine the deviation function This function represents the sum of the deviations of the theoretical points from the experimental points. The value of α minimizing this function we can consider approximately as the search value. We shall note that the identification accuracy depends on the number of experimental points, together with the accuracy of the other parameters of the system. The method of parametric identification that we provide compares various nomographic techniques [39,40] (and references therein) and its advantage consists in the accurate quantitative estimation of choosing the search parameter. It is important that the deviation function (35) can be constructed on the entire range of supposed values of the parameter; this improves the accuracy of the identification. For testing the provided technique, we use the experimental data obtained in [35] (and references therein). The values for samples of polymer concrete based on polyester resin (dian-and dihloangidrid-1,1-dichloro-2,2-diethylene) are presented in Table 1.  Finally, using these data, we construct the deviation function (40) presented in Figure 4. The constructed graph shows that the deviation function has the minimum for α ≈ 1.47, and it allows us to assume that the order of the fractional derivative in Expression (31) is equal to 1.47. Figure 5 presents the experimental points and the theoretical curve. Comparing the experimental data presented in Table 1 with the model allows us to draw a conclusion that the provided model is adequate and our techniques for parametric identification have a high level of accuracy. The knowledge of the parameter α in the model (31) allows us, in particular, to predict various stress-strain characteristics of the material (polymer concrete, asphaltic concrete, etc.) when subjected to loading. Now, having a technique for the parametric identification of the order of the fractional derivative in the Begley-Torvik model, we will proceed to the presentation of the method of separation of variables and its application to find the deformation-strength characteristics of polymer concrete.   Finally, using these data, we construct the deviation function (40) presented in Fig. 4. Constructed 303 graph shows that the deviation function has the minimum for α ≈ 1.47, and it is allow us to assume 304 that the order of the fractional derivative in expression (31) is equal to 1.47. Figure 5 presents the 305 experimental points and the theoretical curve. Comparing the experimental data presented in Table 1 306 with the model allows us to draw a conclusion that provided model is adequate and our technique for

Method of Separation of Variables for Time-Space Fractional Vibration Equations
As it was noted in [41] the fractional calculus has attracted the attention of many authors in recent years. In this regard, we should note the paper [14], as a unique comprehensive review of fractional calculus and its application with the authoritative contribution of leading world experts. As it was noted earlier, we can describe anomalous diffusion or dispersion using the fractional space derivatives, and some processes with 'memory' effects using the fractional time derivatives. In [41] (and see the references therein), the following equation was investigated and was used to describe, in particular, the vibration of a string taking into account friction in a medium with fractal geometry. This equation may be used to model changes in the deformation-strength characteristics of polymer concrete under loading.
In domain D = 0 < x < 1, 0 < t < 1, we consider the first boundary value problem for the equation for the equation of vibration of a string with a fractional derivative of order with respect to the spatial variable Here, 0 < α < 2, c-constant, D α 0x u-constant, D α 0x u-fractional derivative of the Riemann-Liouville type of order α. Fractional derivative of order α for function f (x) in a point x(0 = m − 1 < α < m), m ∈ N) defined by the formula Obtained results are applied [14,41] for modeling changes in the deformation-strength characteristics of polymer concrete under loading. The solution to problem (41)-(44) will be sought by the Fourier method u(x, t) = X(x)T(t).
Then the solution to problem (41)-(44) is written out in the standard way putting in the last expression t − 0, we have To find A m differentiate both parts (12) by t and let t = 0 we obtain, from here [A m Z m (0) + B m Z m (0)](X m (x), X m (x)) = (ψ(x), X m (x)).
Finally, we have, which allows us to write a solution to problem (41)-(44) in the form (50).
In the earlier papers of the author, Equation (46) was used to model the certain deformation-strength characteristics of polymer concrete. In this paper, only transverse vibrations are considered, all movements occur in one plane and the granule moves perpendicular to the axis. Then, to simulate changes in the deformation-strength characteristics of polymer concrete under loading, we have the following first boundary-value problem (here u(x, t)-granule displacement in moment t) u(0, t) = u(1, t) = 0, u(x, 0) = ϕ(x) = 0, C k n λ n−k (−C 1 ) k Γ(2n − kβ + 2) x 2n+1−1.47k , x 2n+1−1.47k , Let us find eigenvalues λ j numerically using the high-level language of technical calculations MATLAB taking α = 1.47, C 1 = 1.8 (according by [41]). Eigenvalues are presented in Table 2. Then, an approximate solution to problem (51)-(54) will take the form Formula (55) allows us to write a solution to the problem (51)-(54) if the functions ψ(x) and ϕ(x) are continuously differentiable. Finally, it remains to determine the parameter β. This parameter can again be determined by the technique developed in [14,41] since the parameter α is already defined.

Remark 3.
We shall note the papers [42,43], where the same problem is solved by numerical methods and compared with the solution (55).
Funding: This research received no external funding.