Exact Solutions of Fractional Order Oscillation Equation with Two Fractional Derivative Terms

The fractional oscillation equation with two fractional derivative terms in the sense of Caputo, where the orders α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document} satisfy 1<α≤2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1<\alpha \le 2$$\end{document} and 0<β≤1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0<\beta \le 1$$\end{document}, is investigated, and the unit step response and two initial value responses are obtained in two forms by using different methods of inverse Laplace transform. The first method yields series solutions with the nonnegative powers of t, which converge fast for small t. The second method is our emphasis, where the complex path integral formula of the inverse Laplace transform is used. In order to determine singularities of integrand we first seek for the roots of the characteristic equation, which is a transcendental equation with four parameters, two coefficients and two noninteger power exponents. The existence conditions and properties of the roots on the principal Riemann surface are given. Based on the results on the characteristic equation, we derive these responses as a sum of a classical exponentially damped oscillation, which vanishes in an indicated case, and an infinite integral of the Laplace type, which converges fast for large t, with a steady component in the unit step response. Asymptotic behaviors of solutions for large t are derived as algebraic decays in negative power laws characterized by the orders α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha$$\end{document} and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta$$\end{document}. The fractional system exhibits a magical transition from oscillation to monotonic decay in negative power law.

A fractional oscillator is a mass adjacent to a set of parallel spring and viscoelastic materials obeying the fractional stress-strain relationship. Then the fractional oscillator equation is a fractional ordinary differential equation. Equations of this type were investigated in [3,5,9,[21][22][23][24][25].
We use the Caputo fractional derivative to set up the model in this paper since classical initial value condition can be equipped with.
The Laplace transform is frequently-used in analytical research for fractional calculus [5,26,27]. Suppose f(t) is locally integrable on (0, +∞) and exponentially bounded for large enough t, the Laplace transform of the fractional integral of Re(s) > 0 , where 0 is the threshold value for exponentially boundedness. Hereafter, by putting a tilde we denote the corresponding Laplace transform. If f(t) has the continuous, locally integrable, exponentially bounded for large enough t, and Laplace transformable derivative of order m on (0, +∞) , then the Laplace transform of the Caputo fractional derivative is [26] Note that if = m , then Eq. (6) degenerates to the formula of the Laplace transform of integer-order derivatives. We will use the formula of the Laplace transform of power functions In [22], the solution of the fractional oscillation equation x �� (t) + b D t x(t) + cx(t) = 0 , where b > 0 , c > 0 and 0 < < 1 , was considered by using the complex path integral formula of the inverse Laplace transforms. More information on the solution can be derived by such method than the solution in the form of fractional power series in t or in terms of a series involving the higher order derivatives of the Mittag-Leffler functions [5,25]. But the complex path integral formula involves the computation of residues of integrand, where the characteristic equation is a transcendental equation with a power of a noninteger degree without an explicit formula for the roots. So additional work and technique are necessary. In [28], the case of 1 < < 2 was considered. In [23,25], the range of was extended to 0 < < 2 and an external excitation was taken into account.
We note that other forms and methods for the fractional oscillator equation including other fractional definitions are existent. In [29], more general fractional integral operator and its boundedness were described. In [30], the considered oscillator equation is a two-term-equation free from the second order derivative term. In [31], different classes of fractional oscillators in the sense of the Weyl fractional calculus were presented and analyzed. In [32], fractional oscillator process was obtained as solution of a stochastic differential equation, where the Weyl

3
fractional calculus and Fourier transform were used. In [33], the oscillation equation was solved by the averaging method to obtain the approximately analytical solution. In [34,35], distributed-order cases were involved. We also note that the method of fractional power series is practical for many fractional ordinary differential equations [36][37][38] and particularly operational for nonlinear problems, such as the fractional logistic equation in [36], where a convergence analysis was also carried out.
In this paper we investigate the oscillation equation with two fractional derivative terms in the sense of Caputo by using two methods of the inverse Laplace transform. The first method yields series solutions with the nonnegative powers of t, which converge fast for small t. The second method is our emphasis, where the complex path integral formula of the inverse Laplace transform is used. The derived solutions involve infinite integrals of the Laplace type, which converge fast for large t, and can provide more information for the solution structures.
The text is organized as follows. In the next Sect. 2, we present the considered model, conduct a preliminary analysis by the Laplace transform and introduce three basic and significant solutions, i.e., the initial unit rate response, the unit step response and the initial unit displacement response. In Sect. 3, we give these solutions in series forms with nonnegative powers of t. In Sect. 4, we consider the roots of the characteristic equation, a transcendental equation with two powers of noninteger degree, and give the existence conditions and properties of the roots on the principal Riemann surface. In Sect. 5, the complex path integral formula of the inverse Laplace transform is used to derive these solutions, which are found to be a sum of a classical exponentially damped oscillation, which vanishes in an indicated case, and an infinite integral of the Laplace type. More details are yielded from the results. Section 6 summarizes our conclusions.

The Fractional Order Oscillation Equation with Two Fractional Derivative Terms
We consider the initial value problem of fractional order oscillation equation with two fractional derivative terms, where and are limited to the intervals 1 < ≤ 2 and 0 < ≤ 1 , f(t) is the prescribed excitation or forcing term, and x(t) is the system response. In order to make the excitation cover the usual cases with discontinuities such as the rectangular pulse or boxcar functions, we require f(t) to be piecewise absolutely continuous on [0, +∞) with at most finitely many jump discontinuities, and have exponentially bounded first order derivative for large enough t.
Operating the Laplace transform to Eq. (8) with the initial condition (9) gives Solving for the Laplace transform of the response we have the expression By introducing the following two functions by taking the inverse Laplace transforms, which satisfy the relationship and denoting so the solution of the initial value problem, (8) and (9), can be expressed as We note that G(t) is the unit step response of the system (8), G 1 (t) is the initial unit rate response, and G 0 (t) is the initial unit displacement response. We checked that although the inverse Laplace transform of 1 s +bs +c exists, the fractional derivative of order of the inverse Laplace transform could encounter a divergent improper integral. Thus, instead of "impulse response", the unit step response is considered in this work.
For the integer-order case, = 2 and = 1 , the inverse Laplace transform in Eq. (10) gives the classical results for the initial unit rate response, which may be obtained by calculating residues of the integrand at its poles, From Eq. (12), the unit step response in the integer-order case is Here superscript asterisks are labelled to denote the integer-order case. We note that the cases b 2 − 4c > 0 , b 2 − 4c = 0 and b 2 − 4c < 0 correspond to overdamping, critical damping and underdamping oscillation systems.
In the sequel, we devote to derivation of the initial unit rate response G 1 (t) and the unit step response G(t) analytically by using two different methods of the inverse Laplace transform, while the initial unit displacement response G 0 (t) is direct from Eq. (13).

Series Solutions with the Nonnegative Powers of t
First, we derive the series expression of the initial unit rate response G 1 (t) comprised of the nonnegative powers of t. For this purpose, the Laplace transform G 1 (s) needs to be decomposed into a series constituted with the negative powers of s. Deforming and expanding the Laplace transform G 1 (s) in Eq. (10) firstly yield Next, using the formula Note that there is a constant such that as long as Re(s) is greater than the constant the above expansions hold. Now taking the inverse Laplace transform term by term using the formula (7), we obtain the series expression of the initial unit rate response with the positive powers of t, .
Note that the method of series expansion of the image function and then termwise inverse transformation is frequently used in fractional differential equations to obtain a series solution convergent for all t ≥ 0 . The effectiveness of the method was investigated and guaranteed [39]. By introducing the index n = j + i instead of i, we have another expression The unit step response can be obtained readily from Eq. (12) by operating integrations for Eq. (18) term by term, These series solutions converge fast for small t, and have the following asymptotic behaviors for small t, We take b = c = 1 and = 0.5 fixed, simulate the two responses G 1 (t) and G(t) in Figs. 1 and 2, respectively, by the sums of n = 0 through n = 75 in Eqs. (18) and (19) for = 2, 1.75, 1.5, 1.25 . From the two figures, the initial behaviors conform with the asymptotic expressions in (20) and (21), and with decreasing of , the oscillation phenomena are gone. More information will be given in Sect. 5. . .

Roots of the Characteristic Equation
This section is a preparation for the next section. We consider the roots of the characteristic equation, a transcendental equation with noninteger powers, where and are not both integers. This is necessary to be considered for calculating residues by using the complex path integral formula of the inverse Laplace transforms in Eqs. (10) and (11). The power functions with noninteger degrees are multi-valued with the branch point s = 0 . In the next section, the negative real axis will be served as a branch cut. So we look for roots of Eq. (23) on the principal Riemann surface − < arg s < . Owing to the fact that s satisfies Eq. (23) if and only if its conjugate s satisfies Eq. (23). Thus we only need to confine to the half complex plane 0 ≤ arg s < to analyze the problem, while will state results on the principal Riemann surface − < arg s < .
Let s = re iθ where r > 0 , and then substituting it into Eq. (23) leads to r e iαθ + br e iβθ + c = 0. Separating the real part and the imaginary part we have the system of equations For a root s = re iθ of Eq. (23), it follows that ≠ 0 from Eq. (24). In addition, Eq. (25) requires sin( ) < 0 . Thus we can contract the range of the arguments of roots on the upper half complex plane to < < . This means that the roots of Eq. (23), if any, have negative real parts.
From Eq. (25), we solve the modulus r as    ◻ In summary, in the case (iiib) in Proposition 2, i.e., that the orders and coefficients in Eq. (23) satisfy the inequalities < 1 , − < 1 and (32), then Eq. (23) has no roots on the principal Riemann surface − < arg s < . In other cases in

Solutions Fast Convergent for Large t
We deduce the solutions in Eqs. (10) and (11) by using the complex path integral formula of the inverse Laplace transform. The two orders satisfy 1 < ≤ 2 and 0 < ≤ 1 , where two equal signs cannot hold simultaneously. First, we consider the initial unit rate response, where Br denotes the Bromwich path, i.e. a straight line from s = − i∞ to s = + i∞ , where is chosen such that all the singularities of G 1 (s) lie to the left of the line. A schematic diagram of the Bromwich path is shown in Fig. 3.
Since s = 0 is a branch point of the integrand, we consider the problem in the principal Riemann surface − < arg s < and make a branch cut along the nonpositive real axis. By the residue theorem on contour integration and Jordan's (33) s 1,2 = e ±i = − ± i ,  Fig. 3. We note that the residue contribution G 1R (t) vanishes if the case (iiib) in Proposition 2 is encountered. After some algebraic operations, the residue contribution is reduced to where are constants unrelated with t. The Hankel path is decomposed to three subpaths, where the arrow → denotes the integral from −∞ to 0, while the arrow ← denotes the reverse direction, and C is the counterclockwise small circle: s = e iθ , > 0, − < < . It is verified that the limit of the integration on the small circle in Eq. (38) vanishes. So the Hankel path integral is simplified as

By some algebraic operations the Hankel path integral is reduced to a real infinite integral,
where Thus the initial unit rate response G 1 (t) is expressed as a sum of a classical damped oscillation in (36) and a real infinite integral in (39). The unit step response G(t) should be obtained from the initial unit rate response G 1 (t) by using Eq. (12). But we could encounter the fractional integrals of e − t cos( t) and e − t sin( t) , which are difficult for analytical computation. Here we continue to use the complex path integral formula of the inverse Laplace transform in Eq. (11) similar to the preceding derivation for G 1 (t) to obtain the unit step response G(t) in a similar format. The unit step response G(t) is decomposed into where the residue contribution is calculated to be and the Hankel path integral is reduced to the infinite integral, where We note that for this solution, the limit of the integration on the small circle C is 1/c in Eq. (44). In addition, P in Eq. (43) and Q(r) in Eq. (45) are the same as in Eqs. (37) and (41).
Q(r) = r 2 + 2br + cos( ( − )) + b 2 r 2 + 2cr cos( ) + 2bcr cos( ) + c 2 . (42) By using Eq. (13), the initial unit displacement response G 0 (t) has the form The obtained three solutions, G 1 (t) , G(t) and G 0 (t) , have similar formats, and they are in the form of a sum of an exponentially damped oscillation, which vanishes in the case (iiib) in Proposition 2, and an infinite integral of the Laplace type, which converges fast for large t, with a steady component 1/c for G H (t) . The infinite integrals of the Laplace type can be considered as a continuous superposition of exponential decays at different rates, and K 1 (r) and K(r) may be regarded as the decay rate spectrums for the initial unit rate response and the unit step response. Thus, we call G 1R (t) , G R (t) and −cG R (t) the damped oscillation components of the three responses G 1 (t) , G(t) and G 0 (t) , while G 1H (t) , G H (t) and 1 − cG H (t) the monotonic components. Just the monotonic components make the fractional case different from the integer-order case. We note that the monotonic components will become monotonic for sufficiently large t by the asymptotic behaviors considered below.
(46) The infinite integrals of the Laplace type can be computed using the command "NIntegrate" in Mathematica. Taking b = c = 1 , = 1.5 and = 0.5 , in Fig. 4, the spectrums K 1 (r) and K(r) are depicted on the interval 0 < r < 8 , and in Fig. 5, the monotonic components G 1H (t) and G H (t) are shown. The two spectrums remain positive and negative, respectively, in a right open neighborhood of r = 0 , and approach to zero as r → +∞ at the rate O r − + −2 and the rate O r − −1 , respectively. The initial values of the monotonic components are determined through the spectrums as G 1H (0) = 1 ∫ ∞ 0 K 1 (r)dr and G H (0) = 1 c + 1 ∫ ∞ 0 K(r)dr . More information on the monotonic components will be discussed in this section.
Next, we check the solution curves by the results in this section with that in Sect. 2. In Figs. 6 and 7 , curves of G 1 (t) in Eq. (35) and G(t) in Eq. (42) are shown, respectively, for the same parameter values as in Figs. 1 and 2 . Here the effective regions of curves can be greatly extended than that in Sect. 2. We note that the case of b = c = 1 , = 0.5 and = 1.25 belongs to the case (iiib) in Proposition 2.
Then, we examine several cases where two fractional orders are close to integers. Take b = 2 and c = 0.25 . When = 2 and = 1 , the initial unit rate response and the unit step response are from Eqs. (15) and (16), which belong to overdamping case. In Fig. 8, curves of G * 1 (t) together with G 1 (t) for = 1 and = 1.9, 1.8 and 1.7 are shown. In Fig. 9, curves of G * (t) and G(t) are displayed for the same parameter values.
Take b = 2 and c = 5 . When = 2 and = 1 , the initial unit rate response and the unit step response are from Eqs. (15) and (16), which belong to underdamping case. In Fig. 10, curves of G * 1 (t) together with G 1 (t) for = 2 and = 0.9, 0.8 and 0.7 are shown. In Fig. 11, curves of G * (t) and G(t) are displayed for the same parameter values. The four Figs. 8, 9, 10 and 11 display expected variation tendencies that the solutions G 1 (t) and G(t) approach to that of the integer-order cases as and are close to 2 and 1, respectively.

Asymptotic Behaviors for Large t
The infinite integrals of the Laplace type facilitate to determine the asymptotic behavior for large t from the asymptotic expressions of the spectrums for small r obtained via Eqs. (40) and (45), First, G 1H (t) in Eq. (39) may be written by the replacement of the integration variable z = rt instead of r as Inserting the asymptotic expression (47) and calculating integrals, we derive the asymptotic behavior of G 1H (t) for large t, In a similar manner, we obtain the asymptotic behavior for G H (t), The monotonic components obey algebraic decays in the form of negative power laws, apart from the constant 1/c for G H (t) , and are indeed monotonic for large enough t. The first terms on the right hand side of (50) is dominant and positive as long as 1 < < 2.
Inasmuch as the residue contributions in Eqs. (36) and (43) weaken in oscillating at negative exponential rates, the initial unit rate response G 1 (t) and the unit step response G(t) have the consistent asymptotic behaviors for large t with their respective monotonic components, i.e., the Hankel integral contributions. The asymptotic behavior of the initial unit displacement response G 0 (t) is derived from Eq. (13). From the above explanations and Eqs. (50) and (51), we summarize the results as follows. The response curve x(t) is shown in Fig. 12, where = 4 . For the sinusoidal excitation f (t) = sin(t) , the response is which is plotted in Fig. 13. In the two figures, the responses compared with the excitations, which are also shown together with their outputs, display delayed reactions. These imply energy dissipation in the oscillation system.

Conclusions
Solutions of the fractional oscillation equation D t x(t) + b D t x(t) + cx(t) = f (t) , b > 0 , c > 0 , 1 < ≤ 2 and 0 < ≤ 1 , are considered by using two different methods of inverse Laplace transform. Particular concerns are the three basic and significant solutions, i.e., the unit step response, G(t), and the two initial value responses, G 1 (t) and G 0 (t) . The first method yields series solutions with the nonnegative powers x(t) = G(t) * cos(t), of t, which converge fast for small t, via the series expansion to negative powers of s and the inverse Laplace transform term by term. The second method, where the complex path integral formula of the inverse Laplace transform is used, is our emphasis since more techniques are required and more information for the solution structures is educed. In order to determine singularities of integrand we first seek for the roots of the characteristic equation, which is a transcendental equation with four parameters, two coefficients and two noninteger power exponents. In Sect. 4, we give the existence conditions of the roots on the principal Riemann surface and prove that when the roots exist, they are a pair of conjugated simple complex roots with a negative real part. In Sect. 5, we derive these responses as a sum of a classical exponentially damped oscillation, which vanishes in an indicated case, and an infinite integral of the Laplace type, which converges fast for large t, with a steady component for the unit step response. Asymptotic behaviors of solutions for large t are derived as algebraic decays in negative power laws characterized by the orders and as long as and are not both integers. We also indicate that the fractional system exhibits a magical transition from oscillation to monotonic decay in negative power law, which is verified and illustrated in this work.
Journal of Nonlinear Mathematical Physics (2023) 30:531-552 not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.