Dynamic Response Analysis of Structures Using Legendre–Galerkin Matrix Method

: The solution of the motion equation for a structural system under prescribed loading and the prediction of the induced accelerations, velocities, and displacements is of special importance in structural engineering applications. In most cases, however, it is impossible to propose an exact analytical solution, as in the case of systems subjected to stochastic input motions or forces. This is also the case of non-linear systems, where numerical approaches shall be taken into account to handle the governing differential equations. The Legendre–Galerkin matrix (LGM) method, in this regard, is one of the basic approaches to solving systems of differential equations. As a spectral method, it estimates the system response as a set of polynomials. Using Legendre’s orthogonal basis and considering Galerkin’s method, this approach transforms the governing differential equation of a system into algebraic polynomials and then solves the acquired equations which eventually yield the problem solution. In this paper, the LGM method is used to solve the motion equations of single-degree (SDOF) and multi-degree-of-freedom (MDOF) structural systems. The obtained outputs are compared with methods of exact solution (when available), or with the numerical step-by-step linear Newmark-β method. The presented results show that the LGM method offers outstanding accuracy.


Introduction and State-of-Art
Most structural systems in civil engineering applications, as known, are either discrete or can be estimated as discrete equivalents.The dynamic equation of motion of a continuous system can be hence handled as a discrete system.The advantage of this assumption is that the solution of fundamental equations of motion of discrete systems-which is of utmost importance for engineering applications-can be efficiently calculated to predict the responses of structural systems.Among others, let us consider the following well-known second-order differential equation with given initial conditions: .. y(t) = f (t, y(t), .y(t)), t 0 < t ≤ T, y(t 0 ) = y 0 , .y(t 0 ) = .y 0 . (1) Equation ( 1) is commonly used to express the governing equation of mechanical vibrations, quantum dynamic calculations, dynamic equilibrium of structures, etc. Employing and developing new numerical methods for approaching the exact solution is of great importance.
In structural dynamics, algorithms of direct time integration are usually used to obtain the solution of discrete temporal equations of motion at selected time steps [1].In the past, to this aim, different integration algorithms in the time domain have been developed using various methods.The same algorithms are widely used to solve the equation of systems under dynamic loading.Several well-known algorithms have been presented by researchers, among which the Newmark-β [2], Wilson-Theta [3], Runge-Kutta [4], etc., can be pointed out.Detailed explanations can be found in structural dynamics textbooks [1,5,6].
Alternatively, Equation (1) can also be solved through spectral methods of discretization [34][35][36], which are strategic for the numerical solution of differential equations.The main advantage of these methods lies in their accuracy for a given number of unknowns.For smooth problems with simple geometries, these methods offer exponential rates of convergence/spectral accuracy.In contrast, methods such as finite difference and finite element yield only algebraic convergence rates.Three spectral methods, namely the Galerkin [37], Collocation [38], and Tau [39] are extensively used in the literature.Spectral methods are preferable in numerical solutions of partial differential equations due to their high-order accuracy [40,41].The Standard Spectral and Galerkin methods have been extensively investigated to handle different types of problems [42][43][44][45].Several numerical methods have been also developed to solve different types of differential equations.Some of these methods can be used to solve more practical equations.These methods include sparse multiscale representation of the Galerkin method [46], spectral collocation method [47], Jacob spectral method [48], and many others.Recently, Erfanian et al. [49] developed a new method for solving two-dimensional nonlinear Volterra integral equations, based on the use of rationalized Haar functions in the complex plane.Bernoulli Galerkin matrix method has been also proposed by Hesameddini and Riahi [50] for solving the system of Volterra-Fredholm integro-differential equations.A new hybrid orthonormal Bernstein and improved block-pulse functions method has been studied by Ramadan and Osheba [51] for solving mathematical physics and engineering problems.
In recent years, researchers have developed a number of numerical methods to find the solution for such equations.In this regard, matrix methods including the Euler method [52], Bernoulli collocation method [53], hybrid Legendre block-pulse function method [54], least squares method [55], Bessel colocation method [56], etc., have received much attention since 2010.Given that the governing equations of natural phenomena are usually nonlinear and complex, so the chosen method of solving must be commensurate with the problem's complexity and its dimensions.A nonlinear differential equation system consists of several nonlinear equations that solving such a system requires numerical methods.The basis of matrix methods is the expression of each of the functions in the problem based on selected polynomials.After selecting the type of polynomials and expressing each of the functions, the differential equation system becomes a system of linear equations with several unknown coefficients, which can be solved by Galerkin and collocation methods.
In this paper, the Legendre-Galerkin matrix (LGM) method is used to solve the equations of motion of different SDOF and multi-degree-of-freedom (MDOF) structural systems.The results are compared with those from the exact solution (when available), or from the numerical step-by-step Newmark-β method with linear acceleration.The accuracy of the LGM formulation is thus highlighted in the discussion of comparative calculations.

