# Frequency-Domain Assessment of Integration Schemes for Earthquake Engineering Problems

Authors: Juana Arias-Trujillo; Rafael Blázquez; Susana López-Querol
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Although numerical integration is a technique commonly employed in many time-dependent problems, usually its accuracy relied on a time interval small enough. However, taking into account that time integration formulae can be considered to be recursive digital filters, in this research a criterion based on transfer functions has been employed to characterize a wide range of integration algorithms from a frequency approach, both in amplitude and in phase. By adopting Nyquist’s criterion to avoid the aliasing phenomena, a total of seven integration schemes have been reviewed in terms of accuracy and distortion effects on the frequency content of the signal. Some of these schemes are very well-known polynomial approximations with different degrees of interpolation, but others have been especially defined for solving earthquake engineering problems or have been extracted from the digital signal processing methodology. Finally, five examples have been developed to validate this frequency approach and to investigate its influence on practical dynamic problems. This research, focused on earthquake and structural engineering, reveals that numerical integration formulae are clearly frequency-dependent, a conclusion that obviously has a relevant interest in all dynamic engineering problems, even when they are formulated and solved in the time-domain.


---

## Body

## 1. Introduction

In earthquake engineering problems, integration algorithms are daily used, especially to approach time-dependent problems. Some of the most employed among them are the well-known closed Newton-Cotes formulae with different degrees of interpolation polynomial functions, which includes the endpoints of the integration domain(a,b) as nodes of the problem [1]. In recent years, many other formulae have been developed “ad hoc” for earthquake engineering problems or have been borrowed from other areas of knowledge, in particular from digital signal processing. Usually, the choice of an integration algorithm depends on its simplicity, its computational cost, and its accuracy. The later factor is difficult to assess and very often it is entrusted to an integration interval small enough. However, knowing the computational error associated with a numerical integration scheme is a crucial task in dynamic problems and has been, and it is nowadays, a topic of growing interest in many researches [2–4].Although the main advantages and disadvantages of the most common methods are well-known, in this study the performance of these schemes is analyzed attending to its applicability to earthquake engineering and dynamic problems by means of a frequency-domain assessment, in spite of the fact that time is the usual integration variable. This research extends a frequency approach typical of the field of digital signal processing to various numerical algorithms, which can be considered to be recursive digital filters [5–9]. This methodology allows us to analyze not only the accuracy of the algorithms but also the distortion introduced by them in the frequency content of the signal being integrated, or how the frequency content of the signal affects the accuracy of each algorithm. This new way of approaching the topic of the accuracy of the algorithms, using the frequency domain, can be of great usefulness in earthquake engineering problems. Several researchers have also employed similar methodologies to examine the performance of other time-dependent numerical algorithms within the area of Soil and structural dynamics problems [10–13].In this paper seven integration formulae have been analyzed by comparing the corresponding numerical transfer functions to the analytical one since all of them are frequency-dependent. By means of this comparison, a rational criterion is proposed to predict in advance the performance of a given numerical scheme, depending on time interval (h) and the frequency content of the input signal (ω). Two examples are shown to check the performance of the different methods and to verify the proposed criterion. An additional example is included to highlight these findings applied to a practical earthquake engineering problem.


## 2. Brief Revision of Integration Formulae

In general, any numerical integration formula tries to approximate the function to be integratedX(t) (1) by an interpolation polynomial function of degree n, through known values of it, Xi=X(ti), which are defined at the integration points ti placed at equal time interval, h=tn+1-tn:(1)Yt=∫abXtdt.In this research seven formulae are included. The first four are used in plenty of disciplines and are defined as the closed Newton-Cotes formulae. They have different degree(n) of interpolation polynomial function [1] and they are as follows: trapezoidal rule (n=1), Simpson’s rule (n=2), Simpson three-eights rule (n=3), and Newton-Cotes scheme with n=4, which is labeled as NC n=4 along this document. The fifth algorithm is the rectangle rule, which is widely used due to its simplicity. The sixth formula is called Tick’s rule [5], which is well-known in the field of digital signal processing. This formula was designed to have a transfer function as close to unity as possible throughout the lower half of the Nyquist interval involving only three consecutive points [5]. Transfer functions and Nyquist’s criterion are also used in this research and they are explained in Section 3. Finally, the seventh method is the digital integration filter defined by Schuessler and Ibler [14]. The mathematical definitions of all the aforementioned algorithms are listed as follows:(i)
Trapezoidal rule:(2)Yn+1=Yn+h2Xn+Xn+1.(ii)
Simpson’s rule:(3)Yn+2=Yn+h3Xn+4Xn+1+Xn+2.(iii)
Simpson three-eights rule:(4)Yn+3=Yn+3h8Xn+3Xn+1+3Xn+2+Xn+3.(iv)
Newton-Cotesn=4:(5)
Y
n
+
4
=
Y
n
+
2
h
45
7
X
n
+
32
X
n
+
1
+
12
X
n
+
2
+
32
X
n
+
3
+
7
X
n
+
4
.(v)
Rectangle rule:(6)Yn+1=Yn+hXn.(vi)
Tick’s rule:(7)Yn+2=Yn+hb0Xn+b1Xn+1+b2Xn+2,whereb0=b2=0.3584 and b1=1.2832. It must be highlighted that if b0=b2=1/3 and b1=4/3, Tick’s rule (7) becomes identical to Simpson’s rule (3).(vii)
Schuessler-Ibler rule:(8)Yn+1=Yn+h3∑k=0k=7bkXn+1-kn=0,1,2,…,whereb0=b7=-0.01330 and b1=b6=0.07760 and b2=b5=-0.28279 and b3=b4=1.71822.


## 3. Methodology

Normally, integration algorithms are defined as a linear combination of the data or sampled points of the function to be integrated. When these are data equally spaced, the integration algorithms can be regarded as digital filters [5]. In fact, integration formulae have been a typical application of filtering techniques in digital signal processing methodology, but only some of the most classical rules have been considered in practice [5–7].The performance of any numerical integration algorithm can be analyzed in the frequency domain by means of its transfer function (TF). The transfer function of a system can be defined as a mathematical operator which connects the input signal to the output response. In this particular case, the integration algorithm is the system itself. To obtain the transfer function, it is assumed a complex harmonic function is input signal,eiωt, where ω is the angular frequency of the signal and t is the time. Since all of the above integration algorithms, (2) to (8), are linear and time-invariant systems, the output signal (after integration) can be defined as H(iω)e(iωt), where H(iω) is the complex transfer function (TF) of the corresponding integration scheme which depends on ω and t [5]. Also, since H(iω) is a complex fuction, it is possible to derive the corresponding expressions for amplitude and phase for the function.In the following sections, it is shown how the exact and numerical transfer functions formulae can be obtained. Subsequently, comparing the corresponding numerical transfer function to the exact one, it is possible to develop a criterion to analyze the performance of a particular integration scheme using a frequency approach.


### 3.1. Exact Transfer Function

