A Closed-form Approximation for Impulse Response of Fractionally Damped Oscillators

TL;DR

This paper introduces a computationally efficient and robust closed-form approximation for the impulse response of fractionally damped oscillators, overcoming numerical instabilities of the infinite series representation.

Contribution

The authors develop a novel closed-form approximation for the impulse response of fractional oscillators, improving computational stability and efficiency over the traditional series approach.

Findings

Approximation reduces computational instability.

Maintains reasonable accuracy in practical applications.

Enhances efficiency for engineering use.

Abstract

We consider a fractionally damped oscillator, where the damping term is expressed by the Caputo fractional derivative of order $\beta\in (0,1).$ The impulse response of this oscillator can be expressed in terms of the bivariate Mittag-Leffler function consisting of a double infinite series. Although this series is uniformly convergent, its numerical implementation suffers from computational instabilities. In this contribution, we propose an approximate closed-form solution that avoids these numerical pitfalls while maintaining a reasonable accuracy. The resulting approximation is computationally efficient and robust, making it suitable for practical engineering applications.

Equations

$$\left(D^{\beta }x\right)(t):=\left({}^C{D}_{0+}^{\beta }x\right)(t)=\frac{1}{\Gamma (1-\beta )}{\int }_0^t(t-s{)}^{-\beta }{x}^{\prime }(s)\mathrm{d}s,$$

$$\ddot{x}(t)+2\zeta {\omega }_n^{2-\beta }{D}^{\beta }x(t)+{\omega }_n^{2}x(t)=h(t),t>0,$$

$$x(0)=0,\dot{x}(0)=0.$$

$$E_{(a_1,a_2),b}(z_1,z_2)=\sum _{k=0}^{\infty }\sum _{\begin{array}{c}{k}_1+k_2=k\\ k_j\ge 0\end{array}}\left(\genfrac{}{}{0ex}{}{k}{k_1,k_2}\right)\frac{{\prod }_{j=1}^2{z}_j^{k_j}}{\Gamma (b+{\sum }_{j=1}^2{a}_j{k}_j)}=\sum _{m=0}^{\infty }\sum _{l=0}^{\infty }\frac{\left(\genfrac{}{}{0ex}{}{m+l}{l}\right)z_1^l{z}_2^m}{\Gamma (b+a_{1}l+a_{2}m)}.$$

$$I_{\beta }(t)$$

$$=t{E}_{(2,2-\beta ),2}(-({\omega }_{n}t{)}^2,-2\zeta ({\omega }_{n}t{)}^{2-\beta }).$$

$$x(t)={\int }_0^{t}h(t-s)I_{\beta }(s)\mathrm{d}s,t>0,$$

$$\ddot{x}(t)+(1+2\zeta ){\omega }_n^{2}x(t)=h(t).$$

$${\omega }_d={\omega }_n\sqrt{1+2\zeta }.$$

$$I_0(t)=\frac{\sin ({\omega }_{d}t)}{{\omega }_d}.$$

$$I_0(t)=\sum _{k=0}^{\infty }(-1{)}^k\frac{t^{2k+1}}{(2k+1)!}\sum _{m=0}^k(\genfrac{}{}{0ex}{}{k}{m}){\left(2\zeta {\omega }_n^2\right)}^m({\omega }_n^2{)}^{k-m},$$

$$\ddot{x}(t)+2\zeta {\omega }_n\dot{x}(t)+{\omega }_n^{2}x(t)=h(t),$$

$$I_1(t)=e^{-\zeta {\omega }_{n}t}\frac{\sin ({\omega }_{d}t)}{{\omega }_d},$$

$${\omega }_d={\omega }_n\sqrt{1-{\zeta }^2}$$

$$I_1(t)=\sum _{m=0}^{\infty }\sum _{l=0}^{\infty }\frac{(l+m)!}{m!l!}\frac{(-1{)}^{l+m}{\left(2\zeta \right)}^m({\omega }_{n}t{)}^{2l+m}t}{(2l+m+1)!},$$

$$x(t)=t\sum _{m=0}^{\infty }{\left(-2\zeta ({\omega }_{n}t{)}^{2-\beta }\right)}^m{E}_{2,2+(2-\beta )m}^{(m+1)}(-{\omega }_n^2{t}^2),$$

$$E_{\alpha ,\beta }^{(\gamma )}(z)=\frac{1}{\Gamma (\gamma )}\sum _{l=0}^{\infty }\frac{\Gamma (l+\gamma )z^l}{l!\Gamma (\alpha l+\beta )}$$

$${\tilde{I}}_{\beta }(t)=\exp \left(-{\zeta }_{\text{eq}}{\omega }_{n}t\right)\frac{\sin \left({\omega }_{d,\text{eq}}t\right)}{{\omega }_{d,\text{eq}}},$$

