We present a novel method to solve the Bagley-Torvik equation by transforming it into ordinary differential equations (ODEs). This method is based on the equivalence between the Caputo-type fractional derivative (FD) of order 3/2 and the solution of a diffusion equation subjected to certain initial and boundary conditions. The key procedure is to approximate the infinite boundary condition by a finite one, so that the diffusion equation can be solved by separation of variables. By this procedure, the Bagley-Torvik and the diffusion equations together are transformed to be a set of ODEs, which can be integrated numerically by the Runge-Kutta scheme. The presented method is tested by various numerical cases including linear, nonlinear, nonsmooth, or multidimensional equations, respectively. Importantly, high computational efficiency is achieved as this method is at the expense of linearly increasing computational cost with the solution domain being enlarged.

## Introduction

Originally, the Bagley-Torvik equation was established by modeling the motion of a rigid plate immersed in Newtonian fluid [1]. It was later widely extended to model the dynamic responses of viscoelastically damped structures in which the viscoelasticity is described by a fractional derivative (FD) [2]. There are different definitions of FD [2–5], with various applications in different academic fields [6–8]; and we use the Caputo-type fractional derivative in this study [3].

Due to the difficulty in solving analytical solutions, a lot of numerical techniques have been developed to solve this equation since its presence. A large category of these approaches belongs to time-marching integrators. For instance, Diethelm and Ford [9] developed a fractional linear multistep method and a predictor–corrector (PC) method of the Adams type, respectively. Later, Diethelm et al. [10] extended the predictor–corrector method to more general fractional differential equations (FDEs). Ray and Bera [11] obtained an analytical solution by Adomian decomposition method. Zolfaghari et al. [12] used the enhanced homotopy perturbation method to find more accurate solution. Çenesiz et al. [13] employed the generalized Taylor collocation method to solving the Bagley–Torvik equation. Wang and Wang [14] solved the equation by changing it into a sequential FDE with constant coefficients.

Other solution approaches are based more or less on collocation technique such as the Haar wavelet method [15,16]. Al-Mdallal et al. [17] proposed a collocation-shooting method, which can be used to solve this equation equipped with boundary conditions. And Raja et al. [18] developed a stochastic technique. Numerical solution can also be obtained by the Bessel collocation method [19].

In addition, some methods make use of the expansion of the integral kernel in FD. Atanackovic and Zorica [20] established a new method based on a derived expansion formula for FDs. Krishnasamy and Razzaghi [21] suggested a numerical method via using fractional Taylor vector approximation. The method established by Gülsu et al. [22] is based on Taylor matrix method. More recently, Arqub and Maayah [23] constructed an iterative reproducing kernel algorithm.

As is well known, the key point to solve a FDE is how to deal with the FD. Generally speaking, time-marching and collocation-type numerical methods require repeated usage of the past responses in processing the convolutions describing FDs [24]. For these reasons, the computational cost increases in proportion to order *n*^{2} with *n* as the number of integration steps. To this end, a large number of contributions have been made to enhance the computation efficiency. By using the short memory principle [25], for instance, one can significantly reduce the computational cost to order *n*log(*n*) without losing substantial accuracy. This is very close to the cost of order *n* for solving ordinary differential equations (ODEs). Actually, for both spectrum based and memory-free algorithms, the computational cost has even been controlled to *n* [26,27].

In this paper, we will present a method with the computational cost of order *n* to solve FDEs of order 3/2 exemplified by the Bagley-Torvik equation. Our method is based on a physical meaning behind the Bagley-Torvik equation. Torvik and Bagley [1] talked about the relationship between this equation and a diffusion equation (partial differential equation (PDE)). Fitt et al. [28] proved that the FD of order 3/2 is the solution of a diffusion equation with certain initial and boundary conditions. Motivated by this, we will incorporate the FDE into the PDE and solve them together by the method of separation of variables, according to which a set of ODEs is deduced and solved by the Runge-Kutta scheme.

## Solution Procedures

## Numerical Examples

In this section, we will present several examples to test the presented method. To demonstrate the wide applicability, the considered numerical cases will include both linear and nonlinear, smooth and nonsmooth, single- and multidimensional equations.

