Fractional derivative boundary control in coupled Euler-Bernoulli beams: stability and discrete energy decay

1.   Introduction

This article examines the solutions to the initial boundary value problem of connected Euler-Bernoulli beam equation, expressed as

accompanied by the following boundary conditions:

where $\gamma$ and $\beta$ are positive constants, and the following variables have specific engineering meanings: $\varphi$ and $\psi$ denote the vertical displacements, $\varphi_t$ and $\psi_t$ represent the velocities, $\varphi_x$ and $\psi_x$ signify the rotations, and $\varphi_{xt}$ and $\psi_{xt}$ indicate the angular velocities.

The system's initial conditions are specified as

where the initial data $(\varphi_{0}, \varphi_{1})$ and $(\psi_{0}, \psi_{1})$ belong to an appropriate function space.

The notation $\partial^{\alpha, \eta}_{t}$ denotes the tempered Caputo time-fractional derivative of order $\alpha$ with respect to time $t$, where $0 < \alpha < 1$, as introduced in ^[1,2,3,4]. It is defined as

where $\Gamma(\cdot)$ represents the Gamma function.

The tempered Caputo time-fractional derivative, a generalization of the classical Caputo derivative, incorporates an exponential tempering factor that allows for a finite memory horizon, making it particularly suitable for processes with decaying memory effects. These exponentially modified fractional integro-differential operators were first proposed by Choi and MacCamy in ^[5], presenting a slightly different approach from the classical Caputo derivative formulated by Michele Caputo in ^[6]. Unlike the standard Caputo derivative, the tempered version is especially suited to applications where the memory of past states diminishes over time, as observed in fields such as anomalous diffusion and viscoelasticity. This choice provides both improved model stability and greater accuracy in representing transient dynamics, offering a more realistic approach to systems with decaying temporal correlations.

A basic model that describes the transverse vibration of a system of non-homogeneous connected Euler-Bernoulli beams, as presented in ^[7], is represented by the following system:

While this model provides a foundational description of transverse vibrations in such systems, it primarily addresses the standard beam dynamics without incorporating advanced effects such as memory and complex dissipation. To enhance the model's applicability to real-world scenarios, we extend the classical Euler-Bernoulli framework by introducing fractional derivatives in the boundary conditions. This extension allows for a more accurate representation of systems exhibiting memory effects and intricate dissipative behavior. By integrating these fractional derivatives, we capture more realistic dissipative dynamics and gain deeper insights into the stability properties of the system.

In recent years, the scientific community has shown growing interest in unraveling the intricate dynamics and practical applications of wave equations. The behavior of waves, whether occurring naturally, such as seismic waves in the Earth's crust, or in engineered systems, such as acoustic waves in materials, captivates both researchers and practitioners alike. Considerable research effort has been devoted to investigating wave equations with diverse damping types and exploring their stability and controllability. These waves emerge when a vibrating source disrupts the surrounding medium. Researchers have shown a keen interest in addressing damping-related challenges, whether local or global, and have illustrated different forms of stability.

In ^[4], B. Mbodje investigated the asymptotic behavior of solutions using the system:

He established that the corresponding semigroup lacks exponential stability, showing only strong asymptotic stability. Additionally, the system's energy diminishes over time, approaching a decay proportional to $t^{-1}$ as time extends to infinity.

In ^[8], Akil and Wehbe investigated a multidimensional wave equation featuring boundary fractional damping applied to a portion of the domain's boundary:

They demonstrated the system's strong stability while establishing that it lacked uniform stability. In addition, they derived a polynomial energy decay for smooth solutions of the form $t^{-\frac{1}{1-\alpha}}$. This analysis assumed specific geometric conditions for the boundary control region and leveraged the exponential decay of the wave equation with standard damping.

Recently, in ^[9], Beniani et al. examined a system comprising coupled wave equations featuring a diffusive internal control of a general nature:

They showed the absence of exponential stability and explored the asymptotic stability of the model, establishing a general decay rate that depends on the density function $\varrho$. See also the references in ^{[7,10,11,12,13]} for further existing results related to the stability and numerical analysis of (1.1).

