Dynamic Response Analysis of Structures Using Legendre–Galerkin Matrix Method

Abstract

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.

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:

$$\left\{\begin{array}{l}\ddot{y}(t)=f(t,y(t),\dot{y}(t)),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{t}_0<t\le T,\\ y(t_0)=y_0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\dot{y}(t_0)={\dot{y}}_0.\end{array}\right.$$

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

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,13,14,15,16,17,18,19,20,21,22,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,25,26,27,28,29,30,31,32,33].

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.

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,58,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_0\sin (\mathrm{\Omega }t)$ with excitation frequency ς (Figure 1).

The governing equation of this system is [1]:

$$m\ddot{u}(t)+c\dot{u}(t)+ku(t)=P_0\sin (\mathrm{\Omega }t)$$

(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_0\sin (\mathrm{\Omega }t)$ with initial conditions $u(0)=u_0$ and $\dot{u}(0)={\dot{u}}_0$ is the sum of the displacements of the transient state of the system ($u_{transient}$) and the steady state ($u_{steady\hspace{0.17em}state}$), and is expressed as [5]:

$$u_{total}=u_{transient}+u_{steady\hspace{0.17em}state},$$

(3)

in which $u_{transient}$ and $u_{steady\hspace{0.17em}state}$ are defined as:

$$\begin{array}{l}{u}_{transient}=e^{-\xi \omega t}\left(u_0\cos {\omega }_{d}t+\frac{{\dot{u}}_0+u_0\omega \xi }{{\omega }_d}\sin {\omega }_{d}t\right)+\\ \frac{p_0}{k}\frac{e^{-\xi \omega t}}{{(1-{\beta }^2)}^2+{(2\xi \beta )}^2}\left[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\right],\end{array}$$

(4)

and

$$u_{steady\hspace{0.17em}state}=\frac{p_0}{k}\frac{1}{{(1-{\beta }^2)}^2+{(2\xi \beta )}^2}\left\{(1-{\beta }^2)\sin \mathrm{\Omega }t-2\xi \beta \cos \mathrm{\Omega }t\right\}$$

(5)

where $\xi$ and $\omega =\sqrt{\raisebox{1ex}{$k$}\!\left/ \!\raisebox{-1ex}{$m$}\right.}$ are the damping ratio and natural frequency of the system, respectively, and ${\omega }_d=\omega \sqrt{1-{\xi }^2}$ indicates the damped frequency. Finally, $\beta =\frac{\mathrm{\Omega }}{\omega }$ 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:

• Determination of the initial conditions, $\{u_0\}$, $\left\{{\dot{u}}_0\right\}$ and $\left\{{\ddot{u}}_0\right\}$, where the value of $\left\{{\ddot{u}}_0\right\}$ is given by:

  $${\ddot{u}}_0=\frac{p_0-C{\dot{u}}_0-K{u}_0}{M}$$

  (6)
• Determination of $\Delta t$ (assumed time step, e.g., 0.005)
• Determination of the effective stiffness matrix:

  $$\widehat{k}=k+\frac{\gamma }{\beta \Delta t}C+\frac{1}{\beta {\left(\Delta t\right)}^2}M$$

  (7)

  where $\beta =1/6$ and $\gamma =1/2$ for linear acceleration method.
• Calculation of a and b factors:

  $$\left\{\begin{array}{l}a=\frac{1}{\beta \Delta t}M+\frac{\gamma }{\beta }C\\ b=\frac{1}{2\beta }M+\Delta t\left(\frac{\gamma }{2\beta }-1\right)C\end{array}\right.$$

  (8)
• Calculation of the effective load:

  $$\Delta {\widehat{p}}_i=\Delta p_i+a{\dot{u}}_i+b{\ddot{u}}_i$$

  (9)
• Calculation of displacement, velocity, and acceleration:

  $$\left\{\begin{array}{l}{u}_{i+1}=u_i+\Delta u_i,\\ {\dot{u}}_{i+1}={\dot{u}}_i+\Delta {\dot{u}}_i,\\ {\ddot{u}}_{i+1}={\ddot{u}}_i+\Delta {\ddot{u}}_i,\end{array}\right.$$

  (10)

  where,

  $$\Delta u_i=\frac{\Delta {\widehat{p}}_i}{\hat{K}}$$

  (11)

  $$\Delta {\dot{u}}_i=\frac{\gamma }{\beta \Delta t}\Delta u_i-\frac{\gamma }{\beta }{\dot{u}}_i+\Delta t\left(1-\frac{\gamma }{2\beta }\right){\ddot{u}}_i$$

  (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.$$

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

$${\displaystyle \sum _{j=1}^f\left\{m_{ij}{\ddot{y}}_j(t)+c_{ij}{\dot{y}}_j(t)+k_{ij}{y}_j^{}(t)\right\}=p_i(t),}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,2,\dots f\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t\in [0,t_1],$$

(14)

under the following initial conditions:

$$\left\{\begin{array}{c}{y}_i\left(0\right)={\lambda }_i,\\ {y^{\prime }}_i\left(0\right)={\gamma }_i,\end{array}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\hspace{0.17em}\hspace{0.17em}2,\hspace{0.17em}\hspace{0.17em}\dots ,f.\right.$$

(15)

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 ${\lambda }_i$ and ${\gamma }_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.

3.1. Approximation of the Function Using Shifted Legendre Polynomials

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

  $$\begin{array}{l}{L}_0(t)=1,\\ L_1(t)=t,\\ L_2(t)=\frac{3}{2}{t}^2-\frac{1}{2},\\ ⋮,\\ L_{m+1}(t)=\frac{2m+1}{m+1}t\hspace{0.17em}{L}_m(t)-\frac{m}{m+1}{L}_{m-1}(t).\end{array}$$

  (16)

For further application of these polynomials, they are shifted to the interval [0, L] with the change of variable $\frac{2}{L}t-1$. Therefore, these shifted polynomials which are shown by $L_m^{*}(t)$ are represented as the following:

$$L_m^{\ast }(t)=L_m(\frac{2}{L}t-1)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}t\in [0,L].$$

(17)

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

  $$〈f(t),g(t)〉={\int }_a^{b}f(t)g(t)w(t)dt.$$

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

  $$〈f(t),g(t)〉=0.$$

  (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^{\ast }(t)L_n^{\ast }(t)dt=\left\{\begin{array}{l}0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}m\ne n,\\ \frac{L}{2m+1}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}m=n.\end{array}\right.$$

(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)={\displaystyle \sum _{m=0}^{\infty }{f}_m{L}_m^{*}(t)},$$

(21)

where f_m can be obtained from:

$$f_m=\frac{〈f(t),L_m^{*}(t)〉}{〈L_m^{*}(t),L_m^{*}(t)〉}=\frac{2m+1}{L}{\int }_0^{L}f(t)L_m^{*}(t)dt.$$

(22)

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

$$f(t)={\displaystyle \sum _{m=0}^N{f}_m{L}_m^{*}(t)={\mathbf{F}}^T{\mathsf{\varphi }}_N(t)},$$

(23)

where $\mathbf{F}$ and ${\mathsf{\varphi }}_N(\mathrm{t})$ are factors of the shifted Legendre vectors expressed as:

$$\mathbf{F}={[f_{0\hspace{0.17em}}\hspace{0.17em}{f}_{1\hspace{0.17em}}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{f}_N]}^T,$$

(24)

$${\mathsf{\varphi }}_N(t)={[L_0^{*}(t)\hspace{0.17em}\hspace{0.17em}{L}_1^{*}(t)\hspace{0.17em}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{L}_N^{*}(t)]}^T.$$

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

$$y_{i,N}^{}(t)={\displaystyle \sum _{m=0}^N{a}_{i,m}{L}_m^{*}(t)}={\mathbf{A}}_{{}_{i,N}}^T{\mathsf{\varphi }}_N(t)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em},$$

(26)

where,

$$\begin{array}{l}{\mathbf{A}}_{i,N}={[a_{i,0}\hspace{0.17em}\hspace{0.17em}{a}_{i,1}\hspace{0.17em}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{a}_{i,N}]}^T,\\ {\mathsf{\varphi }}_N(t)=[L_0^{*}(t)\hspace{0.17em}\hspace{0.17em}{L}_1^{*}(t)\hspace{0.17em}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{L}_N^{*}(t)].\end{array}$$

(27)

The matrix form of the derivative vector will be:

$$\frac{d}{dt}{\mathsf{\varphi }}_N(t)=\mathbf{D}{\mathsf{\varphi }}_N(t),$$

(28)

where $\mathbf{D}$ is a matrix of dimensions (N + 1) × (N + 1) and is obtained as below:

$$\mathbf{D}=[d_{ij}]=\left\{\begin{array}{l}\frac{2(2j+1)}{L},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=i-k,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\begin{array}{l}k=1,3,\dots ,N\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{if} \mathrm{N} \mathit{is} \mathit{odd} \mathit{number},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ k=1,3,\dots ,N-1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{if} \mathrm{N} \mathit{is} \mathit{even} \mathit{number},\end{array}\right.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ 0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}otherwise.\end{array}\right.$$

(29)

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

$$\begin{array}{l}If N=1 \to D=\left[\begin{array}{ll}0\hfill & 0\hfill \\ \frac{2}{L}\hfill & 0\hfill \end{array}\right]\hfill \\ If N=2 \to D=\left[\begin{array}{lll}0\hfill & 0\hfill & 0\hfill \\ \frac{2}{L}\hfill & 0\hfill & 0\hfill \\ 0\hfill & \frac{6}{L}\hfill & 0\hfill \end{array}\right]\hfill \\ If N=3 \to D=\left[\begin{array}{llll}0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \frac{2}{L}\hfill & 0\hfill & 0\hfill & 0\hfill \\ 0\hfill & \frac{6}{L}\hfill & 0\hfill & 0\hfill \\ \frac{2}{L}\hfill & 0\hfill & \frac{10}{L}\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

Using Equation (27), the k-th derivative of ${\mathsf{\varphi }}_N(t)$ vector is:

$$\frac{d^k}{d{t}^k}{\mathsf{\varphi }}_N(t)={\mathbf{D}}^k{\mathsf{\varphi }}_N(t),$$

(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)={\displaystyle \sum _{m=0}^N{a}_{i,m}{L}_m^{*}(t)}={\mathbf{A}}_{{}_{i,N}}^T\mathbf{D}{\mathsf{\varphi }}_N(t),$$

(31)

and

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

(32)

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

$$p_i(t)={\displaystyle \sum _{m=0}^N{p}_{i,m}{L}_m^{*}(t)}={\mathbf{P}}_{i,N}^T{\mathsf{\varphi }}_N(t),$$

(33)

where ${\mathit{P}}_{i,N}^T$ is a (N+1) dimensional vector with the following expression:

$${\mathbf{P}}_{i,N}={[p_{i,0}\hspace{0.17em}\hspace{0.17em}{p}_{i,1}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{p}_{i,N}]}^T,$$

(34)

where,

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

(35)

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

$${\displaystyle \sum _{j=1}^f\left\{m_{ij}{\mathbf{A}}_{{}_{j,N}}^T{\mathbf{D}}^2+c_{ij}{\mathbf{A}}_{{}_{j,N}}^T\mathbf{D}+k_{ij}{\mathbf{A}}_{{}_{j,N}}^T\right\}{\mathsf{\varphi }}_N(t)={\mathbf{P}}_{i,N}^T{\mathsf{\varphi }}_N(t)}$$

(36)

Now, let us assume that:

$$\begin{gathered} {\mathbf{Y}}_N=\mathbf{A}{\mathsf{\varphi }}_N(t)=\left[\begin{array}{c}{y}_{1,N}\\ y_{2,N}\\ ⋮\\ y_{f,N}\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{A}=\left[\begin{array}{c}{\mathbf{A}}_{1,N}^T\\ {\mathbf{A}}_{2,N}^T\\ ⋮\\ {\mathbf{A}}_{f,N}^T\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{P}=\left[\begin{array}{c}{\mathbf{P}}_{1,N}^T\\ {\mathbf{P}}_{2,N}^T\\ ⋮\\ {\mathbf{P}}_{f,N}^T\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} \\ \mathbf{M}=[m_i{}_j],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{C}=[c_{ij}],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{K}=[k_{ij}],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i,j=1,2,\dots ,f, \end{gathered}$$

(37)

where ${\mathbf{Y}}_N$ is an f dimensional vector, A and P are $f\times (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}){\mathsf{\varphi }}_N(t)=\mathbf{P}{\mathsf{\varphi }}_N(t),$$

(38)

or

$$\mathbf{U}{\mathsf{\varphi }}_N(t)=\mathbf{P}{\mathsf{\varphi }}_N(t),$$

(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 ${\mathsf{\varphi }}_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:

$$〈\mathbf{U}{\mathsf{\varphi }}_N(t),L_m^{*}(t)〉=〈\mathbf{P}{\mathsf{\varphi }}_N(t),L_m^{*}(t)〉,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}m=0,1,\dots ,N.$$

(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{array}{l}\mathbf{A}{\mathsf{\varphi }}_N(0)=\mathsf{\Lambda },\\ \mathbf{A}\mathbf{D}{\mathsf{\varphi }}_N(0)=\mathsf{\Pi },\end{array}$$

(41)

where,

$$\mathsf{\Lambda }={[{\lambda }_1\hspace{0.17em}\hspace{0.17em}{\lambda }_2\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{\lambda }_f]}^T,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathsf{\Pi }={[{\gamma }_1\hspace{0.17em}\hspace{0.17em}{\gamma }_2\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em}{\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 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.

3.3. Solution of a Calculation Example

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

$$\left\{\begin{array}{l}0.1y^{″}(t)+4y(t)=4(e^{-t}-e^{-15t}),\\ y(0)=1,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{y}^{\prime }(0)=-1.\end{array}\right.$$

(42)

The exact solution of this system is:

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

(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}){\mathsf{\varphi }}_N(t)=\mathbf{P}{\mathsf{\varphi }}_N(t),$$

(44)

where $\mathbf{M}=0.1$ and $\mathbf{K}=4$. Furthermore:

$$\begin{gathered} {\mathbf{D}}^2=\left[\begin{array}{cccccc}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{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathsf{\varphi }}_N(t)=\left[\begin{array}{c}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{array}\right], \\ \mathbf{A}=\left[\begin{array}{cccc}{a}_0& a_1& \dots & a_5\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathbf{P}=\left[\begin{array}{ccccc}2.2618& -0.5503& -0.6653& 0.7842& -0.6008\end{array}\right], \end{gathered}$$

(45)

where ${\mathsf{\varphi }}_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:

$$\left[\begin{array}{cccccc}4& 0& 1.2& 0& 4& 0\\ 0& 1.33& 0& 2& 0& 5.6\\ 0& 0& 0.8& 0& 2.8& 0\\ 0& 0& 0& 0.5714& 0& 3.6\\ 1& -1& 1& -1& 1& -1\\ 0& 2& -6& 12& -20& 30\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\end{array}\right]=\left[\begin{array}{c}2.2618\\ -0.1834\\ -0.1331\\ 0.1120\\ 1\\ -1\end{array}\right]$$

(46)

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\hspace{0.17em}{t}^5-69.1914\hspace{0.17em}{t}^4+63.2204\hspace{0.17em}{t}^3-17.6372\hspace{0.17em}{t}^2-t+1.$$

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

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 $\mathbf{M}=\left[\begin{array}{cc}1.5& 0\\ 0& 2\end{array}\right]$, $\mathbf{C}=0$ and $\mathbf{K}=\left[\begin{array}{cc}300& -300\\ -300& 800\end{array}\right]$ for mass, damping and stiffness parameters.

The initial conditions are $u(0)=\left[\begin{array}{c}1\\ 1/2\end{array}\right]$ and $\hspace{0.17em}\hspace{0.17em}\dot{u}(0)=\left[\begin{array}{c}0\\ 0\end{array}\right]$.

• Free vibration, with damping: it is $\mathbf{C}=\left[\begin{array}{cc}1.628& -0.256\\ -0.256& 2.512\end{array}\right]$, 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 Figure 3 and Figure 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:

$$p(t)\hspace{0.17em}=\left[\begin{array}{c}\sin 3t\\ 0\end{array}\right]\hspace{0.17em},\hspace{0.17em}u(0)\hspace{0.17em}=\left[\begin{array}{c}0\\ 0\end{array}\right]\hspace{0.17em}\mathrm{and}\hspace{0.17em}\dot{u}(0)\hspace{0.17em}=\left[\begin{array}{c}0\\ 0\end{array}\right]$$

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 Figure 5 and Figure 6, respectively, and compared with the exact solution, proving again a rather acceptable consistency in the examined time interval.

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 $\mathbf{M}=\left[\begin{array}{ccc}2& 0& 0\\ 0& 2& 0\\ 0& 0& 2\end{array}\right]$, $\mathbf{C}=0$ and $\mathbf{K}=\left[\begin{array}{ccc}2& -2& 0\\ -2& 3& -1\\ 0& -1& 1\end{array}\right]$ for mass, damping and stiffness, respectively. The initial conditions are $u(0)=\left[\begin{array}{c}1\\ 2/3\\ 1/3\end{array}\right]$ and $\hspace{0.17em}\hspace{0.17em}\dot{u}(0)=\left[\begin{array}{c}\begin{array}{c}0\\ 0\end{array}\\ 0\end{array}\right]$.
• Free vibration, with damping: it is $\mathbf{C}=\left[\begin{array}{ccc}0.8& 0& 0\\ 0& 0.7& 0\\ 0& 0& 0.6\end{array}\right]$, while the other parameters are the same as without damping.

The results for $u_1(t)$, $u_2(t)$ and $u_3(t)$ are shown in Figure 7, Figure 8 and Figure 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 damping) for $p(t)\hspace{0.17em}=\left[\begin{array}{c}\sin t\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]\hspace{0.17em}$ with initial conditions $u(0)=\left[\begin{array}{c}0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]\hspace{0.17em}$ and $\hspace{0.17em}\dot{u}(0)=\left[\begin{array}{c}0\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right]$ is solved further by using Legendre’s method.

The results obtained in this case for $u_1(t)$, $u_2(t)$ and $u_3(t)$ are proposed in Figure 10, Figure 11 and Figure 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.

4.3. 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 $\mathbf{M}=\left[\begin{array}{ccccc}1& 0& 0& 0& 0\\ 0& 2& 0& 0& 0\\ 0& 0& 2& 0& 0\\ 0& 0& 0& 2& 0\\ 0& 0& 0& 0& 2\end{array}\right]$, $\mathbf{C}=0$ and $\mathbf{K}=\left[\begin{array}{ccccc}4& -4& 0& 0& 0\\ -4& 12& -8& 0& 0\\ 0& -8& 20& -12& 0\\ 0& 0& -12& 24& -12\\ 0& 0& 0& 12& 28\end{array}\right]$, respectively.

The initial conditions are $u(0)=\left[\begin{array}{c}1\\ 4/5\\ 3/5\\ 2/5\\ 1/5\end{array}\right]\hspace{0.17em}\hspace{0.17em}$ and $\dot{u}(0)=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right]$.

• Free vibration, with damping: it is assumed that $\mathbf{C}=\left[\begin{array}{ccccc}0.5& 0& 0& 0& 0\\ 0& 0.5& 0& 0& 0\\ 0& 0& 0.5& 0& 0\\ 0& 0& 0& 0.5& 0\\ 0& 0& 0& 0& 0.5\end{array}\right]$, while the other parameters are the same as in the undamped case.

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 Figure 13, Figure 14, Figure 15, Figure 16 and Figure 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.

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