For the input signalX(t)=eiωt, the exact integration is shown in the first part of (9). Applying to (9) the definition of transfer function given above, the transfer function for the exact integrator (named as H1∗(iω) along this research) can be obtained (10). The symbol ∗ is referring to the exact integrator and the subindex 1 is referring to the single integral. The exact transfer function can be expressed in different forms according to the properties of complex numbers:(9)∫Xtdt=∫eiωtdt=eiωtiω=H1∗iωeiωt⟶(10)H1∗iω=1iω=-1ωi=1ωe-iπ/2.Now if (9) is integrated again, it is possible to obtain the exact transfer function for the double integral, H2∗(iω), as developed in (11):(11)∬Xtdt=∫H1∗iωeiωtdt=H1∗iω∫eiωtdt=H1∗iω2eiωt=H2∗iωeiωt⟶H2∗iω=H1∗iω2.


### 3.2. Numerical Transfer Functions

Similarly to the exact case, it is possible to compute the transfer functions corresponding to each numerical integration method, when the input signal iseiωt [5, 8–10]. For the sake of conciseness, hereafter the transfer functions of the seventh numerical algorithms previously mentioned are listed, while the complete derivation of every formula can be found in the Appendix:(i)
Trapezoidal rule:(12)H1iω=h21+eiωheiωh-1.(ii)
Simpson’s rule:(13)H1iω=h31+4eiωh+eiω2heiω2h-1.(iii)
Simpson three-eights rule:(14)H1iω=3h81+3eiωh+3eiω2h+eiω3heiω3h-1.(iv)
Newton-Cotesn=4:(15)H1iω=2h457+32eiωh+12eiω2h+32eiω3h+7eiω4heiω4h-1.(v)
Rectangle rule:(16)H1iω=heiωh-1.(vi)
Tick’s rule:(17)H1iω=hb0+b1eiωh+b2eiω2heiω2h-1,whereb0=b2=0.3584 and b1=1.2832.(vii)
Schuessler-Ibler rule:(18)H1iω=h/3b0eiωh+b1+b2e-iωh+b3e-2iωh+b4e-3iωh+b5e-4iωh+b6e-5iωh+b7e-6iωheiωh-1,whereb0=b7=-0.01330 and b1=b6=0.07760 and b2=b5=-0.28279 and b3=b4=1.71822.The procedure followed in the exact case to obtain the transfer function for the double integral (11) is also analogous for the numerical case H2(iω) (19):(19)H2iω=H1iω2.


## 3.1. Exact Transfer Function

For the input signalX(t)=eiωt, the exact integration is shown in the first part of (9). Applying to (9) the definition of transfer function given above, the transfer function for the exact integrator (named as H1∗(iω) along this research) can be obtained (10). The symbol ∗ is referring to the exact integrator and the subindex 1 is referring to the single integral. The exact transfer function can be expressed in different forms according to the properties of complex numbers:(9)∫Xtdt=∫eiωtdt=eiωtiω=H1∗iωeiωt⟶(10)H1∗iω=1iω=-1ωi=1ωe-iπ/2.Now if (9) is integrated again, it is possible to obtain the exact transfer function for the double integral, H2∗(iω), as developed in (11):(11)∬Xtdt=∫H1∗iωeiωtdt=H1∗iω∫eiωtdt=H1∗iω2eiωt=H2∗iωeiωt⟶H2∗iω=H1∗iω2.


## 3.2. Numerical Transfer Functions

Similarly to the exact case, it is possible to compute the transfer functions corresponding to each numerical integration method, when the input signal iseiωt [5, 8–10]. For the sake of conciseness, hereafter the transfer functions of the seventh numerical algorithms previously mentioned are listed, while the complete derivation of every formula can be found in the Appendix:(i)
Trapezoidal rule:(12)H1iω=h21+eiωheiωh-1.(ii)
Simpson’s rule:(13)H1iω=h31+4eiωh+eiω2heiω2h-1.(iii)
Simpson three-eights rule:(14)H1iω=3h81+3eiωh+3eiω2h+eiω3heiω3h-1.(iv)
Newton-Cotesn=4:(15)H1iω=2h457+32eiωh+12eiω2h+32eiω3h+7eiω4heiω4h-1.(v)
Rectangle rule:(16)H1iω=heiωh-1.(vi)
Tick’s rule:(17)H1iω=hb0+b1eiωh+b2eiω2heiω2h-1,whereb0=b2=0.3584 and b1=1.2832.(vii)
Schuessler-Ibler rule:(18)H1iω=h/3b0eiωh+b1+b2e-iωh+b3e-2iωh+b4e-3iωh+b5e-4iωh+b6e-5iωh+b7e-6iωheiωh-1,whereb0=b7=-0.01330 and b1=b6=0.07760 and b2=b5=-0.28279 and b3=b4=1.71822.The procedure followed in the exact case to obtain the transfer function for the double integral (11) is also analogous for the numerical case H2(iω) (19):(19)H2iω=H1iω2.


## 4. Frequency Performance of Integration Algorithms