In this investigation, we extend the traditional Euler-Bernoulli beam model, commonly used to describe deformations under external loads, by incorporating fractional derivatives in the boundary conditions. This extension is essential, as fractional derivatives offer more precise models for systems exhibiting memory effects and complex dissipation phenomena frequently observed in practical applications such as the analysis of material microstructures. Incorporating fractional derivatives into the boundary conditions constitutes a novel advancement within the Euler-Bernoulli framework, enabling a more accurate representation of dissipative effects in diverse materials and structures. The findings from this study provide new insights into the stability characteristics of systems with fractional boundary dissipation, enhancing both theoretical understanding and the numerical experiments used to validate these observations.

The organization of the paper is outlined below: In Section 2, we transform the system (1.1) into an augmented system by integrating the wave equations with compatible diffusion equations. In Section 3, we derive global solutions to the problem (1.1) using arguments based on semigroup theory. In Section 4, we establish the strong stability of our system, despite the lack of resolvent compactness, by applying general criteria from Arendt and Batty. In Section 5, we propose a numerical scheme capable of replicating various decay rate profiles and confirm numerically the known stability outcomes related to the energy decay through several examples.

2.   Augmented model of the system

This section focuses on transforming the model (1.1) into an augmented system. The following theorem is required:

Theorem 2.1. ^[4] Consider $\mu$ be the function given by

Then the relationship between the "input" $U$ and the "output" $O$ of the system

and

can be described as follows:

where

Lemma 2.1. ^[4] If $\lambda \in D = \left\{\lambda\in \mathbb{C}/ \Re \lambda+\eta > 0\right\}\cup \left\{\Im \lambda\neq 0\right\}$, then

and

Based on the previous theorem, the system (1.1) can be expressed as the following augmented model:

where $\zeta = \gamma \pi^{-1} \sin\left(\alpha\pi\right)$. The energy associated with the solution to the system (2.1) is expressed by

where $\|\cdot\|_{2}$ refers to the $L^2(0, L)$ norm.

Then, we have the following lemma:

Lemma 2.2. Given that $\left(\varphi, \psi, \phi\right)$ is a solution to the problem (2.1), the energy functional defined by (2.2) fulfills

Proof. Multiplying Eqs $(2.1)_{1}$ and $(2.1)_{2}$ by $\varphi_{t}$ and $\psi_{t}$, respectively, using integration by parts over $(0, L)$ and Eq $(2.1)_{6}$, and combining the two equations, we derive

Scaling Eq $(2.1)_{3}$ by $\zeta \phi$ and integrating over $\mathbb{R}$ yields

By combining (2.4) and (2.5), we obtain (2.3). The proof of the lemma is thereby finished. □

3.   Problem formulation and well-posedness

Within this section, we establish the well-posedness of (2.1). For this purpose, we define the following Hilbert space (the energy space):

where

For $U = \left(\varphi, u, \psi, v, \phi \right)^{T}$ and $\overline{U} = \left(\overline{\varphi}, \overline{u}, \overline{\psi}, \overline{v}, \overline{\phi}\right)^{T}$, we introduce the inner product in the space $\mathcal{H}$ as follows:

We then reformulate (2.1) in the framework of semigroup theory: By introducing the vector function $U = \left(\varphi, u, \psi, v, \phi\right)^{T}$, the system (2.1) can be expressed as

where $U_{0}: = \left(\varphi_{0}, \varphi_{1}, \psi_{0}, \psi_{1}, \phi_{0}\right)^T$. The operator $\mathcal{A}$ is linear and specified as follows:

The domain of $\mathcal{A}$ is then

We state the following theorem on existence and uniqueness:

Theorem 3.1. (1) If $U_{0}\in D\left(\mathcal{A}\right)$, then there exists a unique, strong solution to the system (2.1):

(2) If $U_{0}\in\mathcal{H}$, then system (2.1) admits a unique weak solution $U\in \mathcal{C}^0\left(\mathbb{R}_{+}, \mathcal{H}\right).$

Proof. Initially, we prove that the operator $\mathcal{A}$ exhibits dissipative properties. For any $U \in D(\mathcal{A})$, where $U = \left(\varphi, u, \psi, v, \phi \right)^{T}$, we find

Thus, $\mathcal{A}$ is dissipative.

