Next Article in Journal
Finite Difference Method for Two-Sided Two Dimensional Space Fractional Convection-Diffusion Problem with Source Term
Next Article in Special Issue
A Nonlinear Model Predictive Control with Enlarged Region of Attraction via the Union of Invariant Sets
Previous Article in Journal
Nonlinear Observability Algorithms with Known and Unknown Inputs: Analysis and Implementation
Previous Article in Special Issue
Modeling of Strength Characteristics of Polymer Concrete Via the Wave Equation with a Fractional Derivative
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Temirkhan Aleroev
National Research Moscow State University of Civil Engineering (NRU MGSU), Yaroslavskoe Shosse, 26, 129337 Moscow, Russia
Mathematics 2020, 8(11), 1877; https://doi.org/10.3390/math8111877
Submission received: 16 September 2020 / Revised: 18 October 2020 / Accepted: 22 October 2020 / Published: 29 October 2020
(This article belongs to the Special Issue Dynamical Systems and Optimal Control)

Abstract

:
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.

In memoriam of my Father Sultan
and my son Bibulat

1. Introduction

Let φ ( x ) L 1 ( 0 , 1 ) . Then the function
d α d x α φ ( x ) 1 Γ ( α ) 0 x ( x t ) α 1 φ ( t ) d t L 1 ( 0 , 1 )
is known as a fractional integral of order α > 0 beginning at x = 0 [1]. Here Γ ( α ) is the Euler gamma-function. As is known (see [1]), the function ψ ( x ) L 1 ( 0 , 1 ) is called the fractional derivative of the function φ ( x ) L 1 ( 0 , 1 ) of order α > 0 beginning at x = 0 , if
φ ( x ) = d α d x α ψ ( x ) ,
which is written
ψ ( x ) = d α d x α φ ( x ) .
Then
d α d x α
denote the fractional integral for α < 0 and the fractional derivative for α > 0 .
Let { γ k } 0 n be a set of real numbers satisfying the condition 0 < γ j 1 , ( 0 j n ). We denote
σ k = j = 0 k γ j 1 ;
μ k = σ k + 1 = j = 0 k γ j ( 0 k n ) ,
and assume that
1 ρ = j = 0 n γ j 1 = σ n = μ n 1 > 0 .
Following M. M. Dzhrbashyan [1], we consider the integro-differential operators
D ( σ 0 ) φ ( x ) d ( 1 γ 0 ) d x ( 1 γ 0 ) φ ( x ) ,
D ( σ 1 ) φ ( x ) d ( 1 γ 1 ) d x ( 1 γ 1 ) d γ 0 d x γ 0 φ ( x ) ,
D ( σ 2 ) φ ( x ) d ( 1 γ 2 ) d x ( 1 γ 2 ) d γ 1 d x γ 1 d γ 0 d x γ 0 φ ( x ) ,
D ( σ n ) φ ( x ) d ( 1 γ n ) d x ( 1 γ n ) d γ n 1 d x γ n 1 d γ 0 d x γ 0 φ ( x ) .
We denote D a x α the operator of fractional integro-differentiation of order α beginning at a R and with end at x R of order [ α ] . By definition we have
D a x ( α ) ϕ ( t ) = s i g n ( x a ) Γ ( α ) a x ϕ ( t ) d t ( x t ) α 1 , α < 0 , ϕ ( t ) L 1 [ a , b ] , ϕ ( t ) , α = 0 , ϕ ( t ) L 1 [ a , b ] , s i g n [ α ] + 1 ( x a ) [ α ] + 1 x [ α ] + 1 D a x α [ α ] 1 ϕ ( t ) , α > 0 , ϕ ( t ) L 1 [ a , b ] ,
where [ α ] is the integer part of α , which satisfies [ α ] α < [ α ] + 1 , and x [ a , b ] .

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
2 u t 2 = 2 u x 2 + C 1 D 0 x α u + C 0 D 0 t β u , 0 < α , β < 2 ,
and
u ( x , t ) t = D 0 + α u ( x , t ) .
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
L ( u ; γ 0 , γ 1 , γ 2 , q ( x ) ) = D ( σ 2 ) u [ λ + q ( x ) ] u ( x ) = 0 ,
X ( x ) + C 1 D 0 x α X = λ X ( x ) ,
(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 L2(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.

2. 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 α , 1 , 0 ) = 1 Γ ( α ) d d x 0 x u ( ζ ) ( x ζ ) 1 α d ζ + λ u ( x ) = 0
u ( 0 ) = 0 , u ( 1 ) = 0 ,
and
L ( u ; 1 , 1 , 1 α , 0 ) = 1 Γ ( α ) 0 x ( x ζ ) α 1 u ( ζ ) d ζ + λ u ( x ) = 0 ,
u ( 0 ) = 0 , u ( 1 ) = 0 .
These problems are the focus of many researchers.
First, we note that in [4] (and references therein), the following problem was considered
u + λ d α d x α u = 0 , 0 < α < 1 ,
u ( 0 ) = 0 , u ( 1 ) = 0 .
with studying of the spectrum of the operator
D ( β ) u = d α d x α d 2 d x 2 u = 1 Γ ( α ) 0 x ( x ζ ) α 1 u ( ζ ) d ζ , ( β = 2 α )
(the operator D ( σ 2 ) transforms to the operator D ( β ) if γ 0 = γ 1 = 1 and γ 2 = 1 α ) .
The operator D ( β ) arose great interest after F. Mainardi’s paper [3]. In this paper, the following equation was considered
1 Γ ( 2 γ ) 0 x u ( ζ ) ( t ζ ) γ 1 d ζ + ω γ u ( t ) = 0
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
D ( β ) u = 1 Γ ( α ) 0 x ( x ζ ) α 1 u ( ζ ) d ζ = λ u
coincides with the solution for Equation (4);
2.
Equations (4) and (7) are equal if
lim x 0 d ( 1 α ) d x ( 1 α ) u ( x ) = 0 .
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:
u ( x ) λ Γ ( 2 α ) 0 1 x ( 1 ζ ) 1 α u ( ζ ) d ζ 0 x ( x ζ ) 1 α u ( ζ ) d ζ =
= u ( x ) λ Γ ( 2 α ) 0 1 G 0 ( x , ζ ) u ( ζ ) d ζ = 0 ,
u ( x ) λ Γ ( 2 α ) 0 1 x 1 α ( 1 ζ ) 1 α u ( ζ ) d ζ 0 x ( x ζ ) 1 α u ( ζ ) d ζ =
u ( x ) λ Γ ( 2 α ) 0 1 G 1 ( x , ζ ) u ( ζ ) d ζ = 0
where G 0 ( x , ζ ) is the Green’s function of the problem (4) and (5), which was constructed in [2] and G 1 ( x , ζ ) is the Green function of the problem
L ( u ; 1 , 1 α , 1 , 0 ) = 1 Γ ( α ) d d x 0 x u ( ζ ) ( x ζ ) 1 α d ζ + λ u ( x ) = 0 ,
u ( 0 ) = 0 , u ( 1 ) = 0 ,
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.
Theorem 1.
If
0 < α < 32 π 2 9 + 2 3 1 ,
then the first eigenvalue of the problem (4) and (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 ) ).
Proof. 
That the first eigenvalue of the problem (4) and (5) is positive and simple for
0 < α < 32 π 2 9 + 2 3 1
was proved in [7,8]. We show now that basic (main) tone of the problem (4) and (5) has no nodes.
It is known that a number λ will be an eigenvalue of the problem (4) and (5) [4] if and only if this value λ is the root (zero) of the function E 1 / β ( λ ; 2 ) and the corresponding eigenfunctions of the problem (4) and (5) are
u n ( x ) = x E 1 β ( λ n x β ; 2 ) , n = 1 , 2 , 3
where
λ 1 , λ 2 , , λ n ,
are zeros of the function E β ( λ ; 2 ) , numbered according to the non-decreasing of their modules,
E ρ ( z ; μ ) = k = 0 z k ( μ + k ρ 1 )
is a Mittag–Leffler function. We shall show that the function
u 1 ( x ) = x E 1 β ( λ 1 x β ; 2 )
does not vanish in ( 0 , 1 ) . Let x 0 ( 0 , 1 ) be such that
x 0 E 1 β ( λ 1 x 0 β ; 2 ) = 0 .
Then the number λ 1 x 0 β is a zero of E β ( λ ; 2 ) , moreover λ 1 x 0 β < λ 1 (since x 0 ( 0 , 1 ) ). This contradicts the assumption that λ 1 is the first zero of the function E 1 / β ( λ ; 2 ) . Theorem 1 is proved. □
Since for α > 2 3 , the function E 1 / β ( λ ; 2 ) has no real zeros [9], then the problem (4) and (5) has this main (at least first eigenvalues are real) oscillatory property only for small α .
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) and (9) is real, simple, and satisfies the condition
0 < λ 1 1 < Γ ( 2 + 2 α ) Γ ( 1 + α ) ,
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
0 < λ 1 1 < Γ ( 2 + 2 α ˜ ) Γ ( 1 + α ˜ ) , α ˜ = 1 α
holds. As it was mentioned above, the problem (8) and (9) is equivalent to the integral equation of Fredholm (II kind)
u ( x ) + λ Γ ( 1 + α ˜ ) 0 x ( x ζ ) α ˜ u ( ζ ) d ζ x α ˜ 0 1 ( 1 ξ ) α ˜ u ( ζ ) d ζ = 0 ,
and the value λ is an eigenvalue of the problem (8) and (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
A u = 1 Γ ( 1 + α ˜ ) 0 x ( x ζ ) α ˜ u ( ζ ) d ζ x α ˜ 0 1 ( 1 ζ ) α ˜ u ( ζ ) d ζ
as
A u = A 0 u A 1 u ,
where
A 0 u = 1 Γ ( 1 + α ˜ ) 0 x ( x ζ ) α ˜ u ( ζ ) d ζ ,
and
A 1 u = 1 Γ ( 1 + α ˜ ) 0 1 x α ˜ ( 1 ζ ) α ˜ u ( ζ ) d ζ .
As operators A 0 and A 1 are trace class operators [10], then
s p A = s p ( A 0 A 1 ) = s p ( A 0 ) s p ( A 1 ) .
Since A 0 is Volterra’s operator, then s p ( A 0 ) = 0 , and so
s p ( A ) = s p ( A 1 ) .
It is easy to find the trace of the operator A 1 ( A 1 it is one-dimensional operator). Let us consider the equation
u ( x ) λ Γ ( 1 + α ˜ ) 0 1 x α ˜ ( 1 ζ ) α ˜ u ( ζ ) d ζ = 0 .
The Fredholm determinant of this equation is
d ( λ ) = | 1 λ K 11 | ,
where
K 11 = 1 Γ ( 1 + α ˜ ) 0 1 ζ α ˜ ( 1 ζ ) α ˜ d ζ = Γ ( 1 + α ˜ ) Γ ( 2 + 2 α ˜ ) .
From this we obtain
s p ( A ) = Γ ( 1 + α ˜ ) Γ ( 2 + 2 α ˜ ) .
Thus
λ 1 1 + i = 2 λ i 1 = Γ ( 1 + α ˜ ) Γ ( 2 + 2 α ˜ ) .
Since the kernel of the operator A is non-negative, then λ 1 is a positive number, and i = 2 λ i 1 is positive, so
λ 1 1 < Γ ( 1 + α ˜ ) Γ ( 2 + 2 α ˜ ) .
Theorem 2 is proved. □
Corollary 1.
Since λ is an eigenvalue of the problem (8) and (9) [5,7] if and only if λ is a zero of the function E 1 / ( 1 + α ˜ ) ( λ ; 1 + α ˜ ) , and the corresponding eigenfunctions of the problem (8) and (9) are
u n ( x ) = x α ˜ E 1 / ( 1 + α ˜ ) ( λ n x 1 + α ˜ ; 1 + α ˜ )
n = 1 , 2 , 3 , where λ 1 , λ 2 , , λ n , are zeros of the function E 1 / ( 1 + α ˜ ) ( λ ; 1 + α ˜ ) , numbered by their non-decreased modules, then the function E 1 / ( 1 + α ˜ ) ( λ ; 1 + α ˜ ) has positive and simple first zero, and the function
u 1 ( x ) = x α ˜ E 1 / ( 1 + α ˜ ) ( λ 1 x 1 + α ˜ ; 1 + α ˜ )
does not vanish in ( 0 , 1 ) .
Conclusion. From these theorems follows
(a)
Which of Equations (4) or (8) is the correct choice as an oscillatory;
(b)
How the spectral structure of L ( u ; γ 0 , γ 1 , γ 2 , q ( x ) ) depends on generatrices γ 0 , γ 1 , γ 2 of order of fractional differential equation;
(c)
How the nature of the modeling process helps to indentify generatrices γ 0 , γ 1 , γ 2 of order of operator.

2.2. 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
D ( σ 2 ) u [ λ + q ( x ) ] u ( x ) = 0 ,
where
D ( σ 2 ) u = 1 Γ ( 1 γ ) d d x 0 x u ( ζ ) ( x ζ ) γ d ζ , 0 < γ < 1 , σ 2 = 1 + γ .
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
u ( 0 ) = 0 , u ( 1 ) = 0 ,
for Equation (10) with q ( x ) = 0 is equivalent to the integral equation
1 Γ ( 2 γ ) 0 x ( x t ) 1 γ u ( t ) d t 0 1 x 1 γ ( 1 t ) 1 γ u ( t ) d t = λ u .
We have:
Theorem 3.
Let γ 0 = γ 1 = 1 , q ( x ) 0 . Then the system of eigenfunctions and associated functions of the problem (10) and (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) and (11), is sectorial [12].
Theorem 4.
All eigenvalues of the problem (10) and (11) for q ( x ) 0 are in the angle | a r g z | < π ( 1 γ ) 2 , 0 < γ < 1 .
Proof. 
Consider the expression ( D ( σ 2 ) f , f ) . It is obvious that
( D ( σ 2 ) f , f ) = 1 Γ ( 1 γ ) d d x 0 x f ( t ) ( x t ) γ d t , f ( x )
= ( 1 Γ ( 1 γ ) 0 x f ( t ) ( x t ) γ d t , f ( x ) ) = ( 0 J x α f , f )
where α = 1 γ and 0 J x α is the operator of fractional integration of order α :
( 0 J x α f ) ( x ) = 1 Γ ( α ) 0 x ( t s ) 1 α f ( s ) d s .
By a well-known Matsaev–Palant theorem ([13], p. 481), the values of the form ( 0 J x α f , f ) is in the angle | a r g z | < π α 2 . This proves Theorem 4. □
Since the number λ is an eigenvalue of the problem (10) and (11) if λ is a zero of the function E 1 / μ ( λ ; μ ) ( μ = 1 + γ ) [5], the following proposition is valid.
Corollary 2.
All zeros of the function E 1 / μ ( λ ; μ ) are in the angle | a r g z | < ( π ( 1 γ ) ) / ( 2 ) , 0 < γ < 1 . Here μ = 1 + γ .
Theorem 5.
The problem (10) and (11) for q ( x ) 0 has no eigenvalues inside the circle with radius Γ ( 4 2 γ ) / Γ ( 2 γ ) centered at the coordinate origin.

2.3. 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
A ρ ( u ) = 0 1 G ( x , t ) u ( t ) d t = 1 Γ ( ρ 1 ) 0 x ( x t ) 1 ρ 1 u ( t ) d t 0 1 x 1 ρ 1 ( 1 t ) 1 ρ 1 u ( t ) d t ,
which was for the first time studied in [7]. Here, 0 < ρ < 2 , and
G ( x , t ) = ( 1 t ) 1 ρ 1 x 1 ρ 1 ( x t ) 1 ρ 1 Γ ( ρ 1 ) , 0 t x 1 , ( 1 t ) 1 ρ 1 x 1 ρ 1 Γ ( ρ 1 ) , 0 x t 1 ,
is the Green function of the following problem S (for λ = 0 ):
1 Γ ( n ρ 1 ) d n d x n 0 x ( x s ) n ρ 1 1 u ( s ) d s + λ u = 0 ,
( n 1 ρ 1 < n , n = [ ρ 1 ] + 1 , where [ ρ 1 ] is the integer part of the number ρ 1 )
u ( 0 ) = 0 , u ( 0 ) = 0 , , u ( n 2 ) ( 0 ) = 0 , u ( 1 ) = 0 .
In this case [7], if γ 0 = γ 1 = = γ n = 1 then the problem S takes the form
u ( n ) + λ u = 0 ,
u ( 0 ) = 0 , u ( 0 ) = 0 , , u ( n 2 ) ( 0 ) = 0 , u ( 1 ) = 0 ,
of which the Green function G ( x , t ) (for λ = 0 ) reads
G ( x , t ) = ( 1 t ) n 1 x n 1 ( x t ) n 1 ( n 1 ) ! , 0 t x 1 , ( 1 t ) n 1 x n 1 ( n 1 ) ! , 0 x t 1 . .
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:
  • For ρ > 1 the operator A ρ is completely nonself-adjoint [5,7,15];
  • For ρ 1 the operator A ρ is sectorial [6] (and see the references therein);
  • 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
T ( n ) u c n 1 ( a | | u | | + b | | T u | | ) , u D , n = 1 , 2 , .
Then for | κ | < 1 / c the series
T ( κ ) u = T u + κ T ( 1 ) u + κ 2 T ( 2 ) u + , u D
defines the operator T ( κ ) with the domain of definition D. If | κ | < ( b + c ) 1 , 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
A ( ε ) u = 0 x ( x t ) 1 + ε u ( t ) d t + 0 1 x 1 + ε ( 1 t ) 1 + ε u ( t ) d t
forms a holomorphic family of type ( A ) , i.e.,
A ( ε ) u = A ( 0 ) u + ε A 1 u + ε 2 A 2 u + + ε n A n u +
where
A ( 0 ) u = 0 x ( x t ) u ( t ) d t + 0 1 x ( 1 t ) u ( t ) d t
the unperturbed operator, and
A n u ( x ) = 0 x K ˜ ( x , t ) n K ( x , t ) n u ( t ) d t ,
K ˜ ( x , t ) n = x ( 1 t ) ln n ( 1 t ) x n ! ,
K ( x , t ) n = x ( 1 t ) ln n ( 1 t ) x n ! , t < x , 0 , t x .
Theorem 8.
If | ε | < 3 / 2 , then the operator
B ˜ ( ε ) u = 0 x ( x t ) 1 + ε u ( t ) d t + 0 1 x ( 1 t ) 1 + ε u ( t ) d t
forms a holomorphic family of type ( A ) where
B ( 0 ) u = 0 x ( x t ) u ( t ) d t + 0 1 x ( 1 t ) u ( t ) d t
is the unperturbed operator, and
B n u ( x ) = 0 x K ¯ ( x , t ) n K ( x , t ) n u ( t ) d t ,
where
K ¯ ( x , t ) n = x ( 1 t ) ln n ( 1 t ) x n ! ,
and
K ( x , t ) n = x ( 1 t ) ln n ( 1 t ) x n ! t < x , 0 , t x .
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) and (5). In [8], it was written that “... in the limiting case α = 0 the problem (4) and (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.
Theorem 9.
lim α α 0 + λ n ( α ) = lim α α 0 λ n ( α ) = λ n ( α 0 ) for any α 0 [ 0 , 1 ] .
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 , I m ( μ ) = 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]:
| ε | < m i n ζ Γ ( a | | R ( ζ , B ˜ ( 0 ) ) | | + b | | B ˜ ( 0 ) R ( ζ , B ˜ ( 0 ) ) | | + c ) 1
(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 γ [ α , β ] (and not only operators of the form A γ [ α , β ] ) in the general scheme of perturbation theory.

2.4. Solving the Problem of Finding the Radon Flux Density by Its Concentration at Different Depths of the Earth’s Surface by the Method of Separated Variables

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].
u ( x , t ) t = D 0 + α u ( x , t ) .
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
u ( x , t ) t = D 0 + α u ( x , t ) ,
where D 0 + α u ( x , t ) —the Riemann-Liouville fractional derivative of the order α , with boundary conditions
u ( 0 , t ) = u ( 1 , t ) = 0 ,
u ( x , 0 ) = ϕ ( x ) ,
Using the method of separation of variables [21], we can write out the solution to this problem
u ( x , t ) = n = 1 δ n e x p ( λ n t ) x α 1 E α , α ( λ n x α ) .
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 ,
-
Fourier coefficients ϕ ( x ) , and the system of functions z n n = 1 = ( 1 x ) α 1 E α , α ( λ n ( 1 x ) α ) is the biorthogonal system of eigenfunctions ω n ( x ) = x α 1 E α , α ( λ n x α ) ( z n —is the system of eigenfunctions of the contiguous boundary problem).
For an approximate solution of this problem, one can use the formula
u ( x , t ) n = 1 N δ n e x p ( λ n t ) x α 1 E α , α ( λ n x α ) .
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
u ( x , t ) n = 1 50 δ n e x p ( λ n t ) x α 1 E α , α ( λ n x α ) .
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.

