**1. 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:

$$\begin{cases} \quad \ddot{y}(t) = f(t, y(t), \dot{y}(t)), \qquad t\_0 < t \le T, \\\ y(t\_0) = y\_{0'} \quad \dot{y}(t\_0) = \dot{y}\_0. \end{cases} \tag{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

**Citation:** Momeni, M.; Riahi Beni, M.; Bedon, C.; Najafgholipour, M.A.; Dehghan, S.M.; JavidSharifi, B.; Hadianfard, M.A. Dynamic Response Analysis of Structures Using Legendre–Galerkin Matrix Method. *Appl. Sci.* **2021**, *11*, 9307. https:// doi.org/10.3390/app11199307

Academic Editor: José A. F. O. Correia

Received: 1 August 2021 Accepted: 29 September 2021 Published: 7 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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].

Since 1994, the Chebyshev polynomials [7], Legendre polynomials [8], Bessel polynomials [9], Hermite polynomials [10], Laguerre polynomials [11], and matrix methods have been used in many research studies to solve linear and nonlinear equations with high orders including partial differential equations, hyperbolic partial differential equations, delay equations, integral and integro-differential equations, SDOF and MDOF systems, etc. One of the approaches to find the solution to an initial value problem is taking semi-analytical procedures such as the differential transform method (DTM) [12–23]. The nature of dynamic equations of motion of SDOF systems makes them differential equations with initial values; hence, DTM has been also applied to solve non-linear SDOF problems [24–33].

