Existence of Bounded Solutions to a Modified Version of the Bagley–Torvik Equation

This manuscript reanalyses the Bagley–Torvik equation (BTE). The Riemann–Liouville fractional differential equation (FDE), formulated by R. L. Bagley and P. J. Torvik in 1984, models the vertical motion of a thin plate immersed in a Newtonian fluid, which is held by a spring. From this model, we can derive an FDE for the particular case of lacking the spring. Here, we find conditions for the source term ensuring that the solutions to the equation of the motion are bounded, which has a clear physical meaning.


Introduction
It is presently recognized that fractional calculus (FC) has a vast horizon of possible applications [1][2][3][4][5][6][7][8]. This mathematical tool is useful in many different fields, and we can cite a few examples such as: modeling the motion of a projectile [9], the study of oscillators [10], the behavior of diodes and memristors [11,12], viscoelasticity [13], robotics [14], or even issues regarding special relativity [15]. Several applications to the dynamics of particles can be found in [16]. This interconnection between FC and real models is natural in the framework of fluid dynamics [17]. The aim of this manuscript is to study in detail some aspects regarding the Bagley-Torvik equation (BTE). This fractional differential equation (FDE) models the vertical motion of a thin large plate immersed in a Newtonian fluid and hung by a spring, under the action of an external force (or source term). In particular, this paper addresses the issue of finding conditions on the source term that ensure that the solution is bounded.
The paper is organized into four sections. In Section 2, after introducing some fundamental notions regarding FC and special functions, we revisit the main ideas behind the deduction of the BTE [17,18]. In Section 3, we consider the particular case of lacking the spring. We compute the solutions and study their stability depending on the initial values. Moreover, we describe the results for a particular example of the source term. Finally, in Section 4, we find conditions on the source term ensuring that the motion of the plate is bounded. This result is obtained after describing the solution to our problem as a result of applying a fractional integral operator to the solution of an ordinary differential equation (ODE).

Fundamental Notions about Fractional Calculus and the Bagley-Torvik Equation
In this section, we recall some fundamental notions concerning FC, which will be used in the rest of the document. Besides, we detail some definitions and identities of special functions, since they will be relevant for the calculations in the follow-up. Finally, we reproduce the deduction of the BTE in [17,18], to derive afterwards a model for the particular case of lacking the spring.

Fundamental Notions of Fractional Calculus
In this subsection, we give the essential notions of FC. We recall that a function f is said to be absolutely continuous of order n ∈ N in [0, b] if, and only if, f is n times differentiable in the weak sense of the fundamental theorem of calculus; see [19].
In other words, f can be obtained as the n-iterated integral of a certain function g ∈ L 1 [0, b] plus a polynomial of degree strictly lower than n. Mathematically, we can write:

Definition 1.
We define the left Riemann-Liouville fractional integral of order γ > 0 of a function f ∈ L 1 [0, b] with base point zero as: for almost every t ∈ [0, b]. In the case that γ = 0, we just define I γ 0 + as: We define the Riemann-Liouville fractional derivative of order γ and base point zero of f as: where the last derivative is well defined due to the condition I It is important to have in mind a result that allows the fractional integration of a functional series term by term. In this sense, we can find the following result as Lemma 15.3 in [19]. Lemma 1. Given an analytic function g(t) defined in an open interval (0, b) and expressed as: where any α j > −1, it is possible to obtain the Riemann-Liouville fractional integral of g with base point zero after integrating term by term. Thus, we have:

Short Survey of Special Functions
We describe briefly some relevant special functions for our computations, together with some identities involving them [18][19][20][21]. The functions that we are going to present appear when computing fractional integrals or derivatives of trigonometric or exponential functions. A wide survey of these formulae can be found in [22].
The Mittag-Leffler functions of two parameters α and β are defined as: There are some common identities such as, However, the interest in this family of functions is given by results that link them with the fractional integrals of well-known functions. For instance, we can write ( [19]): We will also use identities that involve the error function, which is defined as ( [19]): In particular, we have the following expression, labeled as (2.7) in [23], Finally, we introduce the Fresnel integrals (see [19] or [20]), which are defined via the identities: These functions are important since they appear in the calculation of the fractional integrals of the trigonometric functions for half-integer orders. We can find in [20], page 49, the following formulae:

Deduction of the Bagley-Torvik Equation
In this subsection, we revisit the deduction of the BTE following the ideas presented in [17,18]. Let us suppose that we have an infinite and thin table that is moving horizontally on the surface of an infinitely wide and infinitely deep Newtonian fluid. Of course, the infinities that we are considering are, in a real-world application, just finite, but large enough. Our ideal two-dimensional table produces a movement in the fluid. Due to symmetry, the motion of any particle of the fluid at time t depends only on the depth z < 0. Moreover, the velocity function v(t, z) can be described through the diffusion equation. Thus, we have the following partial differential equation (PDE): where ρ and µ are the density and the viscosity of the fluid, respectively. We assume that the fluid is initially in equilibrium and that the effect of the movement tends to vanish when going deeper: Moreover, we assume that the velocity of the table v(t, 0), which coincides with the velocity of the upper layer of the fluid, is a bounded and continuous function. Of course, the function v(t, 0) is a given datum. Besides, in a Newtonian fluid, the shear stress σ(t, z) is given by: An argument involving the Laplace transform in the variable t turns the PDE (8) into a solvable ODE. Besides, from the solution to the ODE and Equation (10), we derive the property: where L denotes the Laplace transform andt stands for the variable in the transformed domain. We highlight that, in order to develop these deductions, it is necessary to consider v(0, z) = 0 when transforming the PDE into the ODE. From (11) and the solution to the ODE, we get: Moreover, if we denote by u(t, z) the accumulated horizontal displacement at time t of the layer at depth z, then we know that: Hence, we can use the fundamental theorem of calculus to derive: which allows us to rewrite (12) as: Let us imagine now that: (i) the table is hung vertically by a spring of constant K > 0; (ii) we introduce it into the Newtonian fluid; and (iii) we measure its vertical displacement as z(t).
By virtue of (14), the friction force will be given by 2 · S · √ µ ρ D where S is the area of the table. The motion of the table will obey the third Newton's law, and if we denote the external forces by f (t), then we have: where z(0) = z (0) = 0. This last expression is called the Bagley-Torvik equation (BTE).

Stability for the Motion of a Plate Immersed in a Newtonian Fluid
In the scope of the BTE model introduced in Section 2.3, if there is no spring, we obtain the equation: where the initial values z(0) = 0 and z (0) = 0 have to be trivial, in order to ensure that the deduction of the equation is correct. However, we can ask what happens if this is not the case for such initial conditions. For instance, let us suppose that z(0) = a. In this case, we can make the following remark. It is clear that the motion of the table should not depend on the initial point, since the choice for the origin of the coordinates is arbitrary. Hence, if the table fulfills z(0) = a, then it is clear that the equation must be the following: with initial values z(0) = a and z (0) = 0. Let us observe how the function z(t) − a behaves under the hypotheses ensuring that the deduction of the BTE still holds. In standard calculus, the translations of coordinates do not alter the differential equations. However, in our case, the previous expression can be rewritten as: with initial values z(0) = a and z (0) = 0. We verify that Equations (16) and (17) are very similar, but a "fictitious force"f (t) appears in the right-hand side of (17). Once the parameter b is fixed, Equation (17) has exactly one solution for each value of a. We will study first the associated homogeneous equation: so that, afterwards, we derive the results for the non-homogeneous case.

Assimilable Cases
Initially, we make a small digression, recalling equations similar to (18), for which we already know this limit behavior.

The Classical Pure Harmonic Oscillator
We can consider the following equation, where we substituted the fractional derivative by an identity term: This equation gives the solution to the well-known harmonic oscillator problem. Since this expression is, indeed, an ODE, there is a well-known answer to the problem of the limit behavior.
In this case, the characteristic polynomial of Equation (19), X 2 + b, has no roots of the strictly negative real part for b > 0 (which is the only realistic physical situation). Hence, there is no pair of solutions that converge to each other when time grows. In fact, this could be proven by direct calculation of the solutions, which are linear combinations of sin √ bt and cos √ bt . We highlight that this is coherent with our intuition. The motion of two identical springs that only differ in the prescribed initial values (initial position and velocity) will never converge to be the same. The underlying physical reason is that we need some "dissipative term" to ensure this approximation. However, since all the solutions to the homogeneous problems are bounded, it is also true that a solution to the non-homogeneous problem is bounded if and only if all of the solutions to the non-homogeneous problem are bounded.

The Classical Damped Harmonic Oscillator
It is also standard to check that if we add the classical damping term and we consider the ODE for the difference between two solutions of the damped oscillator: D 2 y(t) + b D 1 y(t) + c y(t) = 0, for b, c > 0, then the solution y(t) goes to zero when time grows, independently of the prescribed initial values. We recall that we are not discarding that both solutions to the damped harmonic oscillator do not tend to zero when time grows (for instance, they can oscillate), since this behavior depends strongly on the source term. Nevertheless, it is sure that one solution will get closer to the other one, as time passes. The mathematical reason is that, now, both roots of the characteristic polynomial are negative (recall that, due to physical arguments, we have b, c > 0). However, when c = 0, we have a strictly negative root and the zero root for the characteristic polynomial. Again, the set of functions solving the homogeneous problem is bounded, but not uniformly bounded.