3. 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
u + ε D 0 x α u = λ u ,
and the boundary conditions of the Sturm–Liouville type
u ( 0 ) = 0 , u ( 1 ) = 0 .
Note, that for 0 < α < 2 , the spectral structure of operator B ˜ generated by problem
u + ε D 0 x α u = λ u ,
u ( 0 ) = 0 , u ( 1 ) = 0 ,
was considered in detail in our paper [12]. In particular, the following theorem was shown:
Theorem 11.
If | ε | < 10 20 , then all eigenvalues of operator B ˜ are simple and real.
From this theorem, it follows that the operator B ˜ generated no associated functions.
Theorem 12.
The number λ is an eigenvalue of problem (19) and (20) if λ is a zero of the function
ω ( λ ) = 1 + n = 1 m = 0 n ( ε ) n C n m λ n m Γ ( 2 n m α + 2 ) .
The eigenfunctions of problem ( A ) take the form
χ i ( x ) = x + n = 1 m = 0 n ( ε ) n C n m λ i n m Γ ( 2 n m α + 2 ) x 2 n m α + 1 ,
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) and (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 0 x α 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 M is a compact operator and N is a positive operator, then by Keldysh’s theorem ([10], p. 318) we have
lim x 0 n ( r , B ) n ( r , N ) = 1
if for the distribution function n ( r , B ) of characteristic values of the operator N we may choose non-decreasing function φ ( r ) ( 0 r ) such that [10]:
  • lim r φ ( r ) = ;
  • lim r ( ln φ ( r ) ) < ;
  • lim r n ( r , B ) φ ( r ) = 1 .
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:
R e ( A f , f ) 0 , ( f D A ) .
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.,
I m ( A f , f ) 0 , ( f D A ) .
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 A ˜ , is densely defined in a Hilbert space H with values of the form ( A ˜ f , f ) in the left half-plane
R e ( A ˜ f , f ) 0 , ( f D A ˜ ) .
As the operator B ˜ describes oscillations of the mechanical system, then it should be dissipative (see [30] and references therein).
In this paragraph we show that the operator B ˜ 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
R e ( D 0 x α u , u ) 0 , 0 < α < 1 ,
and
R e ( D 0 x α u , u ) 0 , 1 < α < 2 .
Theorem 14.
If 0 < α < 1 and ε > 0 , then the operator, generated by the problem
u ε D 0 x α u = λ u
u ( 0 ) = 0 , u ( 1 ) = 0
is dissipative.
Proof. 
This theorem follows from the relation (23) and the fact that the operator
T u = u , u ( 0 ) = 0 , u ( 1 ) = 0 ,
is dissipative. □
Theorem 15.
If 1 < α < 2 and ε < 0 , then the operator, generated by the problem
u ε D 0 x α u = λ u
u ( 0 ) = 0 , u ( 1 ) = 0
is dissipative.
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
2 u t 2 = 2 u x 2 + C 1 D 0 x α u + C 0 D 0 t β u , 0 < α , β < 2 ,
u ( 0 , t ) = u ( 1 , t ) ,
u ( x , 0 ) = φ ( x ) ,
u t ( x , 0 ) = φ ( x ) ,
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
u ε D 0 x α u = λ u ,
u ( 0 ) = 0 , u ( 1 ) = 0 ,
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) and (30) is of fixed sign (we shall note, that Green’s function of problem (29) and (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) and (30), for 1 < α < 2 , has the following complex structure
G 2 ( x , τ ) = G 1 ( x , τ ) ε E 1 / 2 ( ε , 2 ) τ 1 E β [ ε ( η τ ) ] β d η 0 1 G ( x , t ) D 0 t α 1 E β [ ε t β ] d t ,
G 1 ( x , τ ) = ( 1 x ) τ x E β [ ε ( t τ ) ] d t x x 1 E β [ ε ( t τ ) β ] d t , x τ , x τ 1 E β [ ε ( t τ ) β ] d t , x τ .
For 0 < α < 1 , | ε | < 1 / 4 , Green’s function of the problem (29) and (30) was constructed in [25] (and references therein). Let us show how this function was constructed. Since the problem (29) and (30), for 0 < α < 1 , is equivalent to the equation
u ( x ) + ε Γ ( 2 α ) 0 x ( x t ) 1 α u ( t ) d t 0 1 x ( 1 t ) 1 α u ( t ) d t = λ 0 1 G ( x , t ) u ( t ) d t
then
u ( x ) = λ ( I ε K ) 1 0 1 G ( x , t ) u ( t ) d t ,
where
G ( x , t ) = t ( x 1 ) , t x , x ( t 1 ) , t > x ,
K u = 1 Γ ( 2 α ) 0 x ( x t ) 1 α u ( t ) d t + 1 Γ ( 2 α ) 0 1 x ( 1 t ) 1 α u ( t ) d t
=   x J 0 , 1 2 α J 0 , x 2 α   u .
We can show that
K 2 u = K · K u = ( K x ) J 0 , 1 2 α u x J 0 , 1 4 2 α u + J 0 , x 4 2 α u
K 3 u = K · K 2 u = ( K 2 x ) J 0 , 1 2 α u ( K x ) J 0 , 1 4 2 α u + x J 0 , 1 6 3 α u J 0 , x 6 3 α u .
By induction, we have
K n = i = 1 n ( 1 ) i + 1 ( K n i x ) J 0 , 1 ( 2 α ) i u + ( 1 ) n + 2 J x ( 2 α ) n u .
Thus,
( I ε K ) 1 = I + n = 1 ( ε K ) n u
= I + n = 1 ( ε ) n i = 1 n ( 1 ) i + 1 ( K n i x ) J 0 , 1 ( 2 α ) i u + ( 1 ) n + 2 J x ( 2 α ) n u
= I + n = 1 i = 1 n ( ε ) n ( 1 ) i + 1 ( K n i x ) J 0 , 1 ( 2 α ) i u + n = 1 ( 1 ) n + 2 J x ( 2 α ) n u .
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) and (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.

3.1. 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).

