Methods for solving a semi-differential vibration equation

Fractional-order (or, shortly, fractional) derivatives are used in viscoelasticity since the late 1980’s, and they grow more and more popular nowadays. However, their efficient numerical calculation is nontrivial, because, unlike integer-order derivatives, they require evaluation of history integrals in every time step. Several authors tried to overcome this difficulty. In the followings, some of the proposed methods will be examined for a derivative of order 2 (that is sometimes called a semiderivative).


Introduction
It is common to date fractional calculus back to Leibniz, and several textbooks are available on the subject (e.g.[1]).Still, since the definition of fractional derivatives contains a time history integral, their efficient numerical calculation is not solved yet.There are authors trying to optimize these integrals.However, others try to avoid them altogether.
Suarez and Shokooh [4] solved a vibration equation with fractional damping of order 1  2 , taking advantage of the derivative order.Therefore, their solution is restricted to cases where the order of damping is one half.Still, the presented numerical results are very widely used as reference solutions.Moreover, Suarez and Shokooh presented a solution for the initial value problem of the vibration system, which is uncommon in the literature.
Yuan and Agrawal [6] tried a different way.Using the Laguerre integral formula, they have rewritten the definition of the fractional derivative to a simplified formula.Their method has been studied by Trinks and Ruge who claimed extending it to non-zero initial conditions [5].However, the formula in [5] has a singularity at t = 0 for x 0 0 or v 0 0. They have given no example to reveal their interpretation for this.Later, Schmidt and Gaul [3] pointed out that the method of Yuan and Agrawal actually turns the system into a linear viscoelastic one, which does not reproduce the asymptotic behaviour of the fractionally damped system.Moreover, the resulting linear viscoelastic model is less than optimal in reproducing the widely used complex modulus vs. frequency functions.
The present paper aims to compare the above-mentioned methods on a simple example.As reference solution, the steadystate analytical solution and a very basic, direct simulation will be used.In Section 2, after a brief introduction of fractional derivatives, the solution methods are presented.Next, in Section 3, the results are discussed.Finally, Section 4 gives a conclusion of the paper.The differential equation to be solved is the vibration equation with fractional damping, with one degree of freedom: where D = d dt is the differential operator.In the followings, a sinusoidal force will be considered: where ω e is the circular frequency of the excitation.Another common form of Eq. ( 1) is with

Mathematical aspects
For a rational α, D α is usually defined [1] using either the Riemann-Liouville definition (5) In both cases, (n − 1) < α < n, and (x) is the gamma function Actually, the two definitions only differ in the consideration of conditions at the start of the interval: In the literature, it is the Riemann-Liouville definition that seems to be used more often.However, this poses problems, requiring initial conditions for the fractional-order derivatives, which are unavailable at the moment, due to the fact that a proper physical interpretation is still sought [2].Contrarily, the Caputo definition requires only the initial values of the familiar integer-order derivatives.Moreover, in the applications (and, therefore, in Eq. ( 1)), D usually means 0 D t , and, as will be seen, initial conditions for integer-order derivatives are often considered zero.In this case, the difference between Riemann-Liouville (4) and Caputo (5) definitions vanishes.(The other such case is when D is −∞ D t , i.e. when transients are neglected.)

Reference solutions
Steady-state solution The fractional derivative of the sine function, with the lower terminal at −∞, is for λ > 0 and α > −1.From this, the steady-state solution of Eq. ( 1) for the excitation (2) is where and Direct simulation following definition The other reference is a simulation method derived directly from the definition.To maintain simplicity, the explicit Euler formula has been used, with a constant time step h.At t = i h (where n is a positive integer), the fractional derivative has been approximated as (Thus, ẋ is considered constant within the time step, as the average of v j and v j−1 .) The rest is simply calculated using and Here, F i is the force at t = i h.

Suarez and Shokooh: solution in the operator space 2.2.1 The method
In [4], the authors re-write Eq. (1) to a system of equations so that z x(t) and z 4 (t) = x(t), and solve it in the Laplace operator space, for some elementary cases.These solutions are presented in the next paragraphs.
An interesting point is that these authors use the Riemann-Liouville definition of a fractional derivative.