Basics for SDOF and MDOF Systems
In the dynamic analysis of structures, from a practical point of view, a real structure is defined as a system with infinite degrees of freedom that can be modeled as a discrete system with SDOF or MDOF in a finite element approach.In many cases, these simple models include complex information of the real structure and are able to simulate the behavior of real structure with a good level of accuracy and high calculation efficiency.As an example, in blast-resistant design of structures, the basis of current books and codes and also simplified engineering tools have been established based on equivalent SDOF models of real structures [57][58][59].

Governing Equation for SDOF Systems
A structure with one degree of freedom with mass m, stiffness k, and damping c is assumed to be under dynamic load P(t) = P 0 sin(Ωt) with excitation frequency ς (Figure 1).The governing equation of this system is [1]: where Equation ( 2) is written with similar connotations to Equation (1).The exact solution to the motion equation of a SDOF system under harmonic load P 0 sin(Ωt) with initial conditions u(0) = u 0 and .u(0) = .
u 0 is the sum of the displacements of the transient state of the system (u transient ) and the steady state (u steady state ), and is expressed as [5]: in which u transient and u steady state are defined as: and where ξ and ω = k m are the damping ratio and natural frequency of the system, respectively, and ω d = ω 1 − ξ 2 indicates the damped frequency.Finally, β = Ω ω is the ratio of the excitation frequency to the system natural frequency.

Linear Newmark-β Method
One of the most efficient methods of handling equations of motion of SDOF and MDOF systems has been proposed by Newmark [2].In 1959, Newmark also presented a set of step-by-step numerical equations that are entirely known as the Newmark-β method [5].
In the present study, the linear Newmark-β method is taken into account from [5] to verify the results of the newly developed LGM method.
In order to solve the equation of motion of a MDOF system using the linear Newmarkβ method, it is important to remind that the following steps should be generally followed: 1.
Determination of the initial conditions, {u 0 }, .u 0 and .. u 0 , where the value of .. u 0 is given by: ..
Calculation of a and b factors: Calculation of the effective load: Calculation of displacement, velocity, and acceleration: where, .

Legendre-Galerkin Matrix Method
Legendre polynomials are among the most important functions in the function approximation theory.The LGM method is highly successful in approximating different functions which mainly stems from the orthogonality of Legendre basis components.This orthogonality makes it easy to find the unknown coefficients of the problem.Another reason for using these basis components is their weight.The weight function will not be a problem for the calculation of the integrals of the Galerkin method and hence the operational matrices of differentiations and other existing functions in the problem can be easily found.
The Galerkin method is an efficient and easy-to-implement approach for solution of the equations of motion of MDOF systems: under the following initial conditions: In Equation ( 14), m ij , c ij and k ij are, respectively, components of mass, damping and stiffness matrices.p i (t) is the load prescribed at the i-th degree of freedom and λ i and γ i are, respectively, the initial displacement and velocity applied to the i-th degree of freedom.y i is the displacement of the i-th degree of freedom which is the unknown of the problem, and t represents time.

•
Legendre polynomials: introduced by L m (t), the Legendre polynomials are defined in the interval [−1, 1] and can be obtained using the following recursive equations: For further application of these polynomials, they are shifted to the interval [0, L] with the change of variable 2 L t − 1.Therefore, these shifted polynomials which are shown by L * m (t) are represented as the following: • Inner product of two functions: the inner product of two functions f(t) and g(t) is represented by f (t), g(t) .If the functions are known and continuous in the interval [0, b], then: • Orthogonality of two functions: two functions f (t) and g(t) are orthogonal with respect to the weight function w(t) on [a, b] if: The shifted Legendre polynomials are orthogonal with respect to the weight function Orthogonality of the functions has wide application in the function approximation theory.Any continuous function such as f (x) can be approximated in the interval [0, L] using these polynomials as below: where f m can be obtained from: Only the first N + 1 terms of Equation ( 21) are used in practice.Therefore: where F and φ N (t) are factors of the shifted Legendre vectors expressed as: φ