Since the exact and numerical transfer functions have been, respectively, calculated in (10) and in (12) to (18), it becomes possible to assess the performance of each algorithm relative to the exact solution by comparing the corresponding numerical transfer function with the exact one. In order to simplify this task, the ratio R has been defined as the quotient between the module of the numerical and exact transfer functions. This ratio can be obtained both for the single integral R1, as shown in (20), and for the double integral R2, as shown in (21). In both cases, the closer the ratio to unity, the more accurate the numerical solution:(20)R1=H1iωH1∗iω,(21)R2=H2iωH2∗iω=H1iω2H1∗iω2=H1iωH1∗iω2=R12.To define the valid range of frequencies, avoiding the aliasing phenomena, Nyquist’s criterion has been taken into account [5]. This criterion defines the maximum frequency value which can be achieved by an input signal depending on the time interval (h) used to sample the signal. Owing to that, the maximum frequency on the input signal (fmaxinput or ωmaxinput) must be less than Nyquist’s frequency (fN or ωN) (22) which implies that (23) must be guaranteed:(22)fN=12h⟶ωN=2π2h=πh,(23)fmaxinput<fN⟶ωmaxinput<ωN=πh.In this research, the most common values for time intervals employed for strong ground motion discretization (h=0.02, 0.01, and 0.005 sec) have been adopted. The corresponding maximum values of the input frequencies which can be reached to avoid aliasing phenomena are shown in Table 1. As it can be observed, the lower the value of h, the higher the allowed input frequency. For simplicity, the dimensionless parameter λ, defined as the quotient between the input frequency (ω) and Nyquist’s frequency (ωN), is adopted as shown in (24)(24)λ=ωωN=ωhπ.Table 1
Maximum input frequencies for different values ofh in order to avoid aliasing phenomena according to Nyquist’s criterion.
h (sec)
f
N (Hz)
ω (rad/sec)
0.02
25
50π
0.01
50
100π
0.005
100
200πFigure1 shows the variations with respect to λ of the ratios R1 (20) and R2 (21) for the single and double integral, respectively, and the phase diagram (ϕ), for the seven integration algorithms selected previously. As previously mentioned, the ideal numerical algorithm should maintain the ratios R1 and R2 as close to unity as possible for a large range of λ, and also the phase value should coincide with the exact one.Figure 1
Ratios between numerical and exact amplitude transfer functions (R1-first integration, R2-second integration) and phase diagrams for tested algorithms.
(a)
R
1: single integration
(b)
R
2: double integration
(c)
Phase diagramSince the transfer function is defined by a complex expression, it is convenient to analyze not only its amplitude but also its phase diagram. From Figure1, relevant conclusions are extracted which ratify the earlier conclusions reported by other researchers [5–9]. The main conclusions, which are later on addressed by the examples included in Section 5, are listed next:(1)
Although the parameterλ guarantees that the aliasing phenomena are avoided, it can be observed that in general all algorithms show values of R1 and R2 far away from unity for the upper-half of λ. This effect is even more critical for the R2 case.(2)
Trapezoidal and Schuessler-Ibler rules tend to deamplify the output signal (after integration), while the rectangle, Tick, and Simpson algorithms tend to amplify it, although in different quantity according toλ. On the other hand, the Simpson 3-8 and Newton-Cotes n=4 rules display a singular behaviour: depending on λ, they can amplify or deamplify the output signal.(3)
For the single integral case (R1), Tick’s rule slightly deamplifies the output signal for intermediate values of λ, while it tends to amplify from λ>0.5. On the other hand, Simpson’s rule consistently amplifies the output signal especially for λ above 0.3. The rectangle rule amplifies for values of λ higher than 0.1, in other words for almost the whole range of λ. By contrast, the trapezoidal rule progressively attenuates the output signal for almost the entire range of frequencies whereas Schuessler-Ibler strongly deamplifies beyond λ>0.6. For Simpson 3-8 rule, R1 is close to unity until λ≈0.3, tending to amplify the signal for higher values of λ (with a discontinuity around λ=0.66) until reaching values around 0.72, deamplifying significantly after that. Finally, the Newton-Cotes n=4 formula is close to unity up to values of λ around 0.4, deamplifying quickly for λ between 0.4 and 0.5. For λ=0.5 the transfer function of this scheme presents also a discontinuity, tending to amplify the output signal in a nonuniform manner for higher values of λ.(4)
For the double integral case (R2), all methods show the same performance (amplification or deamplification) compared to the single integration (R1), but, in most of the formulae, the ratio R2 is further from unity for lower values of λ than R1, except for Tick and Schuessler-Ibler’s rules. In general, this implies that, given a signal and a value of h, the error introduced into the solution by the numerical procedures grows for upper integration orders, although in different magnitude for each method.(5)
Those algorithms tending to amplify for the upper-half ofλ may introduce undesirable errors, specifically high-frequency terms, into the integrated signal. Obviously, these algorithms should be avoided if the input signal is contaminated by spurious high-frequency noise.(6)
Regarding phase values, the exact transfer function has a phase value equal to-π/2 for the whole range of λ. Only the trapezoidal, Simpson, and Tick rules show the same phase behaviour, while the rest of numerical algorithms have a phase performance absolutely different. The Simpson 3-8 and Newton-Cotes n=4 alternatively display phase values equal to -π/2 or equal to π/2. However the phase of rectangle rule linearly varies from Φ=π/2 for λ=0 to Φ=0 for λ=1. The phase of Schuessler-Ibler’s rule periodically varies along the whole range of λ, and in a linear way between Φ=π/2 and Φ=-π/2 in the intervals of λ between 0–1/3, 1/3–2/3, and 2/3–1. Thus, both rules introduce a frequency-dependent phase error in the integration scheme for any value of λ.


## 5. Numerical Examples

Five common examples from earthquake engineering have been worked out to illustrate the applicability of the frequency criterion derived in this research, as well as the validity of the above listed conclusions concerning the performance of the seven selected integration algorithms. In the first three examples, the analytical solutions are known in advance and therefore the accuracy of the numerical results can be clearly assessed. In the fourth one, the permanent displacement of a slope due to an earthquake is computed by the Newmark sliding block method, and, in the last one, the error in the displacement record of a strong ground motion is obtained, so that the frequency-effects of the integration algorithms can be analyzed in practical engineering problems.


### 5.1. Example  1: Variable Frequency Input

In order to assess the accuracy of the integration schemes using the proposed frequency approach, a sine function with variable frequency is considered input signal (acceleration) to be integrated. The acceleration is defined in (25), where Ac is the amplitude equal to 3 m/s2, ω is the frequency of the input signal, and t denotes time, whereas zero initial conditions have been assumed. The acceleration has a duration of 1.5 s and the time interval (h) is equal to 0.005 s, meaning that the input frequency ranges from 2π rad/s to 200π rad/s reaching a value of λ between 0.01 and 1, respectively. The input signal is integrated twice to calculate the velocity and displacement records, and after that the maximum velocity (vmax) and the maximum displacement (dmax) are calculated and compared with the analytical solutions which can be obtained from (26) and (27). Consider(25)at=Acsin⁡ωt,(26)vt=∫0tAcsin⁡ωtdt=Acω1-cos⁡ωt,(27)dt=∫0tAcω1-cos⁡ωtdt=Acωt-sin⁡ωtω.Figure2 shows the relative error (in percentage) for both |vmax| and |dmax|, for the whole range of values of λ. A positive relative error means that the corresponding numerical algorithm amplifies the numerical solution and vice versa. From Figure 2, it is concluded that the ratios R1 and R2, plotted in Figure 1, properly predict the performance of each numerical integration algorithm, since the error curves shown in Figure 2 follow exactly the trend plotted in Figure 1. It means that the higher the value of λ, the higher the error, except for the rectangle and Schuessler-Ibler formulae. In the case of the rectangle rule, R1 and R2 predict an amplification in the numerical response whereas the error plotted in Figure 2 shows a deamplification tendency, which is even closer to the trapezoidal rule for the displacement case. On the other hand, this rule presents a positive linear variation of the phase depending on λ, while the exact phase has a negative constant value. This contradictory phase performance could explain the abnormal response of this algorithm. As concerns the Schuessler-Ibler formula, it deamplifies the output record according to R1 and R2 but in a nonuniform way, since the error curve shows an evident periodic behaviour according to its phase diagram. It must be highlighted that the highest errors take place when the phase curve is equal to zero and the lowest errors take place when the phase reaches any of the extreme values.Figure 2
Relative error (in percentage) reached by each algorithm along the whole range ofλ for a sinusoidal signal (both velocity and displacement records are shown).
(a)
Relative error (%) for|vmax|
(b)
Relative error (%) for|dmax|Although this example has been computed with a constant value ofh, it can be observed that numerical integration accuracy is frequency-dependent, a matter of paramount interest in earthquake engineering problems. Moreover, if this example is computed for other values of time interval h, with λ ranging between 0 and 1, these algorithms show similar trends.


### 5.2.  Example  2: Influence of Noise Content of a Signal