78
Attila Pálfalvi Free vibration with initial conditions The solution for F(t) = 0 and initial conditions x 0 and v 0 is Here, where is the error function.λ j and j are the solution of the eigenproblem with , and 4 j is the 4th coordinate of the jth eigenvector.R j comes from the initial conditions, by Step function response Another elementary case studied by Suarez and Shokooh is the response to a force described by the Heaviside function u(t): Here, x 0 and v 0 are considered zero.The solution is then General case For an arbitrary function F(t) and arbitrary initial conditions, the proposal of Suarez and Shokooh is to consider the initial conditions according to the above.As for the force, they suggest discretizing the function and handling it as a series of step loads.However, with this method, the time history integral of the fractional derivative is replaced by a sum of displacements from step loads, as Eq. ( 22) shows.

Implementation
In the present paper, this method has been implemented as follows.The force F i is calculated for every time step.Next, its change is calculated, also in every time step.At t = 0, F 0 = F 0 .Then, where x u (t) is the solution (21) for the step force with f 0 = 1, i is the time step number and h is the (constant) time step size.
2. 3 Yuan and Agrawal: Laguerre integral 2.3.1 The method The idea of Yuan and Agrawal is to rewrite the Caputo definition (5) to an infinite integral, using the definition (6) of the gamma function, and use the Laguerre integration formula to calculate it.This leads to a system of ordinary linear differential equations: with the unknown functions x(t), v(t) and i (t).The functions i (t) are a certain sort of internal variables, with no clear physical meaning attributed to them.
The great advantage of this approach is that Eq. ( 23) only contains current values of position, velocity and n internal variables i , instead of the time history of the positon or velocity.As Eq. ( 23) is a linear differential equation, it can be solved easily by classical methods.However, this advantage is also a drawback, as described by Schmidt and Gaul [3] and mentioned in Section 1.

Implementation
The implementation used in the present paper is very simple.The system of Eqs. ( 23) is written in the matrix form where and are used for the ith time step (where h is the time step).Although this method is excessively simple, and there are much better ones available, it will be shown in Section 3 that most of the error comes from elsewhere.
Solving a semi-differential vibration equation 79 2007 51 2 Eq. ( 24) has also a simple steady-state solution, that has been determined using an excitation from Eq. ( 2) and the well-known trial function x(t) = x 0s sin(ω e t) + x 0c cos(ω e t) .