3.1.1. The Bagley–Torvik Equation and the Laplace Transform

We consider the problem
u ( x ) + c D α u ( x ) + λ u ( x ) = 0 , u ( 0 ) = 0 , u ( x ) = 1 ,
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]:
D α u ( x ) = d 2 d x 2 1 Γ ( 2 α ) 0 t u ( τ ( x τ ) α 1 d τ .
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
0 x u ( t ) d t + c 0 x D α u ( t ) d t + λ 0 x u ( t ) d t = 0 ,
0 x d u ( t ) d t + c Γ ( 2 α ) 0 x d 2 d t 2 0 t u ( τ ) ( t τ ) α 1 + λ 0 x u ( t ) d t = 0 ,
u ( x ) u ( 0 ) + c Γ ( 2 α ) d d t 0 t u ( τ ) d τ ( t τ ) α 1 | 0 x + λ 0 x u ( t ) d t = 0 ,
u ( x ) 1 + c Γ ( 2 α ) d d t 0 x u ( t ) d t ( x t ) α 1 + λ 0 x u ( t ) d t = 0 .
Obtained expression (32) we integrate from 0 to x again:
0 x u ( t ) d t 0 x d t + c Γ ( 2 α ) 0 x d d t 0 t u ( τ ) d τ ( t τ ) α 1 + λ 0 x 0 t u ( τ ) d τ d t = 0 ,
u ( x ) x + c Γ ( 2 α ) 0 x ( x t ) 1 α u ( t ) d t + λ 0 x 0 t u ( τ ) d τ d t = 0 .
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),
U ( s ) = 0 e s t u ( t ) d t .
It is clear that the function
0 x ( x t ) 1 α u ( t ) d t
represents the convolution of the functions u ( x ) and x 1 α . For the convolution of functions, there exists the simple formula of images
0 x f 1 ( x ) f 2 ( x t ) d t = F 1 ( s ) F 2 ( s ) ,
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
0 x 0 t u ( τ ) d τ d t = 0 x d t 0 t u ( τ ) d τ = U ( s ) / s 2 .
After some clearly transformations we obtain
U ( s ) 1 s 2 + c U ( s ) s α 2 + λ U ( s ) / s 2 = 0 .
From this we obtain the formula for the image
U ( s ) = 1 s 2 + c s α + λ .
Formula (36) makes it possible to express the solution of problem (31) using Laplace’s integral
u ( x ) = 1 2 π i σ i σ + i e s t U ( s ) d s .