$${\zeta }_{\text{eq}}(\beta )$$

$${\omega }_{d,\text{eq}}({\omega }_n,\zeta ,\beta )$$

$${\omega }_d=\{\begin{array}{rl}{\omega }_n\sqrt{1-{\zeta }^2}& \text{ for }\beta =1,\\ {\omega }_n\sqrt{1+2\zeta }& \text{ for }\beta =0.\end{array}$$

$$s^2+2\zeta {\omega }_n^{2-\beta }{s}^{\beta }+{\omega }_n^2=0.$$

$$Y_{\beta }=\frac{M_0-M_{\beta }}{M_0-M_1},$$

$$Y_{\beta }={\beta }^{(2.24-0.63\beta )}.$$

$${\omega }_{d,\text{eq}}(\beta )={\omega }_n\sqrt{1+2\zeta -\zeta (2+\zeta ){\beta }^{(2.24-0.63\beta )}},\beta \in [0,1].$$

$${\zeta }_{\text{eq}}=\frac{\delta }{\sqrt{{\delta }^2+4{\pi }^2}},$$

$${\zeta }_{\text{eq}}=\zeta {\beta }^{(0.95-0.85\beta )}.$$

$$H(\omega )=\frac{1}{-{\omega }^2+2\zeta {\omega }_n^{2-\beta }(i\omega {)}^{\beta }+{\omega }_n^2},$$

$${\omega }_n^{2}H(\omega ):=h(g)=\frac{1}{1-g^2+2\zeta (ig{)}^{\beta }}.$$

$$\tilde{H}(\omega )={\int }_{-\infty }^{\infty }{\tilde{I}}_{\beta }(\tau )e^{-i\omega \tau }d\tau =\frac{1}{{\omega }_{d,\text{eq}}^2+({\zeta }_{\text{eq}}{\omega }_n+i\omega {)}^2},$$

$${\omega }_n^2\tilde{H}(\omega ):=\tilde{h}(g)=\frac{1}{({\omega }_{d,\text{eq}}/{\omega }_n{)}^2+({\zeta }_{\text{eq}}+ig{)}^2}.$$

Full text

K. Van Bockstal is supported by the FWO Senior Research Grant G083525N

1. 1 Introduction

2. 2 Series Solution

3. 2.1 $I_{0}(t)$: Limiting case of $I_{\beta}(t)$ when $\beta=0$

4. 2.2 $I_{1}(t)$: Limiting case of $I_{\beta}(t)$ when $\beta=1$

5. 2.3 Convergence and blow-up issues with $I_{\beta}(t)$ when $0<\beta<1$

6. 3 Approximate solution

7. 3.1 Equivalent Frequency of Fractional Oscillators

8. 3.2 Equivalent Damping of Fractional Oscillator

9. 4 Error Analysis

10. 4.1 Comparison in Frequency-domain

11. 4.2 Comparison with results of the eigenvector expansion approach

12. 5 Response to Arbitrary Excitation using Convolution

13. 6 Conclusion

14. A Derivation of impulse response

15. B Limiting case $\beta\searrow 0$

16. C Limiting case $\beta\nearrow 1$

17. D Finite Difference Approach

1. Introduction

Modern structural engineering components use polymers, elastomers, viscoelastic dampers, and viscoelastic damping treatment for enhanced vibration resilience. These materials exhibit frequency-dependent and memory-effect (or hereditary) behavior. By now, it is well proven that the fractional (or non-integer) derivative-based models satisfactorily capture the mechanical behavior of systems involving such materials (see [1, 2, 3] and references therein). The fractional derivative models are compact, analytic, causal, and better fit the experimental data with a smaller number of calibration parameters compared to classical viscoelastic models and complex modulus methods [4, 5].

The use of fractional material models leads to the fractional differential equation (FDE) of motion [6]. Comprehensive reviews detailing the analysis of these FDEs are available in [7, 8]. A range of techniques is available to solve fractional systems, for example, the eigenvector expansion approach [9], Galerkin approach [10], numerical methods [11, 12], and the integral transform (Laplace and Fourier transforms) based approaches [13, 14, 15]. Among these methods, the Laplace transform-based approaches are widely used by the engineering communities. This approach requires the inversion of the fractional-order transfer functions. Mostly, the studies employ Cauchy’s residue theorem when evaluating the inverse transforms, which requires the numerical computation of complex integrals. Some methods, such as Oustaloup’s method and the continued fraction approach, approximate fractional-order transfer functions by higher-order polynomials (see [16] and references therein). However, higher-order polynomials lead to the loss of the intuitive understanding and computational simplicity.

Alternatively, the inverse transforms can also be expressed as an infinite series (involving Mittag-Leffler type functions) when differential equations have constant coefficients. However, such series solutions are generally not favoured in practical applications due to slow convergence, potential numerical instability (e.g., blow-up at large time-instants), and implementation difficulties [8, 15]. For example, Duan et al. [15] attempted to investigate the impulse response of fractionally damped oscillators using the series solution. However, due to the slower convergence at higher time instants, they reverted to using the Laplace inversion method with complex integration via the Gauss–Laguerre quadrature scheme.

In this contribution, we consider a fractionally damped oscillator as shown in Figure 1. The fractional derivative is expressed by the Caputo fractional derivative $D^{\beta}$ of order $\beta\in(0,1)$ defined as [17]:

[TABLE]

where $\Gamma(\cdot)$ is the gamma function and $x^{\prime}(s)$ is the first derivative of $x(s)$ with respect to $s$.

In [6], a detailed derivation is presented that shows that the fractional material models lead to the fractional governing equation of motion

[TABLE]

where $x(t)$ is the displacement response (unit: m) of the oscillator as a function of time $t$ (unit: s), $\ddot{x}(t)$ is the double derivative of $x(t)$ with respect to time, $\zeta\in(0,1)$ is the nondimensional damping ratio, and $\omega_{n}$ is the natural frequency of the oscillator (in rad/s). The function $h(t)$ is the applied excitation (acceleration in m/s2) to the oscillator. The oscillator is initially at rest as shown in Figure 1, i.e., the initial conditions are given by

[TABLE]

It may be noted that $\omega_{n}$ has been raised to an exponent equal to $2-\beta$ due to dimensional consistency (as $D^{\beta}x$ has unit m/sβ), see [6]. If $\beta=1$, this exponent becomes $1$, which is compatible with the classical viscously damped case, noting that $D^{1}x(t)=\dot{x}(t)$.

In this study, motivated by the aforementioned limitations of the series solutions, we will derive an approximate closed-form expression for the impulse response of fractionally damped systems governed by Equations 1 and 2. The proposed solution extends the classical impulse response of integer-order systems to the fractional case and can be convolved with arbitrary excitations to approximate the time-domain response of the system, making it suitable for practical engineering applications.

The paper is structured as follows. In Section 2, we first derive the series solution to problem (1-2). We then examine the limiting cases $\beta=0$ and $\beta=1$ in Sections 2.1 and 2.2. We show that the analytical summation of the resulting series coincides with the classical solution for the two extreme cases. Afterwards, in Section 2.3, we discuss the numerical instabilities associated with the series solution when $0<\beta<1$ and propose a stable closed-form approximation in Section 3. The equivalent frequency and damping of this proposed solution are thoroughly discussed in Sections 3.1 and 3.2. Moreover, in Section 4, we perform an error analysis between the proposed approximation and the series solution (or exact solution). Finally, in Section 5, we demonstrate the applicability of the developed closed-form equation for solving a fractionally damped system with arbitrary excitation.

2. Series Solution

To solve the problem (1-2), we will first derive the impulse response of the oscillator and then, using the convolution integral, we shall obtain the desired response for the applied excitation $h(t)$.

The impulse response $I_{\beta}(t)$ is the solution to (1) with $h(t)=0$ satisfying $x(0)=0$ and $\dot{x}(0)=1.$ This solution can be obtained by applying the Laplace transform method. For more details, we refer the reader to Appendix A. It can be expressed in terms of the bivariate Mittag-Leffler function, which is defined by [18, Eq. (3.16)]

[TABLE]

Here, $\binom{k}{k_{1},k_{2}}=\frac{k!}{k_{1}!k_{2}!}$ is the multinomial coefficient and $\binom{m+l}{l}=\frac{(m+l)!}{m!l!}$ is the binomial coefficient. In particular, it holds that the impulse response $I_{\beta}(t)$ is given by

[TABLE]

Therefore, for any arbitrary fixed $T>0$, the solution $x$ to problem (1-2) for given $h\in\operatorname{C}([0,T])$ belongs to $\operatorname{C}^{2}([0,T])$ and is given by

[TABLE]

see, e.g., [19, Corollary 4.9].

2.1. $I_{0}(t)$: Limiting case of $I_{\beta}(t)$ when $\beta=0$

For $\beta=0$, the governing equation becomes

[TABLE]

This equation does not have an explicit damping term. However, the damping ratio $\zeta$ contributes to the stiffness of the oscillator, causing the effective natural frequency $\omega_{d}$ to become

[TABLE]

From classical structural dynamics, we know that the impulse response of such an undamped ordinary oscillator is

[TABLE]

Using the series solution Equation 4 for $\beta=0,$ we have that

[TABLE]

which exactly matches the classical dynamics solution (6). For more details, the reader may refer to Appendix B.

2.2. $I_{1}(t)$: Limiting case of $I_{\beta}(t)$ when $\beta=1$

For $\beta=1$, the governing equation takes the classical form (viscous damping)

[TABLE]

whose impulse response is well-known and given by

[TABLE]

where

[TABLE]

denotes the effective or damped frequency (valid for $0<\zeta<1$). Substituting $\beta=1$ into the series (4) yields

[TABLE]

which recovers (7). A rigorous derivation of this limiting case is provided in Appendix C.

2.3. Convergence and blow-up issues with $I_{\beta}(t)$ when $0<\beta<1$

Although the series (4) is theoretically uniformly convergent (for $t\in[0,T]$ for arbitrary but fixed $T$), its numerical implementation suffers from significant computational instabilities. As illustrated in Figure 2, if we take $\omega_{n}=10$ rad/s (let $\beta=0.7$ and $\zeta=0.05$), then we observe that the solution diverges (blows up) after approximately 3.5 seconds. This pattern scales with $\omega_{n}$, occurring after $t=35$ s if $\omega_{n}=1$ rad/s.

For $\omega_{n}=10$ rad/s, at $t=3.5$ s, we need to sum at least till $m=13$ and $l=51$, in eq. 4, to obtain an acceptable accuracy. For these values of $m$ and $l$, we have in the denominator of eq. 4, $\Gamma(120.9)$ which is of the order of $4.143\times 10^{198}$ and in the numerator, we get $64!$ which is around $1.2689\times 10^{89}$, and $35^{118.9}\approx 3.8877\times 10^{183}$. These values are so extreme that the limits of floating-point arithmetic are pushed. To achieve better convergence, if one tries further to increase the range of $m$ and $l$, the problem takes the numerical form of an indeterminate $\infty/\infty$. All of this results in a catastrophic loss of precision.

These numerical limitations explain why series solutions remain impractical for solving problem (1-2). Our work addresses this issue by developing an approximate closed-form solution for the impulse response of the fractional oscillator that avoids these numerical pitfalls while preserving the fundamental dynamical characteristics of the original system.

Remark 2.1*.*

Note that the impulse response (4) can also be written as

[TABLE]

where

[TABLE]

is the Prabhakar generalisation of the Mittag-Leffler function. In [20], the authors developed an approach for its implementation in MATLAB (see [21]) based on the numerical inversion of its Laplace transform. However, this approach is not directly applicable here as it is restricted for $\gamma\neq 1$ to $\alpha\in(0,1)$.

3. Approximate solution

In Section 2, we have shown that the impulse response can be represented in terms of the double-infinite series (eq. 4). We have also demonstrated that when $\beta=0$ and $\beta=1$ (see Sections 2.1 and 2.2, respectively), the series converges to the classical impulse solutions consisting of an exponentially decaying term (due to damping $\zeta$) and a sinusoidal function with effective frequency $\omega_{d}$. The analytical summation of the series was possible as the Gamma functions convert to the factorial functions if $\beta=0$ and $1$. Moreover, we have demonstrated in Section 2.3 that the series (4) suffers numerical instabilities.

Our main goal is to propose and validate a closed-form approximation of (4) that is consistent with the classical cases $\beta=0$ and $1$. From Figure 3, we can see that $\beta=0$ corresponds to the undamped oscillator with an amplified spring stiffness, whereas $\beta=1$ corresponds to the classical oscillator with viscous damping. Thus, it is intuitive that the oscillator with $\beta$ somewhere in between [math] and $1$ will exhibit the combined behaviour of a spring and a viscous dashpot. Therefore, the solution for a fractional oscillator can also be approximated to consist of:

1. (i)

an exponentially decaying term with an equivalent damping $\zeta_{\text{eq}}$ as a function of $\beta$ and $\zeta$, and 2. (ii)

sinusoidal term with an equivalent damped frequency $\omega_{d,\text{eq}}$ as a function of $\omega_{n}$, $\beta$ and $\zeta$.

For these reasons, we propose the following closed-form approximation

[TABLE]

where the equivalent damping ratio $\zeta_{\text{eq}}$ and equivalent damped frequency ${\omega}_{d,\text{eq}}(\beta)$ are defined as

[TABLE]

Note that $\tilde{I}_{0}(t)=I_{0}(t)$ and $\tilde{I}_{1}(t)=I_{1}(t)$, so the approximate impulse response coincides with the solution for $\beta=0$ and $\beta=1$, as required.

In contrast to the impulse response series (4), the proposed approximation (9) consists of the product of an elementary exponential and sine function. Since these functions are implemented in standard numerical libraries using optimised and stable algorithms, the approximation avoids the pitfalls associated with Gamma function computations. The proposed closed-form approximation is more reliable for computations and also faster.

In the next Subsections (3.1 and 3.2), we will motivate the definitions for $\zeta_{\text{eq}}$ and ${\omega}_{d,\text{eq}}.$ Afterwards, in Section 4, we will perform an error analysis between the impulse response (4) and its approximation (9).

3.1. Equivalent Frequency of Fractional Oscillators

In Sections 2.1 and 2.2, we defined $\omega_{d}$ for the limiting cases as

[TABLE]

For $\beta\in[0,1],$ the apparent frequency $\omega_{d}$ corresponds to the magnitude of the imaginary part of the root of the characteristic equation (see Equation 22 in Appendix A)

[TABLE]

Numerical methods are required for the evaluation of $\omega_{d}$ as the roots of this transcendental equation cannot be determined analytically for $\beta\in(0,1)$. We used the in-built function ‘fsolve’ of MATLAB [22], which employs the trust-region dogleg algorithm to solve the equations.

It is interesting to note that the fractional damper is an element that has the properties of a viscous damper and a spring depending on the value of $\beta$. For example, when $\beta=1$, it is a purely Newtonian viscous dashpot, and in that case the damped frequency $\omega_{d}=\omega_{n}\sqrt{1-\zeta^{2}}\leq\omega_{n}.$ When $\beta=0$, there is no damping in the system, and the fractional dashpot behaves like a spring such that the natural frequency of the system becomes $\omega_{d}=\omega_{n}\sqrt{1+2\zeta}\geq\omega_{n}.$ When $0<\beta<1$, there is a partial spring nature in the fractional dashpot, which would increase the system stiffness (and hence there will be an increase in frequency), and the remaining partial viscous damping component would attempt to damp (reduce) the frequency. Let us define this effective frequency as $\omega_{d,\text{eq}}$. From the above explanation, it is obvious that $\omega_{n}\sqrt{1+2\zeta}\geq\omega_{d,\text{eq}}\geq\omega_{n}\sqrt{1-\zeta^{2}}$ and $\omega_{d,\text{eq}}$ decreases gradually from $\omega_{n}\sqrt{1+2\zeta}$ to $\omega_{n}\sqrt{1-\zeta^{2}}$ with increasing $\beta$ from 0 to 1, as can also be seen from the numerical results shown in Figure 4. Here, we also refer to two interesting papers [23, 24], where variations of $\omega_{d}$ with $\beta$ were studied, and similar trends were noted.

To obtain an equivalent expression ${\omega}_{d,\text{eq}}$, we define $M_{\beta}:=(\omega_{d}/\omega_{n})^{2}$. Note that $M_{0}=1+2\zeta$ and $M_{1}=1-\zeta^{2}.$ To determine $M_{\beta}$ as a function of $\beta$ and $\zeta$, we introduce a variable $Y_{\beta}$ as follows

[TABLE]

which scales the transition between $\beta=0$ and $\beta=1.$ We plot the scatter of $Y_{\beta}$ with $\beta$ for randomly generated 10,000 statistically independent and uniformly distributed samples of $\omega_{n},\zeta,\beta$. The $\beta$ varies between 0 and 1, $\zeta$ varies from 0.001 to 0.15 (covering the practical ranges), and $\omega_{n}$ varies from 1 to 10 rad/s (covering the practical range for a variety of civil engineering structures). After conducting several numerical experiments, we observed that the regression model of the form $Y_{\beta}=\beta^{(A_{0}-A_{1}\beta)}$ fits well to the scatter, where $A_{0}$ and $A_{1}$ are regression parameters. Note that this form captures the power function transition observed in the numerical data. The fractional calculus is based on the generalization of exponential kernels by power-law-type kernels, and, therefore, it appears natural to have a power function here. The regression parameters are obtained using the ‘fitnlm’ command in MATLAB [25], which uses the Levenberg-Marquardt least squares approach (considering $A_{0}=A_{1}=1$ as initial guesses). From these experiments, we conclude that $A_{0}=2.2380$ (with a 95% confidence interval $[2.2320,2.2441]$) and $A_{1}=0.6316$ (with a 95% confidence interval $[0.6217,0.6415]$). Hence, we have

[TABLE]

Consequently, using eqs. 14 and 15, we get that the apparent frequency $\omega_{d}$ is approximated by the equivalent frequency $\omega_{d,\text{eq}}$ as follows

[TABLE]

Figure 4(a) shows that the proposed eq. 15 fits reasonably well with the numerically obtained data for 100 randomly generated test samples. Figure 4(b) shows a reasonable match between numerically obtained $\omega_{d}$ and proposed $\omega_{d,\text{eq}}$ for three different damping levels and $\omega_{n}=1$ rad/s. Note that the expression matches with the limiting cases $\beta=0$ and $\beta=1.$

3.2. Equivalent Damping of Fractional Oscillator

We already know that the damping ratio is $\zeta$ if $\beta=1$ (pure viscous damping) and [math] if $\beta=0$ (undamped system). However, when $\beta\in(0,1)$, the system is fractionally damped. Therefore, we need to define an equivalent viscous damping ratio $\zeta_{\text{eq}}$, as a function of $\zeta$ and $\beta$, to capture the amplitude decay rate of the fractionally damped oscillator in the classical sense.

Using the logarithmic decrement method [26, Section 2.3], we can estimate the equivalent damping $\zeta_{\text{eq}}$ in the system from a given response signal for $\beta\in(0,1)$. Let the first positive peak in the response $I_{\beta}(t)$ be $x_{1}$ and let the $j^{th}$ peak be $x_{j}$ (refer to Figure 5). If the natural period is $\omega_{n}$, then, following the logarithmic decrement method, we can write the equivalent damping $\zeta_{\text{eq}}$ as

[TABLE]

where $\delta=\dfrac{1}{j-1}\ln{\left(\dfrac{x_{1}}{x_{j}}\right)}$. Figure 5 demonstrates the computation of the equivalent damping ratio for an example case with $\beta=0.5$, $\zeta=0.05,$ and $\omega_{n}=5$ rad/s for $j=2,3,4,$ and 5. It is noted that the rate of response decay reduces slightly with increasing $j$, however, there is no significant difference in the values of $\zeta_{\text{eq}}$ (varies between 0.6879 and 0.6998). Therefore, for the sake of computational efficiency, this study considers only the first two cycles to evaluate $\zeta_{\text{eq}}$.

Now, we will plot $\zeta_{\text{eq}}/\zeta$ as a function of $\beta$ and propose a functional form for $\zeta_{\text{eq}}/\zeta$. For randomly generated 10,000 statistically independent and uniformly distributed samples of $\beta\in(0,1)$, $\zeta\in[0.001,0.15]$ and $\omega_{n}\in[1,10]$ rad/s, $\zeta_{\text{eq}}/\zeta$ is calculated as discussed above. We conducted several numerical experiments and fitted the regression model of the form ${\zeta_{\text{eq}}}/{\zeta}=\beta^{(A_{0}-A_{1}\beta)}$, where $A_{0}$ and $A_{1}$ are regression parameters. As in Section 3.1, this form captures the power function observed in the numerical data and satisfies the boundary conditions. From regression analysis, it is noted that $A_{0}=0.9506$ (with a 95% confidence interval $[0.9486,0.9527]$) and $A_{1}=0.8499$ (with a 95% confidence interval $[0.8452,0.8546]$) can be considered as the mean value of the regression parameters. Thus, we have

[TABLE]

It is worth mentioning here that the numerical experiments did not indicate any significant relation between $\zeta_{\text{eq}}$ and $\omega_{n}$. Figure 6 shows that the proposed function (18) fits reasonably well with the numerically obtained data for 100 randomly generated test samples.

4. Error Analysis

In this section, we inspect the performance of the proposed closed-form approximation (9) for the impulse response of the fractional oscillator. The exact series solution (4) and the approximate solutions are compared for the following four typical cases:

1. (i)

$\omega_{n}=1$ rad/sec, $\beta=0.1$, $\zeta=0.01$, 2. (ii)

$\omega_{n}=10$ rad/sec, $\beta=0.9$, $\zeta=0.15$, 3. (iii)

$\omega_{n}=1$ rad/sec, $\beta=0.5$, $\zeta=0.15$, 4. (iv)

$\omega_{n}=5$ rad/sec, $\beta=0.5$, $\zeta=0.05$.

The comparison of responses is shown in Figure 7, which indicates that the proposed solution is reasonably accurate. Furthermore, we also show the residual in the approximate solution (series solution minus approximate solution) in Figure 8. The residuals are small enough compared to the exact solution.

4.1. Comparison in Frequency-domain

In addition to the time-domain comparison, we also carry out a comparison of frequency response functions. The frequency response $H(\omega)$ of the exact series solution is obtained by substituting $s=i\omega$ into eq. 22 (Appendix A), where $i=\sqrt{-1}$ and $\omega=$ frequency. This gives

[TABLE]

which can be further simplified in terms of non-dimensional frequency $g=\omega/\omega_{n}$ as

[TABLE]

To obtain the frequency response $\tilde{H}(\omega)$ of the proposed approximation eq. 9, we use the definition of the Fourier transform, which leads to

[TABLE]

which in terms of $g$ becomes

[TABLE]

The absolute values of the frequency response functions ($\lvert h(g)\rvert$ and $\lvert\tilde{h}(g)\rvert$) are plotted in Figure 9 for the four sets of $\beta$ and $\zeta$ considered previously. A reasonable match is observed between the two solutions. The approximate solution captures the resonance accurately; however, some minor differences are noted in the static response (near $g=0$), especially when the damping $\zeta$ is high and $\beta$ is away from the extremities (0 and 1). This clearly reinforces the robustness of the approximation in the frequency domain as well.

4.2. Comparison with results of the eigenvector expansion approach

Suarez and Shokooh [9] analysed eq. 1 in the special case $\beta=0.5$. They used the eigenvector expansion approach to obtain the impulse response for $\omega_{n}=10$ rad/s and several damping ratios. We digitised the results corresponding to $\zeta=0.05$ and $1$, contained in their paper, and compared them here with our approximation in Figure 10.

A reasonable match is obtained for $\zeta=0.05$. However, for $\zeta=1$, there is an agreeable qualitative match with some noticeable quantitative deviation because our approximation has been developed for practical values of $\zeta$ which are usually far less than 1. Additionally, the approximated impulse response for $\zeta=1$ depicts an oscillatory nature, which is consistent with the observations of [9, 23, 24] where it is discussed that the fractionally damped oscillators do not strictly exhibit the classical overdamped behaviour (i.e., no oscillations when $\zeta=1$).

In a nutshell, all of the above comparisons clearly demonstrate that the proposed approximation is practically acceptable for engineering applications.

5. Response to Arbitrary Excitation using Convolution

In this section, we demonstrate the application of the developed closed-form approximation $\tilde{I}_{\beta}$ (see eq. 9) by using it along with the convolution technique (see eq. 5) to obtain the responses for any arbitrary excitation.

We consider the same example as in [27], which is given by

[TABLE]

Comparing with the governing equation (1), we get $2\zeta\omega_{n}^{2-\beta}=0.4$ (rad/s)1.44, $\beta=0.56$, $\omega_{n}^{2}=2$ (rad/s)2, and $h(t)=30\cos{(6t)}$ m/s2. Thus, $\omega_{n}=\sqrt{2}$ rad/s and $\zeta=0.1214$. We can now convolve excitation $h(t)$ with the proposed approximate impulse response (9) and the exact series solution (4) to obtain the oscillator response denoted by $x_{p}(t)=\int_{0}^{t}h(t-s)\tilde{I}_{\beta}(s)\,\mathrm{d}s$ and $x_{s}(t)=\int_{0}^{t}h(t-s)I_{\beta}(s)\,\mathrm{d}s$, respectively.

The two responses are compared in Figure 11(a) and the residual $x_{s}(t)-x_{p}(t)$ is plotted in Figure 11(b). A reasonable agreement is obtained between the two responses with negligible residuals throughout. Since the series solution blows up beyond $t=22$s (see Figure 11(c)), a quantitative comparison could be carried out till $t=22$s only.

To evaluate the accuracy of the approximated response $x_{p}(t)$ for $t>22s$, we compared it with the finite-difference solution denoted by $x_{f}(t)$. The details of how to obtain the finite-difference (FDM) solution can be found in Appendix D. We used a time step $\Delta t=0.005s$ in the implementation of the FDM solution. As can be seen from Figure 12, the proposed approximation performs well in the region where the series solution is unable to provide a solution (note the residual $x_{f}(t)-x_{p}(t)$).

Thus, we conclude from this section that the proposed impulse approximation leads to reasonably accurate results when convoluted with arbitrary excitation.

6. Conclusion

In this work, we considered a fractionally damped oscillator governed by the Caputo derivative of order $\beta\in(0,1)$. The exact solution, expressed as a series of bivariate Mittag-Leffler functions, was shown to recover classical results in the limiting cases of $\beta=0$ (undamped oscillator) and $\beta=1$ (viscous damping). Although the series converges uniformly, its numerical evaluation suffers from instabilities that limit its practical utility. In addition to this, the classical structural dynamics with integer-order derivatives provides a well-known expression for the impulse response of the oscillators; however, no such simple formula exists for the fractional oscillators.

To address the above limitations, we proposed an approximate closed-form equation for the impulse response of fractional oscillators that avoids the computational challenges associated with the Mittag-Leffler function and extends the classical impulse response to the case of fractional oscillators. The equivalent frequency and damping of this approximation were rigorously analysed. It was shown that the proposed form is accurate, and its effectiveness in solving fractionally damped systems under arbitrary external excitation was demonstrated. This confirms its suitability for engineering applications.

Future work could explore the derivation of error bounds between the approximate and exact solutions, as well as the study of the case $\beta\in(1,2)$.

Appendix A Derivation of impulse response

We consider the fractional-order differential equation

[TABLE]

with initial conditions $x(0)=0$ and $\dot{x}(0)=1$. Taking the Laplace transform, we obtain

[TABLE]

Recall that for $\left\lvert z\right\rvert<1$ the geometric series gives $\frac{1}{1-z}=\sum_{k=0}^{\infty}z^{k}.$ Put $z=-2\zeta\omega_{n}^{2-\beta}s^{\beta-2}-\omega_{n}^{2}s^{-2}$. We have $\left\lvert z\right\rvert<1$ for sufficiently large $\left\lvert s\right\rvert.$ Hence, we obtain

[TABLE]

Next, employing the binomial theorem $(x+y)^{k}=\sum_{m=0}^{k}\binom{k}{m}x^{m}y^{k-m}$, we expand the $k$-th power and get

[TABLE]

As $\mathcal{L}^{-1}\left\{s^{-\gamma}\right\}=\frac{t^{\gamma-1}}{\Gamma(\gamma)}$ for $\Re(\gamma)>0$, we have

[TABLE]

We can easily check that, $x(0)=0$ and $\dot{x}(0)=1.$

Using $\sum_{k=0}^{\infty}\sum_{m=0}^{k}c_{k,m}=\sum_{m=0}^{\infty}\sum_{k=m}^{\infty}c_{k,m}$ and putting $l=k-m$ give

[TABLE]

Then, (23) can be written as

[TABLE]

Using the bivariate Mittag-Leffler function defined in eq. 3, the solution can be compactly written as

[TABLE]

Appendix B Limiting case $\beta\searrow 0$

For $\beta=0,$ we have from (23) that

[TABLE]

Using the binomial theorem, we get

[TABLE]

since $\sin(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)!}x^{2n+1}.$ This function $x(t)$ is the solution to the second-order differential equation

