Sturm-Liouville Problem with Mixed Boundary Conditions for a Differential Equation with a Fractional Derivative and Its Application in Viscoelasticity Models

: In this study, we obtained a system of eigenfunctions and eigenvalues for the mixed homogeneous Sturm-Liouville problem of a second-order differential equation containing a fractional derivative operator. The fractional differentiation operator was considered according to two deﬁ-nitions: Gerasimov-Caputo and Riemann-Liouville-Visualizations of the system of eigenfunctions, the biorthogonal system, and the distribution of eigenvalues on the real axis were presented. The numerical behavior of eigenvalues was studied depending on the order of the fractional derivative for both deﬁnitions of the fractional derivative.


Introduction
Differential equations containing fractional derivatives are actively used in constructing mathematical models for various phenomena and processes, followed by studying the properties and behavior of the solutions of these equations. Such processes include research on continuous media with memory, filtration of fluids in fractal geometries, physical aspects of stochastic transport and diffusion, mathematical models of viscoelastic materials, models of damped oscillations with fractional damping (e.g., vibrations of rock formations during earthquakes or vibrations of nanoscale sensors), models of nonlocal physical processes and fractal phenomena, climate models, and so on. For a long time, the encyclopedia of fractional calculus was the monograph [1]. Currently, the most detailed publication, presenting the fundamental theoretical principles of fractional calculus and various examples of its practical applications, can be considered [2]. Monographs [3][4][5][6][7][8][9] are also devoted to the application of differential equations with fractional derivatives in constructing models of physical processes.
Especially intensively developed at present is the direction related to the use of differential equations containing fractional differentiation operators in constructing viscoelastic models of solid bodies (for example, [10,11]). A detailed review of viscoelastic models with fractional derivatives and their historical development is given in [12,13].
In addition to direct problems-the study of the properties of solutions of differential equations containing fractional derivatives-inverse problems are always relevantdetermining the type of equation that serves as the best mathematical model for the studied process. In particular, the problem of selecting the order of the fractional derivative included in differential equations used as rheological models is addressed in works [14][15][16][17][18][19][20][21][22][23][24].

Materials and Methods
Let us consider a differential equation containing a fractional-order operator used as a model for viscoelasticity (a Linear Fractionally Damped Oscillator). Such an equation 2 of 10 is used [25][26][27][28][29][30] for modeling the changes in deformational and strength characteristics of various viscoelastic materials. This second-order equation contains a fractional derivative of order less than two: Here w(t) represents the displacement of a fixed point (construction material) at a given time t; T-end of the time period under consideration; D α ot w(t)-represents the Gerasimov-Caputo or Riemann-Liouville fractional differentiation operator of order α.
c, λ, and α-constants. In [22], the following physical interpretation is proposed for these constants: c-represents the material's viscosity modulus, λ-represents the material's stiffness modulus, α-represents the viscoelasticity parameter of the material. The fractional differentiation operators are defined as follows [31][32][33]. The left-sided Gerasimov-Caputo fractional derivative: The left-sided Riemann-Liouville fractional derivative: Here Γ represents the Euler gamma function. We will consider the behavior of a fixed point in a viscoelastic material subjected to an impulsive load at the initial moment of time. Therefore, we impose the following mixed, homogeneous boundary conditions: which means that the displacement of the point is zero at the initial moment of time, and it is given an initial velocity upon impact. The velocity of the point becomes zero at the final moment of time (at the moment of cessation of oscillation). The solution to problems (1)-(2) with homogeneous first-type boundary conditions (Dirichlet problem) is provided in [34].
The general solution of Equation (1) can be found, as in [35], by solving a second-order Volterra equation using a sequence of recursive kernels and expressing it as a power series.