The Immersed Plate in a Fluid
After the previous brief digression though the standard models, we go back to Equation (18): We also state again that, if we want the solution to have physical meaning, then it is necessary to impose the condition y (0) = 0. However, as we shall see in the follow-up, the study of the rest of the solutions has interest in itself, since it will be useful for several calculations.
We want to obtain the general solution for (18), but we will split the calculation into two parts: • Solutions fulfilling y(0) = 0. If we find a non-trivial solution to (18), the rest of the solutions fulfilling y(0) = 0 are obtained after a multiplication by a scalar. Initially, we will make a formal guess via a fractional Taylor expansion. Due to Lemma 1, the "guess" will be a solution to our problem. It seems natural to look for a solution of the form: Since we are looking for a particular solution, we can assign the value a 0 = y (0) = 1. Moreover, since D 2 y(t) + b D 3 2 0 + y(t) = 0, we have the recurrence relations a j+1 = −b · a j , for each j ∈ N. Hence, we conclude that a j = (−1) j · b j . Now, the question is if we can we rewrite: in terms of some well-known functions. The first idea is to separate the integer and the fractional powers, that is: On the one hand, the first summand is easy to recognize, since if we multiply it by b 2 and add one, then we get e b 2 ·t . Thus, we can write: In order to calculate the value of the other series in (20), we recall that from (4), it is straightforward to derive: On the other hand, we see that the second summand in Equation (20) can be rewritten as: By virtue of the identity (21), we get that (22) equals: Hence, we rewrite (20) as: Moreover, we can obtain the value of the limit when t → ∞, If we find a particular solution that is non-trivial, then the solution for the generic condition is obtained after multiplication with parameter a. We look for a solution of the type: Analogously to the previous case, we assume that a 0 = y(0) = 1, since we are looking for a particular solution. Moreover, since D 2 y(t) + b D 3 2 0 + y(t) = 0, we have again the same recurrence relations as before, namely a j+1 = −b · a j , for each j ∈ N. Hence, the solution in this case will be almost identical to the previous one, but with two main differences: • Now, we have coefficients for the degrees zero and 1 2 .

•
The rest of the coefficients coincide with those of the previous case, but after a multiplication by b 2 .
Consequently, we multiply (20) by b 2 , and after adding the two remaining terms of orders zero and 1 2 , we get: Finally, we apply L'Hôpital's rule to obtain the limit when t goes to infinity:

An Example
We will represent the previous results for a particular source term, f (t) = cos(t), and the coefficient b = 1. Of course, the first task is to describe the solutions to: We use the identity: Since we also know that: we obtain a particular solution to (25), given by: Moreover, due to the Taylor expansion of the cosine function, it is immediate to check that this solution fulfills the conditions z(0) = − 1 2 and z (0) = +∞. The particular solution to (25) can be written as: where we recall that C and S are defined via the Fresnel integrals. Moreover, we have the local expansion: Because of (23) and (24), the general solution to (25) with the zero source term can be described as: Besides, z(t) = z P (t) + z GH (t). In other words, the general solution to (25) can be described as the sum of (26) and (27). Due to (23), (24), and (27), we derive the following relations: In this case, z (0) = +∞.
Thus, the only solution z(t) that represents our model with physical meaning is the one associated with β = 1 2 and α = 0, since it is the only possibility to have z (0) = 0. The rest of the solutions do not have a clear physical interpretation, since they have a non-zero derivative at the origin (either finite or infinite). We represent some of the solutions in Figure 1, namely for: The feasible solution to Model (25) is the one described with a blue line. In Figure 1, we see how the three solutions associated with β = 0 (red α = 2, yellow α = 0, purple α = −2) get closer to each other with the passing of time. The same phenomena happens when β = 1 2 (green α = 2, blue α = 0, gray α = −2). Finally, we see that any two solutions associated with different values of β get more distant as time passes.

Sufficient Conditions for Bounded Solutions
We have seen (recall the behavior of the blue line in Figure 1) that the source term f (t) = cos(t) does not provide a bounded solution for a realistic case. It is interesting to provide sufficient conditions on f (t), in order to ensure that the solution with physical meaning is bounded. Observe that this property would mean that the motion of the plate is controlled from below and above. To show this property, it is sufficient to provide conditions on f ensuring that the solution z(t) has a finite limit when t → ∞.
Thus, we have the equation: for z(0) = 0, z (0) = 0, and b > 0, and we are looking for solutions fulfilling the property lim t→∞ z(t) = 0. Of course, the answer will depend on the source function f (t). We will provide sufficient hypotheses on f that ensure the previously mentioned behavior. First, we note the following identity: Thus, if we compute the unique solution to: such that x(0) = 0 and x (0) = 0, then we have: Thus, we only need to provide a condition on f ensuring that the unique solution to (29) with trivial initial values tends to zero when time grows and that this convergence to zero is fast enough, in order not to provide unbounded solutions after the fractional integration step in (30).
For instance, it is sufficient to have that x(t) is bounded and |x(t)| < k · t − 1 2 for some k ∈ R + and t big enough. Due to the classical theory of ODE, we know that the unique solution to Equation (29) with trivial initial values is: x(t) = t 0 r 0 e b 2 (r−s) · f (s) ds dr.
Thus, it will be enough to prove the finiteness of the limit: Note that it is mandatory to assume that the limit on the numerator is zero. In this case, we can apply L'Hôpital's rule and study the limit: Again, we assume that the limit on the numerator is zero. After differentiation in the numerator and denominator and canceling out the exponential factors, we get that (32) equals: Hence, it is enough to ensure that | f (t)| ≤ k · t − 3 2 , for some k ∈ R + . Thus, we have the following theorem. 3. | f (t)| ≤ k · t − 3 2 .
Author Contributions: All authors contributed equally to every part of this work. All authors read and agreed to the published version of the manuscript.
Funding: This research received no external funding.

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