[TABLE]

with initial conditions $x(0)=0$ and $\dot{x}(0)=1$, i.e., it solves problem (21) with $\beta=0$ as $D^{0}x(t)=x(t)-x(0)=x(t).$ Note that taking the Laplace transform of (27) gives

[TABLE]

Using $\mathcal{L}^{-1}\left[\frac{1}{s^{2}+a^{2}}\right]=\frac{1}{a}\sin(at),$ we get (26), which can be written in terms of the effective natural frequency $\omega_{d}=\sqrt{1+2\zeta}\omega_{n}$ as

[TABLE]

Note that $\omega_{d}$ equals the magnitude of the imaginary part of the root of the denominator of (28).

Appendix C Limiting case $\beta\nearrow 1$

In this section, we will show that taking $\beta=1$ in (24) or (25) results in the solution

[TABLE]

to the second-order differential equation

[TABLE]

with initial conditions $x(0)=0$ and $\dot{x}(0)=1$, i.e., solving problem (21) with $\beta=1$. First note that

[TABLE]

As $\zeta<1$, the roots of the denominator are $s_{1,2}=-\zeta\omega_{n}\pm i\omega_{d}.$ Rewriting the denominator of (31) in terms of these roots results in $s^{2}+2\zeta\omega_{n}s+\omega_{n}^{2}=(s+\zeta\omega_{n})^{2}+\omega_{d}^{2}.$ Using the standard inverse Laplace transform identity $\mathcal{L}^{-1}\left[\frac{1}{(s+a)^{2}+b^{2}}\right]=\frac{1}{b}e^{-at}\sin(bt)$ gives (30).

