Frequency/Laplace domain methods for computing transient responses of fractional oscillators

Although computing the transient response of fractional oscillators, characterized by second-order differential equations with fractional derivatives for the damping term, to external loadings has been studied, most existing methodologies have dealt with cases with either restricted fractional orders or simple external loadings. In this paper, considering complicated irregular, but deterministic, loadings acting on oscillators with any fractional order between 0 and 1, efficient frequency/Laplace domain methods for getting transient responses are developed. The proposed methods are based on pole-residue operations. In the frequency domain approach, “artificial” poles located along the imaginary axis of complex plane are designated. In the Laplace domain approach, the “true” poles are extracted through two phases: (1) a discrete impulse response function (IRF) is produced by taking the inverse Fourier transform of the corresponding frequency response function (FRF) that is readily obtained from the exact transfer function (TF), and (2) a complex exponential signal decomposition method, i.e., the Prony-SS method, is invoked to extract the poles and residues. Once the poles/residues of the system are known, those of the response can be determined by simple pole-residue operations. Sequentially, the response time history is readily obtained. Two single-degree-of-freedom (SDOF) fractional oscillators with rational and irrational derivatives, respectively, subjected to sinusoidal and complicated earthquake loading are presented to illustrate the procedure and verify the correctness of the proposed method. The verification is conducted by comparing the results from both the Laplace and the frequency domain approaches with those from the numerical Duhamel integral method. Furthermore, the proposed method applied to a 3-DOF system with fractional damping is also demonstrated.


Introduction
A second-order linear system with a fractional derivative damping term, which can capture the frequencydependence of damping properties in viscoelastic materials, is called a fractional oscillator [11]. On computing the responses of fractional oscillators, existing analytical solutions are limited to cases with specific frac-tional order q, such as q = 1/2 and q = 1/3, subjected to simple loadings, including impulse and step loading [6,22,26]. In contrast, this paper develops efficient frequency/Laplace domain methods to compute the response of fractional oscillators with arbitrary fractional order to complicated irregular, but deterministic, loadings.
Roughly, four distinct strategies in the literature have been pursued to compute the response of a fractional oscillator. As a common way to obtain the transient response of a linear system to any loading is through the Duhamel integral, which employs the impulse response function (IRF) of the system [25], a popular strategy has been attempting to obtain the fractional oscillator's IRF first [3,4,16,21,27]. While considering an oscillator with fractional order q = 1/2, Gaul et al. [12] computed the system's IRF by taking the inverse Fourier transform analytically from the system's frequency response function (FRF). Proceeding in the Laplace domain, Suarez and Shokooh [25] obtained the IRF of a fractional oscillator with q = 1/2 or q = 1/3, which involved solving a 4th-order polynomial for finding the poles of an oscillator with q = 1/2, and a 6th-order with q = 1/3. Using the generalized Mittag-Leffler functions, Achar et al. [1,2] expressed the IRF as a summation of infinity terms. Another common strategy for computing the fractional oscillator's response has been transforming the equation of motion into an alternative form so that it allowed to be solved by ordinary calculus [5,8,9,20,29]. When q is limited to a simple rational number, several papers [6,22,26] transformed the fractional differential equation into a set of equations in the so-called state variable domain where a set of decoupled fractional differential equations with the same fractional order was derived. The solution procedure has been operated in the Laplace domain, and each modal response must be obtained through an inverse Laplace transform. As this approach was often confronted with the difficulty of the inverse Laplace transform, the explicit response solutions of fractional oscillators to complicated loadings were difficult to obtain. Because of the difficulties of calculating inversion of Laplace transforms, Sheng et al. [23] applied the numerical inverse Laplace transform algorithm in fractional calculus. This state variable analysis method not only requires q be a rational number but also rapidly increases the number of state variables when the denominator of q becomes large. The third strategy has been entirely numerical, which approxi-mates the fractional derivative term at discrete times using central difference method [17,24] or Grünwald-Letnikov definition [9,20,29]. The fourth strategy proposed in a few articles [11,28] involved the change of variable and integral approximation schemes to convert the fractional oscillator's equation of motion into a set of coupled linear equations with ordinary derivatives so that they could be computed by any time-domain numerical integration scheme.
This paper proposes novel frequency/Laplace domain methods based on pole-residue operations for computing the transient responses of fractional oscillators with arbitrary q between 0 and 1 to complicated loadings. The key of the proposed method is to obtain the poles and residues of the system transfer function (TF) for fractional oscillators, then those of the response. Knowing the poles and residues of the response allows it to be expressed immediately in the time domain [13,14]. The proposed method includes the following steps: (1) obtaining the poles and residues of the TF for the fractional oscillator from the analytical FRF; (2) computing the poles and residues of the response from those of the TF and the excitation by simple pole-residue operations; and (3) expressing the response time history from the poles and residues of the response.
Two single-degree-of-freedom (SDOF) fractional oscillators with rational and irrational derivatives, respectively, subjected to sinusoidal and complicated earthquake loading are provided to illustrate the proposed method. In both examples, the response calculation will be carried out in both the frequency domain and the Laplace domain. The accuracy of the proposed method will be verified through a time domain approach that employs the numerical Duhamel integral. In addition, the proposed method operated in the Laplace domain for a 3-DOF system with fractional damping will be demonstrated in another numerical example.