Next, we prove that the operator $\lambda I - \mathcal{A}$ is surjective for $\lambda > 0$.

Given $F = (f_1, f_2, f_3, f_4, f_5) \in \mathcal{H}$, we verify the existence of $U \in D(\mathcal{A})$ such that

Then, from $(3.3)_{1}$ and $(3.3)_{3}$, we find that

It is straightforward to see that $u\in H^{2}_{L}\left(0, L\right)$. Furthermore, from $(3.3)_{5}$, we can identify $\phi$ as

By inserting $(3.4)_1$ into $(3.3)_2$ and $(3.4)_2$ into $(3.3)_4,$ we obtain

Finding the solution to system (3.6) is the same as identifying $\varphi, \psi\in H^{4}\left(0, L\right)\cap H^{2}_{L}\left(0, L\right)$ such that

for all $w, \chi\in H^{2}_{L}\left(0, L\right)$. By utilizing Eqs (3.5) and (3.7) and then performing integration by parts, we arrive at

where $\overline{\zeta} = \zeta\displaystyle{\int}^{+\infty}_{-\infty}\frac{\mu^{2}\left(\xi\right)}{\xi^{2}+\eta+\lambda}\, d\xi$. By applying (3.4) once more, we infer that

Substituting (3.9) into (3.8), we arrive at

Consequently, the equivalent form of problem (3.10) is given by the following:

