Journal of Sound and Vibration (1991) 144(2), 293-304

QUASI-PERIODIC SOLUTIONS CALCULATED WITH THE

SIMPLE SHOOTING TECHNIQUE
F. H. LlNGt
Department of Engineering Mechanics, Shanghai Jiao Tong University, Shanghai 200 030,
Peoples Republic of China
(Received 24 April 1989, and in revised,form 22 March 1990)
A shooting-type numerical method for calculating quasi-periodic solutions in a multi-
excited system by using derivatives and with an improved interpolation technique is
presented. Three examples of ordinary differential equations and iterative maps show the
efficiency of the new algorithm. This algorithm can also be effectively applied to a periodic
system.
1. INTRODUCTION
The calculation of quasi-periodic solutions is very important from a practical point of
view. A great number of regular motions in a multi-excited system should be regarded
as quasi-periodic rather than as periodic, since the ratio of two or more exciting frequencies
will usually not be a simple rational number, and if the denominator of this number
is not very small, then the period of the solution will be so large that it is difficult to
reveal the periodicity even within a pretty long observation time. In a real computation,
it is also more convenient to treat them as quasi-periodic solutions as stated below.
On the other hand, there is also a need for calculating quasi-periodic solutions of a
non-linear system in chaotic dynamics. Firstly, a route to chaos through quasi-periodic
motions has been intensively studied. Secondly, one often wants to know the distribution
of different motion patterns in the parameter space, but it is rather time-consuming to
identify a chaotic motion. However, since there exist only four kinds of spatially bounded
motion in non-linear systems-equilibrium, periodic, quasi-periodic and chaotic motion-
and the domain of equilibrium in the parameter space is usually of zero measure, one
can perform the task by determining the distribution of periodic and quasi-periodic
motions.
An efficient way for calculating quasi-periodic solutions is to extract a discrete time 7
map from the continuous flow and its derivative starting from an initial point, where T
is a reference period. This map is usually called a PoincarC map, and the iterative points
of the PoincarC map of quasi-periodic solutions are located on a closed curve in the phase
plane. For a single-excited system, i.e., a periodic system with period T, this curve is
invariant, but for a multi-excited system with non-commensurable frequencies, i.e., a
non-periodic system, although this curve is still closed it is no longer invariant. In other
words, there is a time-independent map for the periodic system, and a time-dependent
map for the non-periodic system.
There are quite a number of publications on the topic of the numerical computation
of the invariant curve of a periodic system. In earlier works [l-4] one tried to locate a
t Present address: Department of Physics and Engineering Physics, Stevens Institute of Technology. Castle
Point, Hoboken, New Jersey 07030, U.S.A.
293
0022-460X/91/020293+12$03.00/0 !$Q 1991 Academic Press Limited
294 F. H. LING
single point on the invariant curve and henceforth the corresponding quasi-periodic
solution. In some recent publications [5-71, one searches for many points on the invariant
curve simultaneously with a collocation method. The collocation method is efficient for
a periodic system, but it is not applicable to a non-periodic system: e.g., most multi-excited
systems.
The method suggested by Kaas-Petersen [&IO] can be regarded as an improvement
of the above-mentioned single point method, in which an interpolation technique is
incorporated. This method is emphasized in the present application to a multi-excited
system, and will be further improved in two aspects in this paper: namely, by using
derivatives instead of the difference approximation and by using all iterative points in
the interpolation. The corresponding new algorithm is described in the next section, and
three examples in section 3 illustrate that the new algorithm is more robust and time-saving.
Moreover, the method is also very good at finding invariant curves in the periodic system,
as shown in section 4.
2. IMPROVED SIMPLE SHOOTING METHOD
In this method, the condition that the Poincar6 map of quasi-periodic solutions of a
non-linear ordinary differential equation should be a closed curve is used for determining
a quasi-periodic solution. For example, a bi-periodic solution x(t) with two incommensur-
able periods T, and Tz can be written as x(r) = X(t, , fJ, where t, = t/ T, , tz =t/ T,, and
the quasi-periodicity requires that Z( 1,) f2) = a( f, + 1, f?) = %( t, , tz+ 1). By defining
To = (kT, mod TJ / T,,
(1)
where k is an integer, one obtains the PoincarC map of the quasi-periodic solution as
P( 7k) =x( kT,) =ri( k, kT,/ TJ =a(0, auk),
and therefore
P(1) =X(0,1) = Z(O, 0) =x(O) and P(0) = Z(O, 0) =x(O).
In this way, the calculation of quasi-periodic solutions is converted to a problem
determining a proper x(O), so that the following condition holds:
P( 1) -x(O) = P(0) -x(O) = 0.
I of
(2)
This is simply a problem of finding zeroes of non-linear functions. Since the non-linear
function itself is obtained through numerical integration, the problem has to be solved
iteratively: e.g., with a simple shooting method.
Before considering the shooting method, one can note that from equation (1) one
knows ?k -s will be distributed randomly in the interval TV E (0,l) and never reach Tr, = 1
or 0, so one has to use an extrapolation. Moreover, since rk = 1 and TV = 0 are in fact the
same thing, one can redefine rk = TV + 1 when Tk < I/2, so that then the Tks are in the
interval (l/2, 3/2), which is convenient for estimating P(1) by using an interpolation.
The application of the simple shooting method to this problem is illustrated schemati-
cally in Figure 1. Suppose one has a linear bi-excited system with the natural frequency
w, the steady state solution of which is x*(t) = cos t + cos fi, which is being sought.
Since there is usually a transient process, if one starts from an arbitrary initial point, one
will obtain a solution of the form
x(t) = C eey cOs(wr+e)+x*(t).
The iterative points of the PoincarC map P(Tk) of x*(f) are located on a bold circle in
the middle of the picture, and the successive iterative points of the PoincarC map P(Tk)
of Z(t) starting from an arbitrary initial point, say A, are located on a spiral which extends
MULTI-EXCITED QUASI-PERIODIC RESPONSES 295
.5)---x(kT,)
Figure 1. Poincare series of a quasi-periodic solution; P(r, 1 transient; bold circle, steady state.
towards the bold circle. The task is to adjust the initial point so that it will be located
on the bold circle. In the shooting process, one moves the initial point in a systematic
way, so that it approaches the bold circle rapidly. In order not to be misled, readers
should notice that Figure 1 does not exhaust all possibilities. For example, the spiral may
also start from a point inside approaching the bold circle from within. Moreover, if the
solution is unstable, the spiral starting near the circle will go away. Nevertheless, all these
cases can be treated uniformly with this method. As in a shooting for periodic solutions,
the distance between the point A and an intersection point of the spiral with the radial
?-k = l(0) (e.g., the first one) should be estimated and then gradually eliminated during
the shooting process. Since the iterative points are sparsely distributed on the spiral, it
is impossible to reconstruct the spiral by using an interpolation and henceforth to
determine a real intersection point. However, one may overcome this difficulty by assuming
there is an average value of the distances between the point A and several intersection
points and the average value of the derivatives of these distances with respect to the
initial conditions. These average values can be easily evaluated by using the interpola-
tion with rk as the argument. Fortunately, the shooting process runs smoothly with these
average amounts.
To perform the interpolation, Kaas-Petersen picks up a few points from the calculated
sequence P(rk), k= I, 2,. . . , which meet the following requirements: (1) by arranging
them in the order of an increasing index k, they form an increasing sequence of 7k; (2)
the value of Tk is located in the interval (1 - E, 1 + E), where E is a small positive number
representing the given error bound.
According to these requirements the points numbered 4, 7, 10 and 13 in Figure 1 are
used to interpolate the average intersection point and the corresponding derivatives.
In Figure 1 a curve through these points represents the interpolation curve, and the
intersection point of this curve with the radial T,, = l(0) can then be regarded as an average
value of P(1). It is evident that only a tiny part of the whole sequence is used, and this
may mean a waste of information gathered in the calculation.
296 F. H. LING
In fact, the essential factor for a smooth run of the shooting process is to evaluate a
proper average value of the difference P( 1) -x(O) and its derivatives with respect to
the shooting parameters. The accuracy of the interpolation in a normal sense is not
important. Therefore, it will be better if one uses all calculated P( TV ) valuest and their
derivatives to perform the interpolation, which are rearranged in the order of the increasing
rks. The interpolation curve of this new scheme is illustrated with a dashed line, and that
of the old one with a solid line, in Figure 2. The new interpolation scheme works very
well in practical calculations, as is shown in section 3. The new scheme enables one to
save a great deal of computer time, since, for a given accuracy, far fewer iterative points
are required than those in the old scheme. From intuitive considerations, we have also
tested other interpolation sets by excepting several bad points, such as 3, 8 and 14, or
the first few points, such as 1, 2 and 3, in Figure 2. But all these variants do not work
better. Since the calculation work of a numerical integration is dominant, to make use
of all calculated points means a great saving of computer time.
:
\
\
2 ;
\
/ \
\
/
\
ig
I 1.1
\
45
\
I
\ 012
2
C
Figure 2. Interpolation curve. -, Kaas-Petersen [8-IO]; - - -, this paper.
For ordinary differential equations, this means that the integration interval can be
greatly reduced, and therefore the integration error is smaller. As is well known, the error
in initial conditions grows exponentially with the length of the integration interval; see,
e.g., reference [ 111.
The second point is to calculate the derivatives of P(rk) with respect to the shooting
parameters directly instead of approximating them with a difference quotient. The main
advantages of this substitution are the following. Firstly, the Newton iteration when using
the derivative is of a quadratic convergence, and that when using the difference quotient
is of a (super)linear convergence. Secondly, it is not very easy to choose a proper increment
of the shooting parameters when using a difference approximation, since too large an
t Although this can be regarded as meaning that the error bound F equals 0.5 in the above-mentioned
requirement, nevertheless, since E = 0.5 means there is no error control at all, the idea behind the new algorithm
is different from that of the old one.
MULTI-EXCITED QUASI-PERIODIC RESPONSES 297
increment would cause a bad approximation, and too small an increment would cause a
severe round-off error, also leading to a poor approximation.
The new algorithm works especially well if the quasi-periodic solution is unstable. In
this case the absolute values of the derivatives are often very large, and the difference
quotient usually fails to be a reasonable approximation. This is, however, much less
serious if one calculates the derivative directly. The practical limitation for evaluating an
unstable solution with the improved method is that the instability should not cause an
overflow. This situation will, however, seldom appear in practical calculations.
The above-mentioned derivatives can be calculated as follows. For ordinary differential
equations (see references [ll, 121, for example)
i =f(x, t), .YE 5%,
(3)
the derivative matrix ~P(T~)/~s = [ax( f)/d~],_~~, is evaluated from an initial value problem
of the ordinary differential equations
ax(t)
1-l
afb, t) ax(t)
ax(t)
=~_
as ax as
[ 1
ax0
- =-
as
,=o
as
(4)
where ax, / as is determined according to the meaning of the parameter vector s. Although
one often takes s =x0 (then ax, / as = Z, , - uni t matrix), many other quantities, e.g., the
coefficients in equations, can also be taken as the components of s. In the latter case,
there will be some zero columns corresponding to these quantities in the derivative matrix.
The main calculation work of evaluating derivatives is to integrate n x n first order
ordinary differential equations (4). On the other hand, if one approximates derivatives
with difference quotients, one needs to integrate n first order ordinary differential equations
(3) n times with different initial conditions to find all the difference quotients. Therefore,
the amount of calculation in both schemes is almost the same. This statement is also true
for an iterative map.
The treatment of the iterative map y,,, =f( y,), y, E $?I, i = 1,2, . . . , where the index
denotes the iteration number, is similar but simpler. The derivative matrix af ( ~~) / as =
d,vi+klas is evaluated by using
aYi+l af ay,
-=--
~,
i=1,2 ,..., k-l,
as dyi as
and the starting value ay, / as is determined according to the meaning of s.
Various interpolation formulas can be used. For a bi-excited system, this is performed
by a one-dimensional interpolation. Lagrangian interpolation has been used in this work,
but a spline interpolation could be an alternative. If there are more than two exciting
frequencies, one can no longer use a Lagrangian or spline interpolation due to the random
distribution of the iterative points. A polynomial of two or more arguments has been
tried for this purpose, but the calculation was quite tedious.
After the initial point of a steady state quasi-periodic solution is decided, a time series
is generated for obtaining the Fourier coefficients by using the Fast Fourier transform or
a triangular function fitting. The latter one has been used in the work reported in this paper.
The stability of the solution is decided by the spectral radius of [ aP( ~~) / ax, ] . , =, . The
solution is stable if this radius is smaller than one, and is unstable if it is larger than one.
3. EXAMPLES
Several examples have been calculated. Although the purpose of developing this
algorithm is to apply it to non-linear problems, two linear examples are discussed first
298 F. H. LING
in this section. Since the aim is to show the advantage of the new algorithm, in the linear
case the exact solution and hence a convincing comparison is available. In the third
example results are presented for a non-linear oscillator with both stable and unstable
solutions. The non-linearity of this oscillator seems to be weak; however, the excitation
is strong enough so that the non-linear effect is apparent.
3.1. EXAMPLE 1: LINEAR OSCILLATOR
The oscillator equation is
jr+2a~++x=(/3-2)cosJZt-2JZasinJZt+(~-1)cost-2~sinf,
(5)
with T, =fin and T2 = 27r. For initial conditions x(O) =2,1(O) =O, the exact quasi-
periodic solution of equation (6) is known as
x=cos~t+cos t.
(6)
To illustrate the method equation (6) is rewritten in the form of equation (3) as
x1=x2,
~,=-2~~x~-~x,+(~-2)cos~t-2~~sinJZt+(~-l)cost-2~sint.
With s = {x,~, x*,,}~ the initial value problem (4) for evaluating the derivatives is
ax, ax,
8x2 ax2
ax, dx,
- - --
8x10 3x20 8x10
i I[-
axzo %o 8x20
1 0
ax, ax2 =
--
axlo axzo
_&.33-$
10 IO
-I[ 1
_2q33f5
,
2 ax,
=o [ 1 1
lo axzo t=.
In reference [8], the calculation was carried out for Q = 2, p = 1 in the interval
[0, KT,], K = 13. For E = 0.2, only four points 4,7,10 and 13 were used in the interpolation.
The resulting error is Ax(O) = 7.38E - 03 and A$O) = 3.79E -04. The results obtained by
using the new algorithm are given in Table 1. It can be seen that to reach a similar
TABLE 1
Accuracy of the numerically calculated
quasi-periodic solution of the ordinary
diferential equation (6) with a = 2,
p=1
K
Ax(O)
Ai(0)
20 -1*20E-10
19 -3.76E-10
18 -6*32E-10
17 4.78E-11
16 7.33E-10
15 1.38E-09
14 7.67E-10
13 4.25E-10
12 8*98E-08
11 7*06E-08
10 2.96E-06
9 5.00E-06
8 -2.83E-04
7 1.38E-04
6 1.85E-02
5 3.21 E-02
4.57E-10
1.40E-09
2.34E-09
-1.79E-10
-2.708-09
-5.28E-09
-1.38E-09
l.lSE-08
4.24E-10
-3.66E-07
1.78E-07
-6.37E-05
-8.33E-05
-9*23E-04
2*98E-03
-6*95E-02
MULTI-EXCITED QUASI-PERIODIC RESPONSES 299
accuracy, one needs only seven points; this means that almost 50% of the computer work
is saved. It can also be seen that with a too large K (K > 17 in this example), the accuracy
will be decreased due to the accumulation of the round-off error. This example shows
clearly that this algorithm works very well for a system with large damping.
3.2. EXAMPLE 2: DIFFERENCE EQUATION
The difference equation is
u(t+2)+u(t+1)+0~75u(t)=cos~~t,
(7)
or, rewritten in the iterative map form,
u,(t+ 1) = u*(t),
u*(t+1)=u,(t)-0~75u,(t)+cosGrt,
(8)
and possesses a steady state quasi-periodic solution
u(t)=Ccosfi&+Dsin&!rrt,
(9)
since the frequency of the solution fir is incommensurable with the basic frequency of
the map. After deciding on the constants C and B, one obtains
u,(O) = C = -1.09085 67916 34538,
u,(O) = C cos fin-+ D sin fir = 1.55572 71093 23350.
437 points were required in reference [ 81 for obtaining an accuracy of Au, (0) =1990E - 09
and Au,(O) =2*62E - 09. The present results are shown in Table 2. Only 14 points were
required for reaching the same accuracy-only 3% of the original calculation amount.
3.3. EXAMPLE 3: NON-LINEAR OSCILLATOR
The non-linear oscillator equation is
TABLE 2
Accuracy of the numerically calculated
quasi-periodic solution of the diference
equation (8) with a =2, p =1
K
Ax(O)
AX(O)
20
19
18
17
16
15
14
13
12
11
10
9
8
6
5
-1.35E-14 5.19E-15
-5*47E-15 -2.16E-15
-1.20E-14 7.048-15
8.93E-14 4.788-14
-3*75E-12 5.67E-12
-5.33E-11 -5.44E-12
-2.06E-10 2.96E-10
3.868-09 4.75E-09
-9.01 E-08 9.84E-08
-3.27E-07 -2.36E-07
-1.75E-06 2.88E-06
-9.79E-05 -4.88E-06
-1.228-05 -2.82E-04
- 1.39E-03 -1.99E-04
-2.12E-02 4.57E-02
-1.42E-01 -2.8OE-03
(10)
300 F. H. LING
with T, = 2rr and Tz = 27rlO.115 when 0 = 1. This equation possesses a periodic solution,
but its period is too large to deal with by a standard shooting method for evaluating a
periodic solution [12, 131. In the linear case, the solution has the form
x(t) = (0.15/a) sin t+A cos (0*115t)+ B sin (0.115t),
where A and B are determined by
(11)
l-o.l15z
-0.23
1_~2l315~][~]=[01~
(12)
With cr = 0.025, 44 points were required in reference [9] for obtaining an accuracy of
Ax(O) = 6*70E -05 and A.%(O) = 1.96E -07. The present results are shown in Table 3. Only
ten points are required for reaching the same accuracy-about 80% computational work
is saved. The results with a smaller damping cr = O-00005 and a negative damping (Y = -0.25
are also shown in Table 3. It is clearly seen from these results that this method is also
applicable to a system with very small, even negative, damping: i.e., for calculating
unstable solutions. If there is no damping at all, then one will have a tri-periodic solution
which depends on initial conditions.
TABLE 3
Accuracy of the numerically calculated quasi-periodic solution of the ordinary di$erential
equation ( 11)
(Y = 0.025 (Y = 0*00005 ff = -0.025
K
Ax(O)
A$01 Ax(O) A1(0)
Ax(O)
Ai(0)
20 -7*61E-16
19 -1.48E-16
18 6+51E-15
17 4.39E-14
16 -8.18E-12
15 -1.6lE-11
14 3.93E-10
13 7.22E-10
12 l.O5E-08
11 -2.14E-08
10 -1*08E-06
9 1.79E-06
8 1.27E-03
7 7.48E-03
6 8.26E-01
5 8.86E-01
6.208-16
1.43E-16
6.51E-15
3.11E-14
-4.99E- 12
-2.98E-11
4.13E-10
1.51E-09
9*77E-09
-3.13E-08
-4.09E-07
2.44E-06
-3.52E-04
1.33E-03
3.46E-01
3.27E-01
-6*37E-14
-5*5OE-15
-4.72E-13
9.22E-12
-1.79E-09
-3.49E-09
8.38E-08
1,43E-07
2.30E-06
-4.94E-06
-2.68E-04
5.06E-04
3.00E-01
2.07E-01
1.32E-00
8.52E-01
1.95E-13
2*80E- 14
2.20E- 12
4.84E-12
-1.47E-09
-7.27E-09
9.97E-08
3*45E-07
2.39E-06
-7.78E-06
-l.l2E-04
6.54E-01
-6.86E-02
-1.97E-01
4.79E-01
-6.71E-01
-2*67E-16
-5*OOE-18
-2.52E-16
-4.268-15
1*05E-12
2.558-12
-5.98E-11
-9.36E-11
-1*69E-09
3.85E-09
2.27E-07
-4.84E-07
-3.57E-04
-2.81 E-03
-9.79E-03
-1.63E-01
-l.l4E-16
-l*OOE-16
-5.74E-16
-1.37E-15
1.42E-12
5.97E-12
-8.OlE-11
-2.64E-10
-1.95E-09
6.50E-09
l.O3E-07
-6.03E-07
1.24E-04
9.87E-04
-2.23E-01
-6.99E-01
Next, one can fix c = 0.01, (Y = 0.025 and vary n for obtaining the resonance curve in
the non-linear case. After finding the quasi-periodic solution, a triangular function fitting
of the steady state time series is performed to obtain the resonance curves shown in
Figure 3 with solid lines. Apart from the basic harmonics with frequency 0 and O-1150
and amplitudes A, and A0.,,5 respectively, the two largest difference harmonics with
frequency 1.230 = (1 + 2 * 0.115)0 and 0-77fl= (1 - 2 * 0.115)fl and the corresponding
amplitudes A, .23 and A,.,, respectively are also shown in the figure. A resonance peak
can be seen near R = 1, which has a clear physical meaning.
MULTI-EXCITED QUASI-PERIODIC RESPONSES 301
a
Figure 3. Resonance curves of a non-linear oscillator, i+O~05i+x+O~Olx3 =0.3 cos Rr+ 1.5 cos (O.llSR1).
0, Harmonic balance; -, shooting method. A o ,,5 is the amplitude of the harmonic with frequency 0~115R.
and so on.
For comparison the resonance has been calculated by the harmonic balance method
with the trial function
x=A,.,,5~~~(0~115t+~,)+A,cos(t+~2), (13)
and the corresponding results are plotted in Figure 3 with dots. It is apparent that for
these harmonics the harmonic balance method gives pretty good results. However, the
error will increase when the amplitude becomes larger. Moreover, if one tries to calculate
other harmonics by taking more terms in the trial function, the derivation will be quite
tedious. It also can be noted that in order to determine the coefficients Ao.,15 and A,,
one has to solve some non-linear algebraic equations numerically. In view of the large
amount of derivation work when using the harmonic balance method, even in relatively
simple cases, the harmonic balance method does not have obvious advantages.
4. PERIODIC SYSTEM
The main purpose of developing this simple shooting algorithm is to calculate quasi-
periodic solutions in a multi-excited system, but it turns out that this method is also very
effective for calculating quasi-periodic solutions of a periodic system. In this case one
has an extra unknown (the winding number), and in order to find a physically meaningful
solution of the crucial equation (2), one solves the problem
I~P(O)-X(O)[[~= Min,
where 11 - [I2 denotes the Euclidean norm, by using an optimization algorithm (IMSL
routine UMINF). Our experience shows that if we have a proper guess of the winding
number (e.g., with a deviation less than 10% of the correct value), the shooting process
converges very rapidly.
302
The example
F. H. LING
i-p(l-x2)i+x=Fcoswt+F
(15)
is treated here. The Poincare section (with y = i) of the case with parameter values
CL = O-2, F=0*2, w = 1.1, Fo= 0.5 has been studied [2] and will be referred to as the
reference case shown in the center of Figure 4. The change of parameter values induces
different solutions as plotted in the other eight parts of Figure 4; each of them differs
from the reference case only by one parameter value as is indicated above the correspond-
ing part. Several phase-locking regions such as F0 3 0.75 and F 3 O-3 are found, where
there are only periodic solutions. The case of F = 0 has a trivial solution, x = 0, but it
possesses also a quasi-periodic solutions as shown in Figure 4. Moreover, the p = 0 case
represents a linear oscillator which does not have quasi-periodic solutions. In Figure 4
one sees 11 discrete points representing a P- 11 response as referred to the natural
frequency unity of the system. However, it is in fact a P- 1 response referred to the
excitation frequency w = 1.1.
Fio=oa
Winding = 0 06491
Reference w=I-3
Winding = 0.06761 Wlnding=0~06710
, , . -.. ,,,.
/ ,.
/
..-..
i
.
i
lJ
\
:*
\ c
. .
.
. . .
- .
,
\
I
I \
,
I
\
I
I
\
,
---
Fo=0.7
Winding = 0.09233
-Y
c,
II
-3 3 -3 3.
x
Figure 4. Quasi-periodic solution of the periodic system
i-~(l-u~)i+x=Fcoswr+Fb,
with the reference parameter values of p = 0.2, F = 0.2, w = 1. I and F0 = 0.5.
i
1
3
(16)
MULTI-EXCITED QUASI-PERIODIC RESPONSES 303
The winding number is very important especially in studying the transition from a torus
solution to chaos, and its precise value has been obtained here, being evaluated during
the iterative process. This is an obvious advantage in comparison with other methods
[ l-7,13].
5. CONCLUDING REMARKS
There are two efficient methods-shooting and collocation-for the boundary-value
problem of ordinary differential equations, and they are both useful for evaluating
quasi-periodic solutions. The shooting-type method can be applied to both a periodic
and a non-periodic system, while the collocation method can be used only for a periodic
system.
A new shooting-type algorithm in which all calculated iterative points of the PoincarC
map are used in the interpolation procedure and derivatives of the map with respect to
the shooting parameters are used has been presented in this paper. Usually, a rather
accurate bi-periodic solution is obtained by integrating over an interval of lo-15 shorter
periods, so that a great deal of computer time will be saved. A system with very small or
even negative damping, which results in an unstable quasi-periodic solution, can also be
treated with the algorithm. Moreover, this algorithm can also be successfully adopted for
calculating quasi-periodic solutions of a periodic system.
It is also possible to obtain the quasi-periodic solution boundary in the parameter
space as has been done for periodic solutions [ 141, but some problems can appear, since
it is not very easy to determine the spectral radius of the matrix corresponding to the
quasi-periodic solutions with a high enough accuracy.
There are also other methods in engineering applications for estimating a @asi-periodic
solution of a multi-excited system, e.g., a Galerkin type method [15] and the technique
in which the estimation of the quasi-periodic solution is converted to a boundary-value
problem of a partial differential equation, which is then solved by a finite difference
method [16]. But these methods do not seem to be able to compete with the shooting
method described in the paper. The cell-to-cell mapping method [ 171 is another useful
tool of studying quasi-periodic solutions; see reference [ 181 for an example of a study
of a circle map by using this method.
ACKNOWLEDGMENTS
This work has been partially supported by the Natural Science Foundation of China
and a grant of the Ministry of Education of China. The revised versions were produced
during the authors visits to Institute B of Mechanics of the University of Stuttgart, and
to the Department of Physics and Engineering Physics, Stevens Institute of Technology.
Thanks are due to W. 0. Schiehlen, E. Kreuzer and G. Schmidt for their hospitality. The
author also appreciates many suggestions made by G. Schmidt to improve the presentation
of the paper. The financial support for the visits was provided by the Ministry of Science
and Art of the State, Baden-Wiittemberg, West Germany, and the U.S. Department of
Energy, Contract No. DE-AC02-84ER 13146.
1.
2.
REFERENCES
M. HEAN 1980 International Journal of Non-linear Mechanics 15,367-376. SW la methode des
sections pour la recherche de certaines solutions presque periodiques de systemes forces
periodiquement.
E. THOULOUZE-PRAY 1981 Nonlinear Analytical Theory and Mathematical Applications 5,
195-202. Existence theorem of an invariant torus of solutions to a periodic differential system.
304 F. H. LI NC
3. E. THOULOUZE-PRAT-I- and M. JEAN 1982 International Journal of Non-Linear Mechanics 17,
319-326. Analyse numerique du comportement dune solution presque periodique dune
equation differentielle periodique par une methode des sections.
4. E. THOULOUZE-PRAY-~ 1983 Lecture Notes in Biomathematics 49. Berlin: Springer-Verlag.
Numerical analysis of the behaviour of an almost periodic solution to a periodic differential tori.
5. T. N. CHAN 1983 Master of Computer Science Thesis, Department of Computer Science, Concordia
Unioersity, Montreal, Canada. Numerical bifurcation analysis of simple dynamical systems.
6. I. G. KEVREKIDIS, R. ARIS, L. D. SCHMIDT and S. PELIKAN lo85 Physica D 16, 243-251.
Numerical computation of invariant circles of maps.
7. M. VAN VELDHUIZEN 1987 SIAM Journal of Scientific and Statistic Computations 8, 951-962.
A new algorithm for the numerical approximation of an invariant curve.
8. CHR. KAAS-PETERSEN 1985 Journal of Computational Physics 58, 395-403. Computation of
quasi-periodic solutions of forced dissipative systems.
9. CHR. KAAS-PETERSEN 1985 Journal of Computational Physics 64, 433-442. Computation of
quasi-periodic solutions of forced dissipative systems II.
10. CHR. KAAS-PETERSEN 1987 Physica D 25, 288-306. Computation, continuation, and bifurca-
tion of torus solutions for dissipative maps and ordinary differential equations.
11. J. STOER and R. BULIRSCH 1973 Einfuhrung in die Numerische Mathematik I I . Berlin: Springer-
Verlag.
12. F. H. LING 1981 Dissertation, Uniuersitiit Stuttgart, West Germany. Numerische Berechnung
period&her Liisungen einiger nichtlinearer Schwingungssysteme.
13. F. H. LING 1982 ZeitschrzftfTir angewandte Mathematik und Mechanik 62, TSS-T58. Numerische
Berechnung periodischer Liisungen nichtlinearer Schwingungssysteme.
14. F. H. LING 1989 Zeitschrzft fur Physik B 75, 539-545. A numerical study of the distribution of
different attractors in the parameter space.
15. L. 0. CHUA and A. USHIDA 1981 I nstitute of Electrical and Electronic Engineers Transactions,
Circuits and Systems CAS-28, 953-971. Algorithms for computing almost periodic steady-state
response of nonlinear systems to multiple input frequencies.
16. B. WEYH 1982 VDI -Bericht Nr. 456, 113-120. Berechnung nichtlinearer quasiperiodischer
Schwingungenen in Maschinenskzen mit Schlupfkupplungen.
17. C. S. HSU 1987 Cell-to-Cell Mapping. New York: Springer-Verlag.
18. F. H. LING and Z. R. LIU 1990 (Preprint). Limiting probability density of a quasi-periodic orbit.