First of all, we take a numerical example of Eq. (1) with $m=1$, $c1=0.2$, $c2=0.1$, $k=1$, $f(t)=sin\u2009t$, and the initial conditions $x(0)=1$ and $x\u02d9(0)=0$. We select the parameters $L=100$ and $N=1000$, and solve system (15) by the fourth-order Runge-Kutta method with the step length $\Delta t=0.001$. Since it is quite cumbersome to find an exact analytical solution, in order to verify the obtained results, we employ the famous PC algorithm [9,10] to get the reference solution. Nice agreement can be observed between these results, as shown in Fig. 1. Note that the step length of the PC algorithm is also fixed as $\Delta t=0.001$.

Figure 2 shows the comparison of the CPU running time spent by the PC algorithm and the presented method, respectively. The running time for the PC algorithm increases in proportion to $nstep2$ with $nstep$ as the number of integration time steps. An overall linear increasing computational cost can be observed for the presented method.

Now we talk about the choice of parameters $L$ and $N$. Theoretically, $L$ should be a large constant as it is introduced to approximate infinity. As shown in Fig. 3, however, we have to raise $N$ much more as *L* increases to guarantee a reasonable precision for the obtained results. The reason consists in that, when $L$ increases, it requires more truncated terms of the Fourier series to maintain the accuracy. As $L$ is fixed, differently, higher precision can be achieved with $N$ increasing as it denotes the truncated number of the Fourier series. Anyway, in order to achieve accurate results, it is suggested to select a moderate constant for $L$ and a large enough number for $N$. Notice that there would be a trade-off between the precision and the computational efficiency as the dimension of the deduced ODE increases linearly with respect to *N*.

*g*(

*x*), can be given as either smooth or nonsmooth functions in

*x*(

*t*).

Likewise, we choose the parameters $L=100$ and $N=1000$, and solve system (16) with $g(x)=2x3$ with the step length $\Delta t=0.001$. As shown in Fig. 4, nice consistence can still be observed between these results. Consider a nonsmooth case with $g(x)={2x3+x,x<02x3+x+2,x\u22650$ and the same initial conditions. As shown in Fig. 5, the presented method can provide accurate results compared to the PC solution.

As shown in Fig. 6, nice agreement can also be observed, indicating that the presented method can be straightforwardly applied in higher dimensional problems. In fact, there is no substantial difference between the solution procedures for a single- or multidimensional system, as in the presented method one needs to only approximate the FD by the solution of the PDE, i.e., Eq. (3). In our opinion, this is the key point to guarantee the wide applicability of the presented method.

## Conclusions and Remarks

We have presented an efficient method for solving FDEs with fractional order as 3/2, based on the fact that the FD can be considered as the solution to a certain diffusion equation given by PDE. After approximating an infinite boundary condition by a finite one, the fractional differential equation and the PDE are solved by the method of separation of variables. A set of ODEs is then deduced, which are solved by the Runge–Kutta scheme. A variety of numerical examples are taken to validate the presented method, indicating that it is effective for both linear and nonlinear equations, even when nonsmoothness is included.

Aside from the wide applicability, another significant merit of the presented method consists in its high computational efficiency. As the FDEs are transformed into ODEs, the computational cost is at the order of *n* with *n* as the number of integration steps. In fact, the computational efficiency has been enhanced to just the same as solving ODEs. As we know, such high computational efficiency for solving FDEs has been shared by several other approaches such as spectrum based and memory-free algorithms [26,27], and a newly initiated method named as differential operator multiplication method [29,30]. To the best of our knowledge, an outstanding merit of the presented method is its applicability to various different equations by maintaining such high computational efficiency.

Notice that we are enlightened by the physical meaning for the FD of order 3/2. In fact, the presented method can be extended to other problems of order m/2 with m as an odd number, with the help of the index law $D\alpha D\beta f(t)=D\alpha +\beta f(t)$ for FDs [2]. For example, the FD of order 5/2 can be considered as $D5/2f(t)=D3/2f\u2032(t)$, which can possibly be handled as the FD of 3/2 under the framework of the presented method. This will be among our future work. Frankly, it will be a challenging topic for a possible extension of the presented method in solving problems with general fractional derivatives.

## Funding Data

National Natural Science Foundation of China (Grant Nos. 11702333, 11572356, and 11672337; Funder ID: 10.13039/501100001809).