Results
The general solution of Equation (1) in the case of the Gerasimov-Caputo fractional derivative, in the form of a power series, is given by: Here n k = n! k!(n−k)! represents a binomial coefficient. Using the boundary conditions, it follows: from the first that B = 0 and from the second that the eigenvalues of the problems (1)- (2) are the zeros of the function Λ GC α (λ): (3) Figure 1 shows the graph of the function Λ GC α (λ) for various values of the fractional derivative order with T = 1 and c = 0.2. (3) Figure 1 shows the graph of the function Λ ( ) for various values of the fractional derivative order with T = 1 and c = 0.2.  Table 1 presents the first ten eigenvalues of the problems (1)-(2) for T = 1 and c = 0.2 obtained through numerical simulation using MATLAB. According to the values in Table 1, as the fractional derivative order increases, the eigenvalues increase slightly. The maximum relative deviation is observed in the second eigenvalue: The system of eigenfunctions for the problems (1)-(2) with the Gerasimov-Caputo fractional derivative is given by: represents the eigenvalues with = 1, 2, …  Table 1 presents the first ten eigenvalues of the problems (1)-(2) for T = 1 and c = 0.2 obtained through numerical simulation using MATLAB. Table 1. Eigenvalues of the problems (1)-(2) for T = 1 and c = 0.2 with various values of the fractional derivative order according to Gerasimov-Caputo. According to the values in Table 1, as the fractional derivative order increases, the eigenvalues increase slightly. The maximum relative deviation is observed in the second eigenvalue: max k=1,10 The system of eigenfunctions for the problems (1)-(2) with the Gerasimov-Caputo fractional derivative is given by: Axioms 2023, 12, 779 4 of 10 λ m represents the eigenvalues with m = 1, 2, . . . Thus, solutions to problem (1) satisfying boundary conditions (2) are found. That is, we obtain the equations of motion of a fixed point of a viscoelastic body to which a non-zero velocity is given at the initial moment (rest position) of time, and at the final moment of time T, the velocity is zero. Figure 2 illustrates the graphs of the first 6 eigenfunctions for the problems (1)-(2) with a fractional derivative order of α = 1.2 according to Gerasimov-Caputo, with T = 1 and c = 0.2. The values of functions W GC 1 (t) ÷ W GC 6 (t) were computed using numerical simulation in MATLAB, where a partial sum of the series was taken instead of the full sum.
Axioms 2023, 12, x FOR PEER REVIEW 4 of Thus, solutions to problem (1) satisfying boundary conditions (2) are found. That we obtain the equations of motion of a fixed point of a viscoelastic body to which a no zero velocity is given at the initial moment (rest position) of time, and at the final mome of time T, the velocity is zero. Figure 2 illustrates the graphs of the first 6 eigenfunctions for the problems (1)    The system of eigenfunctions ( ) , ,… is complete in , but not orthogon [34]. Let us construct a biorthogonal system for it using the approach from [36]. To obta  Thus, solutions to problem (1) satisfying boundary conditions (2) are found. That is, we obtain the equations of motion of a fixed point of a viscoelastic body to which a nonzero velocity is given at the initial moment (rest position) of time, and at the final moment of time T, the velocity is zero. Figure 2 illustrates the graphs of the first 6 eigenfunctions for the problems (1)    The system of eigenfunctions ( ) , ,… is complete in , but not orthogonal [34]. Let us construct a biorthogonal system for it using the approach from [36]. To obtain The system of eigenfunctions W GC m (t) m=1,2,... is complete in L 2 , but not orthogonal [34]. Let us construct a biorthogonal system for it using the approach from [36]. To obtain Axioms 2023, 12, 779 5 of 10 a biorthogonal system for the system W GC m (t) ∞ m=1 , we consider the adjoint equation to (1). For this purpose, we examine the operator: Then, the adjoint operator to A is denoted as A * , such that: where the scalar product is defined as: Assuming the functions f and g satisfy the boundary conditions (2), we have: Using the self-adjointness of the repeated differentiation operator, we have: The linearity of the fractional differentiation operator gives: By changing the integration limits in the double integral, we have Let us denote: then, using the boundary conditions, Furthermore, by using the relation between the right-sided fractional derivatives according to Gerasimov-Caputo and Riemann-Liouville derivatives, we obtain: Therefore, Thus, the adjoint operator takes the following form: A * (g) = g + c· GC D α tT . The adjoint boundary value problem on the interval t ∈ [0; T] for the Gerasimov-Caputo fractional derivative definition of (1)-(2) is given by: The general solution of Equation (5), obtained through the solution of the second-kind Volterra equation using a sequence of recurrent kernels, is given by: Substituting t = T and using (6), we have B = 0. Differentiating the obtained series, we have: Next, we use the second boundary condition: Hence, the eigenvalues coincide, and the eigenfunctions of the adjoint problem have the form: To prove the orthogonality of the systems W GC m (t) m=1;2;3;... and W GC m (t)  Figure 4 shows the graphs of ( ) ; -the first four eigenfunctions of the problems (1)-(2) in black, and the graphs of ( ) ; -the first four eigenfunctions of the adjoint problem (5)-(6) in gray for the fractional derivative with order α = 1.2 according to Gerasimov-Caputo with T = 1 and c = 0.2.  The solution of problems (1)-(2) is given by: The coefficients A m can be found using the standard approach, employing the inner product in L 2 : Thus, we have the representation: Similar calculations can be carried out if the fractional differentiation operator is defined according to Riemann-Liouville. In this case, the eigenvalues are zeros of the function: Table 2 presents the first ten eigenvalues of the problem (1)-(2) for T = 1, c = 0.2, and various values of the order of the fractional derivative according to Riemann-Liouville. According to the values in Table 2, as the order of the fractional derivative increases, the eigenvalues slightly increase. The maximum relative deviation is observed in the first eigenvalue: max k=1,10 It should be noted that according to the values from Tables 1 and 2, the eigenvalues corresponding to different fractional differentiation operators with the same order of the fractional derivative (specifically, α = 1.5) do not coincide. The maximum relative deviation is observed in the first eigenvalue and amounts to 11%. Next, the system of eigenfunctions of problems (1)-(2) in the case of defining the fractional differentiation operator according to Riemann-Liouville takes the form: In conclusion, let us consider the classical particular case: when α = 1, problems (1)-(2) transform into Sturm-Liouville problems for a second-order equation containing the first derivative: The eigenvalues λ n of problems (8)-(9) are zeros of the function: The eigenfunctions of problems (8)-(9) are given by: It should be noted that the eigenvalues calculated in the MATLAB software using the function Λ(λ) for the values T = 1, c ∈ [0.1; 1.0] with a step of 0.1 coincide with the eigenvalues calculated using the function Λ RL 1 (λ).

