INVERSE PROBLEM FOR THE TIME-FRACTIONAL EULER-BERNOULLI BEAM EQUATION

In this paper, the classical Euler-Bernoulli beam equation is considered by utilizing fractional calculus. Such an equation is called the time-fractional EulerBernoulli beam equation. The problem of determining the time-dependent coefficient for the fractional Euler-Bernoulli beam equation with homogeneous boundary conditions and an additional measurement is considered, and the existence and uniqueness theorem of the solution is proved by means of the contraction principle on a sufficiently small time interval. Numerical experiments are also provided to verify the theoretical findings.


Introduction
Fractional derivatives and integrals are very important in many areas of mathematics, physics and nano-engineering [11,13,19], and are widely used to capture natural and physical phenomena which cannot be predicted by classical integral and differential models. Consider a partial differential equation (PDE) with a fractional derivative in time t ∂ α t w(x, t) + w xxxx (x, t) = F (x, t; a, w), (x, t) ∈ D T , (1.1) where D T = {(x, t) : 0 < x < , 0 < t < T }, F (x, t; a, w) = a(t)w(x, t)+f (x, t), and ∂ α t is the left sided Caputo fractional derivative of order α satisfying 1 < α < 2. Here the parameter α represents the material damage [4] or the damage order index [22]. The so-called left sided Caputo fractional derivative of order α on the interval (0, t) is defined by provided that Γ (·) is the Gamma function.
It is important to emphasize that the Equation (1.1) becomes known models for the other values of α, i.e. 0 < α ≤ 1 and α = 2. In the case of 0 < α ≤ 1 the Equation (1.1) is a model of fractional sub-diffusion equation [6,8]. In particular, Equation (1.1) is a well-known model of thermal grooving by surface diffusion for α = 1, see [1,3,15]. In the case of α = 2, the Equation (1.1) is the simplest model for the transverse motion of a beam called the Euler-Bernoulli beam model. Vertical motion of each point x along the beam, with vertical displacement at time t is represented by w(x, t). Sound wave distribution problems, vibration, buckling and dynamic behavior of various building elements widely used in nano-technology are modeled with classical Euler-Bernoulli beam [7,16]. Therefore, in this paper we will call Equation (1.1) the fractional Euler-Bernoulli beam equation.
We consider the Equation (1.1) with the following initial and boundary conditions w(x, 0) =ϕ(x), w t (x, 0) = ψ(x), 0 ≤ x ≤ , (1.2) w( , t) =w xx ( , t) = w x (0, t) = w xxx (0, t) = 0, 0 ≤ t ≤ T. (1. 3) The physical interpretations of the boundary conditions are as follows: w( , t) = 0 means that the right end of the beam is held in a fixed position, w x (0, t) = 0 means that the beam centerline is perpendicular to the wall, w xx ( , t) = 0 means that there is no curvature at the right end of the beam, and w xxx (0, t) = 0 means that the curvature at the left end is not changing, i.e., there is no external torque.
We refer to several works on the mathematical and physical treatments for (1.1). A nonlinear fractional Euler-Bernoulli beam model was established in [2] to investigate the size-dependent, geometrically nonlinear free vibration of the fractional viscoelastic nano beams. In [21], the classical Euler-Bernoulli beam theory is reformulated utilizing fractional calculus. In [23], the authors presented a theoretical analysis of the free axial vibrations of rods described in terms of the fractional continuum mechanics. The natural frequencies and modal shapes for clamped rods are obtained, and the effects of the derivative order, as well as the corresponding length-scale parameters in the fractional model are discussed. Exact linking relationships between the elastic Euler-Bernoulli beam response and the fractional viscoelastic Timoshenko beam response is introduced in [17].
It is important to note that the non-local action of the fractional differential operator can operate on different spaces. The non-locality in time is called the short memory [18]. The short memory principle takes account the behaviour of w(x, t) only in the recent past, i.e., in the interval [t − L, t] where L is the memory length. That is, According to the short memory principle, the fractional derivative with the lower limit 0 is approximated by the fractional derivative with moving lower limit t − L. For instance, in [22], a hyper-elastic fractional damage material model with memory is considered with the application of fractional calculus. Based on the great amount of applications, we believe that the inverse problems involving the short memory principle need to be considered in the future work.
Our aim is to find the pair of solution {a(t), w(x, t)} from the problem (1.1)-(1.3) with the additional condition Here, the extra measurement (1.4) gives us the vertical displacement at x = 0. In contrast to the fractional Euler-Bernoulli beam equation, the inverse coefficient/source problems for the classical Euler-Bernoulli beam equation (i.e., α = 2 in Equation (1.1)) have been widely considered in the literature. The following papers are some important studies of the inverse problems for the classical Euler-Bernoulli equation. The inverse problem for the determination of the time-dependent force function in the Euler-Bernoulli beam equation, with the periodic boundary condition and an additional integral condition, was investigated in [9]. The inverse problem for finding the time-dependent potential for the Euler-Bernoulli beam equation was studied in [24]. In [14], a special technique was introduced to identify of the solution and the unknown coefficient of the Euler-Bernoulli equation. For more details, please see [5] and [12].
Some numerical aspects of the inverse problems involving the classical Euler-Bernoulli beam were investigated in [9,12,14]. In particular, the authors in [9] proposed a direct numerical method based on the finite difference method, to solve an inverse problem for the Euler-Bernoulli equation with a periodic boundary condition. In [14], an iterative finite difference scheme was introduced to solve the coefficient identification problem. In [12], numerical results were obtained by a regularization algorithm. In this paper, we focus on the timefractional Euler-Bernoulli beam equation which is a more general model than the classical Euler-Bernoulli beam equation in the aforementioned work.
The objective of this paper is to determine the time-dependent coefficient in the fractional Euler-Bernoulli equation using the additional information in (1.4), and prove the existence and uniqueness theorem for small T by means of the contraction principle. In addition, we would like to design a simple direct numerical method for solving the inverse problem.
The article is organized as following. In Section 2, we present some preliminaries used in Section 3. In Section 3, we transform the inverse problem (1.1)-(1.4) to a fixed-point system and prove the existence and uniqueness of a solution on a sufficiently small time interval by means of the contraction princi-ple. In Section 4, we introduce the numerical methods for the inverse problem and discuss the numerical results.

