Development of a tri-axial primary vibration calibration system ACTA IMEKO ISSN: 2221-870X March 2019, Volume 8, Number 1, 33 - 39 ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 33 Development of a tri-axial primary vibration calibration system Zhihua Liu1, Chenguang Cai1, Ming Yang2, Mei Yu1 1 National Institute of Metrology, Bei San Huan Dong Lu 18, 100029 Beijing, China 2 Beijing University of Chemical Technology, Bei San Huan Dong Lu 15, 100029 Beijing, China Section: RESEARCH PAPER Keywords: Primary calibration; Spatial orbit; Multi-exciter control; Interferometry; Tri-axial vibration exciter Citation: Zhihua Liu, Chenguang Cai, Ming Yang, Mei Yu, Development of a tri-axial primary vibration calibration system, Acta IMEKO, vol. 8, no. 1, article 6, March 2019, identifier: IMEKO-ACTA-08(2019)-01-06 Editor: Maija Ojanen-Saloranta, VTT Technical Research Centre of Finland, National Metrology Institute MIKES, Finland Received August 24, 2018; In final form October 25, 2018; Published March 2019 Copyright: © 2019IMEKO. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: This work was supported by Special Fund for Scientific Research in the Public Interest (201410009) and the National Natural Science Foundation of China (No. 51605461). Corresponding author: Zhihua Liu, e-mail: liuzhihua@nim.ac.cn 1. INTRODUCTION Multi-exciter vibration systems have the ability to accurately replicate the real-word realistic dynamic environment of multi- axis motions, which provides access to the simultaneous multicomponent calibration of motion transducers [1], [2]. Following single-axial vibration exciters and planar testing systems, tri-axial vibration exciters that can generate spatial orbit vibration is an important development in the field of transducer calibration [3], [4], [5]. PTB set up a tri-axial vibration standard measuring device for the calibration of vibration transducers under multi-axial excitation [6], [7], [8]. The tri-axial vibration exciter can generate arbitrary motion, the linear combination of motion, along each of its three axes. ISO 16063-31 presents an approach to the testing of transverse vibration sensitivity using a tri-axial vibration exciter [9]. The vibration direction in the X-Y plane is changed every 30 ° in order to evaluate transverse sensitivity [10]. Umeda proposed the sensitivity matrix calibration method of accelerometers using the tri-axial vibration generator and three laser interferometers. The tri-axial vibration exciter was employed to generate three independent motion vectors in order to calibrate all the sensitivities of a tri-axis accelerometer [11], [12], [13]. Qualified spatial orbits, such as linear orbits and circular orbits, play an important role in transducer calibration. Vibration control on a multi-exciter testing system is an essential prerequisite to generate accurate spatial orbit vibration. Multi- exciter vibration control deals not only with amplitude control (as in single-exciter vibration control), but also cross-coupling compensation and phase control [14], [15]. The linear iterative algorithm involving correction of the driving spectrum with the error spectrum and the system impedance was first introduced in [16], [17]. Subsequently, the adaptive iterative algorithm was presented in [18], [19], according to which the system impedance keeps updating for system nonlinearity, and a variable adjustment gain is used for driving spectrum correction. The current vibration control approach concerns more amplitudes and phases for each individual exciter instead of the spatial orbit for ABSTRACT A tri-axial primary vibration calibration system has been set up at National Institute of Metrology for simultaneous calibration of motion transducers. The system is driven by three electrodynamics exciters that are mounted along the three orthogonal axes. The cross- coupling unit based on air bearing has been developed for force transference and motion guiding. Spatial orbit vibration is composited from sine vibration components of the three orthogonal axes. The relationships between the shapes and orientations of the spatial orbits and the amplitudes and phases of the sine vibration components are discussed. Multi-exciter vibration control for both cross- coupling compensation and the amplitude and phase control of the sine vibration components is also investigated. The tri-axial measuring system can simultaneously measure the three orthogonal vibration quantities based on the band-pass sampling method. The experiments show that a variety of spatial orbits can be generated by efficiently reducing the cross-coupling of the tri-axial vibration exciter and the magnitude, and the phase shift of the sensitivities of a tri-axis accelerometer can be determined. mailto:liuzhihua@nim.ac.cn ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 34 calibration purposes. Issues concerning the nature of the relationship between spatial orbits and sine vibration components and how to achieve a certain spatial orbit still require further study. To achieve the multicomponent calibration of vibration transducers, the National Institute of Metrology in China has developed a tri-axial primary vibration calibration system that consists of a cross-coupling unit, a spatial orbit control system, and a primary tri-axial vibration measuring system. In order to reduce costs and improve performance, some new technologies have been used in this system. The paper is organized as follows: In Section 2, a description of the tri-axial primary vibration calibration system is given. In Section 3, spatial orbit generation by means of the tri-axial vibration exciter and multi-exciter vibration control are investigated. In Section 4, experimental investigation of the performance of the tri-axial primary vibration calibration system and calibration use is presented. In the last section, conclusions of this paper are given. 2. SYSTEM DESCRIPTION The tri-axial primary vibration calibration system was designed to generate and measure spatial orbit vibration in the three orthogonal axes. As illustrated in Figure 1, it consists of the tri-axial vibration exciter, the vibration control system, and the primary vibration measuring system. 2.1. The air-bearing cross-coupling unit The tri-axial vibration exciter is driven by three electrodynamic exciters that are mounted along the orthogonal axes relative to the floor. The air-bearing cross-coupling unit was developed for the force transference and motion guiding of the tri-orthogonal axial vibration [20], [21]. As illustrated in Figure 2, the air-bearing cross coupling unit is composed of five sets of air bearings. The X-, Y-, and Z-axial force-transferring air bearings are used to transfer the vibration generated by the electrodynamic exciters. The X- and Y-axial motion-guiding air bearings help to guide the X- and Y-axial motions, respectively. As the transducer under test needs to be mounted on the top of the coupling unit, there is no space to install a Z-axis motion- guiding air bearing. Besides, each of the air bearings do not prevent the motion along the other two orthogonal axial directions. Commercial tri-axial vibration exciters commonly develop cross-coupling units by means of oil film bearings that need additional auxiliary equipment and may also heat the device under test. The air bearings have the advantages of being cheap, easy to maintain, and not thermic; therefore, they have been used in the tri-axial vibration system. The restrictions of the tri-axial vibration exciter are listed in Table 1. 2.2. The primary vibration measuring system The spatial orbit vibration is measured by tri-orthogonal heterodyne interferometers, which are mounted on vibration isolation platforms. The interference signals of the three interferometers and the outputs of the transducer under test are sampled by a six-channel data acquisition system. Usually, there are two sampling methods of dealing with interference signals whose central frequency is about 40 MHz. To fulfil Nyquist’s theorem, a high-speed acquisition card to deal with the reference interference signal with a central frequency of 40 MHz is necessary. In this case, a large amount of data needs to be sampled and processed. The other way is to use an analogue mixer and low-pass filter to down-convert the signal frequency. However, the phase shift and inconsistency of the analogue devices influence the measuring results. As the frequency band of the reference signal of the interferometer in this application is lower than 2 MHz, the band-pass sampling method can be used to directly sample the reference signals with a sample rate that is less than 10 MHz, without any analogy device. According to the procedures presented in ISO 16063-11 [3], three synchronized interference signals are demodulated to realize the primary measurement of the spatial orbit vibration. 2.3. The vibration control system The vibration control system was established based on the PXI system, and the control program was developed using LabView software. Two NI PXI-4461 acquisition cards were used to acquire the tri-axial acceleration feedback signals and to Figure 1. The tri-axial vibration calibration system. Figure 2. The air-bearing cross coupling unit. Table 1. The restrictions of the tri-axial vibration exciter Max. peak force 1000 N Max. pay load 10 kg Max. peak acceleration 20 m/s2 Frequency range 5 Hz to 1.6 kHz Air pressure required 4 bar to 8 bar Cross motion under control < 1% ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 35 generate three channels’ analogue voltage outputs for the exciters’ amplifiers. Multi-exciter vibration control deals not only with amplitude control (as with single-exciter vibration control), but also with cross-coupling compensation and phase control. The transfer function matrix should be determined by pre-exciting the system with multi-random signals in order to calculate the impedance matrix. The desired motion signal, the response, and the driving multi-exciter sine signals are both transformed into spectrum vectors. The driving spectrum is corrected according to the error spectrum and the impedance matrix during the control loop. The impedance matrix also keeps updated for system nonlinearity, and a variable adjustment gain is used for driving spectrum correction. The vibration control algorithm is implemented using the MathScript RT Module. A detailed deduction of the spatial orbit control algorithm is presented in Section 3. 3. SPATIAL ORBIT VIBRATION CONTROL 3.1. Spatial orbit generation Sinusoidal vibrations of the three axes with the same frequency (regardless of amplitudes and phases) are bound to result in a spatial elliptical orbit as depicted in Figure 3. The shape and orientation of the spatial elliptical orbit depend on the sinusoidal vibration amplitudes and phases of the three axes. The orientation of the spatial ellipse can be expressed by a set of three unit vectors that are orthogonal to each other. As shown in Figure 4(b), u and v are parallel with the major and minor axes of the spatial ellipse respectively. w is perpendicular to the spatial ellipse plane. The three unit vectors are as follows:  T321 uuu=u ,   T 321 vvv=v , vuw = . (1) From the theory of coordinate transformation, any spatial ellipse can be obtained from a standard ellipse in the X-Y plane by the rotation of certain angles. Figure 4(a) depicts a standard ellipse whose major and minor axes are parallel with the coordinate axes. We can establish a local coordinate frame with its origin located at the centre of the spatial ellipse and its coordinate axes being parallel with the unit vectors, as shown in Figure 4(b). The coordinate transformation of the two ellipses can be achieved by the rotation matrix R consisting of the three unit vectors.  wvuR = (2) The sinusoidal vibration components of the spatial elliptical orbit can be described in its local coordinate frame as ( ) ( ) ( )                 =           t t ta ta ta Z Y X     sin cos 00 0 0 ~ ~ ~ 2 1 (3) By left multiplying equation (3) by the rotation matrix R , we can calculate the sinusoidal vibration components in the global coordinate frame. ( ) ( ) ( )                           =           t t wvu wvu wvu ta ta ta Z Y X     sin cos 00 0 0 2 1 333 222 111 (4) Let us rewrite equation (4) as ( ) ( ) ( )    += += += tvtuta tvtuta tvtuta Z Y X    sincos sincos sincos 2313 2212 2111 (5) Then, the amplitudes and phases of the sinusoidal vibration components for a given spatial elliptical orbit are given by ( ) ( )221 2 11 ˆ  vua X += ,         −= − 11 211 tan    u v X (6) ( ) ( )222 2 12 ˆ  vuaY += ,         −= − 12 221 tan    u v Y (7) ( ) ( )223 2 13 ˆ  vuaZ += ,         −= − 13 231 tan    u v Z . (8) 3.2. Multi-exciter vibration control In order to acquire the accurate amplitudes and phases of the sinusoidal vibration components, vibration control is performed on the multi-exciters. The linear model of the multi-input multi- output system is given by ( ) ( ) ( )fff DHC = (9) where ( )fD is the driving spectrum, ( )fC is the response spectrum, and ( )fH is the frequency response function. The error spectrum is given by subtracting the response spectrum ( )fC from the reference spectrum ( )fR , as ( )( ) ( ) ( ) ( ) ( ) ( )ffffff DHRCRDF −=−= . (10) The objective of multi-exciter vibration control is to search for an appropriate driving spectrum with minimal error spectrum, i.e. ( )( ) 0DF =f . Broyden’s method is a highly powerful t t y t z x x y z (a) (b) Figure 3. Spatial orbit generation. (a) Time-domain waveform; (b) Elliptical orbit. x y z 1 2 x y z x y z u v w (a) (b) Figure 4. Coordinate transformation process. (a) The standard ellipse; (b) A spatial ellipse. ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 36 alternative to solving the equations by maintaining approximations of the solution and the Jacobian matrix as the iteration progresses [22]. If ( )fnD and ( )fnA are the current approximate solution and Jacobian matrix respectively, then the next approximate solution is given by ( ) ( ) ( ) ( ) ( )( )fffff nnnn CRADD −−= − + 1 1 . (11) After the computation of ( )fn 1+D , ( )fnA is updated to form ( )fn 1+A , and the construction of ( )fn 1+A is determined by Broyden’s method as ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( )ff ffff ff nn nnnn nn 1 T 1 T 111 1 ++ +++ + − += ss ssAy AA , (12) where ( ) ( ) ( ) ( ) ( )( ) ( )( )fff fff nnn nnn DFDFy DDs −= −= ++ ++ 11 11 (13) In order to avoid the inverse operation in equation (11), the Sherman-Morrison formula is applied in order to update the inverse matrix of the approximate Jacobian matrix, which is also called the system impedance. ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( )fff fffff ff nnn nnnnn nn 1 T 1 T 111 1 ++ +++ + − += yZs ZsyZs ZZ (14) Equation (11) is then rewritten as ( ) ( ) ( ) ( ) ( )( )ftttt nnnn CRZDD −−=+1 . (15) In addition, an adjustment gain  is supplied in the iterative expression in the case of receiving an overshoot response due to the large driving spectrum. ( ) ( ) ( ) ( ) ( )( )ftfff nnnn CRZDD −−=+ 1 (16) The norm of the error spectrum is used to judge how close the response is to the reference. ( )( ) ( ) ( )( ) ( ) ( )( )ffffff nnn 111 +++ −−= CRCRD T (17) Furthermore, the optimal adjustment gain can be obtained when the derivative of ( )( )tf n 1+D with reference to  is zero. ( ) ( ) ( ) ( )( )( ) ( ) ( )( ) ( ) ( ) ( ) ( )( )( ) ( ) ( ) ( ) ( )( )( )ffffffff ffffff nnnn nnn CRZHCRZH CRCRZH T T −− −−− = (18) In equation (18), the system’s frequency response function is unknown and should be estimated in terms of known inputs and responses. If a known value  is used for the adjustment gain adopted in equation (16), the driving spectrum is given by ( ) ( ) ( ) ( ) ( )( )fffff nnnn CRZDD −−=+ 1ˆ . (19) The corresponding response spectrum is measured as ( )fn 1ˆ +C . By left multiplying equation (19) by the system’s frequency response function ( )fH , we can get ( ) ( ) ( ) ( ) ( ) ( )( )ffffff nnnn CRZHCC −−=+ 1ˆ . (20) Then, we rewrite equation (20) into the following form ( ) ( ) ( ) ( )( ) ( ) ( )( )ffffff nnnn CCCRZH −−=− +1ˆ 1  . (21) The unknown expression in equation (18) can be eliminated by substituting equation (21) into equation (18) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( )ffff ffff nnnn nnn CCCC CRCC T T −− −− = ++ + 11 1 ˆˆ ˆ  . (22) The vibration control flow is illustrated by means of a block diagram, shown in Figure 5. The frequency response function matrix of the multi-exciter vibration system is first experimentally identified. Then, the initial system impedance is calculated by inverting the result. The reference spectrum and response spectrum are determined according to the desired waveforms and the measured response respectively. The vibration control consists of three main procedures: calculating the adjustment gain according to equation (22), updating the impedance matrix according to equation (14), and correcting the driving spectrum by means of the iterative algorithm expressed in equation (16). The driving signal is then generated according to the driving spectrum. Reference Spectrum Matrix Driven signal Response signal FRF Identification Initial Impedance Matrix Response Spectrum Matrix Iterative Algorithm Eq. (16) Driven Spectrum Matrix Desired waveformAdjustment gain Eq. (22) Impedance Matrix Eq. (14) Figure 5. Vibration control flow. 4. EXPERIMENTAL INVESTIGATION The proposed control method was applied on the tri-axial vibration exciter. Three accelerometers were fixed on the vibration table in order to observe the generated spatial orbit. The sensing axes of the feedback accelerometers were aligned with the three vibration exciters. The cross-coupling of the tri-axial vibration exciter was evaluated by sequentially applying single-axis vibration along its three axes (X, Y, and Z). At first, no vibration control was performed on the tri-axial vibration exciter, during which time only the exciter along the vibrating direction was active. The ratios of the cross-coupling motion to the principal axis motion are shown in Figure 6. As shown in Figure 6(a), the Y-axial cross- coupling ratio was about 0.1 and was close to 0.2 at certain frequencies (such as 100 Hz and 1600 Hz). The X-axial cross- coupling ratio was close to 0.1 as shown in Figure 6(b). The X- axial cross-coupling ratio in Figure 6(c) exceeded 0.1 and reached 0.2 from frequencies above 100 Hz and the Y-axial cross- coupling ratio was close to 0.1. The cross-coupling motion of the tri-axial vibration exciter was so large that it was unsuitable for calibration purposes. ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 37 Then, the proposed multi-exciter vibration control was performed on the tri-axial vibration exciter. All the three exciters were operational, among which the principal axis exciter provided the stimulus, and the other two were used for cross- coupling compensation. The ratio for the cross-coupling motion to the principal axis motion is shown in Figure 7. As seen in Figure 7, the ratio of the cross-coupling motion to the principal axis motion is less than 0.01, which is very small. Comparing Figure 7 with Figure 6, we can conclude that the proposed multi- exciter vibration control can efficiently reduce the cross-coupling of the tri-axial vibration exciter. In addition, three spatial linear orbits parallel to the coordinate axes are generated, as shown in Figure 8(a). Furthermore, four spatial circular orbits that rotate every 30 ° around the X-axis are generated, as shown in Figure 8(b). The plots in Figure 8 demonstrate that the proposed multi-exciter vibration control can generate spatial orbits for further calibration application. A tri-axis accelerometer (model: Kistler 8766A50, SN: 2118860) was selected for calibration using the tri-axial primary vibration calibration system. The calibration was carried out from 5 Hz to 1600 Hz. The magnitude and phase lag of the sensitivities measured are shown in Figure 9. The X-, Y-, and Z- axial magnitudes are 9.72 mV/(m/s2), 9.80 mV/(m/s2), and 10.13 mV/(m/s2) at the reference frequency 160 Hz. The X-, Y- and Z-axial phase shifts are 0.06 °, -0.29 °, and -0.62 °, respectively. The X- and Z-axial magnitudes of sensitivities vary slightly from 5 Hz to 1600 Hz. The Y-axial magnitude of sensitivity, however, has several peaks and troughs. The X-, Y- and Z-axial phase shifts of sensitivities ranges from -1.67 ° to -1.40 °, from 5 Hz to 1600 Hz. (a) (b) (c) Figure 6. Cross-coupling without vibration control. (a) X-axial main motion; (b) Y-axial main motion; (c) Z-axial main motion. (a) (b) (c) Figure 7. Cross-coupling with vibration control. (a) X-axial main motion; (b) Y- axial main motion; (c) Z-axial main motion. ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 38 5. CONCLUSIONS The National Institute of Metrology in China has set up a tri- axial primary vibration calibration system that consists of an air bearing-based tri-axial vibration exciter, a spatial orbit vibration control system, and a primary tri-axial vibration measuring system. Multi-exciter vibration control for cross-coupling compensation and the amplitude and phase control of sine vibration components is proposed. The experimental results show that the proposed control method can efficiently reduce the cross-coupling of the tri-axial vibration exciter and can generate a variety of spatial orbits for further calibration purposes. ACKNOWLEDGEMENT The Authors would like to thank Thomas Bruns and Leonard Klaus for their valuable guidance and help in the multicomponent vibration calibration. This work has been supported by the Special Fund for Scientific Research in the Public Interest (201410009) and National Natural Science Foundation of China (No. 51605461). REFERENCES [1] C. Harman, M. B. Pickel, Multi-axis vibration reduces test time, Evaluation Engineering 45, 6 (2006), pp. 44-47. [2] W. E. Whiteman, M. S. Berman, Fatigue failure results for multi- axial versus uniaxial stress screen vibration testing, Shock and Vibration 9 (2002), pp. 319-328. [3] ISO 16063-11, Methods for the calibration of vibration and shock pick-ups – Part 11: Primary vibration calibration by laser interferometry. [4] J. J. Dosch, D. M. Lally, Automated testing of accelerometer transverse sensitivity, http://www.modalshop.com/techlibrary/JDosch%20transverse %20calibration.pdf. [5] R. D. Sill, E. J. Seller, Accelerometer transverse sensitivity measurement using planar orbital motion, Proc. of the 77th Shock and Vibration Symposium, Monterey, CA, USA, 2006, pp. 8-12. [6] H.-J. von Martens, C. Weissenborn, Simultaneous multi- component calibration - A new research area in the field of vibration and shock, Proc. of the 1st Meeting of the Consultative Committee for Acoustics, Ultrasound and Vibration. Bureau International de Poids et Mesures, Sèvres, France, 20 – 21 July 1999. [7] C. A. Weissenborn, A new measurement system for simultaneous multicomponent calibration of motion transducers, Proc. of the 9th International Conference for Sensors, Transducers and Systems (Sensor 99), Nürnberg, Germany, 18 – 20 May 1999. [8] C. Hof, M. Kobusch, Comparison of the Calibration of a Heavy Multi Component Vibration Transducer on Different Exciter Systems (Calibration of Heavy Triax Transducer), Proc. of the 20th IMEKO TC3 Conference “Cultivating Metrological Knowledge”, 27 - 30 November 2007, Merida, Mexico. [9] ISO 16063-31, Methods for the calibration of vibration and shock transducers - Part 31: Testing of transverse vibration sensitivity, 2009. (a) (b) Figure 8. Spatial orbits. (a) Linear orbit; (b) Circular orbit. (a) (b) Figure 9. Calibrated sensitivities. (a) Magnitude; (b) Phase shift. http://www.modalshop.com/techlibrary/JDosch%20transverse%20calibration.pdf http://www.modalshop.com/techlibrary/JDosch%20transverse%20calibration.pdf ACTA IMEKO | www.imeko.org March 2019 | Volume 8 | Number 1 | 39 [10] T. Usuda, H.-J. von Martens, C. Weissenborn, Theoretical and experimental investigation of transverse sensitivity of accelerometers under multiaxial excitation, Meas. Sci. Technol. 15 (2004), pp. 896–904. [11] A. Umeda, M. Onoe, K. Sakata, T. Fukushia, K. Kanari, H. Iioka, T. Kobayashi, Calibration of three-axis accelerometers using a three-dimensional vibration generator and three laser interferometers, Sensors and Actuators A 114 (2004) pp. 93–101. [12] IEC 60747-14-4, Semiconductor devices- Discrete devices- Part 14-4: Semiconductor accelerometers, 2011. [13] T. Tsuchiya, O. Tabata, A. Umeda, Dynamic sensitivity matrix measurement for single-mass SOI 3-axis accelerometer, Proc. of the IEEE 25th International Conference on Micro Electro Mechanical Systems (MEMS), 29 January - 2 February 2012, Paris, France. [14] C. Harman, Historical development of high performance multi- axis vibration test systems, Environmental Engineering 17, 1 (2004), pp. 40-41. [15] R. C. Stroud, G. A. Hamma, M. A. Underwood, et al., A review of multiaxis/multiexciter vibration technology, Sound and Vibration 30, 4 (1996), pp. 20-27. [16] G. A. Hamma, S. Smith, Simulation of dynamic loads by multichannel digital control. AIAA 78, 511 (1978), pp. 3-6. [17] M. Chen, D. Wilson. The new tri-axial shock and vibration test system at hill air force base, Journal of IEST 41, 2 (1989), pp. 27- 32. [18] M. A. Underwood, Adaptive control method for multiexciter sine tests. US Patent 5299459, 1994. [19] M. A. Underwood, T. Keller, Recent system developments for multi-actuator vibration control, Sound and Vibration 35, 10 (2001), pp. 16-23. [20] H. Aoki, S. Tsutsumi, T. Fukushima, et al., Vibration testing apparatus with increased rigidity in static pressure bearing. US Patent 5549005[P], 1996. [21] M. P. Bulat, P. V. Bulat, The history of the gas bearings theory development, World Applied Sciences Journal 27, 7 (2013), pp. 893-897. [22] C. T. Kelley, Iterative methods for linear and nonlinear equations. North Carolina State University, Raleigh, NC, USA, 1995.