This manuscript introduces the concept of generalized integrating factor for one dimen- sional linear ordinary differential equations of order n. The procedure is used to address linear second order equations with varying and with constant coefficients, commonly found in many practical problems. The solutions are analytically derived by means of double convolutions. Analytical solutions for the constant coefficient case with different types of continuous and discontinuous excitations are discussed with examples. The con- cept of Heaviside series is introduced to generalize the solutions for discrete excitations. Keywords: ordinary differential equation, second order, analytical solution, integrating factor 1. Introduction Linear and non-homogeneous ordinary differential equations are common in many areas of applied mathematics. In special, second order ordinary differential equations appear in various practical problems in mathematics, physics and engineering, like for example the dynamic equilibrium of mechanical systems. Their solution is well known for constant coefficients and periodic forces [17], but the solution for general problems usually depends on special transformations. Integratingfactorsareusedasatooltomodifyageneralordinarydifferentialequation to an exact form. Once in this form, its solution is just a matter of integration. Its use is well known for first order equations [5] but its use for second order equations is not so well established in the literature [7, 1, 15, 8, 14, 2, 4, 3, 12]. We propose and discuss a simple but powerful variant of the traditional integrating factor that presents many benefits. The proposed method can be used for general orders and results in integrating factors function of the independent variable only, disregarding the coefficients of the ordinary differential equation. ∗Corresponding author Email addresses: matheusj2009@hotmail.com (Matheus Janczkowski Fogaça), eduardo.cardoso@udesc.br (Eduardo Lenz Cardoso ) Preprint submitted to Elsevier March 22, 2023 Thegeneralconceptisshownindetailsandappliedtoathird-orderODEwithvarying coefficients to exemplify the procedure. The method is then used to solve second-order ODEs with varying coefficients and with constant coefficients. The latter case is used as starting point to derive many important analytical solutions for different types of excitations (non-homogeneous term of the ODE). Finally, we discuss an approach to use the proposed method to solve discrete forms of excitation by means of Heaviside series. 2. The generalized integrating factor method Consider a linear first order ordinary differential equation a1(t) ˙y(t) +a0(t)y(t) =f(t), (1) where yis the dependent variable, tis the independent variable and ˙y(t)is the derivative ofywith respect to t. The integrating factor for first order differential equations, introduced by Leibniz, relies on the relation p(t) ˙y(t) + ˙p(t)y(t) =. (p(t)y(t)), (2) where. ( )means the derivative, with respect to t, of all the expression inside the paren- thesis. An integrating factor µ(t)is multiplied to Eq. 1 to force the appearance of Eq. 2, µ(t)a1(t)|{z} p(t)˙y(t) +µ(t)a0(t)|{z} ˙p(t)y(t) =. (p(t)y(t)). (3) It immediately follows that ˙p(t) =µ(t)a0(t) =. µ(t)a1(t) = ˙µ(t)a1(t) +µ(t)˙a1(t) (4) such that ˙µ(t) µ(t)=a0(t)−˙a1(t) a1(t), (5) with solution µ(t) = expZa0(t)−˙a1(t) a1(t)dt , (6) allowing the use direct integration to solve . (µ(t)a1(t)y(t)) =µ(t)f(t), (7) such that y(t) =1 µ(t)a1(t)Z µ(t)f(t) dt+C1, (8) where C1is an integration constant. Relation given by Eq. 2 can be generalized to higher orders pj(t)(j)y(t) + ˙pj(t)(j−1)y(t) = (˙pj(t)(j−1)y(t)), (9) 2 where(j)y(t)is the j-th derivative of ywith respect to t(used for j >3). Now consider a linear and non-homogeneous ordinary differential equation of order n an(t)(n)y(t) +an−1(t)(n−1)y(t) +···+a2(t)¨y(t) +a1(t) ˙y(t) +a0(t)y(t) =f(t).(10) It can be rewritten in pairs of derivatives of y(t)such that one element is a derivative higher than the other. This is achieved by partitioning the coefficients multiplying the intermediate derivatives into two terms an(t)(n)y(t) +fn,n−1,1(t)(n−1)y(t)| {z } πn,n+fn,n−1,2(t)(n−1)y(t) +fn,n−2,1(t)(n−2)y(t)| {z } πn,n−1+. . . +fn,3,2(t)...y(t) +fn,2,1(t)¨y(t)| {z } πn,3+fn,2,2(t)¨y(t) +fn,1,1(t) ˙y(t)| {z } πn,2+ fn,1,2(t) ˙y(t) +a0(t)y(t)| {z } πn,1=f(t),(11) where πn,jisthe j-thpartition ofEq. 10, whichisthekeyideafortheproposedapproach. Coefficients fn,j,irefer to the order of the differential equation, n, partition jandi= 1 ori= 2such that an,j(t) =fn,j,1(t) +fn,j,2(t), (12) and there are Nπn= 2n−2partitions j. Multiplying Eq. 11 by a generalized integrating factor µn(t)results in µn(t)an(t)(n)y(t) +µn(t)fn,n−1,1(t)(n−1)y(t)| {z }.  pn,n(n−1)y(t)+µn(t)fn,n−1,2(t)(n−1)y(t) +µn(t)fn,n−2,1(t)(n−2)y(t)| {z }.  pn,n−1(n−2)y(t)+. . . +µn(t)fn,3,2(t)...y(t) +µn(t)fn,2,1(t)¨y(t)| {z }. (pn,3..y(t))+µn(t)fn,2,2(t)¨y(t) +µn(t)fn,1,1(t) ˙y(t)| {z }. (pn,2.y(t))+ µn(t)fn,1,2(t) ˙y(t) +µn(t)a0(t)y(t)| {z }. (pn,1y(t))=µn(t)f(t), (13) such that . pn,n(n−1)y(t) +.  pn,n−1(n−2)y(t) +. . .+. pn,1y(t) =µn(t)f(t) (14) can be exactly integrated to pn,n(n−1)y(t) +pn,n−1(n−2)y(t) +. . .+pn,1y(t) =Z µn(t)f(t) dt+Cn (15) an ordinary differential equation of order n−1, where Cnis an integration constant. 3 Coefficients fn,j,ican be found by solving ˙µn µn n=˙µn µn n−1=···=˙µn µn 3=˙µn µn 2=˙µn µn 1. (16) where πn,j=⇒˙µn(t) µn(t) j=fn,j−1,1(t)−˙fn,j,2(t) fn,j,2(t);fn,0,1(t) =a0(t), fn,n,2(t) =an(t).(17) It follows that there aren(n−1) 2combinations of these pair-wise equations, along with the n−1equations relating each pair of partitions with their coefficient. Thus, the number of equations available to evaluate coefficients fn,j,i,Neq, is Neqn=n(n−1) 2+n−1 =n2+n−2 2≥Nπn;n≥2, (18) being always bigger or equal than the number of partitions, Nπn, such that there are enough equations to solve the problem (actually, not all equations must be used). After finding all coefficients fn,j,i, it is possible to find µn(t)by using Eq. 17 for just one partition ˆj(any partition can be used) ˙µn(t) µn(t) ˆj=fn,ˆj−1,1(t)−˙fn,ˆj,2(t) fn,ˆj,2(t), (19) such that µn(t) = exp Zfn,ˆj−1,1(t)−˙fn,ˆj,2(t) fn,ˆj,2(t)dt . (20) The same procedure depicted above can be carried out successively until reaching a first order equation, where the traditional integrating factor µ1(t)can be used. The procedure is applied to a linear ordinary differential equation of order three to exemplify the procedure. 2.1. Example: linear ordinary differential equation of order three Consider the ODE a3(t)...y(t) +a2(t)¨y(t) +a1(t) ˙y(t) +a0(t)y(t) =f(t). (21) The explicit dependency on tis dropped in the following equations to make the notation more concise. Partitioning coefficients a2anda1intoa2=f2,2,1+f2,2,2anda1= f2,1,1+f2,1,2results in a3...y+f2,2,1¨y|{z} π3,3+f2,2,2¨y+f2,1,1˙y|{z } π3,2+f2,1,2˙y+a0y|{z} π3,1=f. (22) Multiplying Eq. 23 by a generalized integrating factor µ3 µ3a3|{z} p3,3...y+µ3f2,2,1|{z} ˙p3,3¨y+µ3f2,2,2|{z} p3,2¨y+µ3f2,1,1|{z} ˙p3,2˙y+µ3f2,1,2|{z} p3,1˙y+µ3a0|{z} ˙p3,1y=µ3f,(23) 4 such that . p3,3¨y +. p3,2˙y +. p3,1y =µ3f, (24) can be integrated with respect to t, resulting in p3,3¨y+p3,2˙y+p3,1y=Z µ3fdt+C3 |{z } h2(25) a second order ordinary differential equation. The same procedure can be applied to this second order equation. Splitting p3,2= f1,1,1+f1,1,2results in p3,2¨y+f1,1,1˙y|{z } π2,2+f1,1,2˙y+p3,1y|{z } π2,1=h2 (26) and multiplying by a generalized integrating factor µ2 µ2p3,2|{z} p2,2¨y+µ2f1,1,1|{z} ˙p2,2˙y+µ2f1,1,2|{z} p2,1˙y+µ2p3,1|{z} ˙p2,1y=µ2h2, (27) such that . p2,2˙y +. p2,1y =µ2h (28) and integrating w.r.t time results in p2,2˙y+p2,1y=Z µ2hdt+C2 |{z } h1, (29) a first order ODE. Finally, multiplying by an integrating factor µ1results in µ1p2,2|{z} p1,1˙y+µ1p2,1|{z} ˙p1,1y=µ1h1 (30) or . p1,1y =µ1h1 (31) such that y=1 p1,1Z µ1h1dt+C1, (32) is the solution to Eq. 23. Solution provided by Eq. 32 depends on µ3,µ2andµ1. Using Eq. 30, it is possible to write ˙p1,1=µ1p2,1=. µ1p2,2 = ˙µ1p2,2+µ1˙p2,2 (33) such that˙µ1 µ1=p2,1−˙p2,2 p2,2(34) 5 with solution µ1= exp Zp2,1−˙p2,2 p2,2dt! . (35) The same procedure can be applied to Eq. 27. There are two possible combinations: ˙p2,2=µ2f1,1,1=. µ2p3,2 = ˙µ2p3,2+µ2˙p3,2 (36) and ˙p2,1=µ2p3,1=. µ2f1,1,2 = ˙µ2f1,1,2+µ2˙f1,1,2, (37) such that ˙µ2 µ2=f1,1,1−˙p3,2 p3,2=p3,1−˙f1,1,2 f1,1,2. (38) Since p3,2=f1,1,1+f1,1,2it is possible to rewrite previous equation as f1,1,2(p3,2−f1,1,2)−f1,1,2˙p3,2=p3,2p3,1−p3,2˙f1,1,2 (39) or f2 1,1,2=−p3,2p3,1+p3,2˙f1,1,2+ (p3,2−˙p3,2)f1,1,2 (40) a Riccati differential Equation. Finally, there are three possible combinations for Eq. 23: ˙p3,3=µ3f2,2,1=. (µ3a3)= ˙µ3a3+µ3˙a3, (41) ˙p3,2=µ3f2,1,1=. µ3f2,2,2 = ˙µ3f2,2,2+µ3˙f2,2,2, (42) and ˙p3,1=µ3a0=. µ3f2,1,2 = ˙µ3f2,1,2+µ3˙f2,1,2. (43) Thus, ˙µ3 µ3=f2,2,1−˙a3 a3=f2,1,1−˙f2,2,2 f2,2,2=a0−˙f2,1,2 f2,1,2. (44) Using the splits a1=f2,1,1+f2,1,2anda2=f2,2,1+f2,2,2 (a2−f2,2,2)−˙a3 a3=(a1−f2,1,2)−˙f2,2,2 f2,2,2=a0−˙f2,1,2 f2,1,2, (45) such that f2 2,2,2=a3˙f2,2,2+ (a2−a3)f2,2,2−a3a1+a3f2,1,2 (46) and f2 2,1,2=f2,2,2˙f2,1,2+ (a1−˙f2,2,2)f2,1,2−a0f2,2,2 (47) form a coupled system of Riccati-like differential equations. Previous equations are much simpler when the coefficients are constant, since most of the rates are null. In this case, Eq. 35 simplifies to µ1= exp Zp2,1 p2,2dt! , (48) 6 the Riccati equation, Eq. 40 turns into a second order algebraic equation f2 1,1,2=−p3,2p3,1+p3,2f1,1,2 (49) and the system of Riccati-like differential equations, Eqs. 46 and 47 reduces to f2 2,2,2= (a2−a3)f2,2,2−a3a1+a3f2,1,2 (50) and f2 2,1,2=a1f2,1,2−a0f2,2,2, (51) a coupled system of quadratic algebraic equations. Thus, the main shortcoming of the proposed procedure is the evaluation of Equations 16. The larger the order of the differential equation, the harder is to find the coefficients fn,j,iat each step of the procedure to decrease the order. For general coefficients aj(t), Eqs. 16 result in a system of Riccati-like coupled differential equations and for constant coefficients a system of coupled quadratic equations. Nonetheless, for linear second order equations the procedure can lead to very inter- esting and practical results, as it will be discussed in the rest of this manuscript. 3. Linear second order ordinary differential equations Consider a linear second order ordinary differential equation (ODE) m(t)¨y(t) +c(t) ˙y(t) +k(t)y(t) =f(t), (52) with initial conditions ˙y(0) = v0andy(0) = u0, where y(t)andf(t)are functions of toverR. This equation is of great interest in physical problems, like for example the vibration of a mass-spring-damper subjected to a force. Thus, independent variable tis also referred to time in the rest of this manuscript. For physical reasons, it is assumed that both m(t)andk(t)∈R>0and that c(t)∈R+, ∀t. The explicit dependency on twill be suppressed in the following equations to simplify the notation. We start by splitting cas π2,1=⇒c=f2,1,1+f2,1,2, (53) where f2,1,1andf2,1,2are also function of time, but over C, such that m¨y+f2,1,1˙y|{z} π2,2+f2,1,2˙y+ky|{z} π2,1=f. (54) Multiplying the ODE by integrating factor µ2(t)results in µ2m|{z} p2,2¨y+µ2f2,1,1|{z} ˙p2,2˙y+µ2f2,1,2|{z} p2,1˙y+µ2k|{z} ˙p2,1y=µ2f, (55) or p2,2¨y+ ˙p2,2˙y+p2,1˙y+ ˙p2,1y=µ2f. (56) 7 Analysing the term ˙p2,2=µ2f2,1,1=. (µ2m) = ˙µ2m+µ2˙m, (57) it is possible to state that ˙µ2 µ2=f2,1,1−˙m m. (58) Following the same procedure, ˙p2,1=µ2k=. µ2f2,1,2 = ˙µ2f2,1,2+µ2˙f2,1,2 (59) such that˙µ2 µ2=k−f2,1,2 f2,1,2. (60) Thus, by relating Eqs. 58 and 60, we obtain the particular form of Eq. 16 f2,1,1−˙m m=k−f2,1,2 f2,1,2. (61) Since c=f2,1,1+f2,1,2it is possible to state that f2,1,1=c−f2,1,2. Using Eq. 61  c−f2,1,2−˙m f2,1,2= k−˙f2,1,2 m (62) such that f2 2,1,2=˙f2,1,2m−km+ (c−˙m)f2,1,2, (63) is a Riccati differential equation1. Equation 60 is a first order ODE with known solution µ2= expZk−˙f2,1,2 f2,1,2dt (66) where f2,1,2is obtained from Eq. 63. It is important to stress that only the particular solution of the Riccati equation is needed. By knowing µ2it is possible to re-write Eq. 56 as . p2,2˙y +. p2,1y =µ2f (67) such that integrating with respect to tresults in (p2,2˙y) + (p2,1y) =Z µ2fdt+C2 |{z } h(68) 1Conversely, it is also possible to define f2,1,2=c−f2,1,1such that f2,1,1−˙m m=k−˙c+˙f2,1,1 c−f2,1,1(64) or f2 2,1,1=−˙f2,1,1m+ ( ˙m+c)f2,1,1−˙mc+ (˙c−k)m, (65) also a Riccati differential equation. 8 a first order ODE. Using another integrating factor µ1such that µ1p2,2|{z} p1,1˙y+µ1p2,1|{z} ˙p1,1y=µ1h (69) or p1,1˙y+ ˙p1,1y=µ1h. (70) Following the same procedure ˙p1,1=µ1p2,1=. µ1p2,2 = ˙µ1p2,2+µ1˙p2,2, (71) such that˙µ1 µ1=p2,1−˙p2,2 p2,2, (72) with known solution µ1= expZp2,1−˙p2,2 p2,2dt. (73) Equation 70 can be written as. p1,1y =µ1h (74) integrating with respect to time results in p1,1y=Z µ1hdt+C1 (75) such that y=1 p1,1Z µ1hdt+1 p1,1C1. (76) Thus, by using the definition of both p1,1andh y(t) =1 µ1(t)µ2(t)m(t) Zt 0µ1(t) Zt 0µ2(t)f(t) dt! dt+Zt 0µ1(t)C2dt+C1 (77) is the general solution for the second order ordinary differential equation stated in Eq. 52. This general solution can be split into its particular, yp(t), and homogeneous, yh(t) parts yp(t) =1 µ1(t)µ2(t)m(t)Zt 0µ1(t)Zt 0µ2(t)f(t) dtdt (78) and yh(t) =1 µ1(t)µ2(t)m(t) Zt 0µ1(t)C2dt+C1! , (79) such that y(t) =yp(t) +yh(t). Constants C1andC2can be found by considering the solution at known t/values. Thus, the proposed solution depends on the solution of a Riccati equation and always results in integrating factors µ1(t)andµ2(2)function of tonly, disregarding the form of the coefficients m(t),c(t)andk(t). 9 Solution of Eq. 63 is fundamental forthe success of the proposedformulation. Indeed, the solution of the Riccati equation is not an easy task if we consider general coefficients m(t),c(t)andk(t). Nonetheless, only the particular solution is needed. Various analyt- ical solutions can be found in the literature for specific forms of coefficients [6, 16, 10], as well as numerical methods [9]. Example I - Cauchy-Euler equation Consider the second order ODE t2¨y(t)−2t˙y(t) + 2y(t) = 6 t2+ 4 ln( t). (80) The coefficients are m(t) =t2,c(t) =−2tandk(t) = 2and the excitation is f(t) = 6t2+ 4 ln( t). The Riccati equation, Eq. 63, can be written as f2 2,1,2=t2˙f2,1,2−2t2−4tf2,1,2, (81) with particular solution f2,1,2=−t. Using Eq. 66 µ2= expZ2 + 1 −tdt=1 t3. (82) Thus, p2,2=µ2m=t−1andp2,1=µ2f1,2=−t−2. Using Eq. 73 µ1= expZ−t−1+t−2 t−1dt= 1. (83) Complete solution is given by Eq. 77 y(t) =1 1t−3t2 Zt 01(Zt 0t−3 6t2+ 4 ln( t) dt) dt+Zt 01C2dt+C1 (84) such that y(t) = (6 t2+ 2) ln( t)−6t2+ 3 + C2t2+C1t, (85) the analytical solution. Example II - Bessel equation Consider EDO t2¨y(t) +t˙y(t) + t2−1 4 y(t) =f(t) =t3 2, (86) with known values y(t0) = ˙y(t0) = 0att0= 0.1. It is known that solution to this equation is singular at t= 0. The corresponding Riccati equation, Eq. 63, is f2 2,1,2=t2˙f2,1,2−tf2,1,2+t2 4−t4, (87) 10 whose candidate solution is a second order polynomial ˜f2,1,2=z0+z1t+z2t2. (88) Applying this polynomial into Eq. 87 yields z2 0+2z0z1t+ 2z0z2+z2 1 t2+2z1z2t3+z2 2t4=z1t2+2z2t3−z0t−z1t2−z2t3+t2 4−t4,(89) that simplifies to z2 0+ 2z0z1t+ 2z0z2+z2 1 t2+ 2z1z2t3+z2 2t4=−z0t+t2 4+z2t3−t4,(90) whose solution is z0= 0,z1=1 2andz2=i. Thus, the integrating factor, Eq. 66, is µ2= exp Zt2−1 4−1 2−2it t 2+it2dt! = exp Zt2−2it−3 4 t 2+it2dt! . (91) The polynomials in the integrand can be factored and simplified to µ2= exp Z t−i 2 t−3i 2 it t−i 2 dt = expZ −idt−3 2Z1 tdt ,(92) such that µ2= exp −it−3 2ln|t| =t−3 2exp (−it). (93) The integrating factor to integrate the differential equation from first order to an algebraic equation is evaluated using Eq. 73, µ1=t1 2exp (it). (94) Solution is then given by Eq. 77 y(t) =t−1 2exp (−it)Z exp (2 it)Z t−3 2exp (−t)f(t) dtdt+ C1t−1 2exp (it) +C2t−1 2exp (−it), (95) and replacing the expression for f(t) y(t) =t−1 2+C1t−1 2exp (it) +C2t−1 2exp (−it). (96) Constants, C1andC2, can be obtained by solving a system of linear equations 11  t−1 2 0exp (it0) t−1 2 0exp (−it0) it−1 2 0−t−3 2 0 2! exp (it0)− it−1 2 0+t−3 2 0 2! exp (−it0) C1 C2 =  y(t0)−t−1 2 0 ˙y(t0) +t−3 2 0 2  .(97) Using values at t0= 0.1results in y(t) =t−1 2+ (−0.4975 + 0 .04992 i)t−1 2exp (it)−(0.4975 + 0 .04992 i)t−1 2exp (−it).(98) Figure 1 compares the real part of the solution obtained by the proposed approach, Eq. 95, with the solution obtained by using the Tsitouras 5/4Runge-Kutta method [18] with adaptative time step. The numerical solution starts at t= 0.1(dotted line) but the analytical solution captures the singularity at t→0(blue curve). Figure 1: Analytical response y(t)obtained for the Bessel equation (Eq. 86) and the solution obtained by using a numerical method, yTsit 5(t). 3.1. Constant coefficients Previous equations are much simpler when m,candkare constant since ˙m,˙c,˙f2,1,2 and˙kare zero. As m > 0, it is possible to normalize the equation by msuch that ¨y(t) + ¯c˙y(t) +¯ky(t) =¯f(t), (99) 12 where ¯c=c/m,¯k=k/mand ¯f=f/m. It follows that Eq. 63 reduces to an algebraic quadratic equation f2 2,1,2=−¯k+ ¯cf2,1,2 (100) withdirect solution f2,1,2=¯c±p ¯c2−4¯k 2, (101) a complex number when ¯c2−4¯k < 0(under damped problems). Any one of the two roots can be used. Equation 66 reduces to µ2= exp ˆkt (102) where ˆk=¯k f2,1,2(103) and Eq. 73 to µ1=exp f2,1,2t exp ˆkt= exp f2,1,2−ˆk t . (104) Thus, the complete solution is y(t) =1 µ1(t)µ2(t) Zt 0exp f2,1,2−ˆk t Zt 0exp ˆkt ¯f(t) dt+C2! dt+C1  (105) where the term1 µ1(t)µ2(t)= exp −f2,1,2t (106) such that y(t) = exp −f2,1,2t Zt 0exp f2,1,2−ˆk t Zt 0exp ˆkt ¯f(t) dt+C2! dt+C1 , (107) is the complete solution. Additionally, it is possible to split the solution in its particular yp(t) = exp −f2,1,2t Zt 0exp f2,1,2−ˆk t Zt 0exp ˆkt ¯f(t) dt! dt ,(108) and homogeneous parts yh(t) = exp −f2,1,2t Zt 0exp f2,1,2−ˆk t C2dt+C1! , (109) which can be further simplified to 13 yh(t) = exp −f2,1,2t C1+C2exp −ˆkt . (110) Particular solution for the constant coefficient case, Eq. 108, can be further devel- oped if we assume some particular form for ¯f(t). The following subsections are devoted to investigate two situations: Continuous excitation (periodic and polynomial) and dis- continuous (unitary impulse and Heaviside). 3.2. Continuous excitation functions 3.2.1. Periodic excitation Lets consider a further hypothesis: excitation f(t)is periodic or its Fourier series is convergent. Thus, normalized excitation ¯f(t)can be represented as a series of exponen- tials ¯f(t) =nkX k=1ckexp (βkt+ϕk), (111) where nkis the number of terms, ck∈Ris an amplitude, βk=iωk∈Cis a complex angular frequency and ϕk∈Cis a complex phase. Previously, it was seen that the integrating factor is an exponential by definition, and that the particular solution, i.e., the solution corresponding to the excitation function, naturally appears through the successive integrations. These integrations are convolu- tions over the integrating factor, thus, if the excitation can be expressed in terms of exponentials, these convolutions can be trivially calculated. Applying Eq. 111 into Eq. 108 yields yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk t Zt 0exp ˆktnkX k=1ckexp (βkt+ϕk) dt dt, (112) using the multiplication property between exponentials yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk t Zt 0nkX k=1ckexp βk+ˆk t+ϕk dt dt, (113) whose inner integral is trivial and results in yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk tnkX k=1ck βk+ˆkexp βk+ˆk t+ϕk dt. (114) Again, rearranging the multiplication of exponentials, yields yp(t) = exp −f2,1,2tZt 0nkX k=1ck βk+ˆkexp βk+f2,1,2 t+ϕk dt,(115) 14 which is again trivially integrated to yp(t) = exp −f2,1,2tnkX k=1ck βk+ˆk βk+f2,1,2exp βk+f2,1,2 t+ϕk ,(116) such that yp(t) =nkX k=1ckexp (βkt+ϕk) βk+ˆk βk+f2,1,2. (117) Equation 117 gives a closed-form and analytic particular solution for any periodic excitation or an excitation represented through its Fourier series expansion. It is worth noticing that all constants given by the two integrations at time t0= 0were omitted for they can be coupled to the constants C1andC2, which are by definition constants of these very integrations. Example with periodic excitation Consider the ODE ¨y(t) + 2 ˙y(t) + 10 y(t) =¯f(t), (118) with initial conditions y(0) = 0 .2and˙y(0) = 0 .0, where ¯f(t) =−cos(0.5t) +sin(t) +cos(1.5t−1.5)−2sin(2t) + 2sin(10t).(119) For this example, one gets ¯k= 10, (120) f2,1,2= 1 + 3 i, (121) ˆk= 1−3i, (122) By using well known relations sin ωηt+ϕη =ie−i(ωηt+ϕη)−ei(ωηt+ϕη) 2(123) and cos ωηt+ϕη =e−i(ωηt+ϕη)+ei(ωηt+ϕη) 2, (124) one gets the excitation in a complex form as ¯f=−1 2e0.5it−1 2e−0.5it−i 2eit+i 2e−it+1 2e1.5i(t−1)+1 2e1.5i(1−t)+ ie2it−ie−2it−ie10it+ie−10it=10X j=1cjeβjt+ϕj.(125) 15 Looking at Eq. 125 and the series explicit in the end, an algorithmic approach may be derived by writing these coefficients and powers in a tabular manner, as in Tab. 1. Using such matrix representation, all the parameters needed for each operation of the sum can be summarized in each column. Thus, the summation can be carried out column by column, as the row with the index jin Tab. 1 hints. Table 1: Tabular representation for the excitation. j 1 2 3 4 5 6 7 8 9 10 c−0.5−0.5−0.5i0.5i0.5 0 .5 i−i−i i β0.5i−0.5i i −i 1.5i−1.5i2i−2i10i−10i ϕ 0.0 0 .0 0 .0 0 .0−1.5i1.5i0.0 0 .0 0 .0 0 .0 Following Tab. 1 to evaluate Eq. 117 yields yp=−78 + 8 i 1537e0.5it−78 + 8 i 1537e−0.5it−2 + 9 i 170eit+−2 + 9 i 170eit+ 62−24i 1105e1.5i(t−1)+62 + 24 i 1105e1.5i(1−t)+2 + 3 i 26e2it+ 2−3i 26e−2it+−2 + 9 i 850e10it−2 + 9 i 850e−10it,(126) and for Eq. 110 yh=C2e−¯k f1,2t+C1e−f1,2t, (127) or yh=C2e(−1+3i)t+C1e−(1+3i)t. (128) The homogeneous solution, yh, is straightforward. Using Eq. 126 and Eq. 128 at t= 0, the integration constants came out as complex conjugates, C2≈0.10564−0.11980 i andC1≈0.10564 + 0 .11980 i. Figure 4 shows the real part of the solution for both the homogeneous and the partic- ular solutions (top) as well as the complete solution (bottom). The solution was obtained for a time span of 40swith intervals of 0.1s. It is worth mention that although y(t)is complex, the maximum amplitude of the complex part of the response was 8.3267×10−17 in the analysed time interval, which is zero when compared to the real part. Since it is known that the response to this problem is real it is worth showing that this is indeed just a matter of representation. To this end, Eqs. 126 and 128 can be transformed to a real-valued form by applying Eq. 123 and Eq. 124 in reverse, yp=−156 1537cos(0.5t)−16 1537sin(0.5t)−2 85cos(t) +9 85sin(t) + 124 1105cos 1.5 (t−1) +48 1105sin 1.5 (t−1) +2 13cos(2t)− 3 13sin(2t)−2 425cos(10t)−9 425sin(10t),(129) 16 Figure 2: Particular and homogeneous solution (top) and complete solution (bottom). and yh=C2e−tcos(3t) +C1e−tsin(3t), (130) with C2≈0.2107andC1≈0.0702. Complete solution y(t) =yp(t) +yh(t)is correct and satisfies Eq. 118 for all t. 3.2.2. Polynomial excitation Assuming a polynomial excitation in the form ¯f(t) =npX k=0ck(t−ts)k, (131) where ck∈Rare coefficients, npis the number of terms and ts∈Ra time shift. Applying Eq. 131 into Eq. 108 yields yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk t Zt 0exp ˆktnpX k=0ck(t−ts)kdt dt, (132) as integration is a linear operator, the summation can be transferred to the whole inte- gration and the coefficient of each power of tcan also be put out of the integral, 17 yp(t) = exp −f2,1,2tnpX k=0ckZt 0exp f2,1,2−ˆk t Zt 0exp ˆkt (t−ts)kdt! dt. (133) The convolution of a polynomial over an exponential is recursively evaluated using integrationbyparts. Forit, letthepowerof tbeapositiveinteger α. Thefirstintegration by parts is Zt 0exp (βt) (t−ts)αdt=Zt 0.1 βexp (βt) (t−ts)αdt= 1 βexp (βt) (t−ts)α t 0−1 βαZ exp (βt) (t−ts)α−1dt, (134) such that another convolution with a smaller power appeares in the RHS of Eq. 134. This procedure can be used recursively until the null power, where the integral is over the exponential only, Zt 0exp (βt) (t−ts)αdt=1 β(t−ts)αexp (βt) t 0 −1 β2 α(t−ts)α−1exp (βt) t 0+1 β3 α(α−1) (t−ts)α−2exp (βt) t 0 −1 β4 α(α−1) (α−2) (t−ts)α−3exp (βt) t 0. . . + (−1)α1 βα α!Z exp (βt) (t−ts)α−αdt.(135) Rarranging the terms, Zt 0exp (βt) (t−ts)αdt= αX l=1(−1)l+11 βlα! (α−l+ 1)!(t−ts)α−l+1exp (βt) t 0 + (−1)α1 βα α! exp ( βt)1 β , (136) which can be further simplified to 18 Zt 0exp (βt) (t−ts)αdt= α+1X l=1(−1)l+11 βlα! (α−l+ 1)!(t−ts)α−l+1exp (βt) t 0. (137) Applying this result to Eq. 133 and neglecting the evaluation of this integral at t= 0, since it can be summed with constant C1, yields yp(t) = exp −f2,1,2tnpX k=0ckZt 0exp f2,1,2−ˆk tk+1X l=1(−1)l+1 ˆk−l k! (k−l+ 1)!(t−ts)k−l+1exp ˆkt dt.(138) Rearranging terms, using the linearity of the integration operator and defining r= k−l+ 1 yp(t) = exp −f2,1,2tnpX k=0ckk+1X l=1(−1)l+1 ˆk−lk! r!Zt 0exp f2,1,2t (t−ts)rdt.(139) Using the result from Eq. 137, yp(t) = exp −f2,1,2tnpX k=0ckk+1X l=1(−1)l+1 ˆk−lk! (r)! r+1X p=1(−1)p+1f−p 2,1,2r! (r+ 1−p)!(t−ts)r+1−pexp f2,1,2t , (140) that can be further simplified to yp(t) =npX k=0ckk+1X l=1(−1)l+1 ˆk−lk! r!r+1X p=1(−1)p+1f−p 2,1,2r! (r+ 1−p)!(t−ts)r+1−p,(141) which is itself another polynomial. Example with polynomial excitation A common example of excitation in electric circuitry analysis is the ramp function, a polynomial of first degree. Consider the ode subjected to ¨y(t) + 2 ˙y(t) + 10 y(t) =¯f(t) =t, (142) with ˙y(0) = y(0) = 0,t∈[0,4]s. 19 Steps used to solve Eq. 141, yp(t), are depicted in an algorithm form in Alg. 1, resulting in yp(t) = 0 .1t−0.02. The homogeneous solution is given by Eq. 110, that yields a system of equations for its integration constants, C1andC2, " exp (0) exp (0) −f2,1,2exp (0) −ˆkexp (0)#C1 C2 =y(0)−yp(0) ˙y(0)−˙yp(0) , (143) which simplifies to 1 1 −1−3i−1 + 3 iC1 C2 =0.02 −0.1 , (144) whose solution is C1= 0.0100−0.0133iandC2= 0.0100 + 0 .0133i. Thus, the complete solution is y(t) = 0 .1t−0.02 + (0 .0100−0.0133i) exp −(1 + 3 i)t + (0.0100 + 0 .0133i) exp −(1−3i)t (145) Algorithm 1: Evaluation of Eq. 141 k= 0 c0= 0 k= 1 c1= 1 l= 1 = ⇒(−1)l+1 ˆk−lk! r!= 0.1 + 0 .3i p= 1 = ⇒(−1)p+1f−p 2,1,2r! (r+1−p)!(t−ts)r+1−p= (0.1−0.3i)t p= 2 = ⇒(−1)p+1f−p 2,1,2r! (r+1−p)!(t−ts)r+1−p= 0.08 + 0 .06i l= 2 = ⇒(−1)l+1 ˆk−lk! r!= 0.08−0.06i p= 1 = ⇒(−1)p+1f−p 2,1,2r! (r+1−p)!(t−ts)r+1−p= 0.1−0.3i yp= 1 (0.1 + 0 .3i) (0.1−0.3i)t + 1 (0.1 + 0 .3i) (0.08 + 0 .06i) + 1 (0.08−0.06i) (0.1−0.3i) yp= 0.1t−0.02 The complete solution, Eq. 145, is then compared to the numerical solution obtained by using the Tsitouras 5/4Runge-Kutta method [18], Fig. 3.2.2. Both solutions converge to the same result. 3.3. Discontinuous excitation functions 3.3.1. Particular solution due to unitary impulses - Dirac’s deltas Consider that the normalized excitation ¯f(t)is given by nδDirac’s deltas at times tk ¯f(t) =nδX k=1ckδ(t−tk), (146) 20 Figure 3: Comparison between ramp solutions using Eq. 145 (green) and using the numerical method [18] (blue). with coefficients ck∈R. Particular solution given by Eq. 108 can be written as yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk t Zt 0exp ˆktnδX k=1ckδ(t−tk) dt dt, (147) or, using the linearity of the integral yp(t) =nδX k=1ckexp −f2,1,2tZt 0exp f2,1,2−ˆk t Zt 0exp ˆkt δ(t−tk) dt! dt. (148) The inner integral can be found by using the Filtering property of the Dirac’s delta ( Appendix, Eq. 214), such that Zt 0exp ˆkt δ(t−tk) dt= exp ˆktk H(t−tk) (149) where H(t−tk)is the Heaviside function at tk. Thus, yp(t) =nδX k=1ckexp −f2,1,2tZt 0exp f2,1,2−ˆk t exp ˆktk H(t−tk) dt,(150) 21 or yp(t) =nδX k=1ckexp −f2,1,2t exp ˆktkZt 0exp f2,1,2−ˆk t H(t−tk) dt.(151) ThisconvolutionoverafunctionmultipliedbytheHeavisidefunctioncanbeevaluated splittingtheintegrationdomainaccordingtothestepfunction, thus, resultinginachange of integration limits, Zt 0exp (βt)fk(t)H(t−tk) dt=Ztk 0exp (βt)fk(t)H(t−tk) dt | {z } =0+Zt tkexp (βt)fk(t)H(t−tk) dt, (152) where βis a generic exponent, such that Zt 0exp (βt)fk(t)H(t−tk) dt= Zt tkexp (βt)fk(t) dt! H(t−tk).(153) Thus, Eq. 151 can be written as yp(t) =nδX k=1ckexp −f2,1,2t exp ˆktkZt tkexp f2,1,2−ˆk t dtH(t−tk),(154) or yp(t) =nδX k=1ck f2,1,2−ˆk exp ˆk(tk−t) −exp f2,1,2(tk−t) H(t−tk).(155) Additionally, for under damped problems, f2,1,2is complex and ˆk=¯k f2,1,2=f∗ 2,1,2 (where ∗stands for complex-conjugate). Thus yp(t) =nδX k=1ck 2iℑ(f2,1,2) ef∗ 2,1,2(tk−t)−ef2,1,2(tk−t) H(t−tk), (156) or yp(t) =nδX k=1ck 2iℑ(f2,1,2)eℜ(f2,1,2)(tk−t) e−iℑ(f2,1,2)(tk−t)−eiℑ(f2,1,2)(tk−t) H(t−tk) (157) and user Euler’s identity for sin yp(t) =nδX k=1ck 2ℑ(f2,1,2)eℜ(f2,1,2)(tk−t)sin ℑ(f2,1,2)(t−tk) H(t−tk).(158) Example with unitary impulse Consider the problem 22 2¨y(t) + ˙y(t) + 2y(t) =δ(t−5), (159) with y(0) = ˙ y(0) = 0. The analytical solution is [5] yp(t) =(0 t <5 2√ 15e5−t 4sin√ 15 4(t−5) t≥5. (160) From the data, ¯c= 1/2,¯k= 1,nδ= 1andt1= 5. Constant f2,1,2can be found by solving Eq. 101 such that f2,1,2= 0.25 + 0 .96824 i. Equation 155 reduces to yp(t) =1 2(0.96824) i e(0.25−0.9682i)(5−t)−e(0.25+0 .96824 i)(5−t) H(t−5)(161) and, although complex, has negligible complex values, matching the analytical solution. Alternatively, using Eq. 158 yp(t) =1 2(0.96824)e0.25(5−t)sin 0.96824( t−5) H(t−5) (162) which is also identical to the analytical solution provided by [5]. 3.3.2. Particular solution due to step functions - Heaviside Lettheexcitationbedefinedasasumoffunctions fk(t)multipliedbyastepfunctions, ¯f(t) =nkX k=1fk(t)H(t−tk). (163) Inserting Eq. 163 into Eq. 108 yields yp(t) = exp −f2,1,2tZt 0exp f2,1,2−ˆk tZt 0exp ˆktnkX k=1fk(t)H(t−tk) dtdt, (164) as integration is a linear operator, yp(t) =nkX k=1exp −f2,1,2tZt 0exp f2,1,2−ˆk tZt 0exp ˆkt fk(t)H(t−tk) dtdt. (165) Using Eq. 153 in Eq. 165 for the inner integral results in yp(t) =nkX k=1exp −f2,1,2tZt 0exp f2,1,2−ˆk tZt tkexp ˆkt fk(t) dtH(t−tk) dt, (166) and for the outer integral 23 yp(t) =nkX k=1exp −f2,1,2tZt tkexp f2,1,2−ˆk tZt tkexp ˆkt fk(t) dtdtH(t−tk). (167) Example with Heaviside function Consider the problem from [5] ¨y(t) + 4y(t) =¯f(t) =  0, 0≤t <5 t−5 5,5≤t <10 1, t ≥10 ˙y(0) = y(0) = 0 , (168) such that the excitation function can be written as ¯f(t) =t−5 5H(t−5)−t−10 5H(t−10). (169) Using Eq. 101 f2,1,2=0 +√0−16 2= 2i (170) and Eq. 167 yp(t) =e−2it Zt te4itZt 5e−2itt−5 5 dtdtH(t−5)+ Zt te4itZt 5e−2it10−t 5 dtdtH(t−10)! (171) such that yp(t) = ie4it+ 4e10it−20e10i e2it−ie20i e−2it−10i 80H(t−5)+ − ie4it+ 4e20it−40e20i e2it−ie40i e−2it−20i 80H(t−10).(172) Figure 4 shows the real part of the solution given by Eq. 172 (the maximum complex part has a magnitude of 1×10−18when Eq. 172 was evaluated). This solution matches the result given by [5]. 24 Figure 4: Real part of the solution, Eq. 172, for problem defined in Eq. 168. 4. Further extensions - discrete or hard to convolute excitations The general solution, Eq. 108, can be used for general as well as some particular continuous forms of excitation f(t). Nonetheless, sometimes we are interested in using excitations only known at discrete time steps tkor using some f(t)hard to evaluate by convolution, although the integral of f(t)itself can be easy to evaluate (at least numerically). These cases can be addressed by using Heaviside series. 4.1. Heaviside series Let a function f(t)be approximated as a series of polynomials ck,mtmup to power nmmultiplied by Heaviside step functions at nkdiscrete time steps tk, f(t)≈˜f(t) =nkX k=0nmX m=0ck,mtmH(t−tk). (173) As the function is essentially approximated by a polynomial of order nm, it is fair to say that Eq. 173 provides an approximation of order nmforf(t), much like using Taylor series. Nonetheless, the coefficients of the polynomials are updated according to t, such that there is no need to center the approximation around a point like in Taylor series. Thus, the calculation of coefficients ck,mis carried out differently. The main idea for evaluating the coefficients ck,mis to preserve the integral of the original function, f(t)and its derivatives at each point tk. Examples with zero, fist 25 and second order approximations are presented in the following and provide a glance at the niche of application of this series representation. The application and advantage of representing the excitation function as a series of Heaviside steps is shown in Fig. 5, where the use of the integrating factor (and its convolution) is straightforward even beyond the subspace of excitation functions whose convolution is simply calculated. Application of Heaviside series Ease to integrateDiscrete domain Ease to convolute Space of functions ℝ → ℝ Figure 5: Diagram of the extension of the function space covered by the use of Heaviside series rep- resentation as excitation function, despite the simple convolution of the original excitation with the exponential function. The Heaviside series extends the viability of the integrating factor technique to exci- tation functions that are easy to integrate (inner integration in Eq. 108) but are hard to convolute. These series also allow the consideration of excitation functions (or distribu- tions) defined on discrete domains, as summarized in Fig. 5. 4.1.1. Zero order approximation For zero order approximation, nm= 0, Eq. 173 reduces to ˜f(t) =nkX k=0ckH(t−tk). (174) Using the preservation of the integral of f(t)at each interval t∈[tl, tl+1]and the mean value theorem for integrals ˜f(t) =lX k=0ck=1 ∆tlZtl+1 tlf(t) dt, tl≤t≤tl+1, (175) where ∆tl=tl+1−tl. Thus, the l-th coefficient, cl, is evaluated by cl=1 ∆tlZtl+1 tlf(t) dt−l−1X k=0ck. (176) 26 4.1.2. First order approximation For first order approximation, nm= 1, Eq. 173 reduces to ˜f=nkX k=0 ck,0+ck,1t H(t−tk). (177) The approximation at interval t∈[tl, tl+1]is given by ˜f(t) =lX k=0ck,0+lX k=0ck,1t=a0,l+a1,lt, tl≤t≤tl+1. (178) The slope of this polynomial, a1,l, can be tailored to be equal to a1,l=f(tl+1)−f(tl) ∆tl, (179) which, by the mean value theorem for derivatives, guarantees that the derivative of f(t)coincides with the derivative of ˜fat one point within the given interval, at least. Coefficient a0,lis obtained by Ztl+1 tl˜f(t) dt=a0,l∆tl+a1,l 2 t2 l+1−t2 l =Ztl+1 tlf(t) dt, (180) such that, a0,l=1 ∆tlZtl+1 tlf(t) dt−a1,l 2∆tl t2 l+1−t2 l . (181) Thus, using Eq. 178, the coefficients are evaluated as cl,0=1 ∆tlZtl+1 tlf(t) dt−a1,l 2∆tl t2 l+1−t2 l −l−1X k=0ck,0 (182) and cl,1=f(tl+1)−f(tl) ∆tl−l−1X k=0ck,1. (183) 4.1.3. Second order approximation For second order approximation, nm= 2, Eq. 173 reduces to ˜f(t) =nkX k=0 ck,0+ck,1t+ck,2t2 H(t−tk). (184) At interval t∈[tl, tl+1]the approximation function is given by ˜f(t) =lX k=0ck,0+lX k=0ck,1t+lX k=0ck,2t2=a0,l+a1,lt+a2,lt2, tl≤t≤tl+1.(185) 27 For this order of approximation, the first derivative of the Heaviside series represen- tation is set to be equal to the derivative of f(t)at the vicinity of each time point tk, what gives the advantage that the first derivative is smooth, although not defined at points tk. Thus, the derivative of the representation, ˜f, is continuous in the domain Ω:(−∞,∞)−{tl}, l∈0,1, . . . , n k. It is straightforward to observe that, if fis differen- tiable at tl,tlis an accumulation point and the limit of the derivative of ˜fconverges from both the left and the right sides. These conditions yield a linear system of equations 2a2,ltl+1+a1,l=˙f(tl+1), (186) and 2a2,ltl+a1,l=˙f(tl) (187) whose solutions are the coefficients a1,landa2,l, a2,l=˙f(tl+1)−˙f(tl) 2∆tl, (188) and a1,l=˙f(tl+1)−2a2,ltl+1. (189) Using the mean value theorem for integrals, one gets the coefficient a3to make the integral of the representation equal to the integral of the represented function at each point tk. Hence, the integral of the representation is Ztl+1 tla2t2+a1t+a0dt=a2 3 t3 l+1−t3 l +a1 2 t2 l+1−t2 l +a0(tl+1−tl),(190) which, comparing to the integral of the represented function, reduces to a0=1 ∆tlZtl+1 tlf(t) dt−1 ∆tla2 3 t3 l+1−t3 l −1 ∆tla1 2 t2 l+1−t2 l .(191) Turning these equations in terms of ck,jyields cl,2=˙f(tl+1)−˙f(tl) 2∆tl−l−1X k=0ck,2, (192) cl,1=˙f(tl+1)−2a2,ltl+1−l−1X k=0ck,1 (193) and cl,0=1 ∆tlZtl+1 tlf(t) dt−1 ∆tla2 3 t3 l+1−t3 l −1 ∆tla1 2 t2 l+1−t2 l −l−1X k=0ck,0.(194) 28 4.2. Heaviside series as excitation function Let the approximation function, ˜f, of Subsec. 4.1.1, 4.1.2 and 4.1.3, be used as excitation functions for for ordinary second order differential equations with constant coefficients. Let the order of the approximation be equal to nlk, ˜f=nkX k=0nlkX l=0ck,ltlH(t−tk). (195) Thus, applying Eq. 195 in Eq. 167 and using the linearity of the integration operator yields yp(t) =nkX k=1nlkX l=0exp −f2,1,2tZt tkexp f2,1,2−ˆk tZt tkexp ˆkt ck,ltldtdtH(t−tk). (196) The inner convolution is evaluated by using Eq. 137 yp(t) =nkX k=1nlkX l=0ck,lexp −f2,1,2tZt tkexp f2,1,2−ˆk t l+1X p=1(−1)p+1 ˆk−pl! s!tsexp ˆkt −l+1X p=1(−1)p+1 ˆk−pl! s!ts kexp ˆktk dtH(t−tk),(197) where s=l−p+ 1. The integral of the outer convolution is split among its integrands, yp(t) =nkX k=1nlkX l=0ck,lexp −f2,1,2t Zt tkexp f2,1,2tl+1X p=1(−1)p+1 ˆk−pl! s!tsdt −Zt tkexp f2,1,2t−ˆk(tk−t)l+1X p=1(−1)p+1 ˆk−pl! s!ts kdt H(t−tk).(198) Taking the constant terms out of the integrals yp(t) =nkX k=1nlkX l=0ck,lexp −f2,1,2t l+1X p=1(−1)p+1 ˆk−pl! s!Zt tkexp f2,1,2t tsdt −l+1X p=1(−1)p+1 ˆk−pl! s!ts kZt tkexp f2,1,2t−ˆk(tk−t) dt H(t−tk),(199) which has, then, the first convolution evaluated again using Eq. 137, and the second convolution evaluated by simple integration, 29 yp(t) =nkX k=1H(t−tk)nlkX l=0ck,lexp −f2,1,2t  l+1X p=1(−1)p+1 ˆk−pl! s!  s+1X q=1(−1)q+1 1 f2,1,2!q (s)! (s+ 1−q)!ts+1−qexp f2,1,2t −s+1X q=1(−1)q+1 1 f2,1,2!q s! (s+ 1−q)!ts+1−q kexp f2,1,2tk  − l+1X p=1(−1)p+1 ˆk−pl! s!ts k  exp f2,1,2t+ˆk(tk−t) −exp f2,1,2tk+ˆk(tk−tk)|{z} 0   1 f2,1,2−ˆk  .(200) Rearranging the terms yp(t) =nkX k=1H(t−tk)nlkX l=0ck,l    l+1X p=1(−1)p+1 ˆk−pl! s!s+1X q=1(−1)q+1 1 f2,1,2!q s! (s+ 1−q)!ts+1−q + exp f2,1,2(tk−t)  1 f2,1,2−ˆkl+1X p=1(−1)p+1 ˆk−pl! s!ts k −l+1X p=1(−1)p+1 ˆk−pl! s!s+1X q=1(−1)q+1 1 f2,1,2!q s! (s+ 1−q)!ts+1−q k   −exp ˆk(tk−t) 1 f2,1,2−ˆk!l+1X p=1(−1)p+1 ˆk−pl! s!ts k  ,(201) which, with further simplification, can be written as yp(t) =nkX k=1H(t−tk)nlkX l=0ck,l    l+1X p=1(−1)p+1 ˆk−pl! s!s+1X q=1(−1)q+1 1 f2,1,2!q s! (s+ 1−q)!ts+1−q + exp f2,1,2(tk−t)  l+1X p=1(−1)p+1 ˆk−pl! s!( ts k f2,1,2−ˆk −s+1X q=1(−1)q+1 1 f2,1,2!q s! (s+ 1−q)!ts+1−q k     −exp ˆk(tk−t) 1 f2,1,2−ˆk!l+1X p=1(−1)p+1 ˆk−pl! s!ts k  .(202) 30 4.2.1. Initial conditions and the Heaviside series Let a second order ordinary differential equation have a function ¯f(t)multiplied by a Heaviside step as excitation function, ¨y(t) + ¯c˙y(t) +¯ky(t) =¯f(t)H(t−t0), (203) which is the same as solving two differential equations, ( ¨y1(t) + ¯c˙y1(t) +¯ky1(t) = 0 ,ift≤t0 ¨y2(t) + ¯c˙y2(t) +¯ky2(t) =¯f(t).otherwise, (204) Att0, the initial conditions of y2must be equal to the values of y1and its derivative at this point, i.e.,y1(t0) =y2(t0)and˙y1(t0) = ˙y2(t0). As there is no excitation before t0 and a purely particular solution is sought after, the solution between t= 0andt=t0is y1(t) = 0and, consequently ˙y1(t) = 0. Thus, it follows that y2(t0) = 0and ˙y2(t0) = 0. This phenomenon is illustrated in Fig. 6. Figure 6: Solution of a linear second order differential equation with constant coefficients due to unitary Heaviside step at t0as excitation function (under damped case). It is straightforward that this holds true even when t0→0. Therefore, all solutions given by Eq. 202, i.e.using Heaviside series as excitation function, have yp(t) = 0and ˙yp(t) = 0as fixed initial conditions. Thus, the imposition of non-homogeneous initial conditions, y(tin)and˙y(tin), attin≤t0gets even simpler, through the following system of linear equations, 31  exp −f2,1,2tin exp −ˆktin −f2,1,2exp −f2,1,2tin −ˆkexp −ˆktin C1 C2 =y(tin) ˙y(tin) , (205) which, particularized for tin= 0, simplifies to " 1 1 −f2,1,2−ˆk#C1 C2 =y(0) ˙y(0) , (206) whose solution is C2=˙y(0) + f2,1,2y(0) f2,1,2−ˆk, (207) C1=y(0)−C2. (208) Thus, the constants C1andC2can be evaluated without any knowledge about the derivative of the particular response. 4.2.2. Examples The following second order differential equation with constant coefficients is used to show the use of Heaviside series as excitation function ¨y(t) + 2 ˙y(t) + 10 y(t) =f(t). (209) The periodic excitation was first addressed to compare the analytical solution (Subsub- sec. 3.2.1) with a solution using Heaviside series representation as excitation function. Then, an example using a random (Gaussian-distribution) was evaluated. The objective of the first example, beside comparing the representation with a purely analytical solu- tion, is to show a common case of excitation in Engineering and Applied Mathematics, but without representing the periodic signal as Fourier series. The objective of the second example is to show the usage of Heaviside series representation when there is no explicit formula for the excitation function. Example I : Let the excitation function f(t) =sin(4t)be represented by Heaviside series with four fixed ∆t:0.2s,0.4s,0.6sand0.8s. Figure 7 shows the original excitation function with blue cross symbols and four different zero order approximations, one for each ∆t. Figures 8 and 9 do the same but for first and second order approximations, respectively. One may observe that the zero order approximation represents the function rather poorly and, thus, needs even smaller steps of time. Quadratic approximation takes the best shape, with smoother corners, but gets out of phase and reduces its amplitude prematurely, much like what happens with zero order approximation. The odd behavior observed in the second order approximation is due to Eq. 188, since for large ∆tthe approximation enforces two very different slopes in the same interval. The amplitude 32 of the linear approximation proved to be less vulnerable to the increase of ∆tin this example. Figure 10 shows the solution obtained with Eq. 202 due to excitation of zero order approximation of the sine function. Again, the same four different ∆twere tested, and the solution using Eq. 202 was compared to the solution due to harmonic (Eq. 117) excitation with initial conditions of y(0) = ˙ y(0) = 0. Figure 11 and Figure 12 show the same procedure with first and second order approximations, respectively. It is straightforward to see that the zero order approximation performed worst, losing amplitude even for small ∆t. The second order approximation yielded slightly better results, but out of phase for ∆t= 0.8, as expected. Solutions obtained by using the first order for ∆t0.2sand0.4smatches the result given by Eq. 117. It is also perceptible that the loss of amplitude is not linearly dependent to the increase of ∆t, as the rate of decay of amplitude grows with the increase in ∆t. It is worth noticing that the analytical solutions obtained herein are not dependent on∆t. The dependence only appears when using the Heaviside series representation for ¯f(t). In this case, one can compare the dependence of the solution with respect to ∆t with traditional numerical integration methods, like the Newmark-Beta, for example. In such case, even when using an unconditionally stable integration methods, there are some constraints on ∆tto avoid aliasing. In the example studied in this section, the smaller period associated to the loading is π/2seconds, such that ∆t <0.15s(preferable less) are to be used to obtain a good representation of the response [11]. Figure 7: Comparison between the original excitation function (blue) and the representations using Heaviside series, ˜f(t), using zero order approximation with different discretizations. 33 Figure 8: Comparison between the original excitation function (blue) and the representations using Heaviside series, ˜f(t), using first order approximation. Example II A set of randomly generated points, obtained by using a Gaussian distribution, is used as the excitation function in this example. Figure 13 shows the random points as blue crosses and the zero order approximation. Figure 14 also shows the random points and the linear approximation. The second order approximation cannot be used in this example, since there is no information about the derivative of the original excitation function. Constants afor the approximations are evaluated by using the mean value theorem for integrals. Thereby, the integral for a discrete set of points is taken using trapezoid rule between each adjacent pair of points. Consequently, the step ∆tis equal to the time between two subsequent points. Figure 13 shows the comparison of the solution obtained by using first order approxi- mation by Heaviside series, ˜yH1(t), by using zero first order approximation by Heaviside series, ˜yH0(t)and by using the Tsitouras 5/4Runge-Kutta method [18] with linear in- terpolation for f(t). A total of 91time steps were needed for the Tsitouras method to match the solution obtained by using the first order Heaviside series ( 50points). Solution obtained with the zero order Heaviside series (red) show some discrepancies with respect to the amplitude at most of the peaks. 34 Figure 9: Comparison between the original excitation function (blue) and the representations using Heaviside series, ˜f(t), using second order approximation. 35 Figure 10: Comparison between the solutions of the differential equation using the original excitation function (blue), and the solutions with Heaviside series as excitation function, ˜y, using zero order ap- proximation. 36 Figure 11: Comparison between the solutions of the differential equation using the original excitation function (blue), and the solutions with Heaviside series as excitation function, ˜y, using first order ap- proximation. 37 Figure 12: Comparison between the solutions of the differential equation using the original excitation function (blue), and the solutions with Heaviside series as excitation function, ˜y, using second order approximation. 38 Figure 13: Comparison between the original excitation function (blue) and the representations using Heaviside series, ˜f(t), using zero order approximation with ∆t= 0.2s. 39 Figure 14: Comparison between the original excitation function (blue) and the representations using Heaviside series, ˜f(t), using first order approximation with ∆t= 0.2s. 40 Figure15: SolutionofthedifferentialequationusingfirstorderapproximationHeavisideseriesexcitation, ˜yH1(t)(green line), using zero order approximation Heaviside series excitation, ˜yH0(t)(red line), and using a numerical solver, yTsit 5(t)(blue dotted line). 41 5. Conclusion This work proposes a new technique to solve one dimensional linear ordinary differ- ential equations of order n≥2by a generalization of the Leibniz integrating factor for first order linear differential equations. The methodology was applied to second-order ordinary differential equation, since this type of ODE is of great importance in Applied Mathematics, Physics and Engi- neering. Special care was devoted to the constant coefficient case for different types of excitations. Various analytical excitation functions yielded analytical solutions in closed form, like periodic and polynomial functions. Closed form solutions were also obtained for discontinuous excitation, as for Dirac’s delta impulse and Heaviside step function. The Heaviside case was particularized for steps multiplied by polynomials, enlarging the range of applications with analytical solutions. In contrast to previous and well estab- lished methods, like undetermined coefficients and Laplace transform, no knowledge of a candidate solution nor the calculation of a inverse transformation was needed. The solution was instead analytically derived by means of double convolutions. The Heaviside series representation was derived by using the closed form solution obtained for Heaviside excitation multiplied by polynomials. Examples showed the ef- fectiveness of representing a function as a series of Heaviside steps, thus, augmenting the set of excitation functions solved within the generalized integrating factor method. A discussion about the discretization of the series was performed and shed light about the advantages of the proposed methodology compared to numerical alternatives, due to size of discretization, representation capability and the use of analytical solutions. For differential equations with general coefficients, the integrating factor depends on a particular solution of a Riccati differential equation. It was shown that the coefficients themselves might help finding a particular solution rather easily. Thus, the Riccati differentialequationposesnostrictlydirectbarriertothewideapplicationofthemethod, capable of giving accurate results and requiring no assumptions of solution candidates. Declaration of Competing Interest Theauthorsdeclarethattheyhavenoknowncompetingfinancialinterestsorpersonal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request Acknowledgments Authors acknowledge CNPq (process number 303900/2020-2) and FAPESC (process number 2021TR843 and scholarship CP48/2021). This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001 and process 88887.807649/2023-00. The authors are thankful for the suggestions of Professor Rafael Furlanetto in the early stages of this work. 42 Appendix - Convolution over Dirac’s delta distribution The convolution of a function over a delta of Dirac is usually defined as [13] Z∞ −∞f(t)δ(t−tk)dt=f(tk). (210) The integration limits can be split as Z∞ −∞f(t)δ(t−tk)dt=Z0 −∞f(t)δ(t−tk)dt+Zt 0f(t)δ(t−tk)dt+Z∞ tf(t)δ(t−tk)dt, (211) fortkstrictly positive, the integral from −∞to0is0by definition, thus, the integral from 0totcan be rewritten as Zt 0f(t)δ(t−tk)dt=Z∞ −∞f(t)δ(t−tk)dt−Z∞ tf(t)δ(t−tk)dt. (212) Thefilterorsiftingproperty of the delta of Dirac is due to the shape of this dis- tribution, i.e., it is null everywhere except in its discontinuity, thus, the function that multiplies the Dirac’s delta is constant at this point, for the discontinuity of the delta distribution is infinitely close to the tkpoint. Therefore, the value of the function can be taken out of the integral and the definition of the Dirac’s delta is used to show that Z∞ −∞f(t)δ(t−tk)dt=Ztk+τ tk−τf(t)δ(t−tk)dt=f(tk)Ztk+τ tk−τδ(t−tk)dt=f(tk), (213) hence, Equation 212 can be rewritten to Zt 0f(t)δ(t−tk)dt=f(tk)−f(tk)( 0, t≥tk 1, t < t k=f(tk)H(t−tk).(214)