Discussion
In viscoelasticity models (see, for example, [37]), the quantity c 2 is referred to as the coefficient of relative damping (over one period τ n of the trigonometric component, the natural logarithm of the eigenfunction decreases c 2 τ n times ), c-is the coefficient of losses, and λ n − c 2 2 -is the (circular) frequency of damped eigenoscillations (the quantity divided by 2π is inversely proportional to the period τ n of the trigonometric component of the eigenfunction). When building mathematical models of viscoelasticity using the apparatus of fractional differentiation, it is clear that both the damping coefficient and the frequency of damped eigenoscillations should contain the order of the fractional derivative. In our view, the relative damping depends on two parameters, c and α, and the frequency of damped eigenoscillations on three parameters, λ m (eigenvalues), c and α.Therefore, the physical meaning of the coefficients in Equation (1) used in [22] appears to be justified by the authors.

Conclusions
In this study, a mixed homogeneous Sturm-Liouville problem for a second-order differential equation containing a fractional differentiation operator was analyzed, considering two definitions: Gerasimov-Caputo and Riemann-Liouville.
A system of eigenfunctions and eigenvalues for the considered mixed boundary value problem was obtained.
A conjugate boundary value problem for the fractional derivative according to Gerasimov-Caputo was constructed.
A system biorthogonal to the found system of eigenfunctions was constructed in the case when the fractional derivative is defined according to Gerasimov-Caputo. The numerical behavior of eigenvalues was studied depending on the order of the fractional derivative for both definitions of the fractional derivative.
Visualizations of the system of eigenfunctions, the biorthogonal system, and the distribution of eigenvalues on the real axis were presented.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.