In earthquake engineering problems, it is necessary to employ different types of signal as accelerograms or records measured from the response of a structure. Any kind of record could be contaminated by different noise sources (internal characteristics of the measurement apparatus, the placement where the instrument is located, ambient noise, etc.) [15, 16]. In this example, the sensitivity of the different integration algorithms to the noise content in a signal is analyzed. A sine signal has been modified by three types of noise and then has been integrated twice to compute the corresponding displacement record. Thanks to the simplicity of a sine function and since its analytical solution is known in advance ((26) and (27)), it is possible to verify the influence of the noise components of a signal in the final results for each integration algorithm.The original sine function (labeled as original signal is this research) is defined in (25), where Ac=3 m/s2 and ω=10π rad/s, with a duration of 1.5 s and h=0.005 s (λ=0.05). Three different types of noise have been chosen, high-frequency, low-frequency, and phase-noise. In the three cases, the contaminated signal is defined as a sine function too, which has been added to the original one. The characteristics of the high-frequency noise are Ahigh=0.3 m/s2 and ωhigh=150π rad/s, while for the low-frequency noise are Alow=0.6 m/s2 and ωlow=0.1π rad/s. For the phase-noise case, a constant variation in the phase equal to 0.15π rad/s has been added to the frequency of the original sine function. The four signals considered in this example are plotted in Figure 3; as can be observed, the noise-signals are very close to the original one.Figure 3
Original function and contaminated signals by three types of noise: high-frequency noise, low-frequency noise, and phase-noise.The three contaminated signals from Figure3 have been integrated twice to obtain the velocity and displacement records, Figure 4. In all cases, the analytical result from the original sine function has been included for comparative purposes. Moreover, in all cases, the solutions computed by the Schuessler-Ibler rule show an anomalous behaviour relative to the rest of methods, both in velocity and in displacement. In the cases of high-frequency noise and phase-noise, the computed results are very close to the analytical solution, even when this has been computed for the original signal without any kind of noise. However, in the case of low-frequency noise, there can be observed a drift tendency with respect to the analytical solution, both in velocity and in displacement records, for all numerical integration schemes. This behaviour has been already reported by many authors [2, 15, 17] when a strong ground motion is integrated twice for obtaining the ground displacement. The numerical phenomenon observed in this example is in agreement with the recommendation postulated by Wang et al. [15] of applying low-cut filtering to a record before integrating it to reduce the effects of baseline offsets.Figure 4
Velocity (a) and displacement (b) records obtained for the three types of contaminated signal computed by different integration schemes (several numerical solutions are overlapping).
(a)
(b)For analyzing in detail the sensitivity of each integration algorithm to the noise content of a particular signal, in Figure5 the difference between the result computed by each algorithm for every contaminated signal with respect to the solution obtained by each algorithm for the original signal without noise has been plotted, both in velocity and in displacement. From this figure, the following observations can be drawn:(i)
All integration algorithms are sensitive to high-frequency noise content. NCn=4, Simpson, Tick, and Trapezoidal rules perceive the phase of the noise although with different magnitud between them, whereas Schuessler-Ibler, Rectangule, and especialy Simpson 3-8 detect the noise with a different value of phase. All methods display a drift in the displacement result although the most sensitive ones are NC n=4, which also shows the highest differences in amplitude in the velocity records, and Simpson 3-8, which also shows a phase absolutely different to the rest of schemes in the corresponding velocity records.(ii)
With respect to the case of low-frequency noise, integration algorithms do not exhibit a particular sensitivity to this type of noise, neither in velocity nor in displacement, since all of them are overlapping with the exact solution (Figure5). As mentioned above, for this type of noise an important drift is observed in the computed result, which could probably be induced by the presence of the low-frequency noise but does not depend on the integration schemes employed.(iii)
Finally, as concerns the case of phase-frequency noise, most of the integration algorithms show a very similar trend among them, except Schuessler-Ibler rule. All of them detect the phase-noise included in the signal although with a little variation in the amplitude compared with the exact result. Again, in most of the integration schemes, a slight drift is also induced in the displacement record, lower for the rectangle rule and exaggerated for Schuessler-Ibler method.Figure 5
Differences between computed results for signals with and without noise obtained from each algorithm. Both velocity (a) and displacement (b) results for the analytical and numerical integration are shown (several numerical solutions are overlapping).
(a)
(b)


### 5.3.  Example  3: Bogdanoff’s Accelerogram

In this example, the velocity records computed from Bogdanoff’s artificial accelerogram [18] are analyzed through Fourier amplitude spectra. This accelerogram is defined as the amplitude-modulated sum of n harmonic terms with uniformly distributed random phase angles, which allows us to manage the frequency content of the accelerogram. Moreover, it becomes possible to know the analytical solution of the integration and therefore to calibrate the accuracy of each numerical integration algorithm.Bogdanoff’s accelerogram is defined by (28), where B and α are two constants equal to 0.292 and −0.333, respectively, n is the number of terms of the summation, ωj and ϕj are the frequency and phase angles of each harmonic term, t is the time in seconds, and a is the acceleration, expressed in m/sec2. Consider(28)at=Bteαt∑j=1ncos⁡ωjt+ϕj.For this example, the input record has a duration of 20 s and it is sampled with a time interval equal to 0.01 s. The accelerogram is composed by 22 terms with frequencies ranging from 6 rad/s (0.95 Hz) to 80.25 rad/s (12.77 Hz), and hence the parameterλ varies between 0.0191 and 0.2554. The resulting accelerogram is plotted in Figure 6(a) and its frequency content can be visualized by means of the Fourier amplitude diagram shown in Figure 6(b).Figure 6
Artificially generated Bogdanoff’s accelerogram and corresponding Fourier amplitude spectrum.
(a)
(b)The analytical and numerical integration velocity records of this accelerogram, and also the differences between both, are plotted in Figure7 for each numerical scheme. As can be observed, all numerical solutions are very close to the analytical one, especially those obtained using the Simpson and Tick rules, while Simpson 3-8 and NC n=4 truncate the peaks of the signal and the Schuessler-Ibler record is slightly shifted with respect to the analytical solution. The maximum positive and negative velocity computed by each algorithm are shown in Table 2, and they are compared with their analytical values to calculate the relative error of each method (expressed in percentage). A positive error means that the corresponding numerical algorithm amplifies the numerical solution and vice versa. For this example, while Simpson and Tick show the lowest error, other algorithms like NC n=4 have lower error for vmax- case than for the vmax+ case. By contrast, the trapezoidal rule always deamplifies the output whereas the rectangle and Schuessler-Ibler rules amplify in one case and deamplify in the other. Finally, the Simpson 3-8 gives the highest error values for this example.Table 2
Maximum positive and negative velocities computed analytically and numerically for Bogdanoff’s accelerogram, and corresponding relative errors expressed in percentage.
v
max
⁡
+ (m/s)
Relative error (%)
v
max
⁡
- (m/s)
Relative error (%)
Analytical
0.1501
—
-
0.1811
—
Trapezoidal
0.1493
-
0.5543
-
0.1806
-
0.2823
Simpson
0.1502
0.0141
-
0.1811
0.0058
Simpson 3-8
0.1456
-
3.0392
-
0.1782
-
1.6175
Newton-Cotesn=4
0.1455
-
3.0666
-
0.1811
-
0.0016
Rectangle
0.1513
0.8061
-
0.1809
-
0.1429
Tick
0.15
-
0.0714
-
0.1811
-
0.0375
Schuessler-Ibler
0.1511
0.6599
-
0.1801
-
0.5667Figure 7
Analytical and numerical solutions for the first integration of Bogdanoff’s accelerogram computed by each scheme. Comparison between both types of velocity records (a) and numerical integration errors ((b) for plotting purpose different scales have been used).To further understand the frequency performance of each algorithm after integration of Bogdanoff’s accelerogram, the corresponding Fourier amplitude spectra for the numerical velocity records have been computed, and they have been subtracted from the analytical Fourier amplitude spectrum. The resulting deviation Fourier amplitude spectra are plotted in Figure8, where Fourier amplitude spectra have been previously divided by the number of intervals of each record. According to this figure, the trapezoidal rule shows a clearly increasing deamplification for upper values of frequency, that is, for higher λ. By contrast, rectangle formula has the inverse tendency, although the magnitude of the deviation is smaller than in the trapezoidal case. The strong influence of rectangle-phase in its own numerical solution must be taken into account. On the other hand, Simpson and Tick’s formulae are practically close to zero throughout the whole range of frequency. Newton-Cotes n=4 and Schuesller-Ibler schemes do not display a clear trend since both amplify and deamplify the Fourier amplitude spectra. As for the Simpson 3-8 rule, it has an evident periodic deviation of high magnitude and, as indicated above, this formula gives the highest relative error for this example. Notice that these three algorithms have also periodic phase diagrams.Figure 8
Deviations of numerical from analytical Fourier amplitude spectra for first integration records of Bogdanoff’s accelerogram (Simpson 3-8 curve has been cut off for plotting purposes).