3.1.2. 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 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)
ρ ( α , δ ) = | u ( x , α ) u ( x , α ) | d x ,
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.
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
u ( x ) = x n = 1 m = 0 n ( 1 ) n + 1 C n m c m λ n m x 2 n + 1 m α Γ ( 2 n m α + 2 )
Comparing the graphs of the solutions numerically obtained by Formulas (37) and (39), we can draw a conclusion about their identity (Figure 3).
h ( t , α ) : = t n = 1 50 m = 0 n ( 1 ) n + 1 c o m b i n ( n , m ) c m λ n m t 2 n + 1 m α Γ ( 2 n m α + 2 )

3.1.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.
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
F ( α ) = i = 1 N ( U i u ( x i , α ) ) 2 .
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.

3.2. 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
2 u t 2 = 2 u x 2 + C 0 D 0 t α u + C 1 D 0 x β u + F
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
2 u t 2 = 2 u x 2 + C 1 D 0 x α u + C 0 D 0 t β u , 0 < α , β < 2 ,
u ( 0 , t ) = u ( 1 , t ) ,
u ( x , 0 ) = φ ( x ) ,
u t ( x , 0 ) = φ ( x ) ,
Here, 0 < α < 2 , c-constant, D 0 x α u -constant, D 0 x α 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
D α f ( x ) = d m d x m 1 Γ ( m a ) a x f ( τ ) ( x τ ) α + 1 m d τ .
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 ) .
We substitute (45) into Equation (41), then for an unknown function X ( x ) we obtain the two-point Dirichlet problem
X ( x ) + C 1 D 0 x α X = λ X ( x ) ,
X ( 0 ) = X ( 1 ) = 0 .
The solution to problem (46) and (47) was written out in [1,2], see the references therein. In particular, it was shown that the number λ is an eigenvalue of problem (46) and (47), if and only if λ is the zero of the function
ω ( λ ) = n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 )
and the corresponding eigenfunctions X j ( x ) have the form
X j ( x ) = n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 k α , j = 1 , 2 , 3 ,
(here j -j-th eigenfunction of the problem (46) and (47)). The system of the eigenfunctions (48) is complete [14,41] but not orthogonal, thus we construct the system
X ˜ j ( x ) = ( 1 x ) n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 k α , j = 1 , 2 , 3 ,
which is biorthogonal to the system of eigenfunctions
X j ( x ) = n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 k α , j = 1 , 2 , 3 ,
Next, we find the general solution to the equation
T ( t ) + c 0 D 0 t β T ( t ) = λ T ( t ) .
As in the case of Equation (46), we have
T m ( t ) = A m t + n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 k α + B m 1 + n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 1 ) x 2 n + 1 k α .
Let us designate
Z m ( t ) = t + n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 k α , Z ˜ m ( t ) = 1 + n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 1 ) x 2 n + 1 k α .
Then the solution to problem (41)–(44) is written out in the standard way
u ( x , t ) = m = 1 T m ( t ) X m ( x ) = m = 1 [ A m Z m ( t ) + B m Z ˜ m ( t ) ] X m ,
putting in the last expression t 0 , we have
φ ( x ) = n = 0 B m Z ˜ ( 0 ) Z ( 0 ) X m ( x ) ,
so
B m = 1 Z ˜ ( 0 ) ( φ ( x ) , X ˜ m ( x ) ) ( X m ( x ) , X ˜ m ( x ) ) .
To find A m differentiate both parts (12) by t and let t = 0 we obtain,
m = 1 [ A m Z m ( 0 ) + B m Z ˜ m ( 0 ) ] X m ( x ) = ψ ( x )
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,
A m = 1 Z m ( 0 ) 1 X m ( x ) , X ˜ m ( x ) B m Z ˜ m ( 0 )
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)
2 u t 2 = 2 u x 2 + C 0 D 0 x β u + C 1 D 0 t 1.47 u , 0 < α , β < 2 ,
u ( 0 , t ) = u ( 1 , t ) = 0 ,
u ( x , 0 ) = φ ( x ) = 0 ,
u t ( x , 0 ) = φ ( x ) = 0 ,
the solution of which according to Formula (50) has the form
u ( x , t ) = m = 1 [ A m Z m ( t ) + B m Z ˜ m ( t ) ] X m ,
where
X j ( x ) = n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 1.47 k ,
X ˜ j ( x ) = ( 1 x ) n = 0 k = 0 n C n k λ n k ( C 1 ) k Γ ( 2 n k β + 2 ) x 2 n + 1 1.47 k ,
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
u ( x , t ) = m = 1 5 [ A m Z m ( t ) + B m Z ˜ m ( t ) ] X m ,
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.