Expression of the LGM Method
Assuming that Equation ( 14) has a unique solution under the initial conditions in Equation ( 15), the goal is to find an approximate analytical solution to Equation ( 14) using the discretized Legendre series as below: where, The matrix form of the derivative vector will be: where D is a matrix of dimensions (N + 1) × (N + 1) and is obtained as below: For example, for different values of N, the following equations are obtained: Using Equation ( 27), the k-th derivative of φ N (t) vector is: Therefore, by replacing k = 1 and 2 in Equation ( 29), and then combining the results with Equation ( 25), one can write: . .
Moreover, the functions p i (t) can be shown in matrix form as follows: where P T i,N is a (N+1) dimensional vector with the following expression: where, By introducing Equations ( 25) and ( 30)- (32) in Equation ( 14), it is thus possible to obtain: Now, let us assume that: where Y N is an f dimensional vector, A and P are f × (N + 1) dimensional and M, C and K are the known f × f dimensional matrices, respectively.Using Equations ( 35) and ( 36), the matrix form for approximation of Equation ( 14) can be obtained as: or where, U = MAD 2 + CAD + KA.
Equation ( 38) can be solved using different methods, namely the direct method (i.e., omitting φ N (t) from both sides of the equation and assuming it to vanish), the collocation method, or the Galerkin method.
In the present paper, the Galerkin method is selected to solve the equation.This method tries to minimize the error toward zero with the inner product of the equations in Legendre basis L * m (t).Consequently: The above equation is an algebraic system of linear equations with N + 1 equations and N + 1 unknowns that can be easily solved.It is needed, however, to apply the boundary conditions to the problem, as also in accordance with Equation (15).Following Equations ( 25) and (30), it is: where, Finally, in order to find the response of the system in Equation ( 14) under the initial conditions from Equation (15), by replacing 2f of the rows of Equation (40) with Equation (38), a system of algebraic equations with (N + 1 − 2f ) equations (Equation ( 38)) and 2f equations (Equation ( 40)) is created.Their solution yields the analytic approximation of the original problem.The last rows of Equation ( 38) are usually replaced with those equations from Equation (40); however, this is not always obligatory and the rows which will cause the system of equations to be singular can be alternatively replaced.
By replacing the obtained values in Equation ( 17), the approximate solution for y(t) is: Figure 2 illustrates the results of the exact solution as a function of time, along with its approximation.The accuracy of the LGM method compared to the exact solution can be noticed in the whole time interval.

Worked Examples and Discussion of Results
For a more exhaustive discussion of the developed LGM method and for a reliable assessment of its potential in structural analysis, the motion equation of an MDOF system is solved in accordance with Equation (14), and some numerical examples are presented.These examples are representative of simple models of real structures (by introducing mass, stiffness, and damping matrices) with and without damping.Accordingly, the current study and the presented formulation of the LGM method represents a first step towards its further extension and use for other applications, such as reliability analysis, optimization, and modal analysis of structures.

