- Email: [email protected]
THE EFFECTS OF BILINEAR MAPPING OF NON-LINEAR SYSTEMS
Mechanical Systems and Signal Processing (1998) 12(5), 713–719 Article No. pg980164
THE EFFECTS OF BILINEAR MAPPING OF NON-LINEAR SYSTEMS K. S J. K. H ISVR, University of Southampton, Southampton SO17 1BJ, U.K.
(Received February 1998 , accepted after revisions June 1998) The bilinear mapping method is very important in designing digital filters from analogue filters for linear systems, because it provides a very simple conversion from a continuous system to a discrete system without introducing aliasing problems. The use of the bilinear mapping method for non-linear systems is investigated by applying it to a simple non-linear ordinary differential system. The results show that the bilinear mapping method is very stable compared to the fourth-order Runge–Kutta method and could be used for non-linear digital filter design.
7 1998 Academic Press
1. INTRODUCTION
When a continuous system is mapped to an equivalent discrete system, one may use various numerical methods such as the Euler method, the Runge–Kutta method, and the bilinear method. Recently, ‘spurious’ solutions introduced by the time discretisation of ordinary differential equations have been investigated. Solutions of the equivalent difference equations may be completely different from that of the original non-linear differential equations depending on the numerical approximation (e.g. step size). Iserles and others [1–5] have studied such a problem particularly for an ordinary differential system given by dy = f(y), dt
y(t0 ) = y0 ,
y(t)$Rd
(1)
where Rd represents d-dimensional space and y0 is the initial condition, and the corresponding discrete map which approximates (1) is given by yn + 1 = FT (yn ),
y0 $Rd
(2)
where the T is the fixed step size. For the ‘‘Good Method Guide’’, Iserles suggested a multistep method in the form [1] k
k
i=0
i=0
s ai yn + i = T s bi f(yn + i )
(3)
where ai , bi are coefficients. When bk is not zero, the method is implicit. One can intuitively anticipate that the multistep methods are likely to model dynamics more faithfully since increasing the order in equation (3) may result in increasing the accuracy, i.e. an m-step method uses the m previous points. Typically an implicit method has a larger region of stability than the corresponding explicit method, and this is mainly due to the feedback nature in an implicit method [6]. Iserles also suggested the implicit midpoint rule as the 0888–3270/98/050713 + 07 $30.00/0
7 1998 Academic Press
. . .
714
best method [1]. The value yn + 1 is obtained by finding the midpoint between yn + 1 and the previous value yn , and implicitly depends on yn + 1 itself as defined in equation (4). yn + 1 = yn + Tf(12 (yn + yn + 1 )).
(4)
Iserles also mentioned another example of a multistep method which may be used without numerical artifacts: the trapezoidal integration rule which yields yn + 1 = yn +
T [f(yn ) + f(yn + 1 )]. 2
(5)
This is notable since it is recognised as the bilinear mapping method in signal processing context. In signal processing, the bilinear mapping method is widely used for mapping continuous linear systems to discrete systems, because it is a simple conversion which avoids aliasing problems. Thus, the bilinear mapping method is very important in designing digital filters from analogue filters for linear systems. So far it has not been reported the effects of bilinear mapping method for non-linear systems. Thus, in this paper, the bilinear mapping method is extended to a non-linear system by applying the trapezoidal integration rule to a non-linear differential equation. By varing the step size, the effects of bilinear mapping method are described, and also compared to the Runge–Kutta method. The results show that the bilinear mapping method (trapezoidal integration rule) is very stable compared to the Runge–Kutta method which gives erroneous solutions when the step size is large. It is also found that the fundamental property of the bilinear mapping method for linear systems carried over to non-linear systems with aliasing which will be shown later. Details of numerical simulation results of both linear and non-linear systems can be found in [7, 8] including many other numerical aspects. 2. SIMULATIONS
Consider a non-linear differential system without a forcing term x¨ + v02x + ox 3 = 0.
(6)
Equation (6) can be expressed in the form (1) as x˙ = y y˙ = −v02x − ox 3.
140
10
(a)
10120
0.5
Amplitude
Amplitude
1
0 –0.5 –1
(7)
(b)
10100 1080 1060 1040 1020
0
5
10
15 20 Time (s)
25
30
100 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time (s)
Figure 1. Solutions of the system (6) using the Runge–Kutta method (f0 = 2). Time histories for (a) T = 0.1, (b) T = 0.25.
3
10
0 –0.5 –1
715
(b)
102
0.5 Modulus
Amplitude
1
(a)
101 100 10–1
0
1
2 3 Time (s)
4
10–2 0
5
5
10 15 Frequency (Hz)
20
Figure 2. Solutions of the system (6) using the Runge–Kutta method (fR1 1 2.4, fR3 1 7.2 and fR5 1 12). (a) Time history (T = 0.01), (b) frequency domain.
Equation (7) can be mapped by using the well-known fourth-order Runge–Kutta method and the bilinear mapping method [equation (5)]. The details of method solving equation (5) can be found in [7, 8]. The results obtained by these two methods are shown by setting the initial conditions x0 = 1, y0 = 0, o = 100 and f0 = 2, where v0 = 2pf0 . The Runge–Kutta method is very sensitive to the step size T. Thus when the step size is relatively large compared to the inherent system’s oscillation, this method gives an erroneous result (Fig. 1). Note that the original system does not have damping [Fig. 1(a)], and the solution soon diverges after a few iterations [Fig. 1(b)]. Thus, the step size is set to be sufficiently small for the Runge–Kutta method (Fig. 2) to get close to the original differential system [i.e. it is assumed that the result of the Runge–Kutta method with step size T = 0.01 represents the true system (6) for this simulation]. The frequency (Hz) components from the Runge–Kutta method are fR1 1 2.4, fR3 1 7.2 and fR5 1 12, where fR1 is the fundamental frequency and the others are third and fifth harmonics (only odd harmonics are present in this case, i.e., fR3 = 3 × fR1 and fR5 = 5 × fR1 ). The bilinear mapping method is investigated by varying the step size to see the effect of discretisation. With step size T = 0.01 solutions from the Runge–Kutta method and the bilinear mapping method are almost the same as shown in Fig. 3, i.e. the frequency components by the bilinear method are fd1 1 2.4, fd3 1 7.2 and fd5 1 12, where fd1 is the fundamental frequency and the others are third and fifth harmonics. In the case of T = 0.05, the results of the bilinear mapping method are shown in Fig. 4, and the frequency components are fd1 1 2.3, fd3 1 6.9 and fd5 1 8.5 (aliased version of fd5 1 11.5). From the above two cases (T = 0.01 and T = 0.05), one can see that the (a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1
0
1
2 3 Time (s)
4
5
10–2 0
5
10 15 Frequency (Hz)
20
Figure 3. Solutions of the system (6) using the bilinear mapping method (fd1 1 2.4, fd3 1 7.2, fd5 1 12). (a) Time history (T = 0.01), (b) frequency domain.
. . .
716 (a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1 10–2
0
1
2 3 Time (s)
4
5
0
2
4 6 Frequency (Hz)
8
10
Figure 4. Solutions of the system using the bilinear mapping method (fd1 1 2.3, fd3 1 6.9, fd5 1 8.5). (a) Time history (T = 0.05), (b) frequency domain.
fundamental frequency (fR1 1 2.4) is shifted to fd1 by the following equation which is the relationship between digital and analogue frequencies for a linear system. vd =
0 1
va T 2 tan−1 2 T
(8)
where vd is the digital angular frequency by the bilinear mapping method and va is the analogue angular frequency (original differential system). However, unlike the linear system, it is seen that aliasing does occur in the bilinear mapping method for non-linear systems. This is caused by the fact that the harmonic components do not follow equation (8). The harmonic components in the digital frequency domain (obtained using the bilinear mapping method) remain as the harmonics of the fundamental frequency in the digital frequency domain [which is computed from equation (8)], i.e. fd3 = 3 × fd1 and fd5 = 5 × fd1 . Thus, the fifth harmonic fd5 is 11.5 Hz, which is greater than the folding frequency 10 Hz and so appears at 8.5 Hz as an alias. This aliasing effect is more readily shown when the step size is increased further, and the results of the cases T = 0.1 and T = 0.25 are shown in Figs 5 and 6, and summarised in Table 1 by comparing original frequencies (fRi ) to the results from the bilinear mapping method (fdi ). The following non-rigorous argument helps to explain the frequency components that arise. Equation (6) is often solved approximately using perturbation methods when o is small [9], i.e. the solution of equation (6) can be expressed as a power series
(a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1
0
2
4 6 Time (s)
8
10
10–2 0
1
2 3 Frequency (Hz)
4
5
Figure 5. Solutions of the system using the bilinear mapping method (fd1 1 2.06, fd3 1 3.82, fd5 1 0.3). Compare the time series with Figure 1(a). (a) Time history (T = 0.1), (b) frequency domain.
3
10
0 –0.5 –1
717
(b)
102
0.5 Modulus
Amplitude
1
(a)
0
5
10 Time (s)
15
20
101 100 10–1 0
0.5
1 1.5 Frequency (Hz)
2
Figure 6. Solutions of the system using the bilinear mapping method (fd1 1 1.38, fd3 1 0.14, fd5 1 1.1). Compare the time series with Figure 1(b). (a) Time history (T = 0.25), (b) frequency domain.
x(t) = x0 (t) + ox1 (t) + o 2x2 (t) + · · ·
(9)
and one perturbs the frequency v by letting v02 = v 2 − ob1 − o 2b2 − · · ·
(10)
where bi are the corrective terms for the frequency and are dependent on the amplitude of the motion. This yields the hierarchy of equations x¨0 + v 2x0 = 0
(11a)
x¨1 + v x1 = b1 x0 − x 2
3 0
x¨2 + v 2x2 = b1 x1 + b2 x0 − 3x1 x02 .
(11b) (11c)
This shows that the fundamental frequency component shifts from v0 and the evaluations of the bi (reference [9]) allow the modified fundamental frequency v to be computed. Further the solution to equation (11a) acts as an input to equation (11b); the non-linear cubic term indicates the harmonic at 3v. Assuming that applying the trapezoidal integration rule to equation (5) is equivalent to applying it to the set of equations (11), then equation (11a) will result in a digital oscillation occurring at vd = 2/T tan−1 (vT/2) [i.e., as equation (8)]. Equation (11b) is interpreted as a linear system driven by an input b1 x0 − x03 which results in a digital frequency of vd plus a harmonic at 3vd . The right-hand sides of equations (11a–c) are like forcing terms which are subject to aliasing leading to the results obtained in Figs 4–6. In other words, the input–output relationship can be described by the Laplace transformation which becomes X(s) = H(s)F(s) for zero initial conditions, where H(s) is the transfer function of the system. The continuous system
T 1 Summary of the frequency components of the system Original frequency
Bilinear mapping method (T = 0.01)
Bilinear mapping method (T = 0.05)
Bilinear mapping method (T = 0.1)
Bilinear mapping method (T = 0.25)
fR1 1 2.4 fR3 1 7.2 fR5 1 12
fd1 1 2.4 fd3 1 7.2 fd5 1 12
2.3 6.9 8.5* (11.5)
2.06 3.82* (6.18) 0.3* (10.3)
1.38 0.14* (4.14) 1.1* (6.9)
* The aliased version of ( ).
718
. . .
function H(s) can be transformed to the discrete-time system function H(z) by using the bilinear transform, i.e. H(z) = H(s),
where s =
0
1
2 1 − z−1 . T 1 + z−1
(12)
The difference equation corresponding to H(z) can be obtained by applying the trapezoidal integration rule (5) to the differential equation corresponding to H(s), and if the input–output relationship of the system is described pictorially, then f(nT) 0004
b $ 0 H
1% b
2 1 − z−1 H 1 + z−1
x(nT) 0004
Notice that the system is bilinearly mapped, but the input is not. From the above illustration, it is obvious that if the input is aliased by discretising coarsely, then the output is also aliased. If the harmonic forcing term is included in equation (6), it can be written as x¨ + v02 x = −ox 3 + f(t),
where f(t) = A cos vt
(13)
Thus, it can be interpreted that the whole right-hand side of equation (13) is an input whose effect is much the same as that described above (aliasing problem). 3. DISCUSSION
The results of the previous section show that the bilinear mapping method for non-linear systems may be used in practical situations. Since the fundamental frequency of the non-linear system obeys equation (8), one can understand the original differential system from the results of the bilinear mapping method. Accordingly it may be possible to apply this to non-linear digital filter design with care, i.e. aliasing can occur in non-linear systems, and also the results of the bilinear method will give a good guide for the determination of the step sizes for other numerical methods, because the bilinear method is very stable even when the step size is quite large. For unknown systems one can apply the bilinear method first in order to obtain fd1 (fundamental frequency) and then convert this value to the analogue frequency by using equation (8). Then one can easily determine the step sizes for other numerical methods since the inherent oscillation of the system may be understood from the results. One must be aware of the possibility that application of the bilinear method to a non-linear system does not ensure the accuracy of the true trajectories of the original system, i.e. the original frequency of the continuous system always becomes lower when the bilinear method is applied. REFERENCES 1. A. I 1994 Numerical Analysis Report, University of Cambridge, DAMTP 1994 /NA5 . Dynamics of Numerics. 2. A. T. P 1989 M.Sc. Thesis, University of Bath. Dynamics of Numerical Methods of Initial Value Problems. 3. J. M. S-S 1990 Applied Mathematics and Computation Report, Report 1990 /3 , Universidad de Valladolid. Numerical ordinary differential equations vs. dynamical systems. 4. A. I, A. T. P and A. M. S 1991 Siam Journal of Numerical Analysis 28, 1723–1751. A unified approach to spurious solutions introduced by time discretisation. Part I: basic theory.
719
5. D. F. G, P. K. S and H. C. Y 1992 IMA Journal of Numerical Analysis 12, 319–338. On spurious asymptotic numerical solutions of explicit Runge–Kutta. 6. L. L and J. H. S 1971 Numerical Solution of Ordinary Differential Equations. London: Academic Press. 7. K. S 1996 PhD Thesis, Institute of Sound and Vibration Research, University of Southampton. Characterisation and identification of chaotic dynamical systems. 8. K. S 1997 Technical Report of the Institute of Sound and Vibration Research, University of Southampton. The numerical effect of discretisation of continuous systems and bilinear mapping of non-linear systems. 9. F. S. T, I. E. M and R. T. H 1988 Mechanical Vibrations; Theory and Applications. Prentice Hall. 10. A. V. O and R. W. S 1989 Discrete-Time Signal Processing. Prentice Hall. 11. S. Z. A-H 1994 MSc Thesis, University of Southampton. The numerical effect of discretisation of continuous time nonlinear systems. 12. D. S 1995 Journal of Sound and Vibration 180, 513–518. A direct integration algorithm and the consequence of numerical stability.
THE EFFECTS OF BILINEAR MAPPING OF NON-LINEAR SYSTEMS K. S J. K. H ISVR, University of Southampton, Southampton SO17 1BJ, U.K.
(Received February 1998 , accepted after revisions June 1998) The bilinear mapping method is very important in designing digital filters from analogue filters for linear systems, because it provides a very simple conversion from a continuous system to a discrete system without introducing aliasing problems. The use of the bilinear mapping method for non-linear systems is investigated by applying it to a simple non-linear ordinary differential system. The results show that the bilinear mapping method is very stable compared to the fourth-order Runge–Kutta method and could be used for non-linear digital filter design.
7 1998 Academic Press
1. INTRODUCTION
When a continuous system is mapped to an equivalent discrete system, one may use various numerical methods such as the Euler method, the Runge–Kutta method, and the bilinear method. Recently, ‘spurious’ solutions introduced by the time discretisation of ordinary differential equations have been investigated. Solutions of the equivalent difference equations may be completely different from that of the original non-linear differential equations depending on the numerical approximation (e.g. step size). Iserles and others [1–5] have studied such a problem particularly for an ordinary differential system given by dy = f(y), dt
y(t0 ) = y0 ,
y(t)$Rd
(1)
where Rd represents d-dimensional space and y0 is the initial condition, and the corresponding discrete map which approximates (1) is given by yn + 1 = FT (yn ),
y0 $Rd
(2)
where the T is the fixed step size. For the ‘‘Good Method Guide’’, Iserles suggested a multistep method in the form [1] k
k
i=0
i=0
s ai yn + i = T s bi f(yn + i )
(3)
where ai , bi are coefficients. When bk is not zero, the method is implicit. One can intuitively anticipate that the multistep methods are likely to model dynamics more faithfully since increasing the order in equation (3) may result in increasing the accuracy, i.e. an m-step method uses the m previous points. Typically an implicit method has a larger region of stability than the corresponding explicit method, and this is mainly due to the feedback nature in an implicit method [6]. Iserles also suggested the implicit midpoint rule as the 0888–3270/98/050713 + 07 $30.00/0
7 1998 Academic Press
. . .
714
best method [1]. The value yn + 1 is obtained by finding the midpoint between yn + 1 and the previous value yn , and implicitly depends on yn + 1 itself as defined in equation (4). yn + 1 = yn + Tf(12 (yn + yn + 1 )).
(4)
Iserles also mentioned another example of a multistep method which may be used without numerical artifacts: the trapezoidal integration rule which yields yn + 1 = yn +
T [f(yn ) + f(yn + 1 )]. 2
(5)
This is notable since it is recognised as the bilinear mapping method in signal processing context. In signal processing, the bilinear mapping method is widely used for mapping continuous linear systems to discrete systems, because it is a simple conversion which avoids aliasing problems. Thus, the bilinear mapping method is very important in designing digital filters from analogue filters for linear systems. So far it has not been reported the effects of bilinear mapping method for non-linear systems. Thus, in this paper, the bilinear mapping method is extended to a non-linear system by applying the trapezoidal integration rule to a non-linear differential equation. By varing the step size, the effects of bilinear mapping method are described, and also compared to the Runge–Kutta method. The results show that the bilinear mapping method (trapezoidal integration rule) is very stable compared to the Runge–Kutta method which gives erroneous solutions when the step size is large. It is also found that the fundamental property of the bilinear mapping method for linear systems carried over to non-linear systems with aliasing which will be shown later. Details of numerical simulation results of both linear and non-linear systems can be found in [7, 8] including many other numerical aspects. 2. SIMULATIONS
Consider a non-linear differential system without a forcing term x¨ + v02x + ox 3 = 0.
(6)
Equation (6) can be expressed in the form (1) as x˙ = y y˙ = −v02x − ox 3.
140
10
(a)
10120
0.5
Amplitude
Amplitude
1
0 –0.5 –1
(7)
(b)
10100 1080 1060 1040 1020
0
5
10
15 20 Time (s)
25
30
100 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Time (s)
Figure 1. Solutions of the system (6) using the Runge–Kutta method (f0 = 2). Time histories for (a) T = 0.1, (b) T = 0.25.
3
10
0 –0.5 –1
715
(b)
102
0.5 Modulus
Amplitude
1
(a)
101 100 10–1
0
1
2 3 Time (s)
4
10–2 0
5
5
10 15 Frequency (Hz)
20
Figure 2. Solutions of the system (6) using the Runge–Kutta method (fR1 1 2.4, fR3 1 7.2 and fR5 1 12). (a) Time history (T = 0.01), (b) frequency domain.
Equation (7) can be mapped by using the well-known fourth-order Runge–Kutta method and the bilinear mapping method [equation (5)]. The details of method solving equation (5) can be found in [7, 8]. The results obtained by these two methods are shown by setting the initial conditions x0 = 1, y0 = 0, o = 100 and f0 = 2, where v0 = 2pf0 . The Runge–Kutta method is very sensitive to the step size T. Thus when the step size is relatively large compared to the inherent system’s oscillation, this method gives an erroneous result (Fig. 1). Note that the original system does not have damping [Fig. 1(a)], and the solution soon diverges after a few iterations [Fig. 1(b)]. Thus, the step size is set to be sufficiently small for the Runge–Kutta method (Fig. 2) to get close to the original differential system [i.e. it is assumed that the result of the Runge–Kutta method with step size T = 0.01 represents the true system (6) for this simulation]. The frequency (Hz) components from the Runge–Kutta method are fR1 1 2.4, fR3 1 7.2 and fR5 1 12, where fR1 is the fundamental frequency and the others are third and fifth harmonics (only odd harmonics are present in this case, i.e., fR3 = 3 × fR1 and fR5 = 5 × fR1 ). The bilinear mapping method is investigated by varying the step size to see the effect of discretisation. With step size T = 0.01 solutions from the Runge–Kutta method and the bilinear mapping method are almost the same as shown in Fig. 3, i.e. the frequency components by the bilinear method are fd1 1 2.4, fd3 1 7.2 and fd5 1 12, where fd1 is the fundamental frequency and the others are third and fifth harmonics. In the case of T = 0.05, the results of the bilinear mapping method are shown in Fig. 4, and the frequency components are fd1 1 2.3, fd3 1 6.9 and fd5 1 8.5 (aliased version of fd5 1 11.5). From the above two cases (T = 0.01 and T = 0.05), one can see that the (a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1
0
1
2 3 Time (s)
4
5
10–2 0
5
10 15 Frequency (Hz)
20
Figure 3. Solutions of the system (6) using the bilinear mapping method (fd1 1 2.4, fd3 1 7.2, fd5 1 12). (a) Time history (T = 0.01), (b) frequency domain.
. . .
716 (a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1 10–2
0
1
2 3 Time (s)
4
5
0
2
4 6 Frequency (Hz)
8
10
Figure 4. Solutions of the system using the bilinear mapping method (fd1 1 2.3, fd3 1 6.9, fd5 1 8.5). (a) Time history (T = 0.05), (b) frequency domain.
fundamental frequency (fR1 1 2.4) is shifted to fd1 by the following equation which is the relationship between digital and analogue frequencies for a linear system. vd =
0 1
va T 2 tan−1 2 T
(8)
where vd is the digital angular frequency by the bilinear mapping method and va is the analogue angular frequency (original differential system). However, unlike the linear system, it is seen that aliasing does occur in the bilinear mapping method for non-linear systems. This is caused by the fact that the harmonic components do not follow equation (8). The harmonic components in the digital frequency domain (obtained using the bilinear mapping method) remain as the harmonics of the fundamental frequency in the digital frequency domain [which is computed from equation (8)], i.e. fd3 = 3 × fd1 and fd5 = 5 × fd1 . Thus, the fifth harmonic fd5 is 11.5 Hz, which is greater than the folding frequency 10 Hz and so appears at 8.5 Hz as an alias. This aliasing effect is more readily shown when the step size is increased further, and the results of the cases T = 0.1 and T = 0.25 are shown in Figs 5 and 6, and summarised in Table 1 by comparing original frequencies (fRi ) to the results from the bilinear mapping method (fdi ). The following non-rigorous argument helps to explain the frequency components that arise. Equation (6) is often solved approximately using perturbation methods when o is small [9], i.e. the solution of equation (6) can be expressed as a power series
(a)
3
10
0 –0.5 –1
(b)
102
0.5 Modulus
Amplitude
1
101 100 10–1
0
2
4 6 Time (s)
8
10
10–2 0
1
2 3 Frequency (Hz)
4
5
Figure 5. Solutions of the system using the bilinear mapping method (fd1 1 2.06, fd3 1 3.82, fd5 1 0.3). Compare the time series with Figure 1(a). (a) Time history (T = 0.1), (b) frequency domain.
3
10
0 –0.5 –1
717
(b)
102
0.5 Modulus
Amplitude
1
(a)
0
5
10 Time (s)
15
20
101 100 10–1 0
0.5
1 1.5 Frequency (Hz)
2
Figure 6. Solutions of the system using the bilinear mapping method (fd1 1 1.38, fd3 1 0.14, fd5 1 1.1). Compare the time series with Figure 1(b). (a) Time history (T = 0.25), (b) frequency domain.
x(t) = x0 (t) + ox1 (t) + o 2x2 (t) + · · ·
(9)
and one perturbs the frequency v by letting v02 = v 2 − ob1 − o 2b2 − · · ·
(10)
where bi are the corrective terms for the frequency and are dependent on the amplitude of the motion. This yields the hierarchy of equations x¨0 + v 2x0 = 0
(11a)
x¨1 + v x1 = b1 x0 − x 2
3 0
x¨2 + v 2x2 = b1 x1 + b2 x0 − 3x1 x02 .
(11b) (11c)
This shows that the fundamental frequency component shifts from v0 and the evaluations of the bi (reference [9]) allow the modified fundamental frequency v to be computed. Further the solution to equation (11a) acts as an input to equation (11b); the non-linear cubic term indicates the harmonic at 3v. Assuming that applying the trapezoidal integration rule to equation (5) is equivalent to applying it to the set of equations (11), then equation (11a) will result in a digital oscillation occurring at vd = 2/T tan−1 (vT/2) [i.e., as equation (8)]. Equation (11b) is interpreted as a linear system driven by an input b1 x0 − x03 which results in a digital frequency of vd plus a harmonic at 3vd . The right-hand sides of equations (11a–c) are like forcing terms which are subject to aliasing leading to the results obtained in Figs 4–6. In other words, the input–output relationship can be described by the Laplace transformation which becomes X(s) = H(s)F(s) for zero initial conditions, where H(s) is the transfer function of the system. The continuous system
T 1 Summary of the frequency components of the system Original frequency
Bilinear mapping method (T = 0.01)
Bilinear mapping method (T = 0.05)
Bilinear mapping method (T = 0.1)
Bilinear mapping method (T = 0.25)
fR1 1 2.4 fR3 1 7.2 fR5 1 12
fd1 1 2.4 fd3 1 7.2 fd5 1 12
2.3 6.9 8.5* (11.5)
2.06 3.82* (6.18) 0.3* (10.3)
1.38 0.14* (4.14) 1.1* (6.9)
* The aliased version of ( ).
718
. . .
function H(s) can be transformed to the discrete-time system function H(z) by using the bilinear transform, i.e. H(z) = H(s),
where s =
0
1
2 1 − z−1 . T 1 + z−1
(12)
The difference equation corresponding to H(z) can be obtained by applying the trapezoidal integration rule (5) to the differential equation corresponding to H(s), and if the input–output relationship of the system is described pictorially, then f(nT) 0004
b $ 0 H
1% b
2 1 − z−1 H 1 + z−1
x(nT) 0004
Notice that the system is bilinearly mapped, but the input is not. From the above illustration, it is obvious that if the input is aliased by discretising coarsely, then the output is also aliased. If the harmonic forcing term is included in equation (6), it can be written as x¨ + v02 x = −ox 3 + f(t),
where f(t) = A cos vt
(13)
Thus, it can be interpreted that the whole right-hand side of equation (13) is an input whose effect is much the same as that described above (aliasing problem). 3. DISCUSSION
The results of the previous section show that the bilinear mapping method for non-linear systems may be used in practical situations. Since the fundamental frequency of the non-linear system obeys equation (8), one can understand the original differential system from the results of the bilinear mapping method. Accordingly it may be possible to apply this to non-linear digital filter design with care, i.e. aliasing can occur in non-linear systems, and also the results of the bilinear method will give a good guide for the determination of the step sizes for other numerical methods, because the bilinear method is very stable even when the step size is quite large. For unknown systems one can apply the bilinear method first in order to obtain fd1 (fundamental frequency) and then convert this value to the analogue frequency by using equation (8). Then one can easily determine the step sizes for other numerical methods since the inherent oscillation of the system may be understood from the results. One must be aware of the possibility that application of the bilinear method to a non-linear system does not ensure the accuracy of the true trajectories of the original system, i.e. the original frequency of the continuous system always becomes lower when the bilinear method is applied. REFERENCES 1. A. I 1994 Numerical Analysis Report, University of Cambridge, DAMTP 1994 /NA5 . Dynamics of Numerics. 2. A. T. P 1989 M.Sc. Thesis, University of Bath. Dynamics of Numerical Methods of Initial Value Problems. 3. J. M. S-S 1990 Applied Mathematics and Computation Report, Report 1990 /3 , Universidad de Valladolid. Numerical ordinary differential equations vs. dynamical systems. 4. A. I, A. T. P and A. M. S 1991 Siam Journal of Numerical Analysis 28, 1723–1751. A unified approach to spurious solutions introduced by time discretisation. Part I: basic theory.
719
5. D. F. G, P. K. S and H. C. Y 1992 IMA Journal of Numerical Analysis 12, 319–338. On spurious asymptotic numerical solutions of explicit Runge–Kutta. 6. L. L and J. H. S 1971 Numerical Solution of Ordinary Differential Equations. London: Academic Press. 7. K. S 1996 PhD Thesis, Institute of Sound and Vibration Research, University of Southampton. Characterisation and identification of chaotic dynamical systems. 8. K. S 1997 Technical Report of the Institute of Sound and Vibration Research, University of Southampton. The numerical effect of discretisation of continuous systems and bilinear mapping of non-linear systems. 9. F. S. T, I. E. M and R. T. H 1988 Mechanical Vibrations; Theory and Applications. Prentice Hall. 10. A. V. O and R. W. S 1989 Discrete-Time Signal Processing. Prentice Hall. 11. S. Z. A-H 1994 MSc Thesis, University of Southampton. The numerical effect of discretisation of continuous time nonlinear systems. 12. D. S 1995 Journal of Sound and Vibration 180, 513–518. A direct integration algorithm and the consequence of numerical stability.