Acknowledgments

I am very thankful for reviewers for their very good grades and remarks, which made my manuscript better.

Conflicts of Interest

The author declares 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.

References

  1. Dzhrbashian, M.M. The boundary-value problem for a differential fractional-order operator of the Sturm–Liouville type. Izv. Akad. Nauk ArmSSR Ser. Mat. 1970, 5, 71–96. [Google Scholar]
  2. Aleroev, T.S. The Sturm–Liouville problem for differential fractional-order equations with fractional derivatives in lower terms. Differ. Uravn. 1982, 18, 341–342. [Google Scholar]
  3. Mainardi, F. Fractional relaxation-oscillation and fractional diffusion-wave phenomena chaos. Solut. Fract. 1996, 7, 1461–1477. [Google Scholar] [CrossRef]
  4. Nakhushev, A.M. Fractional Calculus and Its Application; Fizmatlit: Moscow, Russia, 2003. (In Russian) [Google Scholar]
  5. Aleroev, T.S. Boundary-Value Problems for Differential Equations with Fractional Derivatives. Ph.D. Thesis, Mahatma Gandhi University, Kerala, India, 2000. [Google Scholar]
  6. Aleroev, T.S.; Aleroeva, H.T. The Sturm-Liouville problem for a second-order differential equation with fractional derivatives in the lower terms. Russ. Math. (Izvestiya VUZ. Matematika) 2014, 58, 1–9. [Google Scholar]
  7. Aleroev, T.S.; Aleroeva, H.T. On a class of non-self-adjoint operators that accompany differential equations of fractional order. Ukr. Mat. Bull. 2015, 12, 293–310. [Google Scholar]
  8. Popov, A.Y. On the number of real eigenvalues of the boundary-value problem for a second-order equation with fractional derivative. Fundam. Prikl. Mat. 2006, 12, 137–155. [Google Scholar] [CrossRef]
  9. Pskhu, A.V. Partial Differential Equations of Fractional Order; Nauka: Moscow, Russia, 2005; p. 200. ISBN 5-02-033721-8. [Google Scholar]
  10. Gokhberg, I.T.; Krein, M.G. Introduction to the Theory of Linear Nonself-Adjoint Operators; American Mathematical Society: Providence, RI, USA, 1969. [Google Scholar]
  11. Aleroev, T.S.; Aleroeva, H.T. A problem on the zeros of the Mittag-Leffler function and the spectrum of a fractional-order differential operator. Electron. J. Qual. Theory Differ. Equ. 2009, 25, 1–18. [Google Scholar] [CrossRef]
  12. Aleroev, T.S.; Kirane, M.; Tang, Y.-F. Boundary-value problems for differential equations of fractional order. J. Math. Sci. 2013, 10, 158–175. [Google Scholar] [CrossRef] [Green Version]
  13. Gokhberg, I.T.; Krein, M.G. Theory and Applications of Volterra Operators in Hilbert Space; American Mathematical Society: Providence, RI, USA, 1970. [Google Scholar]
  14. Aleroev, T.; Aleroeva, H. Problems of Sturm–Liouville type for differential equations with fractional derivatives. In Fractional Differential Equations; Kochubei, A., Luchko, Y., Eds.; Walter de Gruyter GmbH: Berlin, Germany; Munich, Germany; Boston, MA, USA, 2019; pp. 21–46. [Google Scholar]
  15. Aleroev, T.S. On a class of operators connected with differential equations of fractional-order. Sib. Math. Z. 2005, 46, 1201–1207. [Google Scholar]
  16. Aleroev, T.S. On the completeness of eigenfunctions of a differential fractional-order operator. Differ. Uravn. 2000, 36, 829–830. [Google Scholar] [CrossRef]
  17. Kato, T. Perturbation Theory for Linear Operators; Mir: Moscow, Russia, 1972. [Google Scholar]
  18. Loginov, B.V. To the estimate of the exactness of the method of perturbations. Izv. Akad. Nauk UzbSSR. Ser. Fiz.-Mat. Nauk 1963, 6, 14–19. [Google Scholar]
  19. Kilbas, A.; Srivastava, H.; Trujillo, J. Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 2006. [Google Scholar]
  20. Podlubny, I. Fractional Differential Equations; Academic Press: New York, NY, USA, 1999. [Google Scholar]
  21. Aleroev, T.S.; Aleroeva, H.T.; Huang, J.; Tamm, M.V.; Tang, Y.; Zhao, Y. Boundary value problems of fractional Fokker–Planck equations. Comput. Math. Appl. Adv. Fract. Differ. Equ. (IV) Time-Fract. PDEs 2017, 73, 959–969. [Google Scholar] [CrossRef]
  22. Ali, M.; Aziz, S.; Malik, S. Inverse source problems for a space–time fractional differential equation. Inverse Probl. Sci. Eng. 2020, 28, 47–68. [Google Scholar] [CrossRef]
  23. Aleroev, T.; Erokhin, S. The problem of finding the radon flux density by its concentration at different depths of the earth’s surface. IOP Conf. Ser. J. Phys. 2020, 1425, 012152. [Google Scholar] [CrossRef]
  24. Yan, S.; Cui, M. Finite difference scheme for the time-fractional Fokker–Planck equation with time- and space-dependent forcing. Int. J. Comput. Math. 2019, 96, 379–398. [Google Scholar] [CrossRef]
  25. Aleroev, T.S.; Kekharsaeva, E.R. Boundary value problems for differential equations with fractional derivatives. Integral Transform. Spec. Funct. 2017, 28, 900–908. [Google Scholar] [CrossRef]
  26. Aleroev, T.S.; Aleroeva, H.T. On the Eigenfunctions and Eigenvalues of a Class of Non-Selfadjoint Operators. Lobachevskii J. Math. 2016, 37, 227–230. [Google Scholar] [CrossRef]
  27. Aleroev, T.S.; Aleroeva, H.T. Erratum to: “On the Eigenfunctions and Eigenvalues of a Class of Non-Selfadjoint Operators”. Lobachevskii J. Math. 2016, 37, 815. [Google Scholar] [CrossRef]
  28. Larionov, E.A.; Zveriaev, E.M.; Aleroev, T.S. On Theory of Normal Operators Weak Perturbations; Keldysh Institute Preprints: Moscow, Russia, 2014; Volume 14, p. 31. [Google Scholar]
  29. Larionov, E.A. On stability of bases in Hilbert spaces. Eurasian Math. J. 2020, 11, 65–71. [Google Scholar] [CrossRef]
  30. Tseytlin, A.I. Applied Methods for Boundary Value Problems Solving in Structural Mechanics; Stroyizda Publication: Moscow, Russia, 1984. [Google Scholar]
  31. Erokhin, S.V.; Aleroev, T.S.; Frishter, L.Y. Sturm-Liouville problem for the equationviscoelastic damping oscillator. Int. J. Comput. Civil Struct. Eng. 2015, 11, 77–81. [Google Scholar]
  32. Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Integrals and Derivatives of Fractional Order and Some of Their Applications; Nauka Tekhnika: Minsk, Belarus, 1987. (In Russian) [Google Scholar]
  33. Bagley, R.L.; Torvik, P.J. A thoretical basis for the application of fractional calculus to viscoelasticity. J. Rheol. 1983, 27, 201–203. [Google Scholar] [CrossRef]
  34. Bagley, R.L.; Torvik, P.J. Fractional calculus—A different approach to the analysis of viscoelastically damped structures. J. AIAA 1983, 21, 741–748. [Google Scholar] [CrossRef]
  35. Kekharsaeva, E.R.; Pirozhkov, V.G. Modeling changes in the deformation-strength characteristics of asphalt concrete under loading using fractional calculus. In Proceedings of the 6th All-Russia Conference with International Participation on Mechanics of Composite Materials and Structures, Complex and Heterogeneous Media (IPRIM RAN), Moscow, Russia, 23–27 November 2016; Obraztsov, I.F., Yanovskii, Y.G., Eds.; IPRIM RAN: Moscow, Russia, 2016; pp. 104–109. [Google Scholar]
  36. Vasques, C.M.A.; Moreira, R.A.S.; Dias Rodrigues, J. Experimental identification of GHM and ADF parameters for viscoelastic damping modeling. In Proceedings of the 3rd European Conference on Computational Mechanics Solids, Structures and Coupled Problems in Engineering, Lisbon, Portugal, 5–8 June 2006. [Google Scholar]
  37. Aleroev, T.; Erokhin, S.; Kekharsaeva, E. Modeling of deformation-strength characteristics of polymer concrete using fractional calculus. IOP Conf. Ser. Mater. Sci. Eng. 2018, 365, 032004. [Google Scholar] [CrossRef] [Green Version]
  38. Ditkin, V.A.; Prudnikov, A.P. Integral Transforms and Operational Calculus; Nauka: Moscow, Russia, 1974; Pergamon: Oxford, UK, 1965. [Google Scholar]
  39. Kekharsaeva, E.R.; Aleroev, T.S. Model of deformation-strength characteristics of chlorine-containing polyesters based on fractional derivatives. Plast. Massy 2001, 3, 35. [Google Scholar]
  40. Man’kovskii, V.A.; Sapunov, V.T. Nomographic properties of fractional exponential E-functions when describing linear viscoelasticity. Zavod. Lab. Diagn. Mater. 2000, 66, 45–50. [Google Scholar]
  41. Aleroeva, H.; Aleroev, T. Some applications of fractional calculus. IOP Conf. Ser. Mater. Sci. Eng. 2020. [Google Scholar] [CrossRef]
  42. Zhang, J.; Huang, J.; Aleroev, T.; Tang, Y. A Linearized ADI for Two-Dimensional Time-Space Fractional Nonlinear Vibration Equations. Commun. Nonlinear Sci. Numer. Simul. 2020, in press. [Google Scholar]
  43. Zhang, J.; Aleroev, T.; Tang, Y.; Huang, J. Numerical Schemes for Time-Space Fractional Vibration Equations. Adv. Appl. Math. Mech. 2020, in press. [Google Scholar]