Results with the different methods
In the previous section, two methods have been presented for solving a vibration equation with a fractional damping of order 1 2 .In the next part of the paper, a harmonic vibration will be examined with some of them.
Numerical parameters in Eqs. ( 1) and (2) are: m = 1, ω n = 10, η = 0.5, F 0 = 1 and ω e = 4π.Simulation time is 10 s, and time steps are constant in all cases.For the Yuan-Agrawal method, 5 Laguerre points have been used.All methods have been implemented in the Maple computer algebra system.
for the ith time step (where h is the time step).Although this method is excessively simple, and e much better ones available, it will be shown in Section ?? that most of the error comes from re. ation (??) has also a simple steady-state solution, that has been determined using an excitation from n (??) and the well-known trial function x(t) = x 0s sin(ω e t) + x 0c cos(ω e t) .
esults with the different methods revious section, two methods have been presented for solving a vibration equation with a fractional g of order 1 2 .In the next part of the paper, a harmonic vibration will be examined with some of erical parameters in Equations (??) and (??) are: m = 1, ω n = 10, η = 0.5, F 0 = 1 and ω e = 4π.ion time is 10 s, and time steps are constant in all cases.For the Yuan-Agrawal method, 5 Laguerre ave been used.All methods have been implemented in the Maple computer algebra system.It is immediate to see that the direct method fits well to lytical steady-state solution.For the other methods, some problems are seen.Results are shown in Figs. 1 through 4. It is immediate to see that the direct method fits well to the analytical steady-state solution.For the other methods, some problems are seen.
Fig. 2 shows that the result of the method of Suarez and Shokooh fits well to the analytical solution, and the transient Figure ?? shows that the result of the method of Suarez and Shokooh fits well to the analytic and the transient part is practically the same as for the direct method.However, near the simulated period of time, a decrease of amplitude is observed.Its reason is unknown; the suspe is a numerical error in the complex error function.
As for the method of Yuan and Agrawal, there is a large error in the steady-state amplitud ?? and ??).Although it decreases when using a smaller time step, it does not disappear com mentioned in Section ??, the method substitutes the fractionally damped system by a linearly da Therefore, the amplitude and phase lag of the steady-state response can be calculated directly (S Results for different numbers of Laguerre points are shown in Table ??.It is immediate to se 5 Laguerre points, the further decrease of the time step will barely improve the solution.Mor with 20 Laguerre points, the error of the amplitude is above 4%.(In addition, with an increasing Laguerre points, the matrix B of Equation (??) becomes worse-conditioned, causing a numerical and requiring a smaller time step (or a better time integration method).) 5  It is also worth mentioning computational needs.Table ?? shows the elapsed CPU times for t methods (on the same system).Although none of the methods was highly optimized, trivial iss the pre-calculation of products which are always the same have been done (and included in the CPU time).It is immediate to see that the Yuan-Agrawal method was very fast (and can be faster by e.g.modal reduction of the system).However, curiously, the method of Suarez and Sh been beaten by our simplist direct simulation method (although for the former, the necessary va error function have been calculated before the main loop).

Conclusion
In this paper, some methods for the solution of fractional differential equations of derivative or been examined.The followings have been seen.
• The method of Yuan and Agrawal, based on Laguerre integrals, is attractive at a first gla was the fastest of the tested methods.However, it lacks precision in predicting harmonic which means that further results obtained with it should be verified by other means.
• The method of Suarez and Shokooh, which gives a solution of a 1-DOF system with a part is practically the same as for the direct method.However, near the end of the simulated period of time, a decrease of amplitude is observed.Its reason is unknown; the suspected cause is a numerical error in the complex error function.
As for the method of Yuan and Agrawal, there is a large error in the steady-state amplitude (Figs. 3 and 4).Although it decreases when using a smaller time step, it does not disappear completely.As mentioned in Section 2.3, the method substitutes the fractionally damped system by a linearly damped one.Therefore, the amplitude and phase lag of the steady-state response can be calculated directly (Section 2.3.2).
Results for different numbers of Laguerre points are shown in Table 1.It is immediate to see that with 5 Laguerre points, the further decrease of the time step will barely improve the solution.Moreover, even with 20 Laguerre points, the error of the amplitude is above 4%.(In addition, with an increasing number of Laguerre points, the matrix B of Eq. ( 24) becomes worse-conditioned, causing a numerical instability, and requiring a smaller time step (or a better time integration method).)It is also worth mentioning computational needs.Table 2 shows the elapsed CPU times for the different methods (on the same system).Although none of the methods was highly optimized, trivial issues such as the pre-calculation of products which are always the same have been done (and included in the necessary CPU time).It is immediate to see that the Yuan-Agrawal method was very fast (and can be made still faster by e.g.modal reduction of the system).However, curiously, the method of Suarez and Shokooh has been beaten by our simplist direct simulation method (although for the former, the necessary values of the error function have been calculated before the main loop).

Conclusion
In this paper, some methods for the solution of fractional differential equations of derivative order 1 2 have been examined.The followings have been seen.
• The method of Yuan and Agrawal, based on Laguerre integrals, is attractive at a first glance, and it was the fastest of the tested methods.However, it lacks precision in predicting harmonic behaviour, which means that further results obtained with it should be verified by other means.
• The method of Suarez and Shokooh, which gives a solution of a 1-DOF system with a derivative order of 1 2 , is usable.However, since it contains the error function, numerical difficulties have been experienced.Moreover, it was slower than a very simple and similarly accurate direct method.Based on the above, and considering also the fact that 1 2 is not the only derivative order used in material modeling, a direct simulation method can be proposed for the solution of a fractionally damped vibration equation with an arbitrary function.As mentioned in the introduction, further literature is available on the efficient calculation of the time history integral.

Figure 1 :
Figure 1: Result of the direct method (Section ??), 2000 time steps.lts are shown in Figures ?? through ??.It is immediate to see that the direct method fits well to lytical steady-state solution.For the other methods, some problems are seen.

Table 1 :
Harmonic results for the Yuan-Agrawal method.

Table 2 :
Calculation times for the different methods.
Harmonic results for the Yuan-Agrawal method.