Preliminaries
Throughout this paper, we use the following definition and lemmas: ]. The generalized Mittag-Leffler function is defined by Lemma 1 [ [18]]. Let 0 < α < 2, and β ∈ R be arbitrary. We suppose that µ is such that πα 2 < µ < min {π, πα}. Then there exists a constant C α,β such that where F(s) is the Laplace transform of f (t).

Solution of the inverse problem
In this section, we will examine the existence and uniqueness of the solution of the inverse initial-boundary value problem for the Equation (1.1) with timedependent coefficient.
. Since a is solely time-dependent and the boundary conditions are homogeneous, we attempt to apply the Fourier method of eigenfunction expansion to the problem (1.1)-(1.4). Consider first the auxiliary spectral problem given by This spectral problem has the eigenvalues λ n = µ 4 n and the eigenfunctions X n (x) = 1 √ cos(µ n x), where µ n = (2n+1)π 2 , n = 0, 1, 2, .... In the rest of the paper, we will consider the following spaces to investigate the inverse problem (1.1)-(1.4): Note that E 5 T is also a Banach space. Let us seek the solution of the problem (1. where w n (t) is the solution of the following initial value problem: ... Applying the Laplace transform on both sides of (3.2) leads to By using the inverse Laplace transform, we can obtain the solution of the Cauchy problems (3.2) as For the determination of a(t), we can derive from the Equation (1.1) with the additional data (1.4). Considering the Equation (3.4) at x = 0 in the last equation of a(t), we get where Φ = [φ 0 , φ 1 ] T and φ 1 and φ 0 are equal to the right hand sides of (3.4) and (3.5), respectively. That is, Let us show that Φ maps E 5 T onto itself continuously. In other words, T . We will use the following assumptions on the data of problem (1.1)-(1.4): By using integration by parts under the assumptions (A 0 )-(A 2 ), it easy to show that Under the assumptions (A 0 )-(A 2 ) and using Cauchy-Schwartz and Bessel inequalities, we obtain from (3.7) that where C α,i , i = 1, 2 are constants defined in Lemma 1. Since the right hand side of (3.8) is bounded, we have φ 0 (z) ∈ C[0, T ]. Next, let us show that φ 1 (z) ∈ B 5 T , i.e., we need to show that where φ 1n (t) is equal to the right hand side of w n (t) as in (3.3). After some manipulations on the last equality under the assumptions (A 0 )-(A 2 ), we get (3.9) From the Bessel inequality and [ ∞ n=0 (µ 5 n max 0≤t≤T |w n (t)|) 2 ] 1/2 < +∞, series on the right side of the (3.9) are convergent. Thus J T (φ 1 ) < +∞ and φ 1 (z) belongs to the space B 5 T . Now we show that Φ is a contraction mapping on E 5 T . Let z i = a i (t), w i (x, t) T for i = 1, 2 be any two elements of E 5 T . We know that Under the assumptions (A 0 )-(A 2 ) and considering (3.8)-(3.9), we obtain that 1 µ 2 n 1/2 and C(a 1 , w 2 ) is a constant that depends on the norms of a 1 (t) C[0,T ] and w 2 (x, t) B 5 T . A(T ) has limit zero as T tends to zero. It means that for sufficient small T , the operator Φ is contraction mapping which maps E 5 T onto itself continuously. Then, according to the Banach fixed point theorem there exists unique solution of (3.6). Thus, we proved the following theorem:

Numerical experiments
In this section, we will discuss the numerical results for the inverse initialboundary value problem (1.1)-(1.4) with a time-dependent coefficient.
Throughout this section, we use w(x, t) and W (x, t) to represent the exact and the numerical solution to the inverse problem (1.1)-(1.4), respectively.
Based on the discussion in Section 3, we know that the exact solution w(x, t) can be written as w(x, t) = ∞ n=0 w n (t)X n (x), where X n (x) = cos(µ n x) with µ n = (2n+1)π 2 , and w n (t) for n = 0, 1, 2, . . . are the Fourier coefficients. In order to obtain the numerical solution, we look for solution of the form W (x, t) = N n=0 w n (t)X n (x), which is the truncated Fourier series of w(x, t). Here (N +1) is the number of Fourier modes in the truncated Fourier series. Such a truncated series will converge to w(x, t) as N goes to infinity. However, in the numerical computations, we have to choose the suitable values of N and the time step size to ensure the numerical stability. We will define our numerical scheme before we specify the assumption of N . Equations (1.1)-(1.3) lead to ∂ α t w n (t) + µ 4 n w n (t) = a(t)w n (t) + f n (t), 0 ≤ t ≤ T, w n (0) = ϕ n , w n (0) = ψ n , n = 0, 1, . . . , N. (4.1) To solve the equations (4.1) numerically, we partition the temporal domain [0, T ] into uniform subintervals, each of which has length ∆t. We then let t m = m∆t for all m. Thus at t = t m , the first fractional differential equation in (4.1) can be approximated as for n = 0, 1, . . . , N . Here W i n denotes the numerical approximation of w n (t) at t = t i and d m, Note that the approximation of the time-fractional derivative in equation (4.2) is of first order [20]. The investigation of higher-order time discretization methods will be left as future work. Also, from equation (4.2), if we let α goes to 2, we can get W m+1 for n = 0, 1, . . . , N . To ensure the stability of the difference scheme above, we assume that (∆t) 2 µ 4 N is of order O(1) or smaller. That is, N must satisfy the assumption that ∆t ((2N + 1)π/(2l)) 2 is of order 1 or smaller. For example, if the spatial domain is [0, 1] and we take ∆t = 10 −3 , then we choose N to be less than 10.
Since d m,m−1 = 1, equation (4.2) can be rewritten as The scheme of the inverse initial-boundary value problem (1.1)-(1.4) with a time-dependent coefficient is given as follows. We first update W 0 m and W 1 m using the initial conditions. In particular, we compute W 0 n for n = 0, 1, . . . , N using the discrete cosine transformation of ϕ(x), and we compute W 1 m using the cosine transformation of ϕ(x) + ψ(x)∆t. In addition, a 0 is updated using In order to compute a 1 , we let x = 0 and t = t 1 in Equation (1.1). Then we have ∂ α t w(0, t 1 ) + w xxxx (0, t 1 ) = a 1 w(0, t 1 ) + f (0, t 1 ). Since w(0, t 1 ) = r(t 1 ), there is , which leads to Next, we can use Equation (4.3) with m = 1 to get the formulation of w 2 n . We then let x = 0 and t = t 2 in Equation (1.1), and use w(0, t 2 ) = r(t 2 ) to obtain the formulation of a 2 . For general m ≥ 1, we compute a m and W m+1 n (with n = 0, 1, . . . , N ) in an alternating order as follows: we first compute a m using Equation (4.4), and then compute W m+1 n for n = 0, 1, . . . , N using Equation (4.5), where and Note that the term w xxxx (0, t m ) in Equation (4.4) should be computed using N n=0 µ 4 n W m n . Finally, the numerical solution of w(x, t) at final time T can be computed using the inverse cosine transformation of W Nt n (for n = 0, 1, . . . , N ) where N t is the time level for t = T . Note that the numerical methods presented here are similar to that in [25]. Now, we discuss some numerical results. Without loss of generality, we choose the spatial domain to be [0, 1]. Example 1. In the first numerical example, we consider the inverse initialboundary value problem (1.1)-(1.4) with the following data: The exact solution to the inverse problem is w(x, t) = (1 + t 2 ) cos( πx 2 ) and a(t) = e −t . In our numerical simulations, we choose α = 1.9, ∆t = 1 × 10 −4 , T = 0.1. The numerical solution of a(t) for t ∈ [0, T ] and its error are shown in Figure 1. We can see that the error of a(t) increases as t increases. The numerical results also show that the the absolute maximum error of a is 6.0411×10 −5 , which indicates that our numerical method can recover the function a(t) accurately. The numerical solution and error of w(x, t) are given in Figure 2.
We observe that at each fixed time t, the maximum absolute error of w occurs at the left boundary x = 0. Due to the fact that we use truncated Fourier series, the error of w at the right boundary stays at zero. The numerical results show that the absolute maximum error of w(x, t) is 1.0023 × 10 −5 . Therefore, our numerical method lead to accurate recovery of both a(t) and w(x, t). In addition, when we further increase the value of α to let it approach 2, we observe that both the errors of a(t) and w(x, t) decrease. In particular, when we choose α = 1.99, then the absolute maximum errors of w(x, t) and a(t) are 1.0016 × 10 −5 and 6.0377 × 10 −5 , respectively. Such a result makes sense because when α approaches 2 from the left, the time-fractional Euler-Bernoulli beam equation approaches the classical Euler-Bernoulli beam equation, and our temporal discretization will approach a second-order accuracy scheme. The absolute maximum error of w(x, t) for different values of α can be seen in Figure 3.  We observe that the absolute maximum error of w(x, t) decreases monotonically as α increases. This is also consistent with the fact that the convergence order for the discretization of ∂ α t w n (t) increases as α increases. In order to see the influence of ∆t on the accuracy of the numerical solution, we fix the value of α and N (i.e., α = 1.9 and N = 4), choose various values of ∆t, and compute the absolute maximum errors of w(x, t) and a(t) for each choice of ∆t. When we use ∆t = 10 −3 , the absolute maximum errors of w(x, t) and a(t) are 1.0022 × 10 −4 and 6.0413 × 10 −4 , respectively. If we decrease the value of ∆t by a factor of 10, such errors decrease to 1.0023 × 10 −5 for w(x, t), and 6.0411 × 10 −5 for a(t). When we further decrease the value of ∆t to be 10 −5 , the absolute maximum error of w(x, t) becomes 1.0022 × 10 −6 , and the error of a(t) becomes 6.0378 × 10 −6 . This implies the convergence of our scheme as ∆t approach 0.
Example 2. Next, we consider the another numerical example, where we choose the following given data: The exact solution to the inverse initial-boundary value problem (1.1)-(1.4) with the data above is w(x, t) = cos( πx 2 ) + 1+t 3 100 cos( 3πx 2 ) and a(t) = sin(t). In this example, we compute the numerical solution when α = 1.8 and T = 0.1. The numerical solution and the error of a(t) are given in Figure 4, from which we can observe that the error of a(t) increases as t increases.  Numerical results show that the absolute maximum error of a(t) is 1.2587 × 10 −5 . The numerical solution and the error of w(x, t) are given in Figure 5. Unlike the distribution of error in Figure 2(b), the absolute maximum error of w(x, t) for this example is not located at the boundary. Numerical results show that the absolute maximum error of w(x, t) is 3.3290 × 10 −8 .
To see the influence of time step size, we choose different ∆t and compute the error of w and a. When we use ∆t = 10 −3 , the absolute error of w(x, t) and a(t) are 2.6788 × 10 −7 and 1.0221 × 10 −4 , respectively. When we use the time step size ∆t = 10 −5 , the absolute maximum error of w(x, t) and a(t) become 3.6648 × 10 −9 and 1.3462 × 10 −6 , respectively. Hence, the observation above suggests that if ∆t decreases within a certain range, then the error of w and a will also decrease. Next, we fix the time step size to be ∆t = 10 −4 and change the value of α. In particular, when α = 1.7, the absolute maximum error of w(x, t) and a(t) are 6.1008 × 10 −8 and 1.1300 × 10 −4 , respectively; when α = 1.9, the absolute maximum error of w(x, t) and a(t) are 1.6229 × 10 −8 and 6.5163 × 10 −6 , respectively. When we compare the results for α = 1.7, 1.8 and 1.9, we observe the pattern that larger value of α leads to more accurate results, which is consistent with the discussion in Example 1. When we fix the value of α and compare the results for ∆t = 10 −3 , 10 −4 and 10 −5 , we can also see the monotonically decreasing of the absolute maximum errors of w(x, t) and a(t), see Figure 6. In particular, when we choose α = 1.8, T = 0.1 and ∆t = 10 −5 , the absolute maximum errors of w(x, t) and a(t) are 1.8281 × 10 −9 and 7.8356 × 10 −7 , respectively. This example also verifies that our numerical methods can lead to accuracy solution to the inverse problem for small T .

Conclusions
The paper considers the problem of determining the time-dependent coefficient for the fractional Euler-Bernoulli equation with homogeneous boundary conditions and an additional measurement. The existence and uniqueness of a solution on a sufficiently small time interval are proved by means of the contraction principle. The fixed-point system is presented via Fourier series. Such a form of the system brings along computations that are technically simpler than the system in the case of the usual variational approach. A numerical method using the truncated Fourier series is proposed. The numerical results show that our method lead to accurate solutions to the inverse problem for small time interval, and the absolute maximum error of the solutions decreases as the order of the time-fractional derivative approaches 2 from the left. Higher-order temporal discretization methods for the inverse problem are currently under investigation.