Alternatively, Equation (1) can also be solved through spectral methods of discretization [34–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–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.

#### **2. 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–59].

## *2.1. 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*<sup>0</sup> sin(Ω*t*) with excitation frequency *ς* (Figure 1).

**Figure 1.** SDOF structure subjected to dynamic load P(t).

The governing equation of this system is [1]:

$$m\ddot{u}(t) + c\dot{u}(t) + ku(t) = P\_0 \sin(\Omega t) \tag{2}$$

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*<sup>0</sup> sin(Ω*t*) with initial conditions *<sup>u</sup>*(0) = *<sup>u</sup>*<sup>0</sup> and . *<sup>u</sup>*(0) = . *u*<sup>0</sup> is the sum of the displacements of the transient state of the system (*utransient*) and the steady state (*usteady state*), and is expressed as [5]:

$$
\mu\_{total} = \mu\_{transient} + \mu\_{steady\,state} \tag{3}
$$

in which *utransient* and *usteady state* are defined as:

$$\begin{split} u\_{\text{transient}} &= e^{-\xi\omega t} \Big( u\_0 \cos\omega\_d t + \frac{\dot{u}\_0 + u\_0\omega\_s^x}{\omega\_d} \sin\omega\_d t \Big) + \\ \frac{p\_0}{k} &\frac{e^{-\xi\omega t}}{\left(1 - \beta^2\right)^2 + \left(2\xi\beta\right)^2} \Big[ 2\xi\beta\cos\omega\_d t + \frac{\omega}{\omega\_d} \left\{ 2\beta\xi^2 - \beta(1 - \beta^2) \right\} \sin\omega\_d t \Big], \end{split} \tag{4}$$

and

$$u\_{\text{steady-state}} = \frac{p\_0}{k} \frac{1}{\left(1 - \beta^2\right)^2 + \left(2\xi\beta\right)^2} \left\{ \left(1 - \beta^2\right) \sin\Omega t - 2\xi\beta\cos\Omega t \right\} \tag{5}$$

where *ξ* and *ω* = q *k m* are the damping ratio and natural frequency of the system, respectively, and *ω<sup>d</sup>* = *ω* p 1 − *ξ* 2 indicates the damped frequency. Finally, *β* = <sup>Ω</sup> *ω* is the ratio of the excitation frequency to the system natural frequency.

#### *2.2. 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:

$$
\ddot{u}\_0 = \frac{p\_0 - \mathbb{C}\dot{u}\_0 - \mathbb{K}u\_0}{M} \tag{6}
$$


$$\hat{k} = k + \frac{\gamma}{\beta \Delta t} \mathcal{C} + \frac{1}{\beta (\Delta t)^2} \mathcal{M} \tag{7}$$

where *β* = 1/6 and *γ* = 1/2 for linear acceleration method. 4. Calculation of a and b factors:

$$\begin{cases} \ a = \frac{1}{\beta \Delta t} M + \frac{\gamma}{\beta} \mathcal{C} \\\ b = \frac{1}{2\beta} M + \Delta t \left(\frac{\gamma}{2\beta} - 1\right) \mathcal{C} \end{cases} \tag{8}$$

5. Calculation of the effective load:

$$
\Delta \mathfrak{P}\_i = \Delta p\_i + a\dot{u}\_i + b\ddot{u}\_i \tag{9}
$$

6. Calculation of displacement, velocity, and acceleration:

$$\begin{cases} \begin{aligned} u\_{i+1} &= u\_i + \Delta u\_{i\nu} \\ \dot{u}\_{i+1} &= \dot{u}\_i + \Delta \dot{u}\_{i\nu} \\ \ddot{u}\_{i+1} &= \ddot{u}\_i + \Delta \ddot{u}\_{i\nu} \end{aligned} \end{cases} \tag{10}$$

where,

$$
\Delta u\_i = \frac{\Delta \mathfrak{p}\_i}{\mathfrak{K}} \tag{11}
$$

$$
\Delta \dot{u}\_{\dot{l}} = \frac{\gamma}{\beta \Delta t} \Delta u\_{\dot{l}} - \frac{\gamma}{\beta} \dot{u}\_{\dot{l}} + \Delta t \left(1 - \frac{\gamma}{2\beta}\right) \ddot{u}\_{\dot{l}} \tag{12}
$$

$$
\Delta \ddot{u}\_i = \frac{\gamma}{\beta \left(\Delta t\right)^2} \Delta u\_i - \frac{1}{\beta \Delta t} \dot{u}\_i - \frac{1}{2\beta} \ddot{u}\_i. \tag{13}
$$

#### **3. 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:

$$\sum\_{j=1}^{f} \left\{ m\_{i\bar{j}} \ddot{y}\_j(t) + c\_{i\bar{j}} \dot{y}\_j(t) + k\_{i\bar{j}} y\_j(t) \right\} = p\_i(t), \qquad i = 1, 2, \dots, f \qquad t \in [0, t\_1], \tag{14}$$

under the following initial conditions:

$$\begin{cases} \ y\_i(0) = \lambda\_{i\nu} \\ \ y\_i'(0) = \gamma\_{i\nu} \end{cases} \text{ for } \ i = 1, 2, \dots, f. \tag{15}$$

In Equation (14), *mij*, *cij* and *kij* are, respectively, components of mass, damping and stiffness matrices. *pi*(*t*) is the load prescribed at the *i*-th degree of freedom and *λ<sup>i</sup>* and *γ<sup>i</sup>* are, respectively, the initial displacement and velocity applied to the *i*-th degree of freedom. *yi* is the displacement of the *i*-th degree of freedom which is the unknown of the problem, and *t* represents time.

#### *3.1. Approximation of the Function Using Shifted Legendre Polynomials*

• Legendre polynomials: introduced by *Lm*(*t*), the Legendre polynomials are defined in the interval [−1, 1] and can be obtained using the following recursive equations:

$$\begin{array}{l} L\_0(t) = 1, \\ L\_1(t) = t, \\ L\_2(t) = \frac{3}{2}t^2 - \frac{1}{2}, \\ \vdots \\ L\_{m+1}(t) = \frac{2m+1}{m+1}t \, L\_m(t) - \frac{m}{m+1}L\_{m-1}(t). \end{array} \tag{16}$$

For further application of these polynomials, they are shifted to the interval [0, *L*] with the change of variable <sup>2</sup> *L t* − 1. Therefore, these shifted polynomials which are shown by *L* ∗ *<sup>m</sup>*(*t*) are represented as the following:

$$L\_m^\*(t) = L\_m(\frac{2}{L}t - 1) \qquad t \in [0, L]. \tag{17}$$

• Inner product of two functions: the inner product of two functions *f(t)* and *g(t)* is represented by h*f*(*t*), *g*(*t*)i. If the functions are known and continuous in the interval [0, *b*], then:

$$
\langle f(t), g(t) \rangle = \int\_{a}^{b} f(t)g(t)w(t)dt. \tag{18}
$$

• 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:

$$
\langle f(t), g(t) \rangle = 0. \tag{19}
$$

The shifted Legendre polynomials are orthogonal with respect to the weight function *w*(*t*) = 1 in the interval [0, *L*]. This means that:

$$\int\_{a}^{b} L\_{m}^{\*}(t)L\_{n}^{\*}(t)dt = \begin{cases} 0 & m \neq n, \\ \frac{L}{2m+1} & m = n. \end{cases} \tag{20}$$

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:

$$f(t) = \sum\_{m=0}^{\infty} f\_m L\_m^\*(t),\tag{21}$$

where *f<sup>m</sup>* can be obtained from:

$$f\_m = \frac{\langle f(t), L\_m^\*(t) \rangle}{\langle L\_m^\*(t), L\_m^\*(t) \rangle} = \frac{2m+1}{L} \int\_0^L f(t) L\_m^\*(t) dt. \tag{22}$$

Only the first *N +* 1 terms of Equation (21) are used in practice. Therefore:

$$f(t) = \sum\_{m=0}^{N} f\_m L\_m^\*(t) = \mathbf{F}^T \Phi\_N(t),\tag{23}$$

where **F** and φ*N*(t) are factors of the shifted Legendre vectors expressed as:

$$\mathbf{F} = \begin{bmatrix} f\_0 \ f\_1 \ \dots \ f\_N \end{bmatrix}^T \tag{24}$$

$$\Phi\_N(t) = \begin{bmatrix} L\_0^\*(t) \ L\_1^\*(t) \ \dots \ L\_N^\*(t) \end{bmatrix}^T. \tag{25}$$

## *3.2. 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:

$$\Phi\_{i,N}(t) = \sum\_{m=0}^{N} a\_{i,m} L\_m^\*(t) = \mathbf{A}\_{i,N}^T \Phi\_N(t) \tag{26}$$

where,

$$\begin{array}{c} \mathbf{A}\_{i,N} = [a\_{i,0} \ a\_{i,1} \ \dots \ a\_{i,N}]^T, \\ \boldsymbol{\Phi}\_N(t) = [L\_0^\*(t) \ L\_1^\*(t) \ \dots \ L\_N^\*(t)]. \end{array} \tag{27}$$

The matrix form of the derivative vector will be:

$$\frac{d}{dt}\Phi\_N(t) = \mathbf{D}\Phi\_N(t),\tag{28}$$

(29)

where **D** is a matrix of dimensions (*N* + 1) × (*N* + 1) and is obtained as below:

$$\mathbf{D} = [d\_{ij}] = \begin{cases} \frac{2(2j+1)}{L}, & \text{for } \ j = i - k, \\\ 0, & \text{otherwise.} \end{cases} \text{ if } \mathbf{N} \text{ is odd number,} \\ \begin{cases} k = 1, 3, \dots, N - 1 & \text{if } \mathbf{N} \text{ is even number,} \\\ 0, & \text{otherwise.} \end{cases}$$

For example, for different values of N, the following equations are obtained:

*I f N* = 1 → **D** = 0 0 2 *L* 0 *I f N* = 2 → **D** = 0 0 0 2 *L* 0 0 0 6 *L* 0 *I f N* = 3 → **D** = 0 0 0 0 2 *L* 0 0 0 0 6 *L* 0 0 2 *L* 0 10 *L* 0 

Using Equation (27), the *k*-th derivative of φ*N*(*t*) vector is:

$$\frac{d^k}{dt^k} \Phi\_N(t) = \mathbf{D}^k \Phi\_N(t) \, \tag{30}$$

Therefore, by replacing *k* = 1 and 2 in Equation (29), and then combining the results with Equation (25), one can write:

$$\dot{y}\_{i,N}(t) = \sum\_{m=0}^{N} a\_{i,m} L\_m^\*(t) = \mathbf{A}\_{i,N}^T \mathbf{D} \boldsymbol{\Phi}\_N(t),\tag{31}$$

and

$$\ddot{y}\_{i,N}(t) = \sum\_{m=0}^{N} a\_{i,m} L\_m^\*(t) = \mathbf{A}\_{i,N}^T \mathbf{D}^2 \boldsymbol{\Phi}\_N(t). \tag{32}$$

Moreover, the functions *pi*(*t*) can be shown in matrix form as follows:

$$p\_i(t) = \sum\_{m=0}^{N} p\_{i,m} L\_m^\*(t) = \mathbf{P}\_{i,N}^T \boldsymbol{\Phi}\_N(t),\tag{33}$$

where *P T i*,*N* is a (*N*+1) dimensional vector with the following expression:

$$\mathbf{P}\_{i,N} = \begin{bmatrix} p\_{i,0} \ p\_{i,1} \ \dots \ p\_{i,N} \end{bmatrix}^T \tag{34}$$

where,

$$p\_{i,m} = \frac{2m+1}{L} \int\_0^L p\_i(t)L\_m^\*(t)dt. \tag{35}$$

By introducing Equations (25) and (30)–(32) in Equation (14), it is thus possible to obtain:

$$\sum\_{j=1}^{f} \left\{ m\_{lj} \mathbf{A}\_{j,N}^T \mathbf{D}^2 + c\_{lj} \mathbf{A}\_{j,N}^T \mathbf{D} + k\_{lj} \mathbf{A}\_{j,N}^T \right\} \boldsymbol{\Phi}\_N(t) = \mathbf{P}\_{i,N}^T \boldsymbol{\Phi}\_N(t) \tag{36}$$

Now, let us assume that:

$$\mathbf{Y}\_{N} = \mathbf{A}\boldsymbol{\Phi}\_{N}(t) = \begin{bmatrix} \mathcal{Y}\_{1,N} \\ \mathcal{Y}\_{2,N} \\ \vdots \\ \mathcal{Y}\_{f,N} \end{bmatrix}, \quad \mathbf{A} = \begin{bmatrix} \mathbf{A}\_{1,N}^{T} \\ \mathbf{A}\_{2,N}^{T} \\ \vdots \\ \mathbf{A}\_{f,N}^{T} \end{bmatrix}, \quad \mathbf{P} = \begin{bmatrix} \mathbf{P}\_{1,N}^{T} \\ \mathbf{P}\_{2,N}^{T} \\ \vdots \\ \mathbf{P}\_{f,N}^{T} \end{bmatrix}, \tag{37}$$

$$\mathbf{M} = [m\_{ij}]\_\prime \cdot \mathbf{C} = [c\_{ij}]\_\prime \cdot \mathbf{K} = [k\_{ij}]\_\prime \quad i, j = 1, 2, \dots, f, \mathbf{c}$$

where **Y***<sup>N</sup>* 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:

$$(\mathbf{M}\mathbf{A}\mathbf{D}^2 + \mathbf{C}\mathbf{A}\mathbf{D} + \mathbf{K}\mathbf{A})\phi\_N(t) = \mathbf{P}\phi\_N(t),\tag{38}$$

or

$$\mathbf{U}\boldsymbol{\Phi}\_N(t) = \mathbf{P}\boldsymbol{\Phi}\_N(t),\tag{39}$$

where,

$$\mathbf{U} = \mathbf{M}\mathbf{A}\mathbf{D}^2 + \mathbf{C}\mathbf{A}\mathbf{D} + \mathbf{K}\mathbf{A}.$$

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* ∗ *<sup>m</sup>*(*t*). Consequently:

$$
\langle \mathbf{U}\Phi\_N(t), L\_m^\*(t)\rangle = \langle \mathbf{P}\Phi\_N(t), L\_m^\*(t)\rangle, \qquad m = 0, 1, \ldots, N. \tag{40}
$$

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:

$$\begin{aligned} \mathbf{A}\boldsymbol{\phi}\_N(0) &= \mathbf{A},\\ \mathbf{A}\mathbf{D}\boldsymbol{\phi}\_N(0) &= \mathbf{H}, \end{aligned} \tag{41}$$

where,

$$\mathbf{A} = [\lambda\_1 \ \lambda\_2 \ \dots \ \lambda\_f]^T, \ \mathbf{II} = [\gamma\_1 \ \gamma\_2 \ \dots \ \gamma\_f]^T.$$

Finally, in order to find the response of the system in Equation (14) under the initial conditions from Equation (15), by replacing 2*f* of the rows of Equation (40) with Equation (38), a system of algebraic equations with (*N* + 1 − 2*f*) equations (Equation (38)) and 2*f* 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.

#### *3.3. Solution of a Calculation Example*

Consider the following equation as a SDOF structure without damping under a prescribed load:

$$\begin{cases} \ 0.1y''(t) + 4y(t) = 4(e^{-t} - e^{-15t}), \\ y(0) = 1, \quad y'(0) = -1. \end{cases} \tag{42}$$

The exact solution of this system is:

$$\begin{array}{l} y(t) = -0.3618\sin(6.3245t) + 0.1753\cos(6.3245t) \\ + 0.9756e^{-t} - 0.1509e^{-15t}. \end{array} \tag{43}$$

Choosing *N* = 5, the LGM method is used to solve the system. The matrix system can be obtained as below:

$$(\mathbf{M}\mathbf{A}\mathbf{D}^2 + \mathbf{K}\mathbf{A})\Phi\_N(t) = \mathbf{P}\Phi\_N(t),\tag{44}$$

where **M** = 0.1 and **K** = 4. Furthermore:

$$\mathbf{D}^2 = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 12 & 0 & 0 & 0 & 0 & 0 \\ 0 & 60 & 0 & 0 & 0 & 0 \\ 40 & 0 & 140 & 0 & 0 & 0 \\ 0 & 168 & 0 & 252 & 0 & 0 \end{bmatrix}, \quad \Phi\_N(t) = \begin{bmatrix} 1 \\ 2t - 1 \\ 6t^2 - 6t + 1 \\ 20t^3 - 30t^2 + 12t - 1 \\ 70t^4 - 140t^3 + 90t^2 - 20t + 1 \\ 252t^5 - 630t^4 + 560t^3 - 210t^2 + 30t - 1 \end{bmatrix}, \tag{45}$$

**A** = - *a*<sup>0</sup> *a*<sup>1</sup> . . . *a*<sup>5</sup> , **P** = - 2.2618 −0.5503 −0.6653 0.7842 −0.6008 ,

where φ*N*(*t*) is the shifted Legendre vector on the interval [0,1].

Using the Galerkin matrix method, the following algebraic system of equations is obtained:


Solving the above algebraic equations, the unknown Legendre coefficients are finally obtained as follows:

$$a\_0 = 0.60901, a\_1 = 0.07098, a\_2 = 0.27735, a\_3 = -0.40712, a\_4 = -0.12676, a\_5 = 0.09574.$$

By replacing the obtained values in Equation (17), the approximate solution for *y*(*t*) is:

$$y(t) \cong 24.1271 \, t^5 - 69.1914 \, t^4 + 63.2204 \, t^3 - 17.6372 \, t^2 - t + 1.\tag{47}$$

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.

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

#### **4. 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.

#### *4.1. 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, without damping: it is **M** = 1.5 0 0 2 , **C** = 0 and

**K** = 300 −300 <sup>−</sup>300 800 for mass, damping and stiffness parameters.

The initial conditions are *<sup>u</sup>*(0) = 1 1/2 and . *<sup>u</sup>*(0) = 0 0 .

• Free vibration, with damping: it is **C** = 1.628 −0.256 <sup>−</sup>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.

**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.** 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.

Successively, the examined structure is analyzed under forced vibration. In this case, it is:

$$p(t) = \begin{bmatrix} \sin 3t \\ 0 \end{bmatrix}, \mu(0) = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \\ \text{and } \dot{u}(0) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

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.

**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.** 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 Structure*

As a further validation example, a structure with three degrees of freedom (3DOF) is analyzed in forced and free vibration conditions, with or without damping.

• Free vibration, without damping: it is **M** = 2 0 0 0 2 0 0 0 2 , **<sup>C</sup>** <sup>=</sup> 0 and

$$\begin{aligned} \mathbf{K} &= \begin{bmatrix} 2 & -2 & 0 \\ -2 & 3 & -1 \\ 0 & -1 & 1 \end{bmatrix} \text{ for mass, damping and stiffness, respectively. The initial conditions are } \dot{u}(0) = \begin{bmatrix} 1 \\ 2/3 \\ 1/3 \end{bmatrix} \text{ and } \dot{u}(0) = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}. \\\ \text{Free vibration, with damping; it is } \mathbf{C} &= \begin{bmatrix} 0.8 & 0 & 0 \\ 0 & 0.7 & 0 \\ 0 & 0 & 0.6 \end{bmatrix} \text{, while the other parameters is } \mathbf{K} \end{aligned}$$

are the same as without damping.

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

$$\begin{array}{rcl} \text{(damping) for } p(t) &=& \begin{bmatrix} \sin t \\ 0 \\ 0 \end{bmatrix} \text{ with initial conditions } u(0) &=& \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \text{ and } \\\ \dot{u}(0) &=& \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \text{ is solved further by using Legendre's method.} \end{array}$$

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.

**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.** 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.** 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.** 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.** 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.** 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.
