CoupledsystemsofsecondorderOrdinaryDifferentialEquationsarefoundinmanyareas of engineering, such as vibrations and electric circuits. Thus, accurate and efficient solu- tionsareofgreatinterest. Thismanuscriptextendstheconceptofgeneralizedintegrating factor, recently introduced by the authors for one dimensional problems, to coupled n- dimensional systems of second order ODEs. The general formulation for time-varying coefficients is found and particularized to the constant coefficient case. Closed-form integrating factors are obtained for the constant coefficient case under mild assump- tions about the damping coefficient matrix. Analytical solutions for different types of continuous and discontinuous excitations are presented for problems with constant coeffi- cients and proportional damping, disregarding the damping level. Examples are provided for each type of excitation studied in this manuscript. Results obtained with the pro- posed approach are compared to the numerical solution obtained by using the traditional Newmark-beta method. Results show that the proposed procedure is accurate, not suf- fering with interpolation errors found in traditional numerical methods. Keywords: coupled system, ordinary differential equation, second order, analytical solution, integrating factor 1. Introduction Couple systems of second-order Ordinary Differential Equations (ODEs) are found in many applications in applied mathematics, physics and engineering. One example is the time-dependent equation of motion obtained by using the finite-element method, where, after the spatial discretization, results in a set of ncoupled ODEs M¨y(t) +C˙y(t) +Ky(t) =f(t), (1) ∗Corresponding author Email addresses: matheusj2009@hotmail.com (Matheus Janczkowski Fogaça), eduardo.cardoso@udesc.br (Eduardo Lenz Cardoso ) Preprint submitted to Elsevier April 25, 2023 where y(t)is the vector of dependent variables, tis the independent variable, M,Cand Kare coefficient matrices. There are many different approaches to solve the aforementioned system of coupled ODEs. Inabroadsense, undersomecircumstancesregarding C, itmaybepossibletoconvert the original coupled problem into nuncoupled problems [12]. When coefficient matrix Cis written as a linear combination of MandKone can use a change of basis given by the eigenvectors of the generalized eigenvalue problem (modal problem), resulting in nindependent one-dimensional second order ODEs [17]. Extensions to other general viscous coefficient matrix Cinvolve the solution of a quadratic eigenvalue problem and are discussed by [15]. This general procedure is known as Modal Superposition and it is common to consider just a subset of the neigenvectors to exclude (usually larger) modes intheresponseandtoavoidtheexpensivecomputationofalltheeigenpairs. Nonetheless, each individual equation still needs to be solved either analytically or numerically. Analytical solutions of one dimensional second order ODEs are generally known for simple harmonic excitations [16]. For other types of periodic functions, it is common to use the Fourier Series to obtain an harmonic approximation [5]. A recent proposal to obtain analytical solutions to a broader cases of excitations is presented by [3] and extended to the multidimensional case in this manuscript. Other approach is to re-write Eq. 2 by using state variables, converting a coupled second order equation to an extended system of first order coupled ODEs [4]. Analytical solutions can be derived by using exponential maps [14]. The most used approach to solve Eq. 2 is to use numerical methods. In this regard, there are many alternatives in the literature and a complete survey is out of the scope of this manuscript. Again, there are methods to directly solve Eq. 2 [10] and methods to solve the extended first order coupled system obtained by using state variables [4]. Numerical methods have as the main advantage the possibility to consider general forms of excitation. The objective of the present manuscript is to extend the methodology presented in [3] to analytically solve coupled second order ODE systems. Different forms of excitations are discussed in details. It is shown that for the constant coefficient case and under mild assumptions over C the evaluation of the complete response depends on the solution of a quadratic matrix equationwithclosedform. Theanalyticalsolutionsdonotdependonthelevelofdamping and can be used for any damping ratio. 2. General Formulation Consider a coupled system of nsecond order ODEs M(t)¨y(t) +C(t)˙y(t) +K(t)y(t) =f(t), (2) where MandKare time-dependent, real, positive-definite and symmetric n×nmatrices, Cis a positive semi-definite, time dependent symmetric real n×nmatrix, fis a general vector also depending on time tand vector yis the solution. For easy of notation, explicit dependency on time t,(t), will not be carried out in the following equations. Also, as M 2 is invertible, we multiply Eq. 2 by M−1to the left, such that I¨y+¯C˙y+¯Ky=¯f. (3) To solve this equation we start by splitting ¯Cas ¯C=F2,1,1+F2,1,2 (4) where both F2,1,1andF2,1,2are also n×ntime dependent matrices, not necessarily real. Thus, I¨y+F2,1,1˙y|{z} π2,2+F2,1,2˙y+¯Ky|{z} π2,1=¯f. (5) Next step is to multiply Eq. 6 by an invertible n×nmatrix µ2(generalized integrating factor) such that µ2I|{z} P2,2¨y+µ2F2,1,1|{z} ˙P2,2˙y+µ2F2,1,2|{z} P2,1˙y+µ2¯K|{z} ˙P2,1y=µ2¯f. (6) Using the definition of P2,2and its time derivative ˙P2,2=µ2F2,1,1=. (µ2I) =˙µ2 (7) and multiplying by µ−1 2results in µ−1 2˙µ2=F2,1,1 (8) such that µ2= expZ F2,1,1dt . (9) The same procedure can be applied to P2,1and its time derivative ˙P2,1=µ2¯K=. µ2F2,1,2 =˙µ2F2,1,2+µ2˙F2,1,2 (10) and multiplying by µ−1 2to the left ¯K=µ−1 2˙µ2F2,1,2+˙F2,1,2 (11) such that µ−1 2˙µ2= ¯K−˙F2,1,2 F−1 2,1,2 (12) with solution µ2= expZ ¯K−˙F2,1,2 F−1 2,1,2dt . (13) Equating Eqs. 8 and 12 F2,1,1= ¯K−˙F2,1,2 F−1 2,1,2 (14) such that F2,1,1F2,1,2= ¯K−˙F2,1,2 (15) 3 and using Eq. 4 F2,1,1¯C−F2,1,1 = ¯K−˙¯C+˙F2,1,1 (16) such that F2 2,1,1=˙¯C+F2,1,1¯C−¯K−˙F2,1,1 (17) is a coupled system of Riccati differential equations. Thus, solving Eq. 17 it is possible to find µ2using Eq. 9. As both P2,2andP2,1depend on known µ2, it is possible to re-write Eq. 6 as. P2,2˙y +. P2,1y =µ2¯f (18) and integrating w.r.t time  P2,2˙y + P2,1y =Z µ2¯fdt+C2 (19) where C2is a constant vector. Multiplying by another invertible time dependent matrix µ1(generalized integrating factor) µ1P2,2|{z} P1,1˙y+µ1P2,1|{z} ˙P1,1y=µ1Z µ2¯fdt+C2 | {z } h(20) such that ˙P1,1=µ1P2,1=. µ1P2,2 =˙µ1P2,2+µ1˙P2,2 (21) and by multiplying by µ−1 1to the left µ−1 1˙µ1= P2,1−˙P2,2 P−1 2,2 (22) with solution µ1= expZ P2,1−˙P2,2 P−1 2,2dt . (23) Again, P1,1is known after evaluating µ1and Eq. 20 can be written as . P1,1y =µ1h (24) with solution y=P−1 1,1 Z µ1Z µ2¯fdt+C2 dt+C1! , (25) where C2is a constant vector. Complete solution can be split into its particular and homogeneous counterparts as yp=P−1 1,1Z µ1Z µ2¯fdtdt, (26) and yh=P−1 1,1Z µ1C2dt+C1 . (27) 4 This procedure can be time consuming depending on how matrices M,CandKvary with time. Other issue is the solution of Eq. 17. It is worth to mention that the proposed procedure does not depend on a time dis- cretization and, despite being solved by using a computer, is analytical. For example, one can evaluate the response y(t)at any given time without knowing the solution of previ- ous times. Numerical methods, like the well established Newmark-beta method, work in a total different way, building the solution from time step to time step. Approximation errors associated to the interpolation hypothesis of each particular numerical method are sensitive to "large" time steps, such that approximation errors (as well as numerical errors) are expected when using methods relying on approximations. In the following we address the special case of constant coefficients. 3. Constant coefficients We now consider a discrete system of ncoupled ODEs with constant coefficients. Recalling Eq. 6 µ2I|{z} P2,2¨y+µ2F2,1,1|{z} ˙P2,2˙y+µ2F2,1,2|{z} P2,1˙y+µ2¯K|{z} ˙P2,1y=µ2¯f. (28) For the following steps, the time-derivatives of F2,1,1andF2,1,2are considered null, since these partitions are constant just like matrix ¯C. Thus, ˙P2,2=˙µ2=µ2F2,1,1. (29) Multiplying to the left by the inverse of the integrating factor yields µ−1 2˙µ2=F2,1,1. (30) The same procedure can be performed to P2,1, ˙P2,1=µ2¯K=˙µ2F2,1,2, (31) which, with further simplifications, yields µ−1 2˙µ2=¯KF−1 2,1,2. (32) Equating Eq. 30 to Eq. 34, and multiplying to the right by F2,1,2=¯C−F2,1,1results in F2,1,1¯C−F2,1,1 =¯K, (33) further simplifying this equation yields F2 2,1,1−F2,1,1¯C+¯K=0, (34) where 0is the null n×nmatrix. The analytical solution procedure for this matrix polynomial equation is discussed in Subsec. 3.1. The first generalized integrating factor, µ2, is found by means of Eq. 30 5 µ2= expZ F2,1,1dt = exp F2,1,1t . (35) Equation 28 can be written as . P2,2˙y +. P2,1y =µ2¯f (36) which can be integrated w.r.t time, yielding µ2˙y+µ2F2,1,2y=Z µ2¯fdt+C2, (37) where C2is a constant vector. This is a system of first order ordinary differential equa- tions. Thus, by multiplying the equation to the left by µ−1 2and by another generalized integrating factor µ1, and using F2,1,2=¯C−F2,1,1, yields µ1I|{z} P1,1˙y+µ1¯C−F2,1,1 |{z } ˙P1,1y=µ1µ−1 2Z µ2¯fdt+µ1µ−1 2C2. (38) Repeating the same procedure for the second integrating factor, ˙P1,1=˙µ1=µ1¯C−F2,1,1 , (39) or µ−1 1˙µ1=¯C−F2,1,1, (40) such that µ1= expZ¯C−F2,1,1 dt = exp¯C−F2,1,1 t . (41) Equation 38 can be written as . P1,1y =µ1µ−1 2Z µ2¯fdt+µ1µ−1 2C2, (42) such that by integrating and multiplying to left by µ−1 1yields y=µ−1 1Z µ1µ−1 2Z µ2¯fdtdt+µ−1 1Z µ1µ−1 2C2dt+µ−1 1C1, (43) where C1is another constant vector. As the integrating factors are already known for this case, they can be substituted in previous equation, restulting in y= exp −¯C−F2,1,1 tZ exp¯C−F2,1,1 t exp −F2,1,1tZ exp F2,1,1t¯fdtdt+ exp −¯C−F2,1,1 tZ exp¯C−F2,1,1 t exp −F2,1,1t C2dt+ exp −¯C−F2,1,1 t C1.(44) According to [7], exp (At)is an exponential map and, for such map, the property exp (A) exp ( B) = exp ( A+B)holds true if AB =BA(AandBcommute). Con- sequently, if ¯CandF2,1,1commute the exponential multiplications in Eq. 44 can be grouped together as 6 y= exp F2,1,1−¯C tZ exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt+ exp F2,1,1−¯C tZ exp¯C−2F2,1,1 t C2dtexp F2,1,1−¯C t C1,(45) The time-derivative of the exponential map can be used to evaluate the integral with theC1term in Eq. 45, y= exp F2,1,1−¯C tZ exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt+ exp −F2,1,1t¯C−2F2,1,1−1C2+ exp F2,1,1−¯C t C1,(46) where, again, there is a particular solution, yp= exp F2,1,1−¯C tZ exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt,(47) and a homogeneous solution, yh= exp −F2,1,1t¯C−2F2,1,1−1C2+ exp F2,1,1−¯C t C1,(48) which can be further simplified if C2absorbs the matrices multiplication, i.e., yh= exp −F2,1,1t C2+ exp F2,1,1−¯C t C1, (49) both valid if ¯CF2,1,1=F2,1,1¯C (50) commute. Complete expressions to evaluate constants C1andC2are developed in Appendix A and an efficient approach to evaluate Eq. 49 is discussed in Appendix F. Commutativity of matrices F2,1,1and¯Cis studied in the following. 3.1. Solving the quadratic equation, Eq. 34 Assume that the damping matrix, C, is given by Rayleigh damping model, i.e., proportional damping. Thus C=αM+βK, (51) normalizing by the mass matrix, one gets ¯C=αM−1M+βM−1K=αI+β¯K, (52) therefore, ¯C¯K=αI¯K+β¯K2=α¯KI+β¯K2=¯K¯C, (53) such that ¯Cand¯Kcommute. 7 As stated before, F2,1,1is given by the following matrix quadratic equation F2 2,1,1−F2,1,1¯C+¯K=0. (54) According to [8], a quadratic equation of the form X2+BX+D=0, (55) has a solution given by X=−1 2B+1 2h B2−4Di1 2, (56) ifBandDcommute. When a matrix solves a matrix polynomial equation it is called a solvent. In Eq. 54, the unknown matrix variable is at the left of the coefficient, F2,1,1¯C, hence, it is called a left solvent [9]; while in Eq. 56 the unknown variable is at the right, BX, thus, a right solvent. The question is whether solvent given by Eq. 56 solves Eq. 54, since Eq. 54 has the first coefficient equal to Iand¯Cand ¯Kcommute. To validate it, let the candidate of solvent, ˜F2,1,1, be calculated by Eq. 56 and be substituted into Eq. 54. ˜F2,1,1=1 2¯C+1 2h ¯C2−4¯Ki1 2, (57) applyied to Eq. 54, 1 4¯C2+1 4¯Ch ¯C2−4¯Ki1 2+1 4h ¯C2−4¯Ki1 2¯C+1 4h ¯C2−4¯Ki − (58) 1 2¯C2−1 2h ¯C2−4¯Ki1 2¯C+¯K. (59) The terms without the square roots are cancelled such that 1 4¯Ch ¯C2−4¯Ki1 2+1 4h ¯C2−4¯Ki1 2¯C−1 2h ¯C2−4¯Ki1 2¯C. (60) From Corollary 1.34 of [8], a special case of function defined on the spectra of the mul- tiplication of two matrices is the square root function, that gives [AB]1 2A=A[BA]1 2. (61) LetA=¯CandB=¯C−4¯C−1¯K, thus, h ¯C2−4¯Ki1 2¯C=¯Ch ¯C2−4¯C−1¯K¯Ci1 2; (62) as¯Cand¯Kcommute, h ¯C2−4¯Ki1 2¯C=¯Ch ¯C2−4¯C−1¯C¯Ki1 2=¯Ch ¯C2−4¯Ki1 2, (63) hence, 8 1 4¯Ch ¯C2−4¯Ki1 2+1 4h ¯C2−4¯Ki1 2¯C−1 2h ¯C2−4¯Ki1 2¯C=0, (64) therefore Eq. 57 is a solvent to Eq. 54 and F2,1,1=1 2¯C+1 2h ¯C2−4¯Ki1 2. (65) Another interesting fact is that, as Eq. 57 was derived for a right solvent and satisfies a left solvent as well, such that F2 2,1,1−F2,1,1¯C+¯K=F2 2,1,1−¯CF2,1,1+¯K=0, (66) thus, by inspection, one can realize that ¯CandF2,1,1commute too and Eq. 46 is consistent when Rayleigh damping is used. Another useful relations valid for Rayleigh damping are ¯CF2,1,1= αI+β¯K F2,1,1=αF2,1,1+β¯KF2,1,1 (67) and F2,1,1¯C=F2,1,1 αI+β¯K =αF2,1,1+βF2,1,1¯K, (68) thus, as ¯CandF2,1,1commute, by comparing the last two equations if follows that ¯K andF2,1,1also commute, i.e.,F2,1,1¯K=¯KF2,1,1. Other form of damping leading to commutativity between ¯CandF2,1,2is C=X jβjKj, (69) for example. A numerical method shall be used to find F2,1,1when other models of constant damp- ing failure to keep commutativity between the damping and the stiffness matrices are used [8]. 3.2. Under damped problems Equations 47 and 49 depend on exp F2,1,1−¯C (70) and exp −F2,1,1 . (71) For under-damped problems, it is possible to show that exp F2,1,1−¯C∗= exp −F2,1,1 . (72) First, as exp (A∗) = exp ( A)∗it is possible to write F∗ 2,1,1−¯C∗=−F2,1,1 (73) and as ¯Cis strictly real F∗ 2,1,1=−F2,1,1+¯C (74) 9 such that ℜ F2,1,1 =−ℜ F2,1,1 +¯C=⇒ ℜ F2,1,1 =1 2¯C (75) and −ℑ F2,1,1 =−ℑ F2,1,1 . (76) Thus, for under-damped problems the imaginary part of F2,1,1should be only related to the second term in the R.H.S of Eq 65 ℑ F2,1,1 =1 2h ¯C2−4¯Ki1 2(77) and this is true iff the matrix inside the square root is strictly real and has only real negative eigenvalues. The first condition is always fulfilled since both ¯Cand ¯Kare real matrices. Additionally, the term ¯C−2F2,1,1in Eq. 49 reduces to ¯C−2F2,1,1=−2ℑ F2,1,1 i (78) a purely imaginary matrix. In the following, we investigate the conditions necessary to satisfy that ¯C2−4¯K (79) has only real negative eigenvalues. Structural damping The simpler form of proportional damping is given by C=βK. In this case, the matrix inside the square root is ¯C2−4¯K=β2¯K2−4¯K. (80) It is worth noticing that the traditional eigenvalue problem is (K−λM)x=0 (81) where λis a strictly real positive eigenvalue and xis the associated eigenvector. Pre- multiplying by M−1¯K−λI x=0, (82) and multiplying by ¯K ¯K2x=λ¯Kx=λλx=λ2x. (83) Thus, Eq. 80 can also be written as β2XΛ2X−1−4XΛX−1(84) where Λis a diagonal matrix containing the eigenvalues of ¯KandXa matrix with its corresponding eigenvectors. Arranging X β2Λ2−4Λ X−1(85) such that β2λ2 i−4λi<0∀i= 1..n (86) andβmust satisfy β <r 4 λi∀i= 1..n. (87) 10 Proportional Damping Considering ¯C=αI+β¯K. (88) The square of ¯Cis ¯C2=α2I+ 2βα¯K+β2¯K2, (89) such that the term inside the square root is ¯C2−4¯K=α2I+ 2βα¯K+β2¯K2−4¯K. (90) Using the change of basis X α2I+ 2βαΛ+β2Λ2−4Λ X−1(91) such that α2+ 2βαλ i+β2λ2 i−4λi<0∀i= 1..n. (92) There are many possibilities to find suitable values of αandβin the previous equation. On can notice that the last two terms should be negative or zero if βsatisfies Eq. 87. If β2λ2 i−4λi=−ϵi (93) then ¯β=s −ϵi+ 4λi λ2 i. (94) Thus, α2+ 2¯βλi α−ϵi<0∀i= 1..n, (95) such that ¯α <−¯βλi ±q¯βλi2+ϵi. (96) 4. Particular solutions obtained by considering specific excitation functions So far, conditions for obtaining the solution of systems of coupled second order dif- ferential equations were discussed. To this end, considering constant coefficients, it was shown that permanent solution is given by Eq. 47. Conditions to solve the quadratic equation associated to the integrating factor F2,1,1were also discussed in details for pro- portional damping. Next sections are devoted to discuss further analytical solutions that can be obtained for some particular forms of excitations. A common benchmark problem is proposed to evaluate the different formulations obtained in the next sections. The problem is a 3DOFs system described by M= 2.0 0.0 0.0 0.0 2.0 0.0 0.0 0.0 1.0 , (97) K= 6.0−4.0 0 .0 −4.0 6 .0−2.0 0.0−2.0 6 .0 ×102(98) 11 and C=βK, (99) with y(0) =.y(0) =  0 0 0  (100) where the damping parameter βis specified in each example. This small problem was chosen since it makes easy to visualize the response. Nonetheless, it is a coupled system of ODEs and is able to show all the intended characteristics of the proposed approach. 4.1. Periodic excitations Let the excitation vector be defined as f(t) =g1(t)e1+g2(t)e2+. . .+gn(t)en, (101) where nis the dimension of the problem, i.e., the problem has ndegrees of freedom, andgj(t)is a function of time multiplying unitary vector ej. Normalizing by the mass matrix ¯f=g1(t)M−1e1+. . .+gn(t)M−1en=g1(t)v1+. . .+gn(t)vn.(102) Assuming gj(t) =nkX k=1cjkexp βjkt+ϕjk , (103) where nkis the number of terms, cjk∈Ris an amplitude, βjk=iωjk∈Ca complex angular frequency and ϕjk∈Ca complex phase. Substitutingthisexcitationvectorintotheinnerconvolutionoftheparticularsolution in Eq. 47, one gets Z exp F2,1,1t¯fdt=Z exp F2,1,1tnX j=1nkX k=1cjkexp βjkt+ϕjk vjdt,(104) as the vectors vjare independent of time, they can be left out of the integral to the right side and as the exponentials cjkexp βjkt+ϕjk are not matrices, they can be commuted with the exponential of matrices, Z exp F2,1,1t¯fdt=nX j=1nkX k=1Z cjkexp βjkt+ϕjk exp F2,1,1t dtvj.(105) Each one of these integrals can be evaluated by parts Z cjkexp βjkt+ϕjk exp F2,1,1t dt=Z" cjk βjk.  exp βjkt+ϕjk# exp F2,1,1t dt= cjk βjkexp βjkt+ϕjk exp F2,1,1t t t0−Zcjk βjkexp βjkt+ϕjk exp F2,1,1t F2,1,1dt,(106) 12 where it is possible to neglect the boundary value at t0as it is implicit in the integration constant C2. Grouping common terms and letting F2,1,1out of the integral at the right side for it being constant Z cjkexp βjkt+ϕjk exp F2,1,1t dt" I+1 βjkF2,1,1# =cjk βjkexp βjkt+ϕjk exp F2,1,1t , (107) such that Z cjkexp βjkt+ϕjk exp F2,1,1 tdt=cjk βjkexp βjkt+ϕjk exp F2,1,1t" I+1 βjkF2,1,1#−1 .(108) Substituting Eq. 108 into Eq. 105, Z exp F2,1,1t¯fdt=nX j=1nkX k=1cjk βjkexp βjkt+ϕjk exp F2,1,1t" I+1 βjkF2,1,1#−1 vj.(109) For the second convolution in Eq. 47, Z exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt=Z exp¯C−2F2,1,1 t nX j=1nkX k=1cjk βjkexp βjkt+ϕjk exp F2,1,1t" I+1 βjkF2,1,1#−1 vjdt,(110) asF2,1,1and¯Ccommute, Z exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt= ZnX j=1nkX k=1cjk βjkexp βjkt+ϕjk exp¯C−F2,1,1 t" I+1 βjkF2,1,1#−1 vjdt.(111) This equation can also be integrated by parts such that Z exp¯C−2F2,1,1 tZ exp F2,1,1t¯fdtdt= nX j=1nkX k=1cjk β2 jkexp βjkt+ϕjk exp¯C−F2,1,1 t" I+1 βjk¯C−F2,1,1#−1" I+1 βjkF2,1,1#−1 vj,(112) again, as ¯CandF2,1,1commute, the particular solution, Eq. 47, is given by 13 yp=nX j=1nkX k=1cjk β2 jexp βjkt+ϕjk" I+1 ωji¯C−F2,1,1#−1" I+1 βjkF2,1,1#−1 vj,(113) asA−1B−1= (BA)−1, yp=nX j=1nkX k=1cjk β2 jkexp βjkt+ϕjk " I+1 βjkF2,1,1#" I+1 βjk¯C−F2,1,1# −1 vj,(114) which simplifies to yp=nX j=1nkX k=1cjk β2 jkexp βjkt+ϕjk" I+1 βjk¯C−1 β2 jkF2,1,1¯C−F2,1,1#−1 vj,(115) and, using Eq. 33, reduces to yp=nX j=1nkX k=1cjk β2 jkexp βjkt+ϕjk " I+1 βjk¯C−1 β2 jk¯K#−1 vj , (116) It is worth noticing that previous expression can be further simplified. First, recall thatvj=M−1ej. Also, as A−1B−1= (BA)−1 " I+1 βjk¯C+1 β2 jk¯K#−1 M−1= M" I+1 βjk¯C+1 β2 jk¯K# −1 (117) such that yp=nX j=1nkX k=1cjk β2 jkexp βjkt+ϕjk" M+1 βjkC+1 β2 jkK#−1 ej. (118) It is also possible to manipulate the term 1/β2 j,kby inserting it into the matrix yp=nX j=1nkX k=1cjkexp βjkt+ϕjkh β2 jkM+βjkC+Ki−1 ej (119) and as βjk=iωjk yp=nX j=1nkX k=1cjkexp βjkt+ϕjkh K+iωjkC−ω2 jkMi−1 ej (120) 14 or yp(t) =nX j=1nkX k=1cjkexp βjkt+ϕjk KD−1 jkej, (121) where KDjk=h K+iωjkC−ω2 jkMi (122) is the Dynamic Stiffness Matrix for jk. It is worth noticing that particular solution given by Eq. 121, a superposition of harmonic responses, can be computed directly using M, CandKwith no need to compute F2,1,1. Solution of Eq. 121 can be written as a linear combination of pre processed vectors kjk KDjkkjk=ej, (123) since these operations do not depend on t. Thus, yp(t)can be efficiently evaluated as yp(t) =nX j=1nkX k=1cjkexp βjkt+ϕjk kjk, (124) for a given time t. The first derivative of Eq. 121 with respect to time t ˙yp(t) =nX j=1nkX k=1cjkβjkexp βjkt+ϕjk kjk, (125) is needed to evaluate the constant vectors C1andC2(Appendix A). When there are complex-conjugate pairs in the βcoefficients, an interesting imple- mentation optimization can be made for Eq. 123. Let two matrices KDjξandKDjνbe the dynamic stiffness matrices for the j-th degree of freedom, and let their βcoefficients be complex-conjugate, i.e.,βjξ=β∗ jν; following from Eq. 122, KDjξandKDjνare also complex-conjugate. Thus, kjξcan be written as kjξ=K−1 Djξej= K∗ Djν−1 ej, (126) which, using the property of inverse of complex-conjugate matrices, is simplified to kjξ= K−1 Djν∗ ej= K−1 Djνej∗ =k∗ jν. (127) Thereby, when there are pairs of complex-conjugate βcoefficients, just half of the dy- namic stiffness matrices must be evaluated. Example Consider the 3DOFs problem with excitation f=  0 3sin(4t) 0  , (128) 15 andβ= 1×10−2. This is an under damped harmonic problem with known permanent solution yp= [K+ 4iC−16M]−1  0 3 0  (129) with amplitude |yp|=  9.687 13.756 4.711  ×10−3, (130) such that yp(t) = sin(4 t)|yp|. (131) Using the proposed formulation, f=g2(t)e2= 3 sin(4 t)  0 1 0  (132) where, by using Euller’s identity, g2(t) = 3 sin(4 t) =3i 2exp (−4it)−3i 2exp (4 it) (133) andnk= 2,c21=3i 2,ω21=−4,β21=−4i,c22=−3i 2,ω22= 4,β22= 4i. Particular solution can then be found by using Eq. 121 yp=2X j=22X k=1cjkexp βjkt KD−1 jkej=3i 2exp (−4it)KD−1 21e2+−3i 2exp (4 it)KD−1 22e2 (134) with KD21=K−4iC−16M, (135) and KD22=K+ 4iC−16M. (136) Term z=KD−1 21e2is conjugate to KD−1 22e2. Thus, yp(t) = 3i 2exp(−4it)z−i 2exp(4 it)z∗ (137) or, by splitting zandz∗into their real and imaginary parts yp(t) = 3 i 2exp(−4it)ℜ(z) +i2 2exp(−4it)ℑ(z)−i 2exp(4 it)ℜ(z) +i2 2exp(4 it)ℑ(z)! . (138) Collecting common terms yp(t) = 3i 2exp(−4it)−i 2exp(4 it) | {z } sin(4t)ℜ(z)+3i 2exp(−4it) +i 2exp(4 it) | {z } 0ℑ(z)(139) 16 such that yp(t) = 3i 2 exp (−4it)−exp (4 it) | {z } sin(4t)  3.229 4.585 1.570  ×10−3(140) or yp(t) = sin(4 t)  9.686 13.756 4.711  ×10−3, (141) the expected solution. We can verify some conclusions obtained in previous sections. We start by computing ¯C2−4¯K= −1187 788 2 788 −1185 391 4 782 −2362  (142) which has real negative eigenvalues −2653.11,−1762.11and−318.77. Thus, its square root is a complex matrix, as discussed in previous sections 1 2h ¯C2−4¯Ki1 2= 16.026−6.29−0.410 −6.29 15 .62−2.53 −0.820−5.06 24 .03 i. (143) Matrix ¯Khas eigenvalues Λ= 79.85 0 .0 0 .0 0.0 445 .49 0 .0 0.0 0 .0 674 .66  (144) such that inequalities given by Eq. 87 for β= 1×10−2are also satisfied 1×10−2<(r 4 79.85,r 4 445.49,r 4 674.66) . (145) Matrix F2,1,1=1 2¯C+1 2h ¯C2−4¯Ki1 2= 1.5−1.0 0 .0 −1.0 1 .5−0.5 0.0−1.0 3 .0 + 16.026−6.29−0.410 −6.29 15 .62−2.53 −0.820−5.06 24 .03 i (146) and we can verify that F2,1,1and¯C−F2,1,1are complex-conjugate. Homogeneous solution, Eq. 49, for this problem is yh= exp −F2,1,1t C2+ exp F2,1,1−¯C t C1 (147) where C1=  2.49 3.38 1.18  ×10−4+i  2.36 2.87 1.08  ×10−3(148) 17 andC2=C∗ 1are obtained by solving Eqs. A.8 and A.9. The maximum order of magni- tude of the complex part of the homogeneous solution was 10−19, therefore, the homo- geneous solution is real-valued, as expected. Exponentials in Eq. 147 can be computed only once if a constant time step ∆tis used, as discussed in Appendix F. The real part of each DOFs of the complete solution y(t) =yh(t) +yp(t)is shown in Fig. 1 (solid lines). The solution yis compared to ˜yobtained by using traditional Newmark-beta method with standard parameters and ∆t= 0.001s (dotted lines). It is clear that both solutions match. Figure 1: Complete solution for the under damped problem with sinusoidal excitation. Solutions y1,y2 andy3, (real parts) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2 and˜y3, obtained by using the Newmark-beta method, are shown as dotted lines. Example We study the previous example with β= 10, an over damped problem. The proposed approach does not make any assumption on the level of damping such that the solution procedure used in the previous example does not change. The main differences are the fact that F2,1,1is a real matrix F2,1,1= 3−2 0 −2 3 −1 0−2 6 ×103, (149) 18 as well as the integration constants C1=  1.87 2.81 0.94  ×10−4(150) and C2=−  5.61 6.23 2.40  ×10−9. (151) Solutions obtained by using the proposed approach (solid lines) and the Newmark- beta method with ∆t= 0.01s (dotted lines) are shown in Fig. 2. Figure 2: Complete solution for the over damped problem with sinusoidal excitation. Solutions y1,y2 andy3, (real part) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2 and˜y3, obtained by using the Newmark-beta method, are shown as dotted lines. 4.2. Polynomial excitation Let the normalized excitation vector be defined as ¯f=g1(t)M−1e1+. . .+gn(t)M−1en=g1(t)v1+. . .+gn(t)vn,(152) where gj(t)are polynomial functions gj(t) =nkX k=0cjk t−tjk, (153) 19 cjk∈Rare coefficients, nkthe number of terms and tj∈Rtime shift. Substituting in Eq. 47, yields yp=e[F2,1,1−¯C]tZ exp¯C−2F2,1,1 tZ exp F2,1,1tnX j=1nkX k=0cjk t−tjkvjdtdt, (154) using the linearity of the integral operator and as vjdoes not depend on time results in yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkZ exp¯C−2F2,1,1 tZ exp F2,1,1t t−tjkdtdtvj. (155) The inner convolution can be evaluated by parts (here we use exinstead of exp(x)to shortening the equations), Z eAt t−tjαdt=Z . eAtA−1 t−tjαdt = t−tjαeAtA−1−Z eAt t−tjα−1dtαA−1 = t−tjαeAtA−1−α t−tjα−1eAtA−2+Z eAt t−tjα−2dtα(α−1)A−2 = t−tjαeAtA−1−α t−tjα−1eAtA−2+α(α−1) t−tjα−2eAtA−3+···+ (−1)m−1α(α−1). . .(α−m+ 2) t−tjα−m+1eAtA−m+ (−1)mZ eAt t−tjα−mdtα(α−1). . .(α−m+ 1)A−m m=α= t−tjαeAtA−1−α t−tjα−1eAtA−2+α(α−1) t−tjα−2eAtA−3+···+ (−1)α−1α! t−tj eAtA−α+ (−1)αα!eAtA−α−1,(156) which, by arranging all the terms in a single sum, yields Z exp (At) t−tjαdt= exp ( At)α+1X l=1(−1)l+1 α! (α−l+ 1)! t−tjα−l+1A−l.(157) Using this integral formula with Eq. 155 and the considering that ¯CandF2,1,1 commute, yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkZ exp¯C−F2,1,1 t k+1X l=1(−1)l+1 k! (k−l+ 1)! t−tjk−l+1F−l 2,1,1dtvj. (158) 20 As the integral is a linear operator yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)! Z exp¯C−F2,1,1 t t−tjk−l+1F−l 2,1,1dtvj, (159) and using again Eq. 157, yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)!exp¯C−F2,1,1 t k−l+2X p=1(−1)p+1 (k−l+ 1)! (k−l−p+ 2)! t−tjk−l−p+2¯C−F2,1,1−pF−l 2,1,1vj.(160) Since ¯CandF2,1,1commute for Rayleigh damping, the matrix exponentials can be cancelled. From Eq. 33, it follows that ¯C−F2,1,1=F−1 2,1,1¯K. Substituting this relation and applying the inverse property A−1B−1= (BA)−1results in yp=nX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)! k−l+2X p=1(−1)p+1 (k−l+ 1)! (k−l−p+ 2)! t−tjk−l−p+2h Fl 2,1,1F−p 2,1,1¯Kpi−1 vj,(161) and as vj=M−1ejand using again the inverse property, yp=nX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)! k−l+2X p=1(−1)p+1 (k−l+ 1)! (k−l−p+ 2)! t−tjk−l−p+2h MFl−p 2,1,1¯Kpi−1 ej. (162) As¯KandF2,1,1commute for Rayleigh damping, Eq. 67, and pandlare integers, property in Eq. C.1 holds, yp(t) =nX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)! k−l+2X p=1(−1)p+1 (k−l+ 1)! (k−l−p+ 2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej. (163) 21 Derivative of previous equation with respect tis needed to evaluate constant vectors C1andC2(Appendix A) ˙yp(t) =nX j=1nkX k=0cjkk+1X l=1(−1)l+1 k! (k−l+ 1)! k−l+2X p=1(−1)p+1 (k−l+ 1)! (k−l−p+ 2)!(k−l−p+ 2) t−tjk−l−p+1h M¯KpFl−p 2,1,1i−1 ej,(164) valid when k−l−p+ 2is not zero. Example Consider the 3DOF problem with β= 1×10−6and excitation f=  0 10t−t2 0  , (165) such that g2(t) =c20+c21t+c22t2with c20= 0,c21= 10,c22=−1andt2= 0. Evaluation of particular solution yp(t), Eq. 163, is shown by means of Alg. 1 where j= 2andnk= 2. The terms obtained in the Algorithm result in yp(t) = (2 t−10) M¯K2F−1 2,1,1−1 + M¯KF2,1,1−1 e2+ 10t−t2 M¯K−1e2− 2 M¯K2−1 e2−2 M¯K3F−2 2,1,1−1 e2−2 M¯KF2 2,1,1−1 e2.(166) The derivative of the particular response w.r.t time tis given by .yp= 2 M¯K2F−1 2,1,1−1 + M¯KF2,1,1−1 e2+ (10−2t) M¯K−1e2(167) such that by using Eqs. A.8 and A.9, one gets the integration constants, C1=  8.8792 0.1457 4.6479  ×10−5+i  1.5639 1.9436 0.7237  ×10−3, (168) withC2=C∗ 1. Thus, by using the real part of Eq. 49, one gets the complete solution y(t) =yh(t) + yp(t)for the three DOFs, as shown in Fig. 3 (solid lines). The solution yis compared to the solution ˜yobtained by using traditional Newmark-beta method with standard parameters and ∆t= 0.001s (dotted lines). It is clear that both solutions match. 4.3. Dirac’s delta distribution Let the normalized excitation vector be defined as ¯f=g1(t)M−1e1+. . .+gn(t)M−1en=g1(t)v1+. . .+gn(t)vn,(169) 22 Algorithm 1: Evaluation of Eq. 163 for j= 2andnk= 2. k= 0 c1,0= 0 k= 1 c1,1= 10 l= 1 = ⇒(−1)l+1 k! (k−l+1)!= 1 p= 1 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= t M¯K−1e2 p= 2 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= − M¯K2F−1 2,1,1−1 e2 l= 2 = ⇒(−1)l+1 k! (k−l+1)!=−1 p= 1 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej=  M¯KF2,1,1−1e2 k= 2 c1,2=−1 l= 1 = ⇒(−1)l+1 k! (k−l+1)!= 1 p= 1 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= t2 M¯K−1e2 p= 2 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= −2t M¯K2F−1 2,1,1−1 e2 p= 3 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= 2 M¯K3F−2 2,1,1−1 e2 l= 2 = ⇒(−1)l+1 k! (k−l+1)!=−2 p= 1 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= t M¯KF2,1,1−1e2 p= 2 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej= − M¯K2−1e2 l= 3 = ⇒(−1)l+1 k! (k−l+1)!= 2 p= 1 = ⇒(−1)p+1 (k−l+1)! (k−l−p+2)! t−tjk−l−p+2h M¯KpFl−p 2,1,1i−1 ej=  M¯KF2 2,1,1−1 e2 23 Figure 3: Complete solution for the three DOFs example subjected to a polynomial excitation. Solutions y1,y2andy3, (real part) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2and˜y3, obtained by using the Newmark-beta method are shown as dotted lines. where gj(t)are Dirac’s delta distributions gj(t) =nkX k=0cjkδ t−tjk , (170) cjk∈Rare coefficients, nkthe number of terms and tjk∈Rtime shifts. Substituting in Eq. 47, yields yp= exp F2,1,1−¯C tZ exp¯C−2F2,1,1 tZ exp F2,1,1tnX j=1nkX k=0cjkδ t−tjk vjdtdt, (171) using the linearity of the integral operator and as vjdoes not depend on time results in yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkZ exp¯C−2F2,1,1 tZ exp F2,1,1t δ t−tjk dtdtvj. (172) Equation B.5, Appendix B, is used to evaluate the inner convolution 24 yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkZ exp¯C−2F2,1,1 t exp F2,1,1tk H t−tjk dtvj, (173) and, as ¯C−2F2,1,1andF2,1,1commute, the two exponential maps can be simplified to yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkZ exp¯C−2F2,1,1 t+F2,1,1tk H t−tjk dtvj, (174) hence, by using the relation derived in [3] Zt 0exp (At)fk(t)H(t−tk) dt= Zt tkexp (At)fk(t) dt! H(t−tk),(175) where Ais a matrix. The outer convolution is evaluated to yp= exp F2,1,1−¯C tnX j=1nkX k=0cjkH t−tjkZt tkexp¯C−2F2,1,1 t+F2,1,1tk dtvj, (176) which is yp=nX j=1nkX k=0cjkH t−tjk exp F2,1,1−¯C t exp¯C−2F2,1,1 t+F2,1,1tk −exp¯C−F2,1,1 tk¯C−2F2,1,1−1vj.(177) The commutativity of the first power is checked,  F2,1,1−¯C t¯C−2F2,1,1 t+F2,1,1tk = t2 F2,1,1¯C−2F2 2,1,1−¯C2+ 2¯CF2,1,1 +ttk F2 2,1,1−¯CF2,1,1 (178) ¯C−2F2,1,1 +F2,1,1tk F2,1,1−¯C t= t2 ¯CF2,1,1−2F2 2,1,1−¯C2+ 2F2,1,1¯C +ttk F2 2,1,1−F2,1,1¯C , (179) as¯CandF2,1,1commute, Equation 178 and Equation 179 are equal. Commutativity for the other power is straightforward, since a matrix commutes with itself regardless of different constants multiplying it. Hence, also substituting vj=M−1ej, 25 yp=nX j=1 nkX k=0cjkH t−tjkh exp −F2,1,1t+F2,1,1tk −exp F2,1,1−¯C t+¯C−F2,1,1 tk! ¯C−2F2,1,1−1M−1ej,(180) rearranging the terms and applying the property A−1B−1= (BA)−1results in yp=nX j=1 nkX k=0cjkH t−tjkh exp F2,1,1(tk−t) −exp¯C−F2,1,1 (tk−t)!  C−2MF 2,1,1−1ej. (181) In Eq. 73, it was shown that, for under-damped problems, ¯C−F2,1,1is the complex- conjugate of F2,1,1, thus, according to [7], exp (A∗) = exp ( A)∗, and it follows that yp=nX j=1nkX k=0cjkH t−tjkh exp F2,1,1(tk−t) −exp F2,1,1(tk−t)∗i C−2MF 2,1,1−1ej, (182) which further simplifies to yp= 2inX j=1nkX k=0cjkH t−tjk ℑ exp F2,1,1(tk−t) C−2MF 2,1,1−1ej.(183) Example Consider the 3DOF problem with β= 1×10−2and excitation f=  0 δ(t−1)−δ(t−5) 0  , (184) such that g2(t) =c20δ(t−t0) +c21δ(t−t1)with c20= 1,c21=−1,t0= 1andt0= 5. Using the conjugacy of ¯C−F2,1,1andF2,1,1, one can use Eq. 183 directly, yp= 2i H(t−1)ℑ exp F2,1,1(1−t) − H(t−5)ℑ exp F2,1,1(5−t)¯C−2MF 2,1,1−1e2, and as shown, for homogeneous initial conditions and impulse excitation, y=yp. 26 TheresponsewascomparedtothesolutionobtainedbyusingthetraditionalNewmark- beta method with ∆t= 0.001. Each impulse was approximated by [6] δ(t−t0)≈1 2ϵ 1 + cosπ(t−t0) ϵ! t0−ϵ≤t≤t0+ϵ. (185) with ϵ= ∆tto impose the δs to the numerical method (the proposed methodology does not need such approximation). Solutions are shown in Fig. 4 where the solid lines correspond to the real part of the proposed approach and the dotted lines to the numerical solution (Newmark-beta method with ∆t= 0.001s). Figure 4: Complete solution for the three DOFs example subjected to two opposed Dirac’s deltas at t= 1 andt= 5, respectively. Solutions y1,y2andy3, (real part) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2and˜y3, obtained by using the Newmark-beta method, are shown as dotted lines. 4.4. Heaviside step function Let the normalized excitation vector be defined as ¯f=g1(t)M−1e1+. . .+gn(t)M−1en=g1(t)v1+. . .+gn(t)vn,(186) where gj(t)are functions multiplied by other time-dependent functions, gj(t) =nkX k=0fjk(t)H t−tjk , (187) 27 cjk∈Rare coefficients, nkthe number of terms and tjk∈Rtime shifts. Substituting in Eq. 47, yields yp= exp F2,1,1−¯C tZ exp¯C−2F2,1,1 tZ exp F2,1,1tnX j=1nkX k=0fjk(t)H t−tjk vjdtdt, (188) using the linearity of the integral operator and as vjdoes not depend in time results in yp= exp F2,1,1−¯C tnX j=1nkX k=0Z exp¯C−2F2,1,1 tZ exp F2,1,1t fjk(t)H t−tjk dtdtvj. (189) Equation 175 is used to evaluate the inner convolution, yp= exp F2,1,1−¯C tnX j=1nkX k=0Z exp¯C−2F2,1,1 tZt tjkexp F2,1,1t fjk(t) dtH t−tjk dtvj, (190) and, again, Equation 175 is used to evaluate the outer integral, yp= exp F2,1,1−¯C tnX j=1nkX k=0H t−tjkZt tjkexp¯C−2F2,1,1 tZt tjkexp F2,1,1t fjk(t) dtdtvj. (191) ThissolutiondependsonthefunctionthatmultipliestheHeavisidefunction; however, the cases derived in previous subsections cover an extent set of widely used excitation functions, e.g., harmonic, polynomial and constant (a constant function is a polynomial of order 0). Thus, the integrals in Eq. 191 can be found by using previous results as reference. In the following, we present the result for a second order polynomial. 4.4.1. Particularizing Heaviside excitation for second order polynomial In Eq. 187, fj,kcan be particularized to be a second order polynomial, i.e., gj(t) =nkX k=0 cjk0+cjk1t+cjk2t2 H t−tjk =nkX k=0fjk. (192) Applying this excitation in Eq. 191 yields yp(2)= exp F2,1,1−¯C tnX j=1nkX k=0H t−tjkZt tjkexp¯C−2F2,1,1 t Zt tjkexp F2,1,1t cjk0dt+Zt tjkexp F2,1,1t cjk1tdt+ Zt tjkexp F2,1,1t cjk2t2dt! dtvj.(193) 28 The inner convolutions are analytically integrated using Eq. 157, and, rearranging the terms, one gets yp(2)= exp F2,1,1−¯C tnX j=1nkX k=0H t−tjkZt tjkcjk2t2exp¯C−F2,1,1 t F−1 2,1,1+ texp¯C−F2,1,1 t cjk1F−1 2,1,1−2cjk2F−2 2,1,1 + exp¯C−F2,1,1 t 2cjk2F−3 2,1,1 −cjk1F−2 2,1,1+cjk0F−1 2,1,1 + exp¯C−2F2,1,1 t −cjk2t2 jkexp F2,1,1tjk F−1 2,1,1+ 2cjk2tjkexp F2,1,1tjk F−2 2,1,1−cjk1tjkexp F2,1,1tjk F−1 2,1,1 −2cjk2exp F2,1,1tjk F−3 2,1,1+cjk1exp F2,1,1tjk F−2 2,1,1 −cjk0exp F2,1,1tjk F−1 2,1,1 dtvj,(194) which can be again convoluted by using Eq. 157. Multiplying by exp F2,1,1−¯C t , yields yp(2)=nX j=1nkX k=0H t−tjk cjk2t2¯C−F2,1,1−1F−1 2,1,1+t¯C−F2,1,1−1( cjk1F−1 2,1,1−2cjk2F−2 2,1,1 −2cjk2¯C−F2,1,1−2F−1 2,1,1 + 2cjk2¯C−F2,1,1−3F−1 2,1,1 −¯C−F2,1,1−2 cjk1F−1 2,1,1−2cjk2F−2 2,1,1 +¯C−F2,1,1−1 2cjk2F−3 2,1,1 −cjk1F−2 2,1,1+cjk0F−1 2,1,1 + exp −F2,1,1t¯C−2F2,1,1−1( −cjk2t2 jkexp F2,1,1tjk F−1 2,1,1+ 2cjk2tjkexp F2,1,1tjk F−2 2,1,1 −cjk1tjkexp F2,1,1tjk F−1 2,1,1−2cjk2exp F2,1,1tjk F−3 2,1,1 +cjk1exp F2,1,1tjk F−2 2,1,1−cjk0exp F2,1,1tjk F−1 2,1,1 + exp¯C−F2,1,1 tjk−t ( −cjk2t2 jk¯C−F2,1,1−1F−1 2,1,1+ 2cjk2tjk¯C−F2,1,1−2F−1 2,1,1 −2cjk2¯C−F2,1,1−3F−1 2,1,1−tjk cjk1F−1 2,1,1−2cjk2F−2 2,1,1 +¯C−F2,1,1−2( cjk1F−1−2cjk2F−2 2,1,1 −¯C−F2,1,1−1 2cjk2F−3−cjk1F−2+cjk0F−1 − exp¯C−F2,1,1 tjk−t −F2,1,1tjk¯C−2F2,1,1−1( −cjk2t2 jkexp F2,1,1tjk F−1 2,1,1+ 2cjk2tjkexp F2,1,1tjk F−2 2,1,1− cjk1tjkexp F2,1,1tjk F−1 2,1,1−2cjk2exp F2,1,1tjk F−3 2,1,1+ cjk1exp F2,1,1tjk F−2 2,1,1−cjk0exp F2,1,1tjk F−1 2,1,1 vj.(195) ByusingEqs. G.1andD.1,theterm exp F2,1,1tjk canbecommutedwith¯C−2F2,1,1−1. Thus, the expression can be further simplified to 29 yp(2)=nX j=1nkX k=0H t−tjk cjk2t2¯C−F2,1,1−1F−1 2,1,1+t¯C−F2,1,1−1  cjk1F−1 2,1,1−2cjk2F−2 2,1,1 −2cjk2¯C−F2,1,1−2F−1 2,1,1 + 2cjk2¯C−F2,1,1−3F−1 2,1,1 −¯C−F2,1,1−2 cjk1F−1 2,1,1−2cjk2F−2 2,1,1 +¯C−F2,1,1−1 2cjk2F−3 2,1,1 −cjk1F−2 2,1,1+cjk0F−1 2,1,1 + exp F2,1,1 tjk−t¯C−2F2,1,1−1  − cjk2t2 jk+cjk1tjk+cjk0 F−1 2,1,1+ 2cjk2tjk+cjk1 F−2 2,1,1−2cjk2F−3 2,1,1 + exp¯C−F2,1,1 tjk−t ¯C−F2,1,1−1 F−1 2,1,1 −cjk2t2 jk−cjk1tjk−cjk0 +F−2 2,1,1 2cjk2tjk+cjk1 −2cjk2F−3 2,1,1 +¯C−F2,1,1−2 F−1 2,1,1 2cjk2tjk+cjk1 −2cjk2F−2 2,1,1 −2cjk2¯C−F2,1,1−3F−1 2,1,1−¯C−2F2,1,1−1 −F−1 2,1,1 cjk2t2 jk +cjk1tjk+cjk0 +F−2 2,1,1 2cjk2tjk+cjk1 −2cjk2F−3 2,1,1) M−1ej.(196) Twoexamplesarediscussedtoshowtheflexibilityprovidedbytheproposedapproach, since complicate loading scenarios cab be built by a linear combination of polynomials and Heavisides. Example Consider the 3DOF problem with β= 1×10−2and excitation f2(t) = −30 + 40 t−10t21≤t≤3 0 t∈[0,1]∪[3,inf](197) or f=  0 (−30 + 40 t−10t2)H(t−1) + (30 −40t+ 10t2)H(t−3) 0  ,(198) such that c200=−30,c201= 40,c202=−10,t20= 1,c210= 30,c211=−40,c212= 10 andt21= 3. Using the real part of the solution given by Eq. 196 and also solving the problem by using the standard Newmark-beta numerical method lead to the complete responses shown in Fig. 5. Solution obtained by using Eq. 196 is shown as solid lines and solution obtained by using the Newmark method with ∆t= 0.001s is shown as dotted lines. Example 30 Figure 5: Complete solution for the three DOFs example subjected to a quadratic load between t= 1 andt= 3s. Solutions y1,y2andy3, (real part) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2and˜y3, obtained by using the Newmark-beta method, are shown as dotted lines. Consider the 3DOF problem with β= 1×10−2. Assuming an unitary step at DOF 2between t= 1andt= 5 f=  0 H(t−1)− H(t−5) 0  , (199) such that g2(t) =c200H(t−t20) +c210H(t−t21)with c200= 1,c210=−1,t20= 1and t21= 5. Solution provided by Eq. 196 can be particularized to order zero yp(0)=nX j=1nkX k=0H t−tjkn¯C−F2,1,1−1cjk0F−1 2,1,1 −cjk0exp F2,1,1 tjk−t¯C−2F2,1,1−1F−1 2,1,1+ exp¯C−F2,1,1 tjk−t  −cjk0¯C−F2,1,1−1F−1 2,1,1+cjk0¯C−2F2,1,1−1F−1 2,1,1 M−1ej.(200) Using the data for this example 31 yp(0)=H(t−1)n¯C−F2,1,1−1F−1 2,1,1 −exp F2,1,1(1−t)¯C−2F2,1,1−1F−1 2,1,1+ exp¯C−F2,1,1 (1−t)  −¯C−F2,1,1−1F−1 2,1,1+¯C−2F2,1,1−1F−1 2,1,1 M−1e2+ H(t−5)n¯C−F2,1,1−1(−1)F−1 2,1,1 + exp F2,1,1(5−t)¯C−2F2,1,1−1F−1 2,1,1+ exp¯C−F2,1,1 (5−t) ¯C−F2,1,1−1F−1 2,1,1−¯C−2F2,1,1−1F−1 2,1,1 M−1e2.(201) Using the real part of the solution given by previous equation and also solving the problem using the standard Newmark-beta numerical method lead to the complete re- sponses shown in Fig. 6. Solution obtained by using Eq. 201 is shown as solid lines and solution obtained by using the Newmark method with ∆t= 0.001s is shown as dotted lines. Figure 6: Complete solution for the three DOFs example subjected to an unitary step between t= 1 andt= 5s. Solutions y1,y2andy3, (real part) obtained by using the proposed approach, are shown as solid lines. Solutions ˜y1,˜y2and˜y3, obtained by using the Newmark-beta method are shown as dotted lines. 32 4.4.2. Initial conditions of Heaviside and Dirac’s delta excitation As all previous excitation functions were broken in the canonical basis of the Rn, each integration was evaluated separately for each functional component and, then, multiplied by the respective basis vector. Nonetheless, for both the Heaviside and the Dirac’s delta excitations, the component of the solution that multiply each basis vector has null response and derivative at the initial point. Consequently, the initial conditions for these excitations are zero such that yp(t) =0∀t[3]. This result was observed in previous examples. The consideration of non homogeneous initial conditions is discussed in Appendix A. 4.4.3. Matrix complexity for polynomial particularized Heaviside Heaviside excitation particularized to polynomial fjkyields expensive matrices eval- uations, like power and inverse operations. This is expected, since similar complexity is observed when addressing polynomial excitations. As it can be observed in Eq. 196, this complexity increases with the order of fjk. However, some calculation can be avoided by using the quadratic matrix equation, Eq. 34. For example, by factoring Eq. 34, F2 2,1,1−F2,1,1¯C+¯K=F2,1,1 F2,1,1−¯C +¯K=0, (202) which is then rearranged to F2,1,1¯C−F2,1,1 =¯K. (203) Taking the inverse of both sides, results in ¯C−F2,1,1−1F−1 2,1,1=¯K−1=K−1M, (204) thus, the inverse of F2,1,1can be calculated by F−1 2,1,1=¯C−F2,1,1 K−1M; (205) conversely, the inverse of¯C−F2,1,1 can be evaluated by ¯C−F2,1,1−1=K−1MF 2,1,1. (206) Hence, the two most common inverses in Eq. 196 are evaluated using just one inverse, the inverse of K. These relations can also be used to decrease the computational cost for polynomial excitation. The inverse of Kis particularly interesting in the Finite Element Analysis (FMA) context, since the stiffness matrix is sparse, enabling the use of fast, efficient and tailored algorithms for this kind of matrices [13, 2]. Therefore, the inverse of more complicate matrices, like F2,1,1and¯C−F2,1,1 , can be computed by using an inverse that is cheaper to evaluate. The matrix¯C−2F2,1,1 is trickier to simplify, but it can be done using Eq. 65, ¯C−2F2,1,1 =− ¯C2−4¯K1 2; (207) if¯Cis given by proportional damping, 33 ¯C−2F2,1,1 =− β2¯K2+ (2αβ−4)¯K+α2I1 2. (208) Equation 208 can be approximated if we assume both αandβas very small values. In such case¯C−2F2,1,1 ≈ −2i¯K1 2, (209) and the inverse of¯C−2F2,1,1 can be approximated to ¯C−2F2,1,1−1≈i 2 M−1K−1 2=i 2 K−1M1 2. (210) Again, a matrix inverse can be calculated using the inverse of the stiffness matrix. The approximation presented in Eq. 209, however, can also be used in Eq. 65 to reduce computation costs. We did not use such approximations in the examples discussed in this manuscript. 5. Conclusion This work extended the use of the Generalized Integrating Factor [3] to coupled sys- tems of second order ODEs. Expressions for the time-dependent coefficient case were derived in general form. These forms were later particularized to the constant coefficient case where it was shown that under mild assumptions about coefficient matrix Cthe in- tegrating factor can be found in closed form. Analytical particular solutions were derived for different forms of continuous and discontinuous excitations. Complicate expressions for loading can be generated by using a linear combination of polynomials multiplied by Heavisides and also by combining the analytical solutions derived in this manuscript. Examples showed that the proposed approach is accurate and can be made efficient, not suffering with common issues found in traditional numerical approaches, like stabil- ity associated to time discretization and interpolation errors. Initial conditions can be imposed at any given time t0, not only in the extremes of the interval. Actually, no time span is needed to evaluate the response when using the proposed approach. An optimized way to calculate the matrix exponential in a discrete set of points was proposed to greatly reduce computational costs, which benefits the evaluation of the homogeneous solution and the solutions due to Dirac’s delta and Heaviside series. Thereby, important real-world problems can be analytically solved by using the proposed approach. Declaration of Competing Interest Theauthorsdeclarethattheyhavenoknowncompetingfinancialinterestsorpersonal relationships that could have appeared to influence the work reported in this paper. 34 Data availability The complete computer implementation of the theory discussed in this manuscript can be found at https://github.com/CodeLenz/Giffndof.jl with the scripts used to gen- erate the images shown in this manuscript. The interested reader can freely change the parameters to study the proposed formulation. 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 also thankful for the kind replies of Professor Nicholas Higham regarding exponential maps and Dr. Olavo M. da Silva Jr. for carefully reading the manuscript. Appendix A. EvaluationofintegrationconstantsforcoupledsystemsofODEs with constant coefficients The homogeneous solution with constant coefficients and when ¯CandF2,1,1commute is given by Eq. 49. The complete solution y=yh+yp= exp −F2,1,1t C2+ exp F2,1,1−¯C t C1+yp,(A.1) and its derivative w.r.t. time is given by ˙ y=−F2,1,1exp −F2,1,1t C2+ F2,1,1−¯C exp F2,1,1−¯C t C1+˙ yp.(A.2) Considering initial conditions at a time t0yields y(t0) = exp −F2,1,1t0 C2+ exp F2,1,1−¯C t0 C1+yp(t0),(A.3) and ˙ y(t0) =−F2,1,1exp −F2,1,1t0 C2+ F2,1,1−¯C exp F2,1,1−¯C t0 C1+˙ yp(t0) (A.4) which can be summarized in a linear system  exp F2,1,1−¯C t0 exp −F2,1,1t0  F2,1,1−¯C exp F2,1,1−¯C t0 −F2,1,1exp −F2,1,1t0 C1 C2 =u0−yp(t0) v0−˙yp(t0) .(A.5) Hence, C1andC2can be found by solving the above linear system with standard techniques, muchaliketheevaluationof C1andC2forasingledegreeoffreedomproblem. The most common choice for t0is0such that Eq. A.5 reduces to 35 " I I F2,1,1−¯C −F2,1,1#C1 C2 =u0−yp(0) v0−˙yp(0) . (A.6) Pivoting, it further simplifies to  I I 0¯C−2F2,1,1 C1 C2 =( u0−yp(0) v0−˙yp(0)− F2,1,1−¯C u0−yp(0)) .(A.7) thus, it follows that C2=¯C−2F2,1,1−1 v0−˙yp(0)− F2,1,1−¯C u0−yp(0) ,(A.8) and C1=u0−yp(0)−C2. (A.9) A very important particularization for previous equations is when evaluating the constants C1andC2forexcitationsdescribedbyDirac’sdeltas, HeavisidesandHeaviside series. As discussed in [3], yp(t)and˙y(t)pare zero for t≤tH, where tHis the first time with non null excitation. Therefore, if t0≤tH, both yp(t0)and˙yp(t0)are0in Eqs. A.5, A.8 and A.9 and there is no need to evaluate the derivative of the particular response with respect to time. Appendix B. Convolution over Dirac’s delta distribution The convolution of a function over a Dirac’s delta is usually defined as [11] Z∞ −∞f(t)δ(t−tk)dt=f(tk). (B.1) 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, (B.2) 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. (B.3) Thefilterorsiftingproperty of the delta of Dirac is due to the shape of this distri- bution,i.e., it is null almost everywhere except at its discontinuity. As a consequence, the function multiplying the Dirac’s delta is constant at this point, for the discontinuity of the delta distribution is infinitely close to tk. 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 36 Z∞ −∞f(t)δ(t−tk)dt=Ztk+τ tk−τf(t)δ(t−tk)dt=f(tk)Ztk+τ tk−τδ(t−tk)dt=f(tk), (B.4) hence, Equation B.3 can be rewritten as Zt 0f(t)δ(t−tk)dt=f(tk)−f(tk)( 0, t≥tk 1, t < t k=f(tk)H(t−tk).(B.5) Appendix C. Commutativity of powers of two matrices LetAandBbe two square matrices n×nthat commute. For positive integer powers aandb, the following property holds AaBb=BbAa. (C.1) Proof. AaBb=AAA . . .AAA|{z } a timesBBB . . .BBB|{z } b times =AAA . . .AABABB . . .BBB =AAA . . .ABABAB . . .BBB =AAA . . .BABABA . . .BBB =AAA . . .ABABAB . . .BBB ... BBB . . .BBB|{z } b timesAAA . . .AAA|{z } a times=BbAa. Appendix D. Commutativity of inverse of matrix LetAandBbe two n×nsquare matrices that commute and Bbe invertible. Then AB−1=B−1A. (D.1) Proof. AB−1=D, A=DB. BA=BDB =AB, BD=A, D=B−1A, =⇒AB−1=B−1A. 37 Appendix E. Exponential map An exponential map is defined as [7], exp (At) =I+At+1 2A2t2+...+1 n!Antn=∞X j=01 j!Ajtj, (E.1) with time-derivative given by d dt exp (At) =A+A2t+...+1 (n−1)!Antn−1=∞X j=01 j!Aj+1tj.(E.2) By direct comparison of Eq. E.1 to Eq. E.2, one immediately gets d dt exp (At) =Aexp (At) = exp ( At)A. (E.3) Appendix F. Efficient evaluation of exp (A t) One critical component of the proposed formulation is the efficient evaluation of exp (At)for various values of t. According to [1], the following relation holds, exp (A)α= exp ( αA)∀α∈Z, (F.1) for general matrix A Assuming that the time span t∈[ti, tf]is discretized in ntintervals ∆t, it is possible to write t=k∆tsuch that E(t) = exp ( At) = exp ( k∆tA) = exp (∆ tA)k, (F.2) where k∈Z. Hence, the exponential map, E(t), in a discrete set of time can be evaluated through the following recursion t1:E(t1) = exp (∆ tA) t2:E(t2) =E(t1) exp (∆ tA) t3:E(t3) =E(t2) exp (∆ tA) ... tk:E(tk) =E(tk−1) exp (∆ tA), (F.3) where matrix exp (∆ tA)has to evaluated just once. 38 Appendix G. Commutativity of exponential map and matrix LetAandBbe two n×nsquare matrices. If AandBcommute, then exp (A)B=Bexp (A). (G.1) Proof.By definition exp (A)B= I+At+1 2A2t2+...+1 n!Antn B =IB+ABt+1 2A2Bt2+...+1 n!AnBtn (G.2) such that exp (A)B=BI+BAt+1 2BA2t2+...+1 n!BAntn=Bexp (A).