System functions
The most essential background material to this study is the system functions, which are employed to characterize the relationship between the response (output) and the excitation (input) of a linear time-invariant system, including the IRF in the time domain, the FRF in the frequency domain, and the TF in the Laplace domain. Denoting the input and output signals associated with a linear system as u(t) and y(t), respectively, they are related in the time domain as where h(t) is the system's IRF. In the frequency domain, one computes the steady-state response where H(ω) is system's FRF, which is the Fourier transform of h(t) where F(·) denotes the Fourier transform operator and When the system is causal, i.e., h(t) = 0 for t ≤ 0, one writes the input-output relationship in the Laplace domain as whereh(s) is system's TF, which is the Laplace trans- where L(·) denotes the Laplace transform operator. Furthermore, H(ω) can be also obtained directly from h(s) by substituting s with iω, that is

Transfer functions of fractional oscillators
For a fractional oscillator, the equation of motion is written as: where m, c and k are the mass, damping and stiffness of the oscillator; x(t),ẋ(t) andẍ(t) are the displacement, velocity and acceleration of the structure response, respectively; f (t) is the external loading; and D q 0 + x(t) for 0 < q < 1 follows the definition of Caputo's fractional derivative [6,7,19] (9) in which (·) is the gamma function. Since the Laplace transform of Eq. 9 is [6] taking the Laplace transform of Eq. 8 yieldsx(s) = h(s)f (s) with the corresponding TF Substituting s by iω into Eq. 11 yields

Response calculation in the frequency domain
The concept of the pole-residue method developed in Ref. [14] will be employed to compute x(t). When the response computation is carried out in the frequency domain, many imaginary poles iω n are assigned, wherê ω n = n ω and n =integer. Referring to Eq. 12, and denoting H n = H(ω n ), one writes From Eq. 4, the system's IRF can be approximated as H n e iω n t ω (14) where N n must be sufficiently large so that H N n /2 → 0. Denote thus A n e iω n t (16) Noting that A n e iω n t and A n s−iω n form a Laplace transform pair [18] and performing the Laplace transform of Eq. 16 yields Now,h(s) is in a pole-residue form with poles iω n and corresponding residues A n . Obviously, the inverse Laplace transform of Eq. 17 will go back to Eq. 16. An important fact drawn from Eqs. 16 and 17 is that once iω n and A n are given, one can readily obtain h(t) and h(s).
For a periodic excitation f (t) with period T , it can always be expressed as a complex Fourier series where ω m = m ω and ω = 2π/T , and the corresponding complex Fourier coefficient Because f (t) is a real-valued signal, C −m , the coefficient at −ω m , is equal to the complex conjugate of C m . Theoretically, a Fourier series representation of f (t) contains an infinite number of terms; however, in practice f (t) can generally be approximated with sufficient accuracy by finite terms. For a digital signal sampled with a constant t, the fast Fourier transform (FFT) algorithm has been often invoked to compute C m . From Eq. 18, the correspondingf (s) is This is a pole-residue form for a periodic function with poles iω m and residues C m . According to the operation in the Laplace domaiñ x(s) =h(s)f (s), using Eqs. 17 and 20, one writes First, assuming the chosenω n and ω m do not overlap, as the common denominator in Eq. 21 is (s − iω n )(s − iω m ), the pole-residue form (partial-fraction form) of x(s) is [10] x(s) = where and Then, carrying out the inverse Laplace transform yields In Eq. 25, while the first summation term is the transient response related to the system poles iω n , the second summation term is the steady-state response related to the excitation frequencies ω m , Second, if same poles are shared byh(s) andf (s), then second-order poles exist in Eq. 21. When ω = ω, and N m = N n , Eq. 21 becomes Thus, the pole-residue form of Eq. 26 is Taking the inverse Laplace transform of Eq. 27 gives Note that since the system and excitation share the same poles, one cannot separate the total response into the natural response and the forced response.