where the sesquilinear form $a:\left[H^{2}_{L}\left(0, L\right)\times H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ and the antilinear form $L:\left[H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ are expressed as

and

One can readily verify that the bilinear form $a$ is continuous and coercive, while $L$ is continuous. Therefore, by applying the Lax-Milgram theorem, we deduce that for all $w, \chi \in H^{2}_{L}(0, L)$, the problem (3.11) admits a unique solution $(\varphi, \psi) \in H^{2}_{L}(0, L) \times H^{2}_{L}(0, L)$. Furthermore, by invoking classical elliptic regularity, it follows from (3.10) that $(\varphi, \psi) \in [H^{4}(0, L)]^2$. As a result, the operator $\lambda I - \mathcal{A}$ is surjective for any $\lambda > 0$. As a result, the proof of Theorem 3.1 is completed by applying the Hille-Yosida theorem. □

4.   Strong stability of the system

In this part, we make use of a general criterion of Arendt and Batty (see ^[14] or ^[15]), which states that a $\mathcal{C}_0$-semi-group of contractions $e^{\mathcal{A} t}$ on a Banach space is strongly stable if $\mathcal{A}$ has no purely imaginary eigenvalues and $\sigma\left(\mathcal{A}\right)\cap i\mathbb{R}$ comprises only a countable set of elements. We express our primary result in the form of the following theorem.

Theorem 4.1. ^[16] The $\mathcal{C}_{0}$-semigroup $e^{\mathcal{A} t}$ is strongly stable in $\mathcal{H}$; i.e., for all $U_{0}\in \mathcal{H}$, the solution of (2.1) satisfies

To prove Theorem 4.1, we need the lemmas listed below:

Lemma 4.1. $\mathcal{A}$ does not have eigenvalues on $i\mathbb{R}$.

Proof. We consider two cases: $i\lambda = 0$ and $i\lambda \neq 0$.

Case 1. Solving $\mathcal{A}U = 0$ under the boundary conditions given in (3.1) leads to the conclusion that $U = 0$. Therefore, $i\lambda = 0$ is not an eigenvalue of the operator $\mathcal{A}$.

Case 2. We proceed by contradiction. Assume there exists $\lambda \in \mathbb{R}, \lambda \neq 0$, and $U = \left(\varphi, u, \psi, v, \phi_{1}, \phi_{2}\right) \neq 0$, for which $\mathcal{A}U = i\lambda U$. This gives the following system:

From (3.2), we deduce that $\phi \equiv 0$. Consequently, from $(4.1)_{5}$, we derive $u\left(L\right) = 0.$

Consequently, from $(4.1)_{1}$ and (3.1), we obtain

Inserting $(4.1)_{1}$ into $(4.1)_{3}$ and $(4.1)_{2}$ into $(4.1)_{4}$, we obtain

Then, $\Phi = \varphi+\psi$ and $\Psi = \varphi-\psi$ satisfy

By applying Lemma 5.2 and Corollary 5.1 from ^[7], we obtain the following boundary conditions:

Similarly, the same conditions hold for $\Psi$, specifically:

Now, consider $X = (\Phi, \Phi_x, \Phi_{xx}, \Phi_{xxx}, \Psi, \Psi_x, \Psi_{xx}, \Psi_{xxx})$. We can rewrite (4.2)–(4.4) as the initial value problem:

where

By virtue of the Picard theorem for ordinary differential equations, the system (4.5) has a unique solution $X = 0$. Hence, $\Phi = \Psi \equiv 0$. Therefore, $\varphi = 0$ and $\psi = 0$, i.e., $U = 0$. As a result, $\mathcal{A}$ does not have purely imaginary eigenvalues. □

Lemma 4.2. For $\lambda\neq 0$ or $\lambda = 0$ and $\eta\neq 0$, the operator $i\lambda I-\mathcal{A}$ is surjective.

Proof. We distinguish the following cases:

Case 1: $\lambda \neq 0$.

Our goal is to prove that the operator $i\lambda I - \mathcal{A}$ is surjective for $\lambda \neq 0$. To this end, let $F = (f_{1}, f_{2}, f_{3}, f_{4}, f_{5}) \in \mathcal{H}$. We aim to identify $X = (u, \varphi, v, \psi, \phi) \in D\left(\mathcal{A}\right)$ such that the following equation holds:

Alternatively,

Inserting $(4.6)_{1}$ into $(4.6)_{2}$ and $(4.6)_{3}$ into $(4.6)_{4}$, we obtain

Finding a solution to system (4.7) is the same as finding $\varphi, \psi\in H^{4}\left(0, L\right)\cap H^{2}_{L}\left(0, L\right)$ such that

for all $w, \chi\in H^{2}_{L}\left(0, L\right)$. Using (4.8) and (3.5), followed by integration by parts, yields

where $\overline{\zeta} = \zeta{\int}^{+\infty}_{-\infty}\frac{\mu^{2}\left(\xi\right)}{\xi^{2}+\eta+\lambda}\, d\xi$. Consequently, problem (4.9) can be reduced to the problem

where the sesquilinear form $a:\left[H^{2}_{L}\left(0, L\right)\times H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ and the antilinear form $L:\left[H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ are given by

and

where

It is straightforward to confirm that the bilinear form $a$ is continuous and coercive and that the functional $L$ is continuous. By implementing the Lax-Milgram theorem, we gather that for all $w, \chi \in H^{2}_{L}(0, L)$, the problem (4.10) admits a unique solution $(\varphi, \psi) \in H^{2}_{L}(0, L) \times H^{2}_{L}(0, L)$. Furthermore, using classical elliptic regularity, it results from (4.9) that $(\varphi, \psi) \in [H^{4}(0, L)]^2$. Accordingly, the operator $\lambda I - \mathcal{A}$ is surjective for any $\lambda > 0$. Hence, by invoking the Hille–Yosida theorem, we obtain the desired result.

Case 2: $\lambda = 0$ and $\eta \neq 0$.

The problem in system (4.6) can be expressed as

Then, from (4.9), we obtain

for all $w, \chi \in H^{2}_{L}(0, L)$.

Thus, the system (4.11) can be expressed as the problem:

where the sesquilinear form $a_{\eta}:\left[H^{2}_{L}\left(0, L\right)\times H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ and the antilinear form $L_{\eta}:\left[H^{2}_{L}\left(0, L\right)\right]^{2}\rightarrow \mathbb{C}$ are given by

and

It is straightforward to confirm that $a_{\eta}$ is continuous and coercive, and $L_{\eta}$ is continuous. Then, by Lax-Milgram's theorem, the variational problem (4.12) admits a unique solution $\left(\varphi, \psi\right) \in H^{2}_{L}(0, L) \times H^{2}_{L}(0, L).$ We then deduce from (4.11) that $\left(\varphi, \psi\right) \in \left[H^{4}(0, L)\right]^{2}$. Hence, the operator $\mathcal{A}$ is surjective.

The proof of the Lemma is thus complete. □

Proof of Theorem 4.1. According to Lemma 4.1, the operator $\mathcal{A}$ lacks purely imaginary eigenvalues. Additionally, Lemma 4.2 establishes that $R(i\lambda-\mathcal{A}) = \mathcal{H}$ for all $\lambda \in \mathbb{R}^{*}$ and $R(i\lambda-\mathcal{A}) = \mathcal{H}$ for $\lambda = 0$ when $\eta > 0$. Consequently, by the closed graph theorem of Banach, we conclude that $\sigma(\mathcal{A}) \cap i\mathbb{R} = \emptyset$ for $\eta > 0$, and $\sigma(\mathcal{A}) \cap i\mathbb{R} = \{0\}$ when $\eta = 0$. This completes the proof. □

5.   Discrete energy of the system

We will start making use of finite difference methods (FDM) to derive a discrete representation of the solution of (2.1). More precisely, we use the classical finite difference discretization for the time variable and the implicit compact finite difference method of fourth-order for discretization of the space variable. The comparison between different finite difference methods (implicit, explicit, and semi-implicit) will be considered in a subsequent work.

5.1.   Discrete formulation

Consider the discrete domain $\Omega_h$ of $\Omega = [0, L]$, with a uniform grid given by $x_i = i \Delta x, \quad i = 0, 1, \ldots, M, \, \, (M\geq 5)$ (see Figure 1). The time discretization of the interval $I = [0, T]$ is expressed as $t_j = j \Delta t, \; j = 0, 1, \ldots, N$; $\Delta t = t_{j} - t_{j-1} = \frac{T}{N}$, where $M$ and $N$ are positive integers. Denote by $\varphi(x_i; t_j) = \varphi_i^j$ and $\psi(x_i; t_j) = \psi_i^j$ the value of the functions $\varphi$ and $\psi$ evaluated at the point $x_i$ and the instant $t_j$.

Now, we define the following approximation of the derivatives of $\varphi$ and $\psi$, respectively:

and

Then, using (5.1)–(5.4), we obtain the discrete formulation of the initial condition of the system (1.1) as follows:

for $j = \overline{1, N-1}$, and $i = \overline{2, M-2}$. The discrete formulation of the initial conditions (1.3) reads as follows:

for $i = \overline{0, M}$. Next, we seek the discrete formulation of the boundary conditions.

Lemma 5.1. The discrete formulation of the boundary conditions $(1.2)_3$ reads as follows:

Proof. For $j = \overline{1, N-1}$, we have $\varphi_{xxx}\left(L, t_{j+1}\right) = \gamma \partial^{\alpha, \eta}_{t}\varphi\left(L, t_{j+1}\right).$ By using finite differences, we arrive at

On the other hand, using the change of variable $u = \eta(t_{j+1} - r)$, we obtain

Inserting (5.9) in (5.8), we obtain (5.7). □

Hence, the discrete formulation of boundary condition in the system (1.1) follows as

In the following, we reformulate the system (5.5)-(5.6) in abstract form as follows:

where we denoted by $\mathcal{W} = (\varphi_{0}^{j+1}, \varphi_{1}^{j+1}, \cdots, \varphi_{M}^{j+1}, \psi_{0}^{j+1}, \psi_{1}^{j+1}, \cdots, \psi_{M}^{j+1})^t$,

with

and $\mathcal{K} = \left(\begin{array}{cc} A+ C & B\\ B & A\\ \end{array}\right)\in\mathcal{M}_{2M}(\mathbb{R})$, where $A$ is an $M \times M$ matrix:

with $c_1 = \frac {1}{\Delta x^4}, c_2 = \beta+ \frac {1}{\Delta t^2}+6c_1, c_3 = \frac {1}{\Delta x^2}, c_4 = \frac {1}{\Delta x^3}$ and the matrices $B$ and $C$ are given as follows:

Algorithm 1 provides a summary of all the steps for the calculation of $\varphi^{j}$ and $\psi^{j}$ using matrix decomposition techniques:

Recalling that the energy $E(t)$ is expressed as

Now, we define the following approximations of the derivatives of $\varphi$ and $\psi$, respectively:

The $L^2$ norm of a discretized function is approximated by

Thus, the discrete energy of the system (1.1)-(1.2) at time $t_j$ is as follows:

Algorithm 2 summarizes all the steps for calculating of the discrete Energy $E(t)$ using the $L^2$ norm define by (5.10).

5.2.   Numerical test

To evaluate the asymptotic properties of the solutions to the system (1.1), we consider the following data: $\alpha = 0.5,$ $\beta = \eta = 1,$ $\Delta t = 10^{-2}, \Delta x = 10^{-1}$, $L = 1$, and the initial conditions:

Figures 2–4 illustrate the decay of the discrete energy for the chosen initial data. Specifically, Figure 3 demonstrates that when $T$ is sufficiently large, the decay becomes more rapid. Figure 5 demonstrates the exponential decay of the discrete energy for the chosen initial data, as well as for other initial conditions. The conjecture of exponential decay arises from the observation that the plot of the logarithm of the discrete energy forms a straight line with a negative slope, indicating an exponential decay pattern. This behavior has been consistently observed for various initial data, further reinforcing the conjecture. This observation also poses an open question regarding the exponential decay of the energy in the system (1.1).

Figures 6 and 7 illustrate the exponential decay of the discrete energy for the chosen initial data and for other initial conditions as well, for different values of $\alpha$: $\alpha = 0.75$ and $\alpha = 0.95$. We observe a slight increase in the decay rate as $\alpha$ approaches 1.

Now, we provide an example to test the correctness and effectiveness of the numerical scheme. This example includes a comparison with the exact solution to demonstrate the scheme's convergence rate and validate its accuracy. These results ensure the reliability of the analysis and conclusions presented in Section 5.2.

We consider the exact solutions to the problem as follows:

The problem is governed by the equations:

with the following boundary conditions:

and the initial conditions:

The source terms are given by

Figures 8 and 9 show the comparison between the exact solutions and the numerical solutions for $\varphi$ and $\psi$.

Now, we define the root mean square error (RMSE) as

where $u$ and $\hat{u}$ denote the exact and the numerical solutions, respectively.

The data presented in Table 1 below were obtained by comparing the exact solutions $\varphi$ and $\psi$ with their respective numerical approximations $\hat{\varphi}$ and $\hat{\psi}$ using different discretization values of $\Delta x$ and $\Delta t$.

Finally, to present the rate of convergence, we begin by defining the pointwise rate of convergence, which corresponds to the opposite of the slope of the logarithmic energy curve. This calculation is performed for the example provided at the beginning of Section 5.2.

Define the pointwise rate of convergence as

and the rate of convergence as

In Figure 10, we observe the behavior of the decay rate over time. We note that the curve becomes horizontal, indicating that the rate is constant, which supports the notion that the energy decay is exponential. In Table 2, we provide the values of the energy decay rate for different values of $\Delta x$ and $\Delta t$.

6.   Conclusions

In the present work, we studied the strong stabilization of a coupled Euler-Bernoulli beam system with one boundary dissipation of fractional derivative type. First, we analyzed the strong stability of the system (1.1)-(1.2). Next, we employed a finite difference scheme to compute the numerical solutions and demonstrated the stability of the discrete energy. Although we achieved exponential energy decay in one example, the question of exponential energy decay remains open.

In advancing the finite difference methods (FDM) for solving complex fractional and coupled PDEs, this manuscript makes a considerable impact on numerical analysis and applied mathematics, offering a significant resource for researchers and practitioners working on advanced modeling challenges.

Author contributions

Boumediene Boukhari, Ahmed Bchatnia and Abderrahmane Beniani: Visualization; Foued Mtiri: Investigation; Ahmed Bchatnia: Supervision, Validation; Ahmed Bchatnia and Abderrahmane Beniani: Conceptualization, Methodology; Boumediene Boukhari, Foued Mtiri and Abderrahmane Beniani: Resources; Boumediene Boukhari, Foued Mtiri, Ahmed Bchatnia and Abderrahmane Beniani: Writing–original draft. All authors have read and approved the final version of the manuscript for publication.

Acknowledgements

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khaled University for funding this work through Large Research Project under grant number RGP2/368/45.

Conflict of interest

The authors declare that there are no conflicts of interest in this paper.