### 5.4. Example  4: Newmark Sliding Block Analysis

The influence of each numerical algorithm is directly analyzed over a practical earthquake engineering problem, precisely the well-known Newmark sliding block method [17]. This method gives an estimation for the permanent displacement of a slope under a seismic ground motion. When the ground acceleration exceeds a yield acceleration threshold, the block starts to move and slides until the velocity (first integration of the acceleration) becomes zero. To compute the relative displacement, it is necessary to integrate the corresponding relative velocity record (second integration of the input acceleration). The permanent slope displacement depends then on the relationship between the yield acceleration (ay) and the maximum acceleration of the ground motion (amax): if ay is greater than amax no displacement will occur, whereas if ay decreases the slope displacement increases.In this example, two strong ground motions with different frequency contents have been employed, namely, the E-W components of the Gilroy number 1, recorded at Bedrock, and the Gilroy number 2, recorded at Soil. The acceleration time-histories and the corresponding Fourier amplitude spectra are plotted in Figure9. It can be observed that the Gilroy number 1 (Bedrock) has a high-frequency content higher than the Gilroy number 2 (Soil), the mean periods (Tm) being 0.387 s and 0.735 s, respectively [19]. These records have been sampled with a time interval equal to 0.005 s, and the maximum accelerations are 0.4731 g in Gilroy number 1 (Bedrock) and 0.3223 g in Gilroy number 2 (Soil). Moreover, three values of yield acceleration have been employed for each accelerogram, specifically 1/4, 1/2, and 3/4 of amax. The higher ay, the higher influence of the proper integration of the peaks of the record in the computed permanent displacement.Figure 9
Acceleration records and corresponding Fourier amplitude spectra for the E-W components of the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil). The three values of yield accelerationay are indicated in the plots.The permanent slope displacements for each accelerogram and for each value ofay computed by different algorithms are listed in Table 3. For comparison purposes, though no exact solution is available, the numerical deviation of each integration formula has been calculated with respect to the trapezoidal scheme, since this is the usual algorithm in the normal practice. These results (in percentage) are also included in Table 3. From these results it can be concluded that the deviations between schemes are significantly higher in the case of Gilroy number 1 (Bedrock), which has a richer high-frequency content, compared to Gilroy number 2 (Soil). Moreover, normally deviations are higher for the case ay=3/4amax where the peaks of the records have a greater influence than for the other yield accelerations values.Table 3
Permanent slope displacements computed by Newmark sliding method for the strong ground motion Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil), applying different integration formulae and several values ofay. The corresponding deviations (in percentage) with respect to the trapezoidal rule are also listed for comparison purposes.
a
y
=
1
4
a
max
⁡
a
y
=
1
2
a
max
⁡
a
y
=
3
4
a
max
⁡
d (cm)
Deviation (%)
d (cm)
Deviation (%)
d (cm)
Deviation (%)
Gilroy number 1 (Bedrock)
Trapezoidal
11.5954
—
3.7081
—
0.663
—
Simpson
11.597
0.0143
3.6869
-
0.5725
0.6861
3.4829
Simpson 3-8
11.6022
0.0585
3.7256
0.4719
0.622
-
6.185
Newton-Cotesn=4
11.3222
-
2.3556
3.6338
-
2.0041
0.466
-
29.7149
Rectangle
11.6881
0.7998
3.7432
0.9445
0.6769
2.1014
Tick
11.5924
-
0.0256
3.6832
-
0.6722
0.6837
3.1203
Schuessler-Ibler
11.7432
1.2745
3.7292
0.5693
0.5402
-
18.5223
Gilroy number 2 (Soil)
Trapezoidal
17.0186
—
2.8038
—
0.3117
—
Simpson
17.0514
0.1932
2.8074
0.126
0.312
0.0936
Simpson 3-8
17.0043
-
0.084
2.7929
-
0.3912
0.301
-
3.4428
Newton-Cotesn=4
17.3194
1.768
2.7726
-
1.1153
0.2968
-
4.7808
Rectangle
17.0362
0.1038
2.812
0.2921
0.3142
0.8138
Tick
17.0496
0.1826
2.8063
0.089
0.3119
0.0596
Schuessler-Ibler
17.4455
2.5086
2.8077
0.1364
0.3119
0.0585With regard to integration formulae, it can be concluded that, normally, trapezoidal rule gives lower permanent displacements while Simpson and Tick are always very close ones to each other and give values greater than trapezoidal formula. On the other hand, Newton-Cotesn=4 yields the highest deviation results in any case compared to the rest of algorithms, followed by the Schuessler-Ibler scheme.In spite of the fact that the time interval is constant for all cases and equal to the time discretization of the accelerograms, the influence of the frequency content of the input signal, which is related toh by means of λ, is very relevant in the accuracy of the integration results, as well as the proper selection of the integration algorithm.


### 5.5.  Example  5: Error in Displacement Record from Integration of an Accelerogram

In this example, the results obtained by each integration algorithm are compared for computing the velocity and displacement histories after integration of a particular accelerogram [17]. The Gilroy strong ground motions employed in Example 4 have been also used for this goal. These accelerograms, with different frequency content, have been well characterized both in time and in frequency domain in Figure 9. In both cases the duration of the accelerogram is 40 seconds and the sampling frequency is h=0.005 sec.Velocity records are obtained after one integration, where no differences between schemes have been observed. On the contrary, the most relevant differences appear for the displacement histories, Figure10. A drift can be observed in displacement when an accelerogram is integrated twice, and the reasons why this phenomenon occurs are currently under research [2, 15, 17]. However, the focus of this example is not to clarify the reasons for this phenomenon but to analyze the influence of each integration algorithm in the computed result. The errors in displacement (drift) for the final time t=40 sec obtained for each algorithm and for both accelerograms are summarized in Table 4. In general, better results (closer to zero) are obtained for the case of Gilroy number 2 (Soil) than for the case of Gilroy number 1 (Bedrock), which has a high-frequency content that is higher. For these particular accelerograms, minimal differences can be observed between algorithms, except for the Schuessler-Ibler rule, although the results obtained by Trapezoidal, Simpson, and Tick rules are better than for the rest of methods. As remarked above, this is an issue in continuous development. Therefore a deeper research is neccesary and the above conclusions should be taken with caution.Table 4
Displacement in cm for timet=40 sec, computed after double integration of the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil) using different integration schemes.
Gilroy number 1 (Bedrock)
Gilroy number 2(Soil)
Trapezoidal
0.00167
-
0.00035
Simpson
0.0168
0.00570
Simpson 3-8
-
0.028
0.06352
Newton-Cotesn=4
-
0.408
-
0.00049
Rectangle
-
0.219
-
0.10834
Tick
0.0146
0.00479
Schuessler-Ibler
1.485
1.17585Figure 10
Displacement records obtained after integrating twice the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil) using different integration algorithms (several numerical solutions are overlapping).