Response calculation in the Laplace domain
As the pole-residue form ofh(s) in Eq. 17 has been based on designated pure imaginary poles, the system becomes a "periodic" oscillator with its IRF as a periodic function. This is certainly against the true physics, but the response calculation might be still valid for a period. To get a physically meaningfulh(s), one can search for the "true" poles μ n of the system and the corresponding residues β n . Consequently, the pole-residue form ofh(s) is expressed as For extracting the true poles and residues of the system, a procedure shown schematically in Fig. 1 is proposed. First, taking the inverse Fourier transform of H(ω) in Eq. 12 yields a real valued IRF, h(t), for t ≥ 0. While deriving h(t) analytically is challenging, if not impossible, a straightforward way to compute a discrete IRF, h(t k ), can be achieved by taking the inverse discrete Fourier transform (IDFT)-which can be efficiently carried out by using inverse fast Fourier trans- has been computed, one applies the Prony-SS method [15] to approximate it as where N n , the number of component, is usually a small integer, and β n and μ n are constant coefficients. A summary of the Prony-SS method is given at Appendix I; its first step is to determine N n . This is fulfilled by imposing a truncated singular value decomposition of a Hankel matrix built from the dealing signal h(t k ), which has the similar concept to that of a principal component method. For a smooth signal, just a small number of N n is needed to achieve a good approximation. The method's second step is to compute μ n obtained by the eigenvalue analysis of a realized matrix, and its last step is to compute β n based on solving a set of linear simultaneous equations. Because h(t) is a real-valued function, μ n must be either real numbers or in complex conjugate pairs. In turn, the coefficients β n must be also real or appear in complex conjugate pairs. When x(t) is computed based on Eq. 30, it is referred to as the Laplace domain approach in this paper, no matter how the input f (t) has been decomposed. For the decomposition of f (t), either a Prony series or a Fourier series can be used.

Input decomposition: Prony series
For an irregular excitation signal f (t) of a finite duration of T , it can always be approximated by using the Prony-SS method [15]: where N is the number of components; α and λ are constant coefficients. Taking the Laplace transform of Eq. 32 yields Fromx(s) =h(s)f (s), using Eqs. 30 and 33, one writes Expressing Eq. 34 in a pole-residue form yields and Performing the inverse Laplace transform to Eq. 35 then yields In Eq. 38, the first sum term is the natural response governed by the system poles, μ n ; and the second term is the forced response related to the excitation poles, λ .

Input decomposition: Fourier series
Similar to the formulation of Eq. 34, using Fourier series forf (s), one writes In its pole-residue form, one has: and So, taking the inverse Laplace transform of Eq. 40 gives Note that although the natural responses of Eq. 38 and Eq. 43 share the same poles, they have different residues; consequently, Eqs. 38 and 43 have distinct natural responses. Clearly, the forced response in Eq. 43 is the steady-state response to the periodic loading. Theoretically, the total responses within T for Eqs. 38 and 43 should be identical.

Numerical studies
A simple fractional oscillator with q = 1/2 investigated in Ref. [17] is chosen in the numerical studies so that independent numerical verification for the computed response is possible; both simple sinusoidal and complicated earthquake loadings are considered. In addition, a fractional oscillator containing two irrational derivative terms excited by realistic earthquake motion will also be investigated to demonstrate a broader application of the proposed method.
5.1 Example 1: Fractional oscillator with q = 1/2 The fractional oscillator considered first has q = 1/2, with other values being mass m = 1, stiffness k = 1, and damping coefficient c = 0.1, noting that throughout the numerical studies the unit system is the MKS (meter-kilogram-second) system. According to Eq. 12 and the chosen numerical quantities, the oscillator's FRF is:

