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.
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 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 () 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 integrated (1) by an interpolation polynomial function of degree , through known values of it, , which are defined at the integration points placed at equal time interval, :
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 of interpolation polynomial function [1] and they are as follows: trapezoidal rule , Simpson’s rule , Simpson three-eights rule , and Newton-Cotes scheme with , which is labeled as NC 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: (ii) Simpson’s rule: (iii) Simpson three-eights rule: (iv) Newton-Cotes : (v) Rectangle rule: (vi) Tick’s rule: where and . It must be highlighted that if and , Tick’s rule (7) becomes identical to Simpson’s rule (3).(vii) Schuessler-Ibler rule: where and and and .
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, , where is the angular frequency of the signal and 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 , where is the complex transfer function (TF) of the corresponding integration scheme which depends on and [5]. Also, since 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 signal , 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 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:
Now if (9) is integrated again, it is possible to obtain the exact transfer function for the double integral, , as developed in (11):
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 is [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: (ii) Simpson’s rule: (iii) Simpson three-eights rule: (iv) Newton-Cotes : (v) Rectangle rule: (vi) Tick’s rule: where and .(vii) Schuessler-Ibler rule: where and and and .
The procedure followed in the exact case to obtain the transfer function for the double integral (11) is also analogous for the numerical case (19):
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 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 , as shown in (20), and for the double integral , as shown in (21). In both cases, the closer the ratio to unity, the more accurate the numerical solution:
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 () used to sample the signal. Owing to that, the maximum frequency on the input signal ( or ) must be less than Nyquist’s frequency ( or ) (22) which implies that (23) must be guaranteed:
In this research, the most common values for time intervals employed for strong ground motion discretization (, 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 , the higher the allowed input frequency. For simplicity, the dimensionless parameter , defined as the quotient between the input frequency () and Nyquist’s frequency (), is adopted as shown in (24)
Figure 1 shows the variations with respect to of the ratios (20) and (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 and as close to unity as possible for a large range of , and also the phase value should coincide with the exact one.
(a) : single integration
(b) : double integration
(c) Phase diagram
Since the transfer function is defined by a complex expression, it is convenient to analyze not only its amplitude but also its phase diagram. From Figure 1, 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 and far away from unity for the upper-half of . This effect is even more critical for the 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 rules display a singular behaviour: depending on , they can amplify or deamplify the output signal.(3)For the single integral case (), Tick’s rule slightly deamplifies the output signal for intermediate values of , while it tends to amplify from . On the other hand, Simpson’s rule consistently amplifies the output signal especially for above . 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 . For Simpson 3-8 rule, is close to unity until , tending to amplify the signal for higher values of (with a discontinuity around ) until reaching values around 0.72, deamplifying significantly after that. Finally, the Newton-Cotes formula is close to unity up to values of around 0.4, deamplifying quickly for between 0.4 and 0.5. For 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 (), all methods show the same performance (amplification or deamplification) compared to the single integration (), but, in most of the formulae, the ratio is further from unity for lower values of than , except for Tick and Schuessler-Ibler’s rules. In general, this implies that, given a signal and a value of , 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 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 alternatively display phase values equal to or equal to . However the phase of rectangle rule linearly varies from for to for . The phase of Schuessler-Ibler’s rule periodically varies along the whole range of , and in a linear way between and 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 is the amplitude equal to m/s2, is the frequency of the input signal, and denotes time, whereas zero initial conditions have been assumed. The acceleration has a duration of s and the time interval () is equal to 0.005 s, meaning that the input frequency ranges from rad/s to 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 () and the maximum displacement () are calculated and compared with the analytical solutions which can be obtained from (26) and (27). Consider
Figure 2 shows the relative error (in percentage) for both and , 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 and , 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, and 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 and 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.
(a) Relative error (%) for
(b) Relative error (%) for
Although this example has been computed with a constant value of , 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 , 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 m/s2 and rad/s, with a duration of s and s (). 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 m/s2 and rad/s, while for the low-frequency noise are m/s2 and rad/s. For the phase-noise case, a constant variation in the phase equal to 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.
The three contaminated signals from Figure 3 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.
(a)
(b)
For analyzing in detail the sensitivity of each integration algorithm to the noise content of a particular signal, in Figure 5 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. NC , 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 , 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 (Figure 5). 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.
(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 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 and are two constants equal to 0.292 and −0.333, respectively, is the number of terms of the summation, and are the frequency and phase angles of each harmonic term, is the time in seconds, and is the acceleration, expressed in m/sec2. Consider
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).
(a)
(b)
The analytical and numerical integration velocity records of this accelerogram, and also the differences between both, are plotted in Figure 7 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 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 have lower error for case than for the 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.
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 Figure 8, 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 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.
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 () and the maximum acceleration of the ground motion (): if is greater than no displacement will occur, whereas if 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 Figure 9. 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 () 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 . The higher , the higher influence of the proper integration of the peaks of the record in the computed permanent displacement.
The permanent slope displacements for each accelerogram and for each value of 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 where the peaks of the records have a greater influence than for the other yield accelerations values.
With 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-Cotes 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 to 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 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 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, Figure 10. 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 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.
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 ratios , for single integration, and , for double integration, have been defined for ranging between 0 and 1, where the parameter links the time interval 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 ratio 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 . 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-Cotes , 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 by and 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.
Appendix
If it is assumed that the function to be integrated is , where is the angular frequency of the input signal and is the time, the integral can be numerically calculated following the definition of each quadrature rule. Consequently, the corresponding transfer function of each numerical integration algorithm can be obtained as follows.
(i) Trapezoidal Rule. The trapezoidal formula is defined as
According to the transfer function definition, every term of (A.1) can be expressed as
Substituting the above equations into (A.1), it is obtained that
(ii) Simpson’s Rule. Simpson’s formula is defined as
According to the transfer function definition, every term of (A.4) can be expressed as
Substituting the above equations into (A.4), it is obtained that
(iii) Simpson Three-Eights Points. The Simpson 3-8 formula is defined as
According to the transfer function definition, every term of (A.7) can be expressed as
Substituting the above equations into (A.7), it is obtained that
(iv) Newton-Cotes The Newton-Cotes formula is defined as
According to the transfer function definition, every term of (A.10) can be expressed as
Substituting the above equations into (A.10), it is obtained that
(v) Rectangle Rule. The rectangle method is defined as
According to the transfer function definition, every term of (A.13) can be expressed as
Substituting the above equations into (A.13), it is obtained that
(vi) Tick’s Rule. Tick’s method is defined aswhere the parameters are and .
According to the transfer function definition, every term of (A.16) can be expressed as
Substituting the above equations into (A.16), it is obtained that
(vii) Schuessler-Ibler’s Formula. Schuessler-Ibler’s formula is defined aswhere and and and .
According to the transfer function definition, every term of (A.19) can be expressed as
Substituting the above equations into (A.19), it is obtained that
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This study has been partially funded by Spanish Ministry of Science and Innovation through Project BIA2007-67401 and scholarship BES-2008-002770. The financial support to complete this research is gratefully appreciated by the authors. Finally, the authors would also like to thank the constructive suggestions and comments provided by the anonymous reviewers during the reviewing process.