## 5.1. Example  1: Variable Frequency Input

In order to assess the accuracy of the integration schemes using the proposed frequency approach, a sine function with variable frequency is considered input signal (acceleration) to be integrated. The acceleration is defined in (25), where Ac is the amplitude equal to 3 m/s2, ω is the frequency of the input signal, and t denotes time, whereas zero initial conditions have been assumed. The acceleration has a duration of 1.5 s and the time interval (h) is equal to 0.005 s, meaning that the input frequency ranges from 2π rad/s to 200π rad/s reaching a value of λ between 0.01 and 1, respectively. The input signal is integrated twice to calculate the velocity and displacement records, and after that the maximum velocity (vmax) and the maximum displacement (dmax) are calculated and compared with the analytical solutions which can be obtained from (26) and (27). Consider(25)at=Acsin⁡ωt,(26)vt=∫0tAcsin⁡ωtdt=Acω1-cos⁡ωt,(27)dt=∫0tAcω1-cos⁡ωtdt=Acωt-sin⁡ωtω.Figure2 shows the relative error (in percentage) for both |vmax| and |dmax|, for the whole range of values of λ. A positive relative error means that the corresponding numerical algorithm amplifies the numerical solution and vice versa. From Figure 2, it is concluded that the ratios R1 and R2, plotted in Figure 1, properly predict the performance of each numerical integration algorithm, since the error curves shown in Figure 2 follow exactly the trend plotted in Figure 1. It means that the higher the value of λ, the higher the error, except for the rectangle and Schuessler-Ibler formulae. In the case of the rectangle rule, R1 and R2 predict an amplification in the numerical response whereas the error plotted in Figure 2 shows a deamplification tendency, which is even closer to the trapezoidal rule for the displacement case. On the other hand, this rule presents a positive linear variation of the phase depending on λ, while the exact phase has a negative constant value. This contradictory phase performance could explain the abnormal response of this algorithm. As concerns the Schuessler-Ibler formula, it deamplifies the output record according to R1 and R2 but in a nonuniform way, since the error curve shows an evident periodic behaviour according to its phase diagram. It must be highlighted that the highest errors take place when the phase curve is equal to zero and the lowest errors take place when the phase reaches any of the extreme values.Figure 2
Relative error (in percentage) reached by each algorithm along the whole range ofλ for a sinusoidal signal (both velocity and displacement records are shown).
(a)
Relative error (%) for|vmax|
(b)
Relative error (%) for|dmax|Although this example has been computed with a constant value ofh, it can be observed that numerical integration accuracy is frequency-dependent, a matter of paramount interest in earthquake engineering problems. Moreover, if this example is computed for other values of time interval h, with λ ranging between 0 and 1, these algorithms show similar trends.


## 5.2.  Example  2: Influence of Noise Content of a Signal

In earthquake engineering problems, it is necessary to employ different types of signal as accelerograms or records measured from the response of a structure. Any kind of record could be contaminated by different noise sources (internal characteristics of the measurement apparatus, the placement where the instrument is located, ambient noise, etc.) [15, 16]. In this example, the sensitivity of the different integration algorithms to the noise content in a signal is analyzed. A sine signal has been modified by three types of noise and then has been integrated twice to compute the corresponding displacement record. Thanks to the simplicity of a sine function and since its analytical solution is known in advance ((26) and (27)), it is possible to verify the influence of the noise components of a signal in the final results for each integration algorithm.The original sine function (labeled as original signal is this research) is defined in (25), where Ac=3 m/s2 and ω=10π rad/s, with a duration of 1.5 s and h=0.005 s (λ=0.05). Three different types of noise have been chosen, high-frequency, low-frequency, and phase-noise. In the three cases, the contaminated signal is defined as a sine function too, which has been added to the original one. The characteristics of the high-frequency noise are Ahigh=0.3 m/s2 and ωhigh=150π rad/s, while for the low-frequency noise are Alow=0.6 m/s2 and ωlow=0.1π rad/s. For the phase-noise case, a constant variation in the phase equal to 0.15π rad/s has been added to the frequency of the original sine function. The four signals considered in this example are plotted in Figure 3; as can be observed, the noise-signals are very close to the original one.Figure 3
Original function and contaminated signals by three types of noise: high-frequency noise, low-frequency noise, and phase-noise.The three contaminated signals from Figure3 have been integrated twice to obtain the velocity and displacement records, Figure 4. In all cases, the analytical result from the original sine function has been included for comparative purposes. Moreover, in all cases, the solutions computed by the Schuessler-Ibler rule show an anomalous behaviour relative to the rest of methods, both in velocity and in displacement. In the cases of high-frequency noise and phase-noise, the computed results are very close to the analytical solution, even when this has been computed for the original signal without any kind of noise. However, in the case of low-frequency noise, there can be observed a drift tendency with respect to the analytical solution, both in velocity and in displacement records, for all numerical integration schemes. This behaviour has been already reported by many authors [2, 15, 17] when a strong ground motion is integrated twice for obtaining the ground displacement. The numerical phenomenon observed in this example is in agreement with the recommendation postulated by Wang et al. [15] of applying low-cut filtering to a record before integrating it to reduce the effects of baseline offsets.Figure 4
Velocity (a) and displacement (b) records obtained for the three types of contaminated signal computed by different integration schemes (several numerical solutions are overlapping).
(a)
(b)For analyzing in detail the sensitivity of each integration algorithm to the noise content of a particular signal, in Figure5 the difference between the result computed by each algorithm for every contaminated signal with respect to the solution obtained by each algorithm for the original signal without noise has been plotted, both in velocity and in displacement. From this figure, the following observations can be drawn:(i)
All integration algorithms are sensitive to high-frequency noise content. NCn=4, Simpson, Tick, and Trapezoidal rules perceive the phase of the noise although with different magnitud between them, whereas Schuessler-Ibler, Rectangule, and especialy Simpson 3-8 detect the noise with a different value of phase. All methods display a drift in the displacement result although the most sensitive ones are NC n=4, which also shows the highest differences in amplitude in the velocity records, and Simpson 3-8, which also shows a phase absolutely different to the rest of schemes in the corresponding velocity records.(ii)
With respect to the case of low-frequency noise, integration algorithms do not exhibit a particular sensitivity to this type of noise, neither in velocity nor in displacement, since all of them are overlapping with the exact solution (Figure5). As mentioned above, for this type of noise an important drift is observed in the computed result, which could probably be induced by the presence of the low-frequency noise but does not depend on the integration schemes employed.(iii)
Finally, as concerns the case of phase-frequency noise, most of the integration algorithms show a very similar trend among them, except Schuessler-Ibler rule. All of them detect the phase-noise included in the signal although with a little variation in the amplitude compared with the exact result. Again, in most of the integration schemes, a slight drift is also induced in the displacement record, lower for the rectangle rule and exaggerated for Schuessler-Ibler method.Figure 5
Differences between computed results for signals with and without noise obtained from each algorithm. Both velocity (a) and displacement (b) results for the analytical and numerical integration are shown (several numerical solutions are overlapping).
(a)
(b)