Now, we will show that

[TABLE]

From eq. 3, we have that

[TABLE]

Using the Taylor series $\sin(x)=\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)!}x^{2n+1}$ and $e^{x}=\sum_{n=0}^{\infty}\frac{x^{n}}{n!}$, and substituting $\omega_{d}=\omega_{n}\sqrt{1-\zeta^{2}}$, we get

[TABLE]

Let $2l+m+1=N$ with $N\in{\mathbb{N}}$. This defines a relationship between $l$ and $m$ given by $m=N-2l-1.$ Since $m\geq 0$ and $l\geq 0$, the valid values of $l$ are $0\leq l\leq\left\lfloor\frac{N-1}{2}\right\rfloor,$ where $\left\lfloor\cdot\right\rfloor$ is the floor function. Hence, we have

[TABLE]

and

[TABLE]

We see that

[TABLE]

and (using the binomial of Newton)

[TABLE]

Now, from [28, Theorem A] or [29, Eq. (3.121)], we have

[TABLE]

Substituting this expression into eq. 34, and replacing $k$ by $l$, yields eq. 33. Therefore, we conclude that (32) is valid.

Appendix D Finite Difference Approach

For $x\in\operatorname{C}^{1}([0,T])$ satisfying $x(0)=0$, the Grünwald-Letnikov approximation for the Caputo fractional derivative $D^{\beta}x(t)$ is given by (see [30, Section 7.1 and 7.5])

[TABLE]

where $W_{\beta j}$ are weights defined as

[TABLE]

To determine the values of $W_{\beta j}$, one can use the following recurrence relation

[TABLE]

Hence, we get the following first-order finite-difference form for eq. 1

[TABLE]

Rearranging the terms leads to the following recursive relation

[TABLE]

where

[TABLE]

Acknowledgements The suggestions of Professor Arnaud Deraemaeker (Université Libre de Bruxelles, Belgium) on frequency domain comparisons are well acknowledged.