Sinusoidal loading
Let the simple sinusoidal loading be f (t) = sin t, which can be expressed as f (t) = αe λt +ᾱeλ t where α = −0.5i and λ = i. In other words, the chosen excitation has a pair of complex conjugate poles λ,λ = ±i and the corresponding residues α,ᾱ = ∓0.5i.
Response calculation in the frequency domain While applying the proposed method in the frequency domain, the formulation ofh(s) is given in Eq. 17, whose poles are assigned at iω n and residues A n = H n ω/(2π) are computed directly from a discretized FRF (see Eq. 13).
In the numerical implementation, Eq. 17 has been utilized with ω = π/256 and N n = 2048. For the simple excitation f (t) = sin t, while using Eq. 21, it has N m = 2, ω 1 = 1, and C 1 = −0.5i. Plotted in Fig. 2 is the x(t) computed based on Eq. 25; it agrees well with that obtained from the time domain method. As indicated in Eq. 25, the response x(t) is composed of two terms, the natural response and the forced response. Figure 3 plots these two terms separately. Note that the frequency domain method's natural response will also become periodic after 2π/ ω = 512 seconds. Response calculation in the Laplace domain To apply the proposed method in the Laplace domain, one needs to follow the steps shown in Fig. 1 to obtain the poleresidue form ofh(s). To begin, a discretized FRF is generated from Eq. 44 with the frequency interval ω = π/64, number of frequency components N n = 4096, and cut-off frequency ω N = 64π . Then, performing the IFFT of this FRF to obtain a discrete IRF, h(t k ), plotted in Fig. 4, which is in a good agreement with that reported in Ref. [17]. Processing h(t k ) by the Prony-SS method results in an approximated IRF, denoted h a (t), to be where μ 1 = −0.0353 − 1.0353i, β 1 = −0.0055 + 0.4957i, μ 2 = −0.2992 and β 2 = 0.0119. Eq. 45 is plotted in Fig. 5, which indicates that the IRF is dominated by an oscillatory term (complex conjugate poles), together with a negligibly small non-oscillatory term (real pole). This feature concurs with the derivation in Refs. [25,26]. In addition, the obtained μ 1 value is nearly identical to that reported in Ref. [25]. To evaluate the accuracy of Eq. 45, the relative error h (t k ) between h(t k ) and h a (t k ), defined as is computed and plotted in Fig. 6, which indicates the maximum of h occurring at t = 0 is only about 0.8%. Furthermore, from Eq. 45, the corresponding FRF is Shown in Fig. 7 is the comparison between Eqs. 44 and 47, including the magnitude in Fig. 7a and phase angle in Fig. 7b. Indeed, they are in excellent agreement. In this example, when the response functionsx(s) and x(t) are expressed as Eqs. 35 and 38, respectively, there are 5 response poles, including two from the excitation poles and three system poles; the corresponding residues to the excitation and system poles are computed from Eqs. 36 and 37, respectively. The numerical values for all poles and residues are listed in Table 1. From them, the computed response x(t) is shown in Fig. 8, which agrees well with the x(t) that has been obtained by performing the numerical Duhamel integral based on the original h(t). The natural response and forced response obtained from Eq. 38 are plotted in Fig. 9. As expected theoretically, they are completely different from those in Fig. 3. While the natural response diminishes with time due to the presence of the damping, the forced response remains sinusoidal with the frequency identical to the excitation frequency.
The proposed method is much more efficient in computational time than the Duhamel integral method, as shown in Fig. 10. In particular, the efficiency of the proposed method increases with the number of the time steps. As the system poles and residues are computed prior to the dynamic response calculation, the computational time for them are excluded in Fig. 10.