## 5.3.  Example  3: Bogdanoff’s Accelerogram

In this example, the velocity records computed from Bogdanoff’s artificial accelerogram [18] are analyzed through Fourier amplitude spectra. This accelerogram is defined as the amplitude-modulated sum of n harmonic terms with uniformly distributed random phase angles, which allows us to manage the frequency content of the accelerogram. Moreover, it becomes possible to know the analytical solution of the integration and therefore to calibrate the accuracy of each numerical integration algorithm.Bogdanoff’s accelerogram is defined by (28), where B and α are two constants equal to 0.292 and −0.333, respectively, n is the number of terms of the summation, ωj and ϕj are the frequency and phase angles of each harmonic term, t is the time in seconds, and a is the acceleration, expressed in m/sec2. Consider(28)at=Bteαt∑j=1ncos⁡ωjt+ϕj.For this example, the input record has a duration of 20 s and it is sampled with a time interval equal to 0.01 s. The accelerogram is composed by 22 terms with frequencies ranging from 6 rad/s (0.95 Hz) to 80.25 rad/s (12.77 Hz), and hence the parameterλ varies between 0.0191 and 0.2554. The resulting accelerogram is plotted in Figure 6(a) and its frequency content can be visualized by means of the Fourier amplitude diagram shown in Figure 6(b).Figure 6
Artificially generated Bogdanoff’s accelerogram and corresponding Fourier amplitude spectrum.
(a)
(b)The analytical and numerical integration velocity records of this accelerogram, and also the differences between both, are plotted in Figure7 for each numerical scheme. As can be observed, all numerical solutions are very close to the analytical one, especially those obtained using the Simpson and Tick rules, while Simpson 3-8 and NC n=4 truncate the peaks of the signal and the Schuessler-Ibler record is slightly shifted with respect to the analytical solution. The maximum positive and negative velocity computed by each algorithm are shown in Table 2, and they are compared with their analytical values to calculate the relative error of each method (expressed in percentage). A positive error means that the corresponding numerical algorithm amplifies the numerical solution and vice versa. For this example, while Simpson and Tick show the lowest error, other algorithms like NC n=4 have lower error for vmax- case than for the vmax+ case. By contrast, the trapezoidal rule always deamplifies the output whereas the rectangle and Schuessler-Ibler rules amplify in one case and deamplify in the other. Finally, the Simpson 3-8 gives the highest error values for this example.Table 2
Maximum positive and negative velocities computed analytically and numerically for Bogdanoff’s accelerogram, and corresponding relative errors expressed in percentage.
v
max
⁡
+ (m/s)
Relative error (%)
v
max
⁡
- (m/s)
Relative error (%)
Analytical
0.1501
—
-
0.1811
—
Trapezoidal
0.1493
-
0.5543
-
0.1806
-
0.2823
Simpson
0.1502
0.0141
-
0.1811
0.0058
Simpson 3-8
0.1456
-
3.0392
-
0.1782
-
1.6175
Newton-Cotesn=4
0.1455
-
3.0666
-
0.1811
-
0.0016
Rectangle
0.1513
0.8061
-
0.1809
-
0.1429
Tick
0.15
-
0.0714
-
0.1811
-
0.0375
Schuessler-Ibler
0.1511
0.6599
-
0.1801
-
0.5667Figure 7
Analytical and numerical solutions for the first integration of Bogdanoff’s accelerogram computed by each scheme. Comparison between both types of velocity records (a) and numerical integration errors ((b) for plotting purpose different scales have been used).To further understand the frequency performance of each algorithm after integration of Bogdanoff’s accelerogram, the corresponding Fourier amplitude spectra for the numerical velocity records have been computed, and they have been subtracted from the analytical Fourier amplitude spectrum. The resulting deviation Fourier amplitude spectra are plotted in Figure8, where Fourier amplitude spectra have been previously divided by the number of intervals of each record. According to this figure, the trapezoidal rule shows a clearly increasing deamplification for upper values of frequency, that is, for higher λ. By contrast, rectangle formula has the inverse tendency, although the magnitude of the deviation is smaller than in the trapezoidal case. The strong influence of rectangle-phase in its own numerical solution must be taken into account. On the other hand, Simpson and Tick’s formulae are practically close to zero throughout the whole range of frequency. Newton-Cotes n=4 and Schuesller-Ibler schemes do not display a clear trend since both amplify and deamplify the Fourier amplitude spectra. As for the Simpson 3-8 rule, it has an evident periodic deviation of high magnitude and, as indicated above, this formula gives the highest relative error for this example. Notice that these three algorithms have also periodic phase diagrams.Figure 8
Deviations of numerical from analytical Fourier amplitude spectra for first integration records of Bogdanoff’s accelerogram (Simpson 3-8 curve has been cut off for plotting purposes).


## 5.4. Example  4: Newmark Sliding Block Analysis

The influence of each numerical algorithm is directly analyzed over a practical earthquake engineering problem, precisely the well-known Newmark sliding block method [17]. This method gives an estimation for the permanent displacement of a slope under a seismic ground motion. When the ground acceleration exceeds a yield acceleration threshold, the block starts to move and slides until the velocity (first integration of the acceleration) becomes zero. To compute the relative displacement, it is necessary to integrate the corresponding relative velocity record (second integration of the input acceleration). The permanent slope displacement depends then on the relationship between the yield acceleration (ay) and the maximum acceleration of the ground motion (amax): if ay is greater than amax no displacement will occur, whereas if ay decreases the slope displacement increases.In this example, two strong ground motions with different frequency contents have been employed, namely, the E-W components of the Gilroy number 1, recorded at Bedrock, and the Gilroy number 2, recorded at Soil. The acceleration time-histories and the corresponding Fourier amplitude spectra are plotted in Figure9. It can be observed that the Gilroy number 1 (Bedrock) has a high-frequency content higher than the Gilroy number 2 (Soil), the mean periods (Tm) being 0.387 s and 0.735 s, respectively [19]. These records have been sampled with a time interval equal to 0.005 s, and the maximum accelerations are 0.4731 g in Gilroy number 1 (Bedrock) and 0.3223 g in Gilroy number 2 (Soil). Moreover, three values of yield acceleration have been employed for each accelerogram, specifically 1/4, 1/2, and 3/4 of amax. The higher ay, the higher influence of the proper integration of the peaks of the record in the computed permanent displacement.Figure 9
Acceleration records and corresponding Fourier amplitude spectra for the E-W components of the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil). The three values of yield accelerationay are indicated in the plots.The permanent slope displacements for each accelerogram and for each value ofay computed by different algorithms are listed in Table 3. For comparison purposes, though no exact solution is available, the numerical deviation of each integration formula has been calculated with respect to the trapezoidal scheme, since this is the usual algorithm in the normal practice. These results (in percentage) are also included in Table 3. From these results it can be concluded that the deviations between schemes are significantly higher in the case of Gilroy number 1 (Bedrock), which has a richer high-frequency content, compared to Gilroy number 2 (Soil). Moreover, normally deviations are higher for the case ay=3/4amax where the peaks of the records have a greater influence than for the other yield accelerations values.Table 3
Permanent slope displacements computed by Newmark sliding method for the strong ground motion Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil), applying different integration formulae and several values ofay. The corresponding deviations (in percentage) with respect to the trapezoidal rule are also listed for comparison purposes.
a
y
=
1
4
a
max
⁡
a
y
=
1
2
a
max
⁡
a
y
=
3
4
a
max
⁡
d (cm)
Deviation (%)
d (cm)
Deviation (%)
d (cm)
Deviation (%)
Gilroy number 1 (Bedrock)
Trapezoidal
11.5954
—
3.7081
—
0.663
—
Simpson
11.597
0.0143
3.6869
-
0.5725
0.6861
3.4829
Simpson 3-8
11.6022
0.0585
3.7256
0.4719
0.622
-
6.185
Newton-Cotesn=4
11.3222
-
2.3556
3.6338
-
2.0041
0.466
-
29.7149
Rectangle
11.6881
0.7998
3.7432
0.9445
0.6769
2.1014
Tick
11.5924
-
0.0256
3.6832
-
0.6722
0.6837
3.1203
Schuessler-Ibler
11.7432
1.2745
3.7292
0.5693
0.5402
-
18.5223
Gilroy number 2 (Soil)
Trapezoidal
17.0186
—
2.8038
—
0.3117
—
Simpson
17.0514
0.1932
2.8074
0.126
0.312
0.0936
Simpson 3-8
17.0043
-
0.084
2.7929
-
0.3912
0.301
-
3.4428
Newton-Cotesn=4
17.3194
1.768
2.7726
-
1.1153
0.2968
-
4.7808
Rectangle
17.0362
0.1038
2.812
0.2921
0.3142
0.8138
Tick
17.0496
0.1826
2.8063
0.089
0.3119
0.0596
Schuessler-Ibler
17.4455
2.5086
2.8077
0.1364
0.3119
0.0585With regard to integration formulae, it can be concluded that, normally, trapezoidal rule gives lower permanent displacements while Simpson and Tick are always very close ones to each other and give values greater than trapezoidal formula. On the other hand, Newton-Cotesn=4 yields the highest deviation results in any case compared to the rest of algorithms, followed by the Schuessler-Ibler scheme.In spite of the fact that the time interval is constant for all cases and equal to the time discretization of the accelerograms, the influence of the frequency content of the input signal, which is related toh by means of λ, is very relevant in the accuracy of the integration results, as well as the proper selection of the integration algorithm.