Two-Degree-of-Freedom (2DOF) Structure
A structure with two degrees of freedom (2DOF) is analyzed in different dynamic conditions.The structure is first analyzed under free vibration with no effective damping, and then with additional damping.Successively, the analysis is further performed with the assumption of forced vibration, with and without damping.Basic input data and parameters are summarized as follows:
• Free vibration, with damping: it is C = 1.628 −0.256 −0.256 2.512 , while all the other parameters are equal to the undamped case.
The solution of the motion equation using Legendre's method for u 1 (t) and u 2 (t) are respectively shown in Figures 3 and 4.Moreover, the LGM results are compared with the exact method, giving evidence of acceptable consistency.It is worth noting that the exact solution is obtained using the Solver of Systems of Equation in the Maple software.Comparisons are then made with the approximate solution.Successively, the examined structure is analyzed under forced vibration.In this case, it is: Once again, the structure is analyzed with and without damping, using the developed LGM method.The analytical results obtained from the LGM method for u 1 (t) and u 2 (t) are shown in Figures 5 and 6, respectively, and compared with the exact solution, proving again a rather acceptable consistency in the examined time interval.The results for u 1 (t), u 2 (t) and u 3 (t) are shown in Figures 7-9, respectively, where it is possible to see the comparison of the LGM method and linear Newmark-β method, with rather good consistency.
The equation of motion of the structure under forced vibration (with and without The results obtained in this case for u 1 (t), u 2 (t) and u 3 (t) are proposed in Figures 10-12, respectively, and compared with the results by the linear Newmark-β method.Moreover, in this case, the comparative data show the consistency of the LGM with the linear Newmark-β method.

Five-Degree-of-Freedom Structure
In conclusion, a structure with five degrees of freedom (5DOF) is analyzed with the LGM method, both under free vibration with and without damping.

•
Free vibration, without damping: it is assumed that M = The solutions of the equation of motion of the structure using the LGM for damped and undamped states, in terms of u 1 (t) to u 5 (t), are shown in Figures 13-17 and compared with those by the linear Newmark-β method.As for the previously discussed calculation examples, the comparative plots prove an appropriate consistency and a high level of accuracy of the LGM procedure.It should be noted that, in the case of stochastic analysis that is interpreted classically as a repetition of solution for different input values with defined probability density functions (taking into account uncertainties associated with input parameters), the methods such as Monte Carlo simulation (MCS) [60] can be easily used along with the results of LGM method to investigate the problem in a probabilistic manner.This is due to the fact that, in the LGM method, the solution is approximated by discretized Legendre series that can be used as a state function in reliability analysis where no mathematical closed-form state function can be found (e.g., through finite element method).Furthermore, such a solution obtained from the LGM method can also be utilized with analytical reliability methods (e.g., jointly distributed random variables method [61]) or approximate ones (e.g., first and second-order reliability methods (FORM and SORM), point estimate method (PEM), etc. [62]), with less computational efforts in comparison to the MCS method.

Conclusions
The Legendre-Galerkin matrix (LGM) method was developed in this study to solve systems of differential equations of motion.As shown, the advantage of this spectral method is that it converts the governing differential equation of a given problem to a system of algebraic equations, based on a set of orthogonal Legendre polynomials.The final solution leads to a good estimate of the solution for a system of differential equations.
In the present research study, the selected differential equations were typical of single degree (SDOF) and multi-degree-of-freedom (MDOF) structural systems.
In order to prove the accuracy of the proposed method in the response calculation of SDOF and MDOF structures, a number of numerical examples for damped and undamped structural systems under free or forced vibrations were developed and discussed.When available, exact solutions were taken into account for the comparative analysis of LGM predictions, otherwise, the results of the numerical linear Newmark-β method were used to verify the estimates from the developed LGM method.The overall comparative data showed that the LGM method is of high accuracy in estimating the response of SDOF and MDOF systems and can be thus effectively employed in the solution of fundamental motion equations of structures.

Figure 2 .
Figure 2. Comparison of LGM and exact solutions for an SDOF system.

Figure 3 .
Figure 3. Free vibration analysis for a 2DOF system.Comparison of u 1 (t) as a function of time, as obtained from the LGM method or the exact solution.

Figure 4 .
Figure 4. Free vibration analysis for a 2DOF system.Comparison of u 2 (t) as a function of time, as obtained from the LGM method or the exact solution.

Figure 5 .
Figure 5. Forced vibration analysis of a 2DOF system.Comparison of u 1 (t) as a function of time, as obtained from the LGM method or the linear Newmark-β method.

Figure 6 . 1  6 
Figure 6.Forced vibration analysis of a 2DOF system.Comparison of u 2 (t) as a function of time, as obtained from the LGM method or the linear Newmark-β method.4.2.Three-Degree-of-Freedom StructureAs a further validation example, a structure with three degrees of freedom (3DOF) is analyzed in forced and free vibration conditions, with or without damping.


is solved further by using Legendre's method.

Figure 7 .
Figure 7. Free vibration analysis of a 3DOF system.Comparison of u 1 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 8 .
Figure 8. Free vibration analysis of a 3DOF system.Comparison of u 2 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 9 .
Figure 9. Free vibration analysis of a 3DOF system.Comparison of u 3 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 10 .
Figure10.Forced vibration analysis of a 3DOF system.Comparison of u 1 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 11 .
Figure 11.Forced vibration analysis of a 3DOF system.Comparison of u 2 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 12 .
Figure 12.Forced vibration analysis of a 3DOF system.Comparison of u 3 (t) as function of time, as obtained using the LGM method or the linear Newmark-β method.

•
Free vibration, with damping: it is assumed that C = parameters are the same as in the undamped case.

Figure 13 .
Figure13.Free vibration analysis of a 5DOF system.Comparison of u 1 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 14 .
Figure 14.Free vibration analysis of a 5DOF system.Comparison of u 2 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 15 .
Figure 15.Free vibration analysis of a 5DOF system.Comparison of u 3 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 16 .
Figure 16.Free vibration analysis of a 5DOF system.Comparison of u 4 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.

Figure 17 .
Figure 17.Free vibration analysis of a 5DOF system.Comparison of u 5 (t) as a function of time, as obtained using the LGM method or the linear Newmark-β method.