Earthquake loading
This example considers the oscillator characterized by Eq. 44 to earthquake-induced excitation f (t) = −mẍ(t), where the chosenẍ(t) is the measured E1 Centro earthquake acceleration signal in the EW direction, sampled with t = 0.02 for N s = 2048 steps (see Fig. 11). In view of the complexity of this f (t), the FFT algorithm is utilized to decompose it into a pole-residue form. The outcome of the FFT includes 1023 pairs of complex Fourier coefficients C m at ω m = (2π m)/(N s t), m = ±1, · · · , ±1023, along with two real Fourier coefficients at zero and the Nyquist frequency, respectively. For the pole-residue form off (s), iω m are the preselected excitation poles and C m are the corresponding residues. Response calculation in the frequency domain While conducting the proposed method in the frequency domain,ω n in Eq. 17, or ω, can be selected freely. Here, choosing ω = ω is for demonstrating the accuracy of the proposed method when it involves second-order poles. Using FFT automatically imposes the frequency resolution ω = 2π/T where T is the signal duration. For increasing the frequency resolution, padding sufficient zeros at the end of the original signal has been a common technique. After padding zeros to triple the length of the excitation, the poles for both the excitation and system have been designated at iω n = i(2π n)/ (3N s t), for n = −3072, · · ·, 3071. The required system residues A n = H n ω/(2π) are calculated accordingly from the discretized FRF (see Eq. 13). From Eq. 29, the response x(t) is computed,  as shown in Fig. 12, which is in excellent agreement with that obtained by the Duhamel integral approach. For completeness, the two components of Eq. 29 are plotted in Fig. 13. Note that the computed x(t) is only meaningful for t < T = 122.88 seconds, which is 3 times the original earthquake signal length, due to the usages of "periodic" oscillator and loading. Actually, when t becomes large, the computed x(t) from Eq. 29 will diverge because the first term of Eq. 29 is linearly proportional to t. Response calculation in the Laplace domain As the pole-residue forms of the excitation and system have been expressed as a Fourier series and Eq. 45, respec-  Figure 14 shows the obtained x(t), which is in excellent agreement with that computed by the Duhamel integral in the time domain. For completeness, the natural response and the steadystate response are also plotted in Fig. 15.

Example 2:
Fractional oscillator with two irrational fractional derivative terms: q 1 = π/30 and q 2 = 2π/7 Unlike most traditional methods, the proposed frequency/Laplace methods is not limited to oscilla-  Fig. 16 Response of the SDOF fractional oscillator with q 1 = π/30 and q 2 = 2π/7 to the earthquake loading based on the frequency domain and Duhamel integral methods tors with rational fractional orders. To demonstrate its broader application, the oscillator considered in this example contains two fractional derivative terms: q 1 = π/30 and q 2 = 2π/7. Let other parameters be m = 5/2, k = 8/5, c 1 = 1/2, and c 2 = 6/5; as a result, the oscillator's FRF is: In the frequency domain approach, the poles assigned at iω n and residues A n = H n ω/(2π) are computed based on ω = π/256 and N n = 2048. The same numerical verifications conducted in previous examples are repeated here. Figure 16, which plots the response caused by the same earthquake loading utilized in the previous example, shows a good agreement between the proposed method and the Duhamel method. For completeness, the two components of Eq. 29 are also plotted in Fig. 17.
To apply the Laplace domain approach, the poleresidue form ofh(s) with true poles must be determined first. As this method approximates a pole-residue form for h a (t) by the data-driven Prony-SS method, its implementation will not be affected by whether the fractional derivatives are rational or irrational. Conducting the same procedure yields the approximated IRF with μ 1 = −0.2899 ∓ 1.0081i, β 1 = −0.0188 ± 0.2230i, μ 2 = −0.1352 and β 2 = 0.0246; the two components of h a (t) are plotted in Fig. 18. Further- x(t) First sum term Second sum term Fig. 17 Two response components of the SDOF fractional oscillator with q 1 = π/30 and q 2 = 2π/7 to the earthquake loading by the frequency domain method more, Fig. 19 shows a good agreement between the original IRF and the reconstructed h a (t). Lastly, the response due to the earthquake loading computed by the Laplace domain approach is shown in Fig. 20, which has been verified by the Duhamel integral approach. The computation efficiency of the proposed method compared to the Duhamel integral method is demon- where the mass matrix M is a diagonal matrix with its diagonal elements m 1 , m 2 and m 3 ; the stiffness matrix K is formulated as: and the damping matrix C shares the same form as K, but replacing k n with the corresponding c n . In the numerical investigation, the following values are used: m 1 = m 2 = 3.50 × 10 5 , m 3 = m 1 /2, k 1 = k 2 = 2.10 × 10 8 , k 3 = k 1 /2, and c n = 3.00 × 10 5 for n = 1, 2 and 3. The proposed method can apply to this 3-DOF system as long as its matrix frequency response function H(ω), which can be obtained numerically from H(ω n ) = [−ω 2 n M+(iω n ) 0.7 C+K] −1 for each discrete ω n , is available. For illustration purpose, consider the output of interest x 3,1 is the displacement at coordinate 3 due to the excitation at coordinate 1. The (3, 1) entry of H(ω) denoted H 31 (ω), i.e., the frequency response function corresponding to the output coordinate 3 and input coordinate 1, is shown in Fig. 22. Following the same procedure as that for a SDOF system, one obtains the pole-residue form for the corresponding transfer functionh 31 (s) = N n n=1 β 31,n s−μ n where the three pairs of dominant system poles μ n and their corresponding residues β 31,n are listed in Table 2. As the pole-residuẽ h 31 (s) is available, the calculation for x 3,1 (t) can follow that has been developed for SDOF systems. Assuming the system is at rest initially and referring to Eqs. 36-38, one writes or where N n = 6. Consider the loading is the same as that in Example 1, i.e., N = 2, λ 1,2 = ±i and α 1,2 = ∓0.5i. The resulting x 3,1 (t) from Eq. 52 and the Duhamel integral method are compared in Fig. 23, which shows an excellent agreement between them.