## 5.5.  Example  5: Error in Displacement Record from Integration of an Accelerogram

In this example, the results obtained by each integration algorithm are compared for computing the velocity and displacement histories after integration of a particular accelerogram [17]. The Gilroy strong ground motions employed in Example 4 have been also used for this goal. These accelerograms, with different frequency content, have been well characterized both in time and in frequency domain in Figure 9. In both cases the duration of the accelerogram is 40 seconds and the sampling frequency is h=0.005 sec.Velocity records are obtained after one integration, where no differences between schemes have been observed. On the contrary, the most relevant differences appear for the displacement histories, Figure10. A drift can be observed in displacement when an accelerogram is integrated twice, and the reasons why this phenomenon occurs are currently under research [2, 15, 17]. However, the focus of this example is not to clarify the reasons for this phenomenon but to analyze the influence of each integration algorithm in the computed result. The errors in displacement (drift) for the final time t=40 sec obtained for each algorithm and for both accelerograms are summarized in Table 4. In general, better results (closer to zero) are obtained for the case of Gilroy number 2 (Soil) than for the case of Gilroy number 1 (Bedrock), which has a high-frequency content that is higher. For these particular accelerograms, minimal differences can be observed between algorithms, except for the Schuessler-Ibler rule, although the results obtained by Trapezoidal, Simpson, and Tick rules are better than for the rest of methods. As remarked above, this is an issue in continuous development. Therefore a deeper research is neccesary and the above conclusions should be taken with caution.Table 4
Displacement in cm for timet=40 sec, computed after double integration of the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil) using different integration schemes.
Gilroy number 1 (Bedrock)
Gilroy number 2(Soil)
Trapezoidal
0.00167
-
0.00035
Simpson
0.0168
0.00570
Simpson 3-8
-
0.028
0.06352
Newton-Cotesn=4
-
0.408
-
0.00049
Rectangle
-
0.219
-
0.10834
Tick
0.0146
0.00479
Schuessler-Ibler
1.485
1.17585Figure 10
Displacement records obtained after integrating twice the strong ground motions Gilroy number 1 (Bedrock) and Gilroy number 2 (Soil) using different integration algorithms (several numerical solutions are overlapping).


## 6. Conclusions

In this research the frequency performance of seven numerical integration algorithms which behave as recursive digital filters has been explored. The methodology employed is based on comparing, for each algorithm, the analytical and numerical transfer functions. For the sake of simplicity, the ratiosR1, for single integration, and R2, for double integration, have been defined for λ ranging between 0 and 1, where the parameter λ links the time interval h with the input frequency ω through Nyquist’s criterion. In general, the error of any algorithm increases for higher values of λ, but it must be highlighted that λ depends on both the time interval and the frequency content of the input signal.In this frequency-domain approach, both amplitude and phase components have been taken into account, as opposed to the normal practice where phase is scarcely investigated. This methodology has been validated by means of five examples from the field of earthquake engineering. From this research, the following general conclusions can be drawn:(i)
The numerical integration methods investigated in this research are clearly frequency-dependent, since no one of them maintains the ratioR1 close to one for the whole range of λ, displaying in all cases an abnormal behaviour for λ higher than 1/2, whereas the accuracy of each scheme improves for lower values of λ.(ii)
In general, the trapezoidal rule tends to deamplify the solution along the whole range ofλ. On the contrary, Simpson’s formula tends to amplify the results, while Tick’s formula slightly deamplifies for intermediate values of λ and tends to amplify for λ>1/2. Hence, both algorithms should be used with caution if the input signal has a relevant high-frequency content, especially if it is due to noise, and they should be avoided if high-frequencies are wanted to be damped out. Normally the results computed by Simpson’s formula are very close to Tick’s algorithm, although the latter is more accurate than the former. No one of these algorithms exhibits phase problems.(iii)
Algorithms like Simpson 3-8 and Newton-Cotesn=4, in spite of being defined with a higher degree of interpolation function, give worse results than the previous ones and show an evident undesirable periodic phase dependency.(iv)
The rest of algorithms, rectangle and Schuessler-Ibler, not always follow the corresponding trends predicted byR1 and R2 and sometimes give acceptable results while others do not. Probably, these unexpected behaviours are due to their particular phase diagrams. It must be highlighted that the accuracy of Schuesller-Ibler’s formula oscillates according to its own phase behaviour.Unfortunately there is no general algorithm which can be used successfully for all cases, at least from the list of schemes analyzed in this research. From the cases analyzed, the Trapezoidal, Simpson, and Tick rules do not exhibit anomalous behaviours and generally give acceptable results, whereas the Schuessler-Ibler rule shows strange tendencies in some of the examples proposed, so it is not recommended to use it. The clues to choose a particular integration formula depend on the frequency content of the input signal, the computational cost due to the size of the interval of integration, the required accuracy, and mainly the further use of the results. Careful attention should be paid to the proper selection of the integration algorithm, not only from the accuracy point of view but from frequency-dependent behaviour of the algorithm.



---

*Source: 259530-2015-07-16.xml*