Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 53, 4, pp. 959-968, Warsaw 2015 DOI: 10.15632/jtam-pl.53.4.959 NUMERICAL SOLUTION OF NON-HOMOGENOUS FRACTIONAL OSCILLATOR EQUATION IN INTEGRAL FORM Mariusz Ciesielski Czestochowa University of Technology, Institute of Computer and Information Sciences, Częstochowa, Poland e-mail: mariusz.ciesielski@icis.pcz.pl Tomasz Błaszczyk Czestochowa University of Technology, Institute of Mathematics, Częstochowa, Poland e-mail: tomasz.blaszczyk@im.pcz.pl In this paper, a non-homogenous fractional oscillator equation in finite time interval is considered.The fractional equationwithderivatives of orderα ∈ (0,1] is transformed into its corresponding integral form.Next, a numerical solutionof the integral formof the considered equation is presented. In the final part of this paper, some examples of numerical solutions of the considered equation are shown. Keywords: fractional oscillator equation, fractional integral equation, numerical solution 1. Introduction Fractional differential equations have numerous applications in physics, engineering and biology. There is a number of monographs (Atanackovic et al., 2014; Baleanu et al., 2012; Hilfer, 2000; Kilbas et al., 2006; Klimek, 2009; Leszczynski, 2011;Magin, 2006;Malinowska andTorres, 2012; Podlubny, 1999) and a huge number of papers (Agrawal, 2002; Ciesielski and Leszczynski, 2006; Katsikadelis, 2012; Klimek, 2001; Riewe, 1996; Sumelka, 2014; Sumelka and Blaszczyk, 2014; Zhang et al., 2014) that cover various problems in fractional calculus. The list is large and is growing rapidly. An important issue is that the derivative of fractional order at any point of the domain has a local property only when the order is an integer number. For non-integer cases, the fractional derivative is a nonlocal operator and depends on the past values of the function (left derivative) or future ones (right derivative). There are twodifferentways to formulate differential equations containing derivatives of frac- tional order. The first way is simply to replace integer order derivatives in differential equations by fractional ones. The second one relies onmodifying the variational principle by replacing the integer order derivative by a fractional one. This leads to fractional differential equations which are known in the literature as the fractional Euler-Lagrange equations. It involves different types of Lagrangians, e.g., depending on Riemann-Liouville or Caputo fractional derivatives, fractio- nal integrals, and mixed integer-fractional order operators (Agrawal et al., 2011; Almeida and Malinowska, 2012; Baleanu et al., 2014; Klimek et al., 2014; Odzijewicz et al., 2012). This second approach, in recent years, seems to become increasingly important. The main feature of the fractional Euler-Lagrange equations is that these equations invo- lve simultaneously the left and right derivatives of fractional order. This is also a fundamental problem in finding solutions to equations of a variational type (Baleanu et al., 2012; Klimek, 2009). Consequently, numerous studies have been devoted to numerical schemes for the fractio- nal Euler-Lagrange equations (Blaszczyk et al., 2011; Blaszczyk and Ciesielski, 2014, 2015a,b; Bourdin et al., 2013; Pooseh et al., 2013; Xu and Agrawal, 2014). 960 M. Ciesielski, T. Błaszczyk In this paper, we present a numerical solution of the non-homogenous fractional oscillator equation in a finite time interval. 2. Preliminaries Weconsider avariational differential equation containing fractional derivatives of orderα ∈ (0,1] in the finite time interval t ∈ [0,b] CDαb−D α 0+x(t)−ω 2αx(t)= f(t) (2.1) where ω > 0, x(t) is an unknown function and f(t) is a given function. Equation (2.1) is supplemented by the following boundary conditions x(0)= 0 x(b)= L (2.2) In this type of equation (derived from the Lagrangian with fractional derivative and the pro- perties of the considered variational problem (Klimek, 2009) we assume x(0) = 0 on the left boundary. We recall some definitions and properties of the fractional operators (Kilbas et al., 2006; Oldham and Spanier, 1974; Podlubny, 1999). The left and right Riemann-Liouville fractional derivatives of order α ∈ (0,1) are defined by Dα0+x(t) := DI 1−α 0+ x(t) Dαb−x(t) :=−DI 1−α b− x(t) (2.3) and the left and right Caputo derivatives are defined as follows CDα0+x(t) := I 1−α 0+ Dx(t) CDαb−x(t) :=−I 1−α b− Dx(t) (2.4) whereD is an operator of thefirst order derivative and the operators Iα 0+ and Iα b− are respectively the left and right fractional integrals of order α > 0 defined by Iα0+x(t) := 1 Γ(α) t∫ 0 x(τ) (t− τ)1−α dτ (t > 0) Iαb−x(t) := 1 Γ(α) b∫ t x(τ) (τ − t)1−α dτ (t < b) (2.5) If α =1 then D1 0+ x = CD1 0+ x = x′ and D1 b− x = CD1 b− x =−x′. The relations between theCaputo and Riemann-Liouville derivatives are of the form CDα0+x(t)= D α 0+x(t)− t−α Γ(1−α) x(0) CDαb−x(t)= D α b−x(t)− (b− t)−α Γ(1−α) x(b) (2.6) The composition rules of the fractional operators (for α ∈ (0,1]) looks as follows Iα0+ CDα0+x(t)= x(t)−x(0) Iαb− CDαb−x(t)= x(t)−x(b) (2.7) Numerical solution of non-homogenous fractional oscillator equation... 961 and the fractional integral of a constant C Iα0+C = C tα Γ(1+α) (2.8) In particular, when α =1, then CD1 b− D1 0+ =−D2 and Eq. (2.1) becomes −D2f(t)−ω2x(t)= f(t) (2.9) and its analytical solution with respect to boundary conditions (2.2) has the form x(t)= sin(ωt) sin(ωb) [L+F(b)]−F(t) (2.10) where F(t)= 1 ω t∫ 0 f(τ)sin(ω(t− τ)) dτ (2.11) This solution we will apply to verifiction of numerical results obtained by our method for the case α =1. 3. Equivalent integral form The first step of our method is to transform a differential equation into an integral equation. We integrate Eq. (2.1) twice: the first time by using right fractional integral operator (2.5)2 and the second time by using left fractional integral operator (2.5)1. Finally, we obtain Iα0+I α b− CDαb−D α 0+x(t)−ω 2αIα0+I α b−x(t)= I α 0+I α b−f(t) (3.1) Next, we use the composition rule of operators (2.7)2 and get Iα0+ ( Dα0+x(t)−D α 0+x(t) ∣∣∣ t=b ) −ω2αIα0+I α b−x(t)= I α 0+I α b−f(t) (3.2) In the above equation there occurs the value Dα 0+ x(t) ∣∣ t=b , and it is unknown at this stage. Here, we treat it as a constant. By using property (2.7)1 and the assumption x(0) = 0 in the boundary condition, and the fractional integral of constant (2.8) we obtain the following form of the considered equation x(t)−Dα0+x(t) ∣∣∣ t=b tα Γ(α+1) −ω2αIα0+I α b−x(t)= I α 0+I α b−f(t) (3.3) The unknown value Dα 0+ x(t) ∣∣ t=b in the above equation can be determined on the basis of the boundary condition. Substituting the value t = b into Eq. (3.3) x(b)−Dα0+x(t) ∣∣∣ t=b bα Γ(α+1) −ω2αIα0+I α b−x(t) ∣∣∣ t=b = Iα0+I α b−f(t) ∣∣∣ t=b (3.4) we obtain Dα0+x(t) ∣∣∣ t=b = Γ(α+1) bα [ x(b)−ω2αIα0+I α b−x(t) ∣∣∣ t=b − Iα0+I α b−f(t) ∣∣∣ t=b ] (3.5) Next, we put the right-hand side of formula (3.5) into Eq. (3.3) and get the integral form of Eq. (2.1) x(t)−ω2α [ Iα0+I α b−x(t)− (t b )α Iα0+I α b−x(t) ∣∣∣ t=b ] = (t b )α[ L−Iα0+I α b−f(t) ∣∣∣ t=b ] +Iα0+I α b−f(t) (3.6) 962 M. Ciesielski, T. Błaszczyk 4. Numerical solution In this Section,wepresent a numerical scheme forEq. (3.6).We introduce ahomogeneous grid of n+1nodes (with the constant time step ∆t = b/n): 0= t0 < t1 < ... < ti < ti+1 < ... < tn = b, and ti = i∆t, i =0,1, . . . ,n. For every grid node ti, we write the following equation x(ti)−ω2α [ Iα0+I α b−x(t) ∣∣∣ t=ti − ( ti tn )α Iα0+I α b−x(t) ∣∣∣ t=tn ] = ( ti tn )α[ L− Iα0+I α b−f(t) ∣∣∣ t=tn ] + Iα0+I α b−f(t) ∣∣∣ t=ti (4.1) In order to simplify notation, we denote the values of functions x(t) and f(t) at the node ti by xi = x(ti) and fi = f(ti), respectively. The numerical problem is to approximate the values of Iα 0+ Iα b− x(t) ∣∣ t=ti and Iα 0+ Iα b− f(t) ∣∣ t=ti , for i =0,1, . . . ,n, properly. In our previousworks (Blaszczyk and Ciesielski, 2015a,b) we have determined the discrete form of the composition of operators. On the basis of these results, we present the final discrete forms Iα0+I α b−x(t) ∣∣∣ t=ti ≈ i∑ j=0 u (α) i,j n∑ k=j v (α) j,k xk Iα0+I α b−f(t) ∣∣∣ t=ti ≈ i∑ j=0 u (α) i,j n∑ k=j v (α) j,k fk (4.2) where the coefficients u(α)i,j and v (α) j,k are as follows u (α) i,j = (∆t)α Γ(α+2)    0 for i =0 and j =0 (i−1)α+1− iα+1+ iα(α+1) for i > 0 and j =0 (i−j+1)α+1−2(i−j)α+1 +(i−j−1)α+1 for i > 0 and 0 < j < i 1 for i > 0 and j = i v (α) j,k = (∆t)α Γ(α+2)    0 for j = n and k = n (n−j−1)α+1− (n−j)α+1+(n−j)α(α+1) for j < n and k = n (k−j+1)α+1−2(k−j)α+1 +(k−j−1)α+1 for j < n and j < k < n 1 for j < n and k = j (4.3) One can observe that in order to compute the values of operators (4.2) at every node ti, we need to use the values of functions at all nodes of the domain. Now we present the numerical scheme of integral equation (3.6). If we substitute (4.2) into (4.1), then the solution can be written as a system of n+1 linear equations xi −ω2α [ i∑ j=0 u (α) i,j n∑ k=j v (α) j,k xk − ( i n )α n∑ j=0 u (α) n,j n∑ k=j v (α) j,k xk ] = ( i n )α [ L− n∑ j=0 u (α) n,j n∑ k=j v (α) j,k fk ] + i∑ j=0 u (α) i,j n∑ k=j v (α) j,k fk i =0, . . . ,n (4.4) Numerical solution of non-homogenous fractional oscillator equation... 963 The above system can also be written in the matrix form A ·x=b (4.5) with the coefficients of matrices A and b Ai,j = { 1 if i = j 0 if i 6= j −ω 2α [min(i,j)∑ k=0 u (α) i,k v (α) k,j − ( i n )α j∑ k=0 u (α) n,kv (α) k,j ] i =0, . . . ,n j =0, . . . ,n bi = ( i n )α [ L− n∑ j=0 u (α) n,j n∑ k=j v (α) j,k fk ] + i∑ j=0 u (α) i,j n∑ k=j v (α) j,k fk i =0, . . . ,n (4.6) respectively, and x= [x0,x1, . . . ,xn]T. Fig. 1. Numerical solution of Eq. (2.1) for α ∈{0.5,0.6,0.8,1}, b =1, and different values of parameters ω, L and functions f(t); left-side: different ω and f(t)=10; right-side: constant ω =5 and different f(t) 964 M. Ciesielski, T. Błaszczyk Fig. 2. Numerical solution of Eq. (2.1) for α ∈{0.5,0.6,0.8,1}, f(t)= 10sin(2πt), b =1, different values of parameter ω, and L =0 (left-side), L =1 (right-side) 5. Example of computations We present the results of computations obtained by our numerical scheme to the forced fractio- nal oscillator equation. The system of linear equations (4.5) has been solved by the Gaussian elimination algorithm (the LUP decomposition (Press et al., 2007)). We present several exam- Numerical solution of non-homogenous fractional oscillator equation... 965 ples of calculations for different values of the parameters α, ω and different functions f. In the examples we assume: b =1 and L ∈{0,1} in the right boundary condition. In Figs. 1 and 2 we present results of the numerical solution of Eq. (2.1). The values of the parameters used in the solution of equation (2.1) are given in plot legends. The time domain t ∈ [0,1] has been divided into n =1000 subintervals in all examples. 5.1. Numerical error analysis The discretization of differential or integral equations always entails certain errors. The knowledge of an analytical solution xanal(t) allows us to evaluate the errors of our numerical solution xnum(t).We define the average error as Erravg = 1 m+1 m∑ i=0 |xnum(ti)−xanal(ti)| for i =0, . . . ,m (5.1) where m is the number of measurement points. Let us considerEq. (2.1)with theparametersω =0,L =0and the functionf(t)= (b−t)2−α. In this case Iα0+I α b−f(t)= I α 0+I α b−(b− t) 2−α = tαb2Γ(3−α) 2∑ j=0 (−1)j Γ(3− j)Γ(α+1+ j) (t b )j (5.2) and the analytical solution (by using Eq. (3.6) and further simplifications) has form xanal(t)= t αb2Γ(3−α) 2∑ j=0 (−1)j[(t/b)j −1] Γ(3− j)Γ(j +α+1) (5.3) In Tables 1 and 2, we present the values of the analytical (exact) solution xanal(t) for α ∈ {0.5,0.75} at the selected values of nodes ti, i = 0, . . . ,10 and the corresponding nu- merical values xnum(t) calculated for different numbers of grid nodes (n ∈{100,200,400,800}). The errors generated by numerical schema are also included in Tables 1 and 2. One can note that the error decreases with a decrease in the time step of the grid. Table 1. Values of analytical and numerical solutions of Eq. (2.1) for α = 0.5, ω = 0, f(t)= (b− t)2−α, b =1, L =0 ti xanal(ti) xnum(ti) xnum(ti) xnum(ti) xnum(ti) n =100 n =200 n =400 n =800 0 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.1 0.159378794072 0.159379032826 0.159378842522 0.159378803745 0.159378795979 0.2 0.186040855728 0.186041292916 0.186040941860 0.186040872522 0.186040858974 0.3 0.184034779322 0.184035363697 0.184034893461 0.184034801419 0.184034783568 0.4 0.166968260457 0.166968963299 0.166968397209 0.166968286848 0.166968265515 0.5 0.141421356237 0.141422156701 0.141421511713 0.141421386199 0.141421361972 0.6 0.111541920371 0.111542800274 0.111542091205 0.111541953283 0.111541926669 0.7 0.080319362547 0.080320302735 0.080319545231 0.080319397768 0.080319369292 0.8 0.050087922696 0.050088897461 0.050088112561 0.050087959382 0.050087929735 0.9 0.022768399153 0.022769357797 0.022768587015 0.022768435655 0.022768406192 1 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 Erravg – 5.92466 ·10−7 1.15654 ·10−7 2.2385 ·10−8 4.30118 ·10−9 966 M. Ciesielski, T. Błaszczyk Table 2. Values of analytical and numerical solutions of Eq. (2.1) for α = 0.75, ω = 0, f(t)= (b− t)2−α, b =1, L =0 ti xanal(ti) xnum(ti) xnum(ti) xnum(ti) xnum(ti) n=100 n =200 n =400 n =800 0 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.1 0.067645942013 0.067646050586 0.067645965523 0.067645947053 0.067645943086 0.2 0.094996909603 0.094997100277 0.094996949929 0.094996918109 0.094996911394 0.3 0.105395550994 0.105395806812 0.105395604628 0.105395562237 0.105395553351 0.4 0.104362700328 0.104363008356 0.104362764589 0.104362713751 0.104362703135 0.5 0.095196902130 0.095197250210 0.095196974490 0.095196917206 0.095196905277 0.6 0.080331627323 0.080332002341 0.080331705051 0.080331643482 0.080331630691 0.7 0.061751765560 0.061752151496 0.061751845314 0.061751782104 0.061751769002 0.8 0.041170523241 0.041170897406 0.041170600282 0.041170539178 0.041170526550 0.9 0.020119449943 0.020119771392 0.020119515708 0.020119463484 0.020119452744 1 0.000000000000 0.000000000000 0.000000000000 0.000000000000 0.000000000000 Erravg – 2.42522 ·10−7 5.03981 ·10−8 1.04972 ·10−8 2.19045 ·10−9 6. Conclusions In this work, the method of solving of the non-homogenous fractional oscillator equation with derivatives of order α ∈ (0,1] has been presented. First, the equation has been transformed into integral form. Next, the numerical scheme for the integral form of the equation has been proposed.Analysing the solutions of the equation for different values of parameters α, ω, various boundary conditions and different functions f(t) one can understand the influence of these parameters on the solution. One can note that if the value of α decreases then the amplitude of oscillations increases, and if the value of ω increases then the oscillation frequency increases as well. The analytical solution of this type of fractional differential equation (except for the case α =1) isnotyetknown.Ourproposednumericalmethodof solutions forα → 1 is consistentwith the analytical solution for α =1. This type of considered equation can be used inmathematical modelling of the oscillatory behavior of systemswhere the use of the classical oscillator equation is insufficient. References 1. AgrawalO.P., 2002,FormulationofEuler-Lagrangeequations for fractionalvariationalproblems, Journal of Mathematical Analysis and Applications, 272, 368-379 2. Agrawal O.P., Muslih S.I., Baleanu D., 2011, Generalized variational calculus in terms of multi-parameters fractional derivatives,Communications in Nonlinear Science and Numerical Si- mulation, 16, 4756-4767 3. AlmeidaR., Malinowska A.B., 2012,Generalized transversality conditions in fractional calcu- lus of variations,Communications in Nonlinear Science and Numerical Simulation, 18, 3, 443-452 4. Atanackovic T.M., Pilipovic S., Stankovic B., Zorica D., 2014, Fractional Calculus with Applications in Mechanics: Vibrations and Diffusion Processes, ISTE-Wiley, London 5. Baleanu D., Trujillo J.J., 2008, On exact solutions of a class of fractional Euler-Lagrange equations,Nonlinear Dynamics, 52, 331-335 6. BaleanuD., Asad J.H., Petras I., 2014, Fractional Bateman-FeshbachTikochinsky oscillator, Communications in Theoretical Physics, 61, 221-225 7. Baleanu D., Diethelm K., Scalas E., Trujillo J.J., 2012, Fractional Calculus Models and Numerical Methods, World Scientific, Singapore Numerical solution of non-homogenous fractional oscillator equation... 967 8. Blaszczyk T., Ciesielski M., 2014, Numerical solution of fractional Sturm-Liouville equation in integral form, Fractional Calculus and Applied Analysis, 17, 307-320 9. Blaszczyk T., Ciesielski M., 2015a, Fractional oscillator equation: analytical solution and algorithm for its approximate computation, Journal of Vibration and Control (in print), http://dx.doi.org/10.1177/1077546314566836 10. BlaszczykT.,CiesielskiM., 2015b,Fractionaloscillator equation– transformation into integral equation and numerical solution,Applied Mathematics and Computation, 257, 428-435 11. Blaszczyk T., Ciesielski M., Klimek M., Leszczynski J., 2011, Numerical solution of frac- tional oscillator equation,Applied Mathematics and Computation, 218, 2480-2488 12. Bourdin L., Cresson J., Greff I., InizanP., 2013,Variational integrator for fractional Euler- Lagrange equations,Applied Numerical Mathematics, 71, 14-23 13. Ciesielski M., Leszczynski J., 2006, Numerical solutions to boundary value problem for ano- malous diffusion equationwithRiesz-Feller fractional operator, Journal of Theoretical and Applied Mechanics, 44, 2, 393-403 14. Hilfer R., 2000,Applications of Fractional Calculus in Physics, World Scientific, Singapore 15. Katsikadelis J.T., 2012, Nonlinear dynamic analysis of viscoelastic membranes described with fractional differential models, Journal of Theoretical and Applied Mechanics, 50, 3, 743-753 16. Kilbas A.A., Srivastava H.M., Trujillo J.J., 2006, Theory and Applications of Fractional Differential Equations, Elsevier, Amsterdam 17. Klimek M., 2001, Fractional sequential mechanics –models with symmetric fractional derivative, Czechoslovak Journal of Physics, 51, 1348-1354 18. Klimek M., 2009,On Solutions of Linear Fractional Differential Equations of a Variational Type, Publishing Office of CzestochowaUniversity of Technology, Czestochowa 19. Klimek M., Odzijewicz T., Malinowska A.B., 2014, Variational methods for the fractional Sturm-Liouville problem, Journal of Mathematical Analysis and Applications., 416, 402-426 20. Leszczynski J.S., 2011, An Introduction to Fractional Mechanics, Publishing Office of Czesto- chowaUniversity of Technology, Czestochowa 21. Magin R.L., 2006,Fractional Calculus in Bioengineering, Begell House Inc, Redding 22. MalinowskaA.B.,TorresD.F.M., 2012, Introduction to the Fractional Calculus of Variations, Imperial College Press, London 23. Odzijewicz T.,MalinowskaA.B., TorresD.F.M., 2012,Fractional variational calculuswith classical and combinedCaputo derivatives,Nonlinear Analysis, Theory, Methods and Applications, 75, 1507-1515 24. Oldham K.B., Spanier J., 1974,The Fractional Calculus: Theory and Applications of Differen- tiation and Integration to Arbitrary Order, Academic Press, San Diego 25. Podlubny I., 1999,Fractional Differential Equations, Academic Press, San Diego 26. Pooseh S., Almeida R., Torres D.F.M., 2013, Discrete direct methods in the fractional cal- culus of variations,Computational and Applied Mathematics, 66, 668-676 27. PressW.H., Teukolsky S.A., VetterlingW.T., FlanneryB.P., 2007,Numerical Recipes: The Art of Scientific Computing (3rd ed.), Cambridge University Press, NewYork 28. Riewe F., 1996,Nonconservative Lagrangian andHamiltonianmechanics,Physical Review E, 53, 1890-1899 29. Sumelka W., 2014, A note on non-associated Drucker-Prager plastic flow in terms of fractional calculus, Journal of Theoretical and Applied Mechanics, 52, 2, 571-574 30. Sumelka W., Blaszczyk T., 2014, Fractional continua for linear elasticity, Archives of Mecha- nics, 66, 3, 147-172 968 M. Ciesielski, T. Błaszczyk 31. XuY., Agrawal O.P., 2014,Models and numerical solutions of generalized oscillator equations, Journal of Vibration and Acoustics, 136, 051005, 7pp 32. Zhang Y., Chen L., Reeves D.M., Sun H.G., 2014, A fractional-order tempered- stable continuity model to capture surface water runoff, Journal of Vibration and Control, http://dx.doi.org//10.1177/1077546314557554 Manuscript received January 29, 2015; accepted for print May 20, 2015