Concluding remarks
A pole-residue method, operated in either the frequency or the Laplace domain, to efficiently compute the response of a fractional oscillator (FO) to complicated irregular loadings was developed and tested. Unlike many existing approaches that had strict limitations on the fractional order q, the proposed method was applicable to any q between 0 and 1. Despite the fact that the exact transfer function (TF) and frequency response function (FRF) of a FO had been derived analytically, using them directly to compute the FO's transient response was not possible. This paper developed two distinct approaches to convert the exact TF into its pole-residue form. The first approach was a frequency domain method that designated "artificial" poles along the imaginary axis. The second approach was a Laplace domain method that extracted the "true" poles through two steps: (1) a discrete IRF was produced by taking the inverse Fourier transform of the corresponding FRF that was readily obtained from the Code availability All code that support the findings of this study are available from the corresponding author by email (jameshu@uri.edu).

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Ethical approval Not applicable.

Consent to participate Not applicable.
Consent for publication Not applicable.

Appendix I: Summary of the Prony-SS method
Given the IRF signal h(t) that are sampled with interval t, h(t k ) = h(k t), k = 0, 1, · · · , N − 1, using the Prony-SS method is to approximate h(t) as N n n=1 β n exp(μ n t). Summarized here is the three sequential steps of the Prony-SS method to determine N n , μ n and β n , respectively.
where ξ + η = N , and a better choice is both ξ and η are close to N /2. Applying the singular value decomposition of H(0) gives Where U 1 ∈ R ξ ×N n , S 1 ∈ R N n ×N n and V 1 ∈ R η×N n ; theoretically, N n is the rank of H(0), i.e., the number of nonzero singular values in Eq. 54. Mathematically, the singular values go to zero when the rank of the matrix is exceeded; however, for data involving random errors or inconsistencies, some singular values will be very small, but not exactly zero. In this situation, N n could be determined based on the magnitudes of the singular values having been ordered sequentially from the largest to the smallest; a conventional way to determine N n is based on a significant drop of the normalized singular values.
Step 2 for determining μ n : Construct the Hankel matrix H(1) ∈ R ξ ×η from the sampled signal h(t k ): After computing the eigenvalues z n , n = 1, · · · , N n , of A, one can get μ n = ln(z n )/ t.
Step 3 for determining β n : Compute the complex coefficients β n by solving a set of linear equations using a least-square procedure, based on the obtained z n and the sampled h(t k ): ⎡