Figure 1. Graphs of solutions when 1 < α < 2 .
Figure 1. Graphs of solutions when 1 < α < 2 .
Mathematics 08 01877 g001
Figure 2. On the question of the sensitivity of the solutions of Problem (31) to errors of the parameter.
Figure 2. On the question of the sensitivity of the solutions of Problem (31) to errors of the parameter.
Mathematics 08 01877 g002
Figure 3. Determination of function by Formula (39) in Mathcad and comparison of this function with solution obtained using Laplace transform.
Figure 3. Determination of function by Formula (39) in Mathcad and comparison of this function with solution obtained using Laplace transform.
Mathematics 08 01877 g003
Figure 4. Deviation function (11) with parametric identification.
Figure 4. Deviation function (11) with parametric identification.
Mathematics 08 01877 g004
Figure 5. Comparison of experimental data with theoretical curve.
Figure 5. Comparison of experimental data with theoretical curve.
Mathematics 08 01877 g005
Table 1. Experimental points for polymer concrete samples.
Table 1. Experimental points for polymer concrete samples.
x i ( c ) 0.250.50.7511.251.5
U i 0.05−0.04−0.010.02−0.01−0.01
Table 2. Numerical results for eigenvalues of the problem (41)–(44).
Table 2. Numerical results for eigenvalues of the problem (41)–(44).
λ 1 λ 2 λ 3 λ 4 λ 5
16.659.4125.0213.4323.4
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Aleroev, T. Solving the Boundary Value Problems for Differential Equations with Fractional Derivatives by the Method of Separation of Variables. Mathematics 2020, 8, 1877. https://doi.org/10.3390/math8111877

AMA Style

Aleroev T. Solving the Boundary Value Problems for Differential Equations with Fractional Derivatives by the Method of Separation of Variables. Mathematics. 2020; 8(11):1877. https://doi.org/10.3390/math8111877

Chicago/Turabian Style

Aleroev, Temirkhan. 2020. "Solving the Boundary Value Problems for Differential Equations with Fractional Derivatives by the Method of Separation of Variables" Mathematics 8, no. 11: 1877. https://doi.org/10.3390/math8111877

APA Style

Aleroev, T. (2020). Solving the Boundary Value Problems for Differential Equations with Fractional Derivatives by the Method of Separation of Variables. Mathematics, 8(11), 1877. https://doi.org/10.3390/math8111877

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop