- Email: [email protected]
Generalized Fourier series and limit cycles of generalized van der Pol oscillators
Journal
ofSound nnd
Vibration (1990) W(3),
GENERALIZED
FOURIER
GENERALIZED
453-466
SERIES AND
VAN
DER
J. GARCIA-MARGALLO Departamento
LIMIT
OF
POL OSCILLATORS
AND
J.
D.
BEJARANO
de Fi’sica, Facultad de Ciencias, Universidad de Extremadura,
(Received
CYCLES
5 December 1988, and in revised form 2
06071 Badajoz, Spain
June 1989)
A very simple generalization, using Jacobian elliptic functions, of the usual Fourier series, appropriate for non-linear systems, is used to study the first order approximate solution of the generalized van der Pol oscillators X +AX +2BX3 + F( z3 + z,X’ + z,X4)A = 0. A generalized harmonic balance method is used to determine the limit cycles. The cases A > 0, B > 0 and A < 0, B > 0 are considered in some detail. For given values of the parameters z,, the values of A and B for which limit cycles exist are found as functions of m. Numerical values for the radius, frequency, and energy of the limit cycles are given. The presence of zero, one, two, three or four limit cycles depends on the value of the parameters of the equation. Stability, bifurcations of fixed points, and the limit cycle stability are studied qualitatively. 1. INTRODUCTION From elementary physics to quantum field theory, the fundamental functions most used to approximate experimental facts (e.g., the “fundamental oscillators”) are solutions of the equation *+AX=O.
(1.1)
In applied physics most approximation methods, such as Fourier series, depend heavily on the oscillatory solutions and on the Laplace and Fourier transforms of the exponential solutions. According to Legendre [l], Euler [2] first proposed, more than 200 years ago, the equation that generalizes equation (l.l), namely %+AX+2BX3=0.
(1.2)
To our knowledge, our group was the first [3] to put forward the exponential solutions of equation (1.2) for all signs and values of A, B, E, and all X regions of possible interest in classical or quantum mechanics, where E is given by the more convenient “energy equation” *2+AX2+
BX~=
E,
x =
dx/dt.
(1.3)
On the basis of these solutions we have been able to develop the corresponding Fourier series and transforms [4]. We have used [5-71 these generalized Fourier series with the harmonic balance method to find an approximate solution of equations of the type X+x’=E(l-X2)X. The purpose of this paper is to determine the first order solutions by the method of harmonic balance using the generalized Fourier series and the Jacobian elliptic functions 453 0022-460X/90/030453
+ 14 $03.00/O
0
1990 Academic PressLimited
454
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
in equations of the type X +g(x) To avoid mathematical case of g(X)=AX+2BX3 der Pol equation,
= &j-(X,X).
(1.4)
difficulties as much as possible, we limit our study to the simple andf(X, X) =(z3+z2X2+z,X4)X: i.e., to a generalized van
X+AX+2BX3+e(z3+z2X2+z,X4)X=0.
(1.5)
Equations (1.4) or (1.5) with E =0 are called generating equations, and their solutions generating solutions. With B = 0 the approximate solutions for the many different cases are given and discussed in the well known textbook by Minorsky [8]. When z, = 0, solutions of these equations are also well known, and we compare our results and limit cycles with those of other authors [9-111 for this very special case and with numerical integration results for other quite complicated and interesting cases. Due to the X4X term there are zero, one, two, three or four limit cycles depending on the value of the parameters. Following references [lo-121, we have studied the stability of the limit cycles and the bifurcations. 2. THE GENERALIZED FOURIER
SERIES
METHOD
To avoid computational difficulties we limit our discussion to a first order approximation. Due to the fact that the usual Fourier series for the elliptic functions contain higher harmonics, this very simple approximation has proved to be quite good in the cases we have studied. In the potentials we are considering, V(X) = AX’+ BX4 for E < V,,, = A2/4B. We assume a solution of the form X(t) = a pq (wt; k2),
(2.1)
where a and w are constants to be determined, and pq denotes a convenient Jacobian elliptic function determined by the generating solutions of equations of the type (1.4). Substituting equation (2.1) into equation (1.4) gives F,(a, w, k2, E, a) cos cp+ F2(a, w, k2, E, a) sin cp+ (higher order harmonics) = 0,
(2.2) where cp is the generalized circular function of the case being considered [4], and (Y collectively denotes any parameter which appears in the non-linear function f( X, X). We first take F, = 0 and F2 = 0, using the method of harmonic balance [ 131. We shall study the different cases of the quartic anharmonic asymmetric oscillator (AS0 and AAO) with the form of the potential given above: i.e., with E = 0 for X = 0. We shall distinguish two cases of the quartic oscillator: (I) two types of potential: A > 0 and B>O, and A-CO, B>O, E>O; (II) A<0 and B>O, and E
OF TWO TYPES
OF QUARTIC
OSCILLATOR
3.1. OSCILLATOR TYPE 1 For this type of oscillator, following references [6,7], we take a function X = a cn (wt; k2),
(3.1)
GENERALIZED
VAN
DER
POL
where a, o and k* = m are constants determined to first order approximation. expression (3.1) into equation (1.5) gives (aA+2amw2-uw2)
455
OSCILLATORS
cn+(2a3B-2umw*)
Substituting
cn3- EUWZ~ sn dn
- cu3wz2 cn* sn dn - &u5wz, cn4 sn dn = 0. Our generalized [14, 151, gives
Fourier expansion,
[u(A+2mw’-w2)+l~5u(Bu2-mw2)] - [ ~uoz,B,(
if the expansion
(3.2)
is limited to the first harmonic
cos cp3
m) + ~u~wz,B,( m) + EU~CIJZ, B,( m)] sin 40~+ (higher harmonics) = 0,
(3.3) where B,(m), B,(m) and B3( m) are given in Appendix I. Setting the coefficient of sin (p3 to zero, one obtains -[.ww(u4z,B,(m)+u2z2B2(m)+z3B,(m))]=0,
and, for u2 = y, y2z,B,(m)+yz,Bz(m)+z3B3(m)=0.
(3.4)
The coefficient of cos rp, can be most simply made zero,
w2= Bu2/m
and
w’=A/(l-2m),
(3.5)
and from these two equations
A/B= (1-2m)u2/m. 3.2.
OSCILLATOR
TYPE
(3.6)
11
For this type of oscillator the solutions corresponding [6,7] are
to a first order approximation
X = a dn (wt; m).
(3.7)
Substituting expression (3.7) into equation (1.5) gives (uA+2uw2-uw2m)dn+(2u3B-2uw2)dn3-Euwmz,sncn - .a3wmz2 dn* sn cn - .w5wmz, dn4 sn dn = 0.
Our generalized Fourier expansion, if the expansion is limited to the first harmonic, gives [u(A+~w~-~w~)+~~~u(Bu~-~~)]
cos (p2
-[[~u~m(z~~~(m)+u*z~93~(rn)+u~z,9?,(m))]
sin cp,+(higher harmonics) =O, (3.8)
where a,(m), L%,(m) and a,(m) sin (p2to zero gives
are given in Appendix
II. Setting the coefficient oc
-[~u~m(z3%13(m)+u2z2~2(m)+u4z,93~(m))]=0,
and, for u* = y, y2z,CB,(m)+yz2B~(m)+z3W3(m)=0.
(3.9)
Setting the coefficient of cos rp2to zero gives w2 = Bu*
and
w’=A/(m-2),
(3.10)
and from these two equations A/B = (m -2)~~.
(3.11)
456
.I. GARCIA-MARGALLO
AND
4. ANALYTICAL
J. D.
BEJARANO
RESULTS
We assume that E -0 always. Consider the parameters z, , z2, and z3 as co-ordinates axes of the three-dimensional space in which we shall study the solution of equations (3.4) and (3.9) in the eight octants, and the following six special cases: (1) z, = z2 = 0; (2) z, = z3 = 0; (3) z2 = z3 = 0; (4) z, = 0; (5) z2 = 0; (6) z3 = 0. As before, we study all cases in first order approximation. 4.1. THE SIX SPECIAL CASES Case (l), z, =z2=0. In this case equation (1.5) becomes X+AX+2BX3+~(z,)X=0.
(4.1)
C.1.1.: oscillator type I. (a) For A > 0,B > 0,with 0 G m S l/2. There is one fixed point (f.p.) at X = 0, X = 0: i.e., (0,O). Let us study its stability. If z3 > 0, the positive coefficient of X means physically dissipation of energy, so that the f.p. is asymptotically stable (S), a sink. The topological configuration (t.c.) is S. If z, 0,with l/2 6 m s 1, and E > 0, equation (4.1) has three f.p. at (0, O), (*w, 0). If z3 > 0, the f.p. are a saddle at (0,O) and two sinks at (km, 0). If z3 ~0, they are one saddle and two sources. Here there also exists one bifurcation of equilibria as z3 passes through zero. Equation (4.1) undergoes a pitchfork bifurcation of equilibria as A passes through zero. C.1.11.: oscillator type II. (a) For A < 0,B > 0,and OS m s 1, the situation is analogous to (C.l.I.b) in the dnoidal oscillator. Case (2), z, = z3 = 0, and case (3), z2 = z3 = 0, are analogous to case (1). Consequently from the physical point of view, cases (l), (2) and (3), having only one coefficient of X, are simply damped if the coefficients are positive and excited if negative, with no possibility of limit cycles. Case (4), z, =O. As noted in the Introduction, [ll]. Equation (1.5) becomes
this case is well known in the literature on the subject
X+AX+2BX3+e(z3+z2X2)X=0.
(4.2)
C.4.1.: oscillator type I. (a) One f.p. at (0,O). Equation (3.4) that gives the stationary state is yt2B2 + z3B3 = 0, where y= a 2 =-z3B3/z2B2.
(4.3)
The limit cycles (l.c.) can exist only if sign (z2) # sign (z3). Due to the fact that B3/B2 is a positive function of m only, a2 depends only on the scale factor -~3/t2 and on this function. In other words it is only necessary to calculate as we do in Figure 1 for one value of z3/zZ to obtain the general behaviour. For given z2 and z3, from equations (3.6) and (4.3) we obtain a2 and A/B us.m. In Figure 1, A/B is shown us. m (for z2 = 1.0, z3 = -3.O), giving the l.c. of equation (4.2). Observe that this function has extremes: i.e., points in a region of (A/B) with two l.c. that coalesce at the extreme of the function. Let us now study the stability and the bifurcations, that is, the t.c. for different values of the parameters. If z2> 0, and z3 < 0, the f.p. is unstable, U. From Figure 1, there exists one limit cycle for each A/B which is stable. The t.c. is US. If z2 < 0, and z3> 0, the point
GENERALIZED
VAN
DER
457
POL OSCILLATORS
2374-!----l vg5F---=-T 18 -
18-
17 -
1716-
1615 -
15-
I
14 -
14-
13 -
13-
A>0
12 -
12-
B>O
11 -
11 lo-
lo9-
9’0
8-
? q
9-
A>0
8-
7-
7-
6-
6? q
= 4-
54-
Cfl
3-
32-
\
1
2l-
l0 -1 -2-3 -4-5-6 -7 -8-
_;:L
O-1 0,2 0.3 0.40.5 060.7 0.8 0.9 1.0 m
Figure 1. The ratio A/B versus m, case (4),for the solution cn, dn. z, = 0, q = 1 and z3 = -3.
-,0-
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 nl Figure 2. The ratio A/B versus m, case (6), for the solution cm, dn. .zI = 1, z2 = -3 and z, = 0.
of rest is stable, S, and the l.c. is U. The t.c. is SU. Equation (4.2) undergoes a bifurcation of equilibria as .zJ passes through zero. (b) Equation (4.2) has three f.p. at (0,0), (*w, 0); if z2 > 0, and z3 < 0, it has one saddle point and two sinks; if z2< 0, and z3> 0, it has one saddle point and two sources. From Figure 1, the l.c. have four regions in the oscillator: (1st) for each value of A/B one l.c. only exists in the cnodial oscillator; (2nd) for each A/B two l.c. exist, one in the cnodial and another in the dnoidal oscillator; (3rd) for each value of A/B there are two l.c. in the cnoidal oscillator with different m; (4th) for the minimum value of A/B there is one limit cycle in the cnodial oscillator, from coalescence of the two in the (3rd) region. The stability of the l.c. is as follows: (1st) if z2 > 0, z3 < 0, S; if z2 < 0, z3 > 0, U; (2nd) as for the 1st region; (3rd.a) if z2 > 0, z3 < 0, S; if z2 < 0, z3 > 0, U; (3rd.b) if z2> 0, z3 ~0, U; if z2 CO, z3 > 0, S; (4th) the l.c. are semi-stable, US. Equation (4.2)
458
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
undergoes a bifurcation of equilibria as z3 passes through zero. When A passes through zero, one saddle-node coalescence occurs. The frequency (w), from equation (3.5), and the total energy (E) [3] of the l.c. are given by 02= Ba2/m and E = rn,a2w2. C.4.11.: oscillator type II. (a) Equation (4.2) has three f.p., at (0, 0), (+dm, 0). Equation (3.9) for the l.c. becomes yz2g1+ %13z,= 0, so that y=a2=-z,~3/z2532.
(4.4)
The stationary amplitude can exist only if sign (z2) # sign ( z3). The equation is the same as for oscillator I but now the %33k functions are those of oscillator II given in Appendix II. For given z2 and z3, from equations (3.11) and (4.4) we obtain A/B us. m. The l.c. in this case is given in Figure 1 only for certain values of A/B. For these values l.c. also exist in the cnoidal oscillator (C.I.b.l), with stabilities. If z2 > 0, and z3 < 0, one saddle point, two sinks, and one dnoidal l.c. that is unstable. The dnoidal t.c. is SU. If z2 < 0, and z3 > 0, one saddle point, two sources, and one dnoidal l.c. that is stable. The dnoidal t.c. is US. For this case there also exists a bifurcation of equilibria as z3 passes through zero. If m + 1 in (C.4.1.b) and (C.4.11.a) one homoclinic cycle exists. The frequency w, from equation (3.10), and E are given by o2 = Ba* and E = -m,a202. Case (5), z2 = 0. Equation (3.4) reduces to y2z, B, + z,B, = 0, so that the radius of the l.c. is now given by a4= -z,B,/z, B,. Also here B, and B2 are always positive, so that the stationary amplitude can exist only if sign (z,) # sign (z~). The t.c. are the same as in case (4). Case (6), z3 = 0. Equation (1.5) becomes X+AX+2BX2+s(z2X2+z,X4)X=0.
(4.5)
C.6.1.: oscillator type I. (a) There is one f.p. in (0,O). Equation (3.4) gives y’z,B, + yz2B2 = 0, or y(yz, B, + z2B2) = 0, so that y = 0 and y=a
2
=-z2B2/z,B,.
(4.6)
The stationary
amplitude can exist only if sign (zi) # sign ( z2). We shall explain the y = 0 root later when studying the second octant of the general case. For given z, and z2, from equations (3.6) and (4.6) we obtain a2 and A/B us. m. In Figure 2, A/B is shown us. m: i.e. the l.c. of equation (4.5). Observe that these functions have no extremes. If zI > 0, and z2 < 0, the point (0,O) is U, and the l.c. is S. The cnoidal t.c. is US. If z1 CO, and z2> 0, the cnodial t.c. is SU. Equation (4.5) undergoes a bifurcation of equilibria as z2 passes through zero. (b) Equation (4.5) has three f.p. at (0,O) and (fw, 0). From Figure 2 we have one l.c. region in the cnodial oscillator. In this region there exists only one l.c. for each A/B. If z, > 0, z, < 0, one saddle and two sources. In the cnoidal oscillator, one l.c. that is stable. The cnoidal t.c. are US. If z, < 0, z2 > 0, one saddle and two sinks. In the cnoidal oscillator, one l.c. that is unstable. The cnoidal t.c. are SU. When A passes through zero, one saddle-node coalescence occurs. The o and E are given by w2 = Ba’/m and E = m 1a2w2 .
C.6.11.: oscillator type II. (a) From Figure 2, limit cycles exist only of A/B. If z, > 0, and z2 < 0, one saddle, two sources, and one l.c. that t.c. are US. If z, < 0, and z2 > 0, one saddle, two sinks, and one l.c. that t.c. are SU. For z2 passing through zero there is one bifurcation. If m +
for certain Values is S. The dnoidal is U. The dnoidal 1 in (C.6.1.b) and
GENERALIZED
VAN
DER
POL
OSCILLATORS
459
(C.6.II.a), the point (0,O) is a saddle, the l.c. coalesce and one homoclinic cycle appears. The w and E are, respectively, o*= Ba* and E = -m,a*w*. 4.2. GENERAL CASES Octant (1): z,>O, z,>O, z,>O. Equation (1.5) becomes ~+AX+2BX3+~(z3+z,X2+z,X4)~=0.
(4.7)
Oc.l.1.: oscillator type I. (a) Equation (3.4) becomes y’z,B, +yz2B2+ z,B, = 0, and has no positive roots, the only stationary state is the f.p. which is S. (b) The f.p. are one saddle and two sinks. Oc.l.11.: oscillator type II. (a) This is analogous to (0c.l.I.b). Octant (2): z, > 0, z2 < 0 and z3 > 0. Oc.2.1.: oscillator type I. (a) The f.p. is stable, S. Equation (3.4) becomes y*z, B, + yz,B, t z,B, = 0, and this polynomial has two roots, (4.8) a2 = *Y = (-z2B2+sqfl ((z2B2)*-44(zIB1)(z3B3)))I(2zIB1), where there are three cases, as follows. (1) If (-z2B2)‘>4(z,B1)(z3B3). For given zl, z2, and z3, from equations (3.6) and (4.8), we obtain A/B vs. m. The limit cycles for this case are given in Figure 3. Hence the two roots are positive. Of these roots the smaller one, a*-, corresponds to an unstable, and the larger, a*+, to a stable limit cycle. The cnoidal t.c. is SUS, bicyclic. If z3 = 0 the smaller l.c. has radius 0.0. It then coincides with the fixed points and we have only the larger l.c. with radius y = -z2B2/zI B, that coincides with special case (6). (2) If (-z,B,)’ < 4(z, B,)(z,B,), the polynomial has no positive roots, and the only stationary state is the point of rest, S; the cnoidal t.c. is acyclic S. (3) If (-z,B,)* = 4(z, B,)(z,B,) there appears a surface, 5, separating these configurations. The t.c. passes from SUS (or bicyclic) to S (acyclic). The surface 4 is thus the locus of the bifurcation points. Topologically this is also a locus of semi-stable cycles US. (b) Here the f.p. are one saddle and two sinks. Equation (3.4) has also three cases, as follows. If (-z2B2)* > 4(z,B,)(z,B,). From Figure 3, we have five l.c. regions: (1st) one small, U, one large, S, root of the cnoidal oscillator for each A/B; (2nd) one small, U, one large, S, root of the cnoidal oscillator and one small root of the dnoidal oscillator for each A/B; (3rd) one large root, S, two small roots, U, but different m, of the cnoidal oscillator for each A/B; (4th) one large root, S, and coalescence of two l.c. of small roots, SU, of the cnoidal oscillator for each A/B; (5th) one large root only, S, of the dnoidal oscillator for each A/B. The other cases are analogous to (Oc.4.I.a). We obtain bifurcations when A/B passes from one to another of these five regions. When A passes through zero a saddle-node coalescence occurs. The w and E are o*= B(a**)/m and E = m,(a**)w*. Oc.2.11: oscillator type II. (a) The f.p. are one saddle, and two sources, U. Equation (3.9) becomes y*z,%, +yz2%z+ z3B3 = 0. This polynomial has two roots:
a2 = *.v =
(-z222*sW
((~~~2)*-4(2,~~)(~~~~)))/(2~,98,).
(4.9)
Consequently, we have the same three cases as before: (1) If ( -z,L%~)*> 4(~,L%,)(z#~), algebraically there are two roots and both are positive. For given zl, z2 and z3, from equations (3.11) and (4.9) we obtain A/B vs. m. The l.c. for this case is given in Figure 3 and has two regions: (1st) the (2nd) case described previously; (2nd) one large root, S, in the dnoidal oscillator for each A/B.
460
J. GARCIA-MARGALLO
AND
J. D.
BEJARANO
”
”
”
1’
‘-
+
L.
1017 16 15-
A'0 8~0
1413 12 ll10 98769 q
5432l0 -1
-
-2
-
-3
-
-4
-
A<0 B>O
-5 -6
-6 A<0 B>O
-7
-7 -6
-6
_;:L
04 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1 .O m
Figure 3. The ratio A/B versus m, Ott.(2), for the solution cn, dn. 2, = 1.2, z2 = -3.11 and z3 = 1.
(2) If (-z2W,)* <4(z,.93,)(z3~~),
-9 -10
dn A 0.4 0.2 0.3 04 0.5 0.6 0.7 0.0 O-9 I.0 In
Figure 4. The ratio A/B versus m, Ott.(3), solution cn, dn. z, = 1, z2 = -4 and zx = -4.
for the
the polynomial has no real roots, the only stationary state is (km, 0); the dnoidal t.c. is acyclic S. there appears a surface, 0, separating these configura(3) If (-z,93J2 = 4(z,9,)(zj98J tions. Algebraically, this is a locus of double roots: u* = -z29B2/2z,B,. Topologically this is also a locus or semi-stable US. The surface, 0, separates the two possible configurations SUS and S. Equation (4.7) undergoes a bifurcation in the surface, 0, as (z&%)* passes through 4( ~,a,)( z3B3). In (Oct.2.I.b) and (Oct.2.II.a) when m + 1, the limit cycles coalesce for the small and large roots in the saddle point (0, O), in two homoclinic cycles. The w and E are given by W*= B(a’*) and E = -m,(uz*)~*. Octant (3): z,>O, z2C0, z,
GENERALIZED
VAN
DER
POL
OSCILLATORS
461
(a) The f.p. is unstable, U. Equation (3.4) is y’z, B, + yzZB2+ z3B3 = 0. This polynomial has only one positive root: a2= ?? Y = (-zJ&+sqrt
((z~B~)~+~(z,B~)(z~B~)))/(~z,B~).
(4.10)
For given z, , z2, and z3, from equations (3.6) and (4.10), we obtain A/B us. m. Figure 4 gives the limit cycle. The cnoidal limit cycle is stable. The t.c. is US. As was the case for the second octant, when z3+0 the l.c. radius has the minimum value for this octant of u2 = -z2B2/z,B,, analogous to the special case (6). (b) The fixed points are one saddle and two sources. From Figure 4, we have one l.c. for each A/B, which is S. When A passes through zero, a saddle-node coalescence occurs. The w and E are given by w2 = Ba’/m and E = m,a2w2. Oc.3.11.: oscillator type II. (a) The three f.p. are one saddle and two sources. From Figure 4 we have one l.c. for each A/B, which is S. The A/B (dnoidal) < A/B (cnoidal);
““V-----l Y.
-7%
1917 16 15
-
14 1312
-
A'0
11 -
B>O
10 9676-
m G
54321 0 -1
-
-2
-
-3
-
-4
-
-5 -6
-
-7
-
-6
-
-9
-
-lO-
Figure
5. The ratio A/B
dn
A< 0 BP0
’ ’ ’ ’ ’ ’ ’ ’ 0.1 0.2 0.3 0.40.5 0.6 0.7 0.0 09 1.0 m
versus m, Ott.(4),
for the solution
cn, dn. z, = 1, z2 = 1 and I) = -3.
462
J. GARCIA-MARGALLO
AND
J. D.
BEJARANO
they are equal only for m = 1. For m = 1, in (Oct.3.I.b) and (Oct.3.II.a), there exists one homoclinic cycle. The w and E are given by w2= Ba2 and E = -m,a2w2. Octant (4): z, > 0, z2 > 0, z3 < 0. Oct.4.1.: oscillator type I. (a) This is analogous to (Oct.3.I.a). (b) the f.p. are one saddle and two sinks. From Figure 5, there are four regions: (1st) one l.c. for each A/B, which is S; (2nd) two l.c. with one in the cnoidal oscillator, S, and the other in the dnoidal oscillator, which is US, for each A/B; (3rd) two l.c. in the cnoidal oscillator for each A/B, but with different m, both S; (4th) coalescence of the two previous l.c. for the minimum value of A/B. A saddle-node coalescence occurs when A passes through zero. The w and E are given by w2 = Ba2/ m and E = m,a202. Oct.4.11: oscillator type II. (a) The f.p. is analogous to (Oct.4.I.b) and the l.c. is analogous to (Oct.4.I.b.11). Homoclinic cycles exist when m + 1 in (Oct.4.I.b) and (Oct.4.II.a.). The w and E are given by w2 = Ba2 and E = -m1a202. For the octants of the lower half-space (zi 0, z2> 0; (6) z, CO, z3> 0, z, 0. The conditions of the stability are reversed both for the point of rest and for the stationary motion, whereas the polynomial remains the same, and this leads to inverted topological configurations: i.e., octant (5) = (3), (6) = (4), (7) = (1) and (8) = (2). The figures of A/B versus m and ZJZ, are the same.
5. COMPARISON WITH NUMERICAL INTEGRATION It is instructive to compare these simple analytical approximations with a numerical integration of the equations. We used the fourth order Runge-Kutta method and illustrate the more important results with the value E = O-1. The results for the limit cycles are shown in Figure 6 for the special case (C.4.1.a) of the oscillator type I, when A/B = 23.125, z1 = 0, z2 = 1 and z3 = -3, and from Figure 1 one has rn = O-25, and from equation (4.3), a = 3.40. The (a) curve is the analytical
:5
Figure 6. The (a) curve is the analytical solution X = 3.40 cn (6.80r; m = 0.25). The (b) curve is the numerical integration with initial conditions X( t = 0) = 3.40 and A( t = 0) = 0.00.
GENERALIZED
VAN
DER
POL
OSCILLATORS
463
solution X = 3.40 cn (6.80t; m =0.25). The (b) curve represents the results of the numerical integration with initial conditions X( t = 0) = 3.40 and X( 1= 0) = 0.00. In Figure 7 are shown (C.4.1.b) and (C.4.11.a) for A/B = -6.70, z, = 0, z2 = 1, zj = -3; from Figure 1 for the cnoidal oscillator, M = 0,755, and from equation (4.3) a = 3.15. For the dnoidal oscillator M = 0.905, and from equation (4.4) a = 2.47. The (a,) curve is the analytical solution X = 3.15 cn (3.62t; WI= O-755). The (az) curve is the analytical solution X = 2.47 dn (2.47t; m = O-905). The (b,) curve represents the results of the numerical integration with initial conditions X(t = 0) = 3.15 and X(t =0) =O*OO. The (b,) curve represents the results of the numerical integration with initial conditions X( t = 0) = 2.47 and X(t=O)=O*OO.
Figure 7. The (a,) curve is the analytical solution X =3.15 cn (3.621; m =0.755). The (a,) curve is the analytical solution X = 2.47 dn (2.47r; m = 0.905). The (b,) curve is the numerical integration with X( t = 0) = 3.15 and X(? = 0) = 0.00. The (b,) curve is the numerical integration with X( t = 0) = 2.47 and X( I = 0) = 0.00.
In Figure 8 is shown the case (Oc.2.I.a.l) for A/B = 5.05, and z, = 1.2, z2 = -3.11, z3 = 1; from Figure 3, m = 0.23 for a-, and m = 0.27 for a+. The (a,) curve is the analytical solution X = 1.47 cn (3*06t; m =0*23). The (b,) curve represents the results of the numerical integration with initial conditions X( t = 0) = 1.47, X( t = 0) = 0.00. The (a,) curve is the analytical solution X = 1.72cn (3.31 t; m = 0.27). The (b2) curve, the numerical integration with initial conditions X( t = 0) = 1.72, X( t = 0) = 0.00. In Figure 9 are shown the cases (Oc.2.I.b) and (Oc.Z.II.a), for A/B = -0.952, and z, = 1.2, z2 = -3.11, z3 = 1; from Figure 3, m = 0.59 for a+, and m = 0.70 for a-, in the cnoidal oscillator, and m = O-965 for a- in the dnoidal oscillator. The (a,) curve is the analytical solution X = 1.76 cn (2.30t; m = 0.59). The (az) curve is x= 1.29 cn (1*54?; m = 0.70). The (aj) curve is X = 0.95 dn (0.95t; m = O-965). The (b,) curve is the numerical integration with initial conditions X( t = 0) = 1.76, X(t = 0) = 0.00, the (b2) curve is that with X( t = 0) = 1.29, X( t = 0) = 0.00, and the (b,) curve is that with X( t = 0) = 0.95, X( t = 0) = 0.00. All the analytical results seem to be very good for E s 0.1.
464
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
Figure 8. The (a,) curve is the analytical solution X = 1.47 cn (3.061; m =0.23). The (b,) curve is the numerical integration with X( t = 0) = 1.47, X( f = 0) = 0.00. The (a2) curve is the analytjcal solution X = 1.72 cn (3.3 1 t; m = 0.27). The (b,) curve is the numerical integration with X( r = 0) = 1.72, X( r = 0) = 0.00.
Figure 9. The (a,) curve is the analytical solution X = 1.76 cn (2.30r; m = 0.59). The (a,) curve is X = 1.29 cn (1.54r; m =0.70). The (a,) curve is X = 0.95 dn (0.95r; m = 0.96). The (b,) curve is the numerical integration with X(r=O)= 1.76, k(r=O)=O.OO. The (b2) curve is with X(r=O)= 1.29, X(r=O)=O.OO. The (b,) curve is with X( r = 0) = 0.95, X( r = 0) = 0.00.
GENERALIZED
VAN
DER POL OSCILLATORS
465
6. CONCLUSIONS We have described how to generalize the method of harmonic balance to obtain the first order approximations to the periodic solutions of differential equations of the type (1.4), using elliptic functions, and we have also studied the stability of the fixed points, and of the limit cycles and bifurcations, from a qualitative viewpoint. We have applied the method to the study of the limit cycles of equation (1.5). The presence of zero, one, two, three or four limit cycles depends on the value of the parameters z,, as shown in Figures l-5. We have compared the results with numerical integration results in Figures 6-9 and this comparison shows the analytic approximations to be very good.
ACKNOWLEDGMENTS Thanks are due to the work Dr R. Chatwin has contributed to improvements of this and previous papers [6,7] on non-linear mechanics. We gratefully acknowledge the Comisi6n Asesora de Investigacibn Cientifica y Ttcnica (CAYCIT), Spain, project no. PB87-0007, for financial support and J.D.B. would like to acknowledge the Plan Movilizador, and the E.P. Division at CERN for its hospitality.
REFERENCES 1825 Traits’des functions elliptiques et des intkgrales eule’riennes. Paris: 1. A. M. LEGENDRE Huzard-Courcie. 2. L. EULER 1763 Nov. Corn. Petropotan X pg 4. Id toms VI, VII, 1761. 3. J. D. BEJARANO and A. M. SANCHEZ 1988 Journal of Mathematical Physics 29, 1847-1853. Generalized exponential, circular and hyperbolic functions for non-linear quantum mechanics. 4. J. D. BEJARANO and A. M. SANCHEZ 1988 Journal of Mathematical Physics 30, 1871-1876. Generalised Fourier transforms for nonlinear systems. 5. S. B. YUSTE and J.D. BEJARANO 1986Journal of Sound and Vibration 110,347-350. Construction of approximate analytical solutions to a new class of non-linear equations. and J. D. BEJARANO 1987JournalofSoundand Vibration 116,591-595. 6. J. GARCIA-MARGALLO A generalization of the method of harmonic balance. 7. J. GARCIA-MARGALLO,J.D.BEJARANO~~~S.B.YUSTE~~~~ Journalofsoundand Vibration 125, 13-21. Generalized Fourier series for the study of limit cycles. 8. N. MINORSKY 1962 Nonlinear Oscillations. Malabar: Litton Educational Publishing, Inc., U.S.A. 9. P. HOLMES and D. RAND 1980 Journal of Non-Linear Mechanics 15,449-458. Phase portraits and bifurcations of the non-linear oscillator: jl + (a + yx*),? + /3x + Sx3 = 0. 10. J. GUCKENHEIMER and P. HOLMES 1986 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. 11. E. KNOBLOCH and M. R. E. PROCTOR 1981 JournalofFluid Mechanics 108,291-316. Nonlinear periodic convection in double-diffusive systems. 12. D. W. JORDAN and P. SMITH 1987 Nonlinear Ordinary Differential Equations. New York: Oxford University Press. 13. R. MICKENS 1987Journal of Sound and Vibration 111,515-518. A generalization of the method of harmonic balance. 14. M. ABRAMOWITZ and J. A. STEGUN (editors) 1972 Handbook of Mathematical Functions. New York: Dover. 15. P. D. BYRD and M. D. FRIEDMAN 1971 Handbook of Elliptic Integrals for Engineers and Scientists. Berlin: Springer-Verlag. 16. S. B. YUSTE. Submittedfor publication. 17. A. M. SANCHEZ and J.D. BEJARANO 1986 Journal of Physics A: Mathematical and General 19, 887. 18. J. D. BEJARANO, A. M. SANCHEZ 85, 5128.
and C. M. RODRIGUEZ
1986 JournalofChemical
Physics
J.
466
GARCIA-MARGALLO
AND
J. D. BEJARANO
APPENDIX I: OSCILLATOR TYPE I-DERIVATION k = 1,2,3
J
237
B,(m)=i
OF THE FUNCTIONS B,(m), 4K
cos4
z
sin z dn u sin (p3dp, = -!-
lr 0 lr I 0 = (4/105~rm~){(9mm,(2m - l)+[lOmm,+4(2m -l)-[lOmm,+4(2m
-(9mm:(2m
cn4 u sn2 u dn* u du - 1)‘](4m - 1)) E(m)
-1)'](2m,-3mm,))K(m)},
271
B,(m)=l
cos'
z
sin z dn u sin p3 dq, = -
cn* u sn2 u dn* u du
r I0 =(4/15~m2)[2(m2+m,)E(m)+m,(m-2)K(m)], sn2 u dn2 u du
B,(m)=:
= (4/3rrm)[(2m
-l)E(m)+
m,K(m)].
Here K(m) and E(m) are the complete elliptic integrals of the first and second kind respectively, and ml = k’* = 1 - m is the complementary parameter. In our treatment in terms of elliptic functions, the argument of sin rp, and cos p3 in the Fourier expansion is, for example, the usual amplitude function cp3 = am (of; m), so that cos 'p3 = cn (wt; m) = cn u and sin cp3 = sn (of; m) = sn u; i.e., cp3 = I dn u du. Observe that in the first order approximation of this paper our method is similar to the principle of energy balance as used in reference [ 1l] and to all orders of approximation to the “cubication” of reference [ 161. The foregoing papers being limited only to the case z, = 0, they need only B2 and s3 with different normalizations and/or notations and factors. The values of these integrals if B = 0, i.e. for m = 0, are given by Minorsky [8], and with our notation and normalizations B, = 5/8, B2 = l/4 and B3 = 1. If A = 0, B # 0 it is well known that for all energies E one has m = l/2 and consequently at all energies of the limit cycles the integrals are also constant: for oscillator type I, B, = 0.0550, B, = 0.2157, B3 = 0.7869. Observe that B3 is also the Q.M. “action integral” 4 p dq for the ASO. In consequence this integral is also found, except for normalization factors, in papers on the Q.M. AS0 oscillators [ 171. B2 is the corresponding integral for the AA0 [ 181. APPENDIX II: OSCILLATOR II-DERIVATION k = 1,2,3
OF THE FUNCTIONS a,(m),
dn4 u sn2 u cn2 u du
-[2m,(3-3m+2m2)K(m)]}, 4K
277
932(m)=L
r
J 0
dn’usnucnusinQ2dyz=~
J
dn’ u sn2 u cn2 u du
0
sn2 u cn2 u du =(4/3nm)[(2-m)E(m)-2m,K(m)].
Fourier expansion, in our treatment in terms of sin Q and cos Q is now calculated with the Q appropriate for this case, i.e., (p2 = k j cn u du, which gives cos cp2 = dn u and sin (p2 = k sn u. Except for the factor m, B2 is the same as for oscillator I.
The
ofSound nnd
Vibration (1990) W(3),
GENERALIZED
FOURIER
GENERALIZED
453-466
SERIES AND
VAN
DER
J. GARCIA-MARGALLO Departamento
LIMIT
OF
POL OSCILLATORS
AND
J.
D.
BEJARANO
de Fi’sica, Facultad de Ciencias, Universidad de Extremadura,
(Received
CYCLES
5 December 1988, and in revised form 2
06071 Badajoz, Spain
June 1989)
A very simple generalization, using Jacobian elliptic functions, of the usual Fourier series, appropriate for non-linear systems, is used to study the first order approximate solution of the generalized van der Pol oscillators X +AX +2BX3 + F( z3 + z,X’ + z,X4)A = 0. A generalized harmonic balance method is used to determine the limit cycles. The cases A > 0, B > 0 and A < 0, B > 0 are considered in some detail. For given values of the parameters z,, the values of A and B for which limit cycles exist are found as functions of m. Numerical values for the radius, frequency, and energy of the limit cycles are given. The presence of zero, one, two, three or four limit cycles depends on the value of the parameters of the equation. Stability, bifurcations of fixed points, and the limit cycle stability are studied qualitatively. 1. INTRODUCTION From elementary physics to quantum field theory, the fundamental functions most used to approximate experimental facts (e.g., the “fundamental oscillators”) are solutions of the equation *+AX=O.
(1.1)
In applied physics most approximation methods, such as Fourier series, depend heavily on the oscillatory solutions and on the Laplace and Fourier transforms of the exponential solutions. According to Legendre [l], Euler [2] first proposed, more than 200 years ago, the equation that generalizes equation (l.l), namely %+AX+2BX3=0.
(1.2)
To our knowledge, our group was the first [3] to put forward the exponential solutions of equation (1.2) for all signs and values of A, B, E, and all X regions of possible interest in classical or quantum mechanics, where E is given by the more convenient “energy equation” *2+AX2+
BX~=
E,
x =
dx/dt.
(1.3)
On the basis of these solutions we have been able to develop the corresponding Fourier series and transforms [4]. We have used [5-71 these generalized Fourier series with the harmonic balance method to find an approximate solution of equations of the type X+x’=E(l-X2)X. The purpose of this paper is to determine the first order solutions by the method of harmonic balance using the generalized Fourier series and the Jacobian elliptic functions 453 0022-460X/90/030453
+ 14 $03.00/O
0
1990 Academic PressLimited
454
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
in equations of the type X +g(x) To avoid mathematical case of g(X)=AX+2BX3 der Pol equation,
= &j-(X,X).
(1.4)
difficulties as much as possible, we limit our study to the simple andf(X, X) =(z3+z2X2+z,X4)X: i.e., to a generalized van
X+AX+2BX3+e(z3+z2X2+z,X4)X=0.
(1.5)
Equations (1.4) or (1.5) with E =0 are called generating equations, and their solutions generating solutions. With B = 0 the approximate solutions for the many different cases are given and discussed in the well known textbook by Minorsky [8]. When z, = 0, solutions of these equations are also well known, and we compare our results and limit cycles with those of other authors [9-111 for this very special case and with numerical integration results for other quite complicated and interesting cases. Due to the X4X term there are zero, one, two, three or four limit cycles depending on the value of the parameters. Following references [lo-121, we have studied the stability of the limit cycles and the bifurcations. 2. THE GENERALIZED FOURIER
SERIES
METHOD
To avoid computational difficulties we limit our discussion to a first order approximation. Due to the fact that the usual Fourier series for the elliptic functions contain higher harmonics, this very simple approximation has proved to be quite good in the cases we have studied. In the potentials we are considering, V(X) = AX’+ BX4 for E < V,,, = A2/4B. We assume a solution of the form X(t) = a pq (wt; k2),
(2.1)
where a and w are constants to be determined, and pq denotes a convenient Jacobian elliptic function determined by the generating solutions of equations of the type (1.4). Substituting equation (2.1) into equation (1.4) gives F,(a, w, k2, E, a) cos cp+ F2(a, w, k2, E, a) sin cp+ (higher order harmonics) = 0,
(2.2) where cp is the generalized circular function of the case being considered [4], and (Y collectively denotes any parameter which appears in the non-linear function f( X, X). We first take F, = 0 and F2 = 0, using the method of harmonic balance [ 131. We shall study the different cases of the quartic anharmonic asymmetric oscillator (AS0 and AAO) with the form of the potential given above: i.e., with E = 0 for X = 0. We shall distinguish two cases of the quartic oscillator: (I) two types of potential: A > 0 and B>O, and A-CO, B>O, E>O; (II) A<0 and B>O, and E
OF TWO TYPES
OF QUARTIC
OSCILLATOR
3.1. OSCILLATOR TYPE 1 For this type of oscillator, following references [6,7], we take a function X = a cn (wt; k2),
(3.1)
GENERALIZED
VAN
DER
POL
where a, o and k* = m are constants determined to first order approximation. expression (3.1) into equation (1.5) gives (aA+2amw2-uw2)
455
OSCILLATORS
cn+(2a3B-2umw*)
Substituting
cn3- EUWZ~ sn dn
- cu3wz2 cn* sn dn - &u5wz, cn4 sn dn = 0. Our generalized [14, 151, gives
Fourier expansion,
[u(A+2mw’-w2)+l~5u(Bu2-mw2)] - [ ~uoz,B,(
if the expansion
(3.2)
is limited to the first harmonic
cos cp3
m) + ~u~wz,B,( m) + EU~CIJZ, B,( m)] sin 40~+ (higher harmonics) = 0,
(3.3) where B,(m), B,(m) and B3( m) are given in Appendix I. Setting the coefficient of sin (p3 to zero, one obtains -[.ww(u4z,B,(m)+u2z2B2(m)+z3B,(m))]=0,
and, for u2 = y, y2z,B,(m)+yz,Bz(m)+z3B3(m)=0.
(3.4)
The coefficient of cos rp, can be most simply made zero,
w2= Bu2/m
and
w’=A/(l-2m),
(3.5)
and from these two equations
A/B= (1-2m)u2/m. 3.2.
OSCILLATOR
TYPE
(3.6)
11
For this type of oscillator the solutions corresponding [6,7] are
to a first order approximation
X = a dn (wt; m).
(3.7)
Substituting expression (3.7) into equation (1.5) gives (uA+2uw2-uw2m)dn+(2u3B-2uw2)dn3-Euwmz,sncn - .a3wmz2 dn* sn cn - .w5wmz, dn4 sn dn = 0.
Our generalized Fourier expansion, if the expansion is limited to the first harmonic, gives [u(A+~w~-~w~)+~~~u(Bu~-~~)]
cos (p2
-[[~u~m(z~~~(m)+u*z~93~(rn)+u~z,9?,(m))]
sin cp,+(higher harmonics) =O, (3.8)
where a,(m), L%,(m) and a,(m) sin (p2to zero gives
are given in Appendix
II. Setting the coefficient oc
-[~u~m(z3%13(m)+u2z2~2(m)+u4z,93~(m))]=0,
and, for u* = y, y2z,CB,(m)+yz2B~(m)+z3W3(m)=0.
(3.9)
Setting the coefficient of cos rp2to zero gives w2 = Bu*
and
w’=A/(m-2),
(3.10)
and from these two equations A/B = (m -2)~~.
(3.11)
456
.I. GARCIA-MARGALLO
AND
4. ANALYTICAL
J. D.
BEJARANO
RESULTS
We assume that E -0 always. Consider the parameters z, , z2, and z3 as co-ordinates axes of the three-dimensional space in which we shall study the solution of equations (3.4) and (3.9) in the eight octants, and the following six special cases: (1) z, = z2 = 0; (2) z, = z3 = 0; (3) z2 = z3 = 0; (4) z, = 0; (5) z2 = 0; (6) z3 = 0. As before, we study all cases in first order approximation. 4.1. THE SIX SPECIAL CASES Case (l), z, =z2=0. In this case equation (1.5) becomes X+AX+2BX3+~(z,)X=0.
(4.1)
C.1.1.: oscillator type I. (a) For A > 0,B > 0,with 0 G m S l/2. There is one fixed point (f.p.) at X = 0, X = 0: i.e., (0,O). Let us study its stability. If z3 > 0, the positive coefficient of X means physically dissipation of energy, so that the f.p. is asymptotically stable (S), a sink. The topological configuration (t.c.) is S. If z,
this case is well known in the literature on the subject
X+AX+2BX3+e(z3+z2X2)X=0.
(4.2)
C.4.1.: oscillator type I. (a) One f.p. at (0,O). Equation (3.4) that gives the stationary state is yt2B2 + z3B3 = 0, where y= a 2 =-z3B3/z2B2.
(4.3)
The limit cycles (l.c.) can exist only if sign (z2) # sign (z3). Due to the fact that B3/B2 is a positive function of m only, a2 depends only on the scale factor -~3/t2 and on this function. In other words it is only necessary to calculate as we do in Figure 1 for one value of z3/zZ to obtain the general behaviour. For given z2 and z3, from equations (3.6) and (4.3) we obtain a2 and A/B us.m. In Figure 1, A/B is shown us. m (for z2 = 1.0, z3 = -3.O), giving the l.c. of equation (4.2). Observe that this function has extremes: i.e., points in a region of (A/B) with two l.c. that coalesce at the extreme of the function. Let us now study the stability and the bifurcations, that is, the t.c. for different values of the parameters. If z2> 0, and z3 < 0, the f.p. is unstable, U. From Figure 1, there exists one limit cycle for each A/B which is stable. The t.c. is US. If z2 < 0, and z3> 0, the point
GENERALIZED
VAN
DER
457
POL OSCILLATORS
2374-!----l vg5F---=-T 18 -
18-
17 -
1716-
1615 -
15-
I
14 -
14-
13 -
13-
A>0
12 -
12-
B>O
11 -
11 lo-
lo9-
9’0
8-
? q
9-
A>0
8-
7-
7-
6-
6? q
= 4-
54-
Cfl
3-
32-
\
1
2l-
l0 -1 -2-3 -4-5-6 -7 -8-
_;:L
O-1 0,2 0.3 0.40.5 060.7 0.8 0.9 1.0 m
Figure 1. The ratio A/B versus m, case (4),for the solution cn, dn. z, = 0, q = 1 and z3 = -3.
-,0-
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 nl Figure 2. The ratio A/B versus m, case (6), for the solution cm, dn. .zI = 1, z2 = -3 and z, = 0.
of rest is stable, S, and the l.c. is U. The t.c. is SU. Equation (4.2) undergoes a bifurcation of equilibria as .zJ passes through zero. (b) Equation (4.2) has three f.p. at (0,0), (*w, 0); if z2 > 0, and z3 < 0, it has one saddle point and two sinks; if z2< 0, and z3> 0, it has one saddle point and two sources. From Figure 1, the l.c. have four regions in the oscillator: (1st) for each value of A/B one l.c. only exists in the cnodial oscillator; (2nd) for each A/B two l.c. exist, one in the cnodial and another in the dnoidal oscillator; (3rd) for each value of A/B there are two l.c. in the cnoidal oscillator with different m; (4th) for the minimum value of A/B there is one limit cycle in the cnodial oscillator, from coalescence of the two in the (3rd) region. The stability of the l.c. is as follows: (1st) if z2 > 0, z3 < 0, S; if z2 < 0, z3 > 0, U; (2nd) as for the 1st region; (3rd.a) if z2 > 0, z3 < 0, S; if z2 < 0, z3 > 0, U; (3rd.b) if z2> 0, z3 ~0, U; if z2 CO, z3 > 0, S; (4th) the l.c. are semi-stable, US. Equation (4.2)
458
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
undergoes a bifurcation of equilibria as z3 passes through zero. When A passes through zero, one saddle-node coalescence occurs. The frequency (w), from equation (3.5), and the total energy (E) [3] of the l.c. are given by 02= Ba2/m and E = rn,a2w2. C.4.11.: oscillator type II. (a) Equation (4.2) has three f.p., at (0, 0), (+dm, 0). Equation (3.9) for the l.c. becomes yz2g1+ %13z,= 0, so that y=a2=-z,~3/z2532.
(4.4)
The stationary amplitude can exist only if sign (z2) # sign ( z3). The equation is the same as for oscillator I but now the %33k functions are those of oscillator II given in Appendix II. For given z2 and z3, from equations (3.11) and (4.4) we obtain A/B us. m. The l.c. in this case is given in Figure 1 only for certain values of A/B. For these values l.c. also exist in the cnoidal oscillator (C.I.b.l), with stabilities. If z2 > 0, and z3 < 0, one saddle point, two sinks, and one dnoidal l.c. that is unstable. The dnoidal t.c. is SU. If z2 < 0, and z3 > 0, one saddle point, two sources, and one dnoidal l.c. that is stable. The dnoidal t.c. is US. For this case there also exists a bifurcation of equilibria as z3 passes through zero. If m + 1 in (C.4.1.b) and (C.4.11.a) one homoclinic cycle exists. The frequency w, from equation (3.10), and E are given by o2 = Ba* and E = -m,a202. Case (5), z2 = 0. Equation (3.4) reduces to y2z, B, + z,B, = 0, so that the radius of the l.c. is now given by a4= -z,B,/z, B,. Also here B, and B2 are always positive, so that the stationary amplitude can exist only if sign (z,) # sign (z~). The t.c. are the same as in case (4). Case (6), z3 = 0. Equation (1.5) becomes X+AX+2BX2+s(z2X2+z,X4)X=0.
(4.5)
C.6.1.: oscillator type I. (a) There is one f.p. in (0,O). Equation (3.4) gives y’z,B, + yz2B2 = 0, or y(yz, B, + z2B2) = 0, so that y = 0 and y=a
2
=-z2B2/z,B,.
(4.6)
The stationary
amplitude can exist only if sign (zi) # sign ( z2). We shall explain the y = 0 root later when studying the second octant of the general case. For given z, and z2, from equations (3.6) and (4.6) we obtain a2 and A/B us. m. In Figure 2, A/B is shown us. m: i.e. the l.c. of equation (4.5). Observe that these functions have no extremes. If zI > 0, and z2 < 0, the point (0,O) is U, and the l.c. is S. The cnoidal t.c. is US. If z1 CO, and z2> 0, the cnodial t.c. is SU. Equation (4.5) undergoes a bifurcation of equilibria as z2 passes through zero. (b) Equation (4.5) has three f.p. at (0,O) and (fw, 0). From Figure 2 we have one l.c. region in the cnodial oscillator. In this region there exists only one l.c. for each A/B. If z, > 0, z, < 0, one saddle and two sources. In the cnoidal oscillator, one l.c. that is stable. The cnoidal t.c. are US. If z, < 0, z2 > 0, one saddle and two sinks. In the cnoidal oscillator, one l.c. that is unstable. The cnoidal t.c. are SU. When A passes through zero, one saddle-node coalescence occurs. The o and E are given by w2 = Ba’/m and E = m 1a2w2 .
C.6.11.: oscillator type II. (a) From Figure 2, limit cycles exist only of A/B. If z, > 0, and z2 < 0, one saddle, two sources, and one l.c. that t.c. are US. If z, < 0, and z2 > 0, one saddle, two sinks, and one l.c. that t.c. are SU. For z2 passing through zero there is one bifurcation. If m +
for certain Values is S. The dnoidal is U. The dnoidal 1 in (C.6.1.b) and
GENERALIZED
VAN
DER
POL
OSCILLATORS
459
(C.6.II.a), the point (0,O) is a saddle, the l.c. coalesce and one homoclinic cycle appears. The w and E are, respectively, o*= Ba* and E = -m,a*w*. 4.2. GENERAL CASES Octant (1): z,>O, z,>O, z,>O. Equation (1.5) becomes ~+AX+2BX3+~(z3+z,X2+z,X4)~=0.
(4.7)
Oc.l.1.: oscillator type I. (a) Equation (3.4) becomes y’z,B, +yz2B2+ z,B, = 0, and has no positive roots, the only stationary state is the f.p. which is S. (b) The f.p. are one saddle and two sinks. Oc.l.11.: oscillator type II. (a) This is analogous to (0c.l.I.b). Octant (2): z, > 0, z2 < 0 and z3 > 0. Oc.2.1.: oscillator type I. (a) The f.p. is stable, S. Equation (3.4) becomes y*z, B, + yz,B, t z,B, = 0, and this polynomial has two roots, (4.8) a2 = *Y = (-z2B2+sqfl ((z2B2)*-44(zIB1)(z3B3)))I(2zIB1), where there are three cases, as follows. (1) If (-z2B2)‘>4(z,B1)(z3B3). For given zl, z2, and z3, from equations (3.6) and (4.8), we obtain A/B vs. m. The limit cycles for this case are given in Figure 3. Hence the two roots are positive. Of these roots the smaller one, a*-, corresponds to an unstable, and the larger, a*+, to a stable limit cycle. The cnoidal t.c. is SUS, bicyclic. If z3 = 0 the smaller l.c. has radius 0.0. It then coincides with the fixed points and we have only the larger l.c. with radius y = -z2B2/zI B, that coincides with special case (6). (2) If (-z,B,)’ < 4(z, B,)(z,B,), the polynomial has no positive roots, and the only stationary state is the point of rest, S; the cnoidal t.c. is acyclic S. (3) If (-z,B,)* = 4(z, B,)(z,B,) there appears a surface, 5, separating these configurations. The t.c. passes from SUS (or bicyclic) to S (acyclic). The surface 4 is thus the locus of the bifurcation points. Topologically this is also a locus of semi-stable cycles US. (b) Here the f.p. are one saddle and two sinks. Equation (3.4) has also three cases, as follows. If (-z2B2)* > 4(z,B,)(z,B,). From Figure 3, we have five l.c. regions: (1st) one small, U, one large, S, root of the cnoidal oscillator for each A/B; (2nd) one small, U, one large, S, root of the cnoidal oscillator and one small root of the dnoidal oscillator for each A/B; (3rd) one large root, S, two small roots, U, but different m, of the cnoidal oscillator for each A/B; (4th) one large root, S, and coalescence of two l.c. of small roots, SU, of the cnoidal oscillator for each A/B; (5th) one large root only, S, of the dnoidal oscillator for each A/B. The other cases are analogous to (Oc.4.I.a). We obtain bifurcations when A/B passes from one to another of these five regions. When A passes through zero a saddle-node coalescence occurs. The w and E are o*= B(a**)/m and E = m,(a**)w*. Oc.2.11: oscillator type II. (a) The f.p. are one saddle, and two sources, U. Equation (3.9) becomes y*z,%, +yz2%z+ z3B3 = 0. This polynomial has two roots:
a2 = *.v =
(-z222*sW
((~~~2)*-4(2,~~)(~~~~)))/(2~,98,).
(4.9)
Consequently, we have the same three cases as before: (1) If ( -z,L%~)*> 4(~,L%,)(z#~), algebraically there are two roots and both are positive. For given zl, z2 and z3, from equations (3.11) and (4.9) we obtain A/B vs. m. The l.c. for this case is given in Figure 3 and has two regions: (1st) the (2nd) case described previously; (2nd) one large root, S, in the dnoidal oscillator for each A/B.
460
J. GARCIA-MARGALLO
AND
J. D.
BEJARANO
”
”
”
1’
‘-
+
L.
1017 16 15-
A'0 8~0
1413 12 ll10 98769 q
5432l0 -1
-
-2
-
-3
-
-4
-
A<0 B>O
-5 -6
-6 A<0 B>O
-7
-7 -6
-6
_;:L
04 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1 .O m
Figure 3. The ratio A/B versus m, Ott.(2), for the solution cn, dn. 2, = 1.2, z2 = -3.11 and z3 = 1.
(2) If (-z2W,)* <4(z,.93,)(z3~~),
-9 -10
dn A 0.4 0.2 0.3 04 0.5 0.6 0.7 0.0 O-9 I.0 In
Figure 4. The ratio A/B versus m, Ott.(3), solution cn, dn. z, = 1, z2 = -4 and zx = -4.
for the
the polynomial has no real roots, the only stationary state is (km, 0); the dnoidal t.c. is acyclic S. there appears a surface, 0, separating these configura(3) If (-z,93J2 = 4(z,9,)(zj98J tions. Algebraically, this is a locus of double roots: u* = -z29B2/2z,B,. Topologically this is also a locus or semi-stable US. The surface, 0, separates the two possible configurations SUS and S. Equation (4.7) undergoes a bifurcation in the surface, 0, as (z&%)* passes through 4( ~,a,)( z3B3). In (Oct.2.I.b) and (Oct.2.II.a) when m + 1, the limit cycles coalesce for the small and large roots in the saddle point (0, O), in two homoclinic cycles. The w and E are given by W*= B(a’*) and E = -m,(uz*)~*. Octant (3): z,>O, z2C0, z,
GENERALIZED
VAN
DER
POL
OSCILLATORS
461
(a) The f.p. is unstable, U. Equation (3.4) is y’z, B, + yzZB2+ z3B3 = 0. This polynomial has only one positive root: a2= ?? Y = (-zJ&+sqrt
((z~B~)~+~(z,B~)(z~B~)))/(~z,B~).
(4.10)
For given z, , z2, and z3, from equations (3.6) and (4.10), we obtain A/B us. m. Figure 4 gives the limit cycle. The cnoidal limit cycle is stable. The t.c. is US. As was the case for the second octant, when z3+0 the l.c. radius has the minimum value for this octant of u2 = -z2B2/z,B,, analogous to the special case (6). (b) The fixed points are one saddle and two sources. From Figure 4, we have one l.c. for each A/B, which is S. When A passes through zero, a saddle-node coalescence occurs. The w and E are given by w2 = Ba’/m and E = m,a2w2. Oc.3.11.: oscillator type II. (a) The three f.p. are one saddle and two sources. From Figure 4 we have one l.c. for each A/B, which is S. The A/B (dnoidal) < A/B (cnoidal);
““V-----l Y.
-7%
1917 16 15
-
14 1312
-
A'0
11 -
B>O
10 9676-
m G
54321 0 -1
-
-2
-
-3
-
-4
-
-5 -6
-
-7
-
-6
-
-9
-
-lO-
Figure
5. The ratio A/B
dn
A< 0 BP0
’ ’ ’ ’ ’ ’ ’ ’ 0.1 0.2 0.3 0.40.5 0.6 0.7 0.0 09 1.0 m
versus m, Ott.(4),
for the solution
cn, dn. z, = 1, z2 = 1 and I) = -3.
462
J. GARCIA-MARGALLO
AND
J. D.
BEJARANO
they are equal only for m = 1. For m = 1, in (Oct.3.I.b) and (Oct.3.II.a), there exists one homoclinic cycle. The w and E are given by w2= Ba2 and E = -m,a2w2. Octant (4): z, > 0, z2 > 0, z3 < 0. Oct.4.1.: oscillator type I. (a) This is analogous to (Oct.3.I.a). (b) the f.p. are one saddle and two sinks. From Figure 5, there are four regions: (1st) one l.c. for each A/B, which is S; (2nd) two l.c. with one in the cnoidal oscillator, S, and the other in the dnoidal oscillator, which is US, for each A/B; (3rd) two l.c. in the cnoidal oscillator for each A/B, but with different m, both S; (4th) coalescence of the two previous l.c. for the minimum value of A/B. A saddle-node coalescence occurs when A passes through zero. The w and E are given by w2 = Ba2/ m and E = m,a202. Oct.4.11: oscillator type II. (a) The f.p. is analogous to (Oct.4.I.b) and the l.c. is analogous to (Oct.4.I.b.11). Homoclinic cycles exist when m + 1 in (Oct.4.I.b) and (Oct.4.II.a.). The w and E are given by w2 = Ba2 and E = -m1a202. For the octants of the lower half-space (zi
5. COMPARISON WITH NUMERICAL INTEGRATION It is instructive to compare these simple analytical approximations with a numerical integration of the equations. We used the fourth order Runge-Kutta method and illustrate the more important results with the value E = O-1. The results for the limit cycles are shown in Figure 6 for the special case (C.4.1.a) of the oscillator type I, when A/B = 23.125, z1 = 0, z2 = 1 and z3 = -3, and from Figure 1 one has rn = O-25, and from equation (4.3), a = 3.40. The (a) curve is the analytical
:5
Figure 6. The (a) curve is the analytical solution X = 3.40 cn (6.80r; m = 0.25). The (b) curve is the numerical integration with initial conditions X( t = 0) = 3.40 and A( t = 0) = 0.00.
GENERALIZED
VAN
DER
POL
OSCILLATORS
463
solution X = 3.40 cn (6.80t; m =0.25). The (b) curve represents the results of the numerical integration with initial conditions X( t = 0) = 3.40 and X( 1= 0) = 0.00. In Figure 7 are shown (C.4.1.b) and (C.4.11.a) for A/B = -6.70, z, = 0, z2 = 1, zj = -3; from Figure 1 for the cnoidal oscillator, M = 0,755, and from equation (4.3) a = 3.15. For the dnoidal oscillator M = 0.905, and from equation (4.4) a = 2.47. The (a,) curve is the analytical solution X = 3.15 cn (3.62t; WI= O-755). The (az) curve is the analytical solution X = 2.47 dn (2.47t; m = O-905). The (b,) curve represents the results of the numerical integration with initial conditions X(t = 0) = 3.15 and X(t =0) =O*OO. The (b,) curve represents the results of the numerical integration with initial conditions X( t = 0) = 2.47 and X(t=O)=O*OO.
Figure 7. The (a,) curve is the analytical solution X =3.15 cn (3.621; m =0.755). The (a,) curve is the analytical solution X = 2.47 dn (2.47r; m = 0.905). The (b,) curve is the numerical integration with X( t = 0) = 3.15 and X(? = 0) = 0.00. The (b,) curve is the numerical integration with X( t = 0) = 2.47 and X( I = 0) = 0.00.
In Figure 8 is shown the case (Oc.2.I.a.l) for A/B = 5.05, and z, = 1.2, z2 = -3.11, z3 = 1; from Figure 3, m = 0.23 for a-, and m = 0.27 for a+. The (a,) curve is the analytical solution X = 1.47 cn (3*06t; m =0*23). The (b,) curve represents the results of the numerical integration with initial conditions X( t = 0) = 1.47, X( t = 0) = 0.00. The (a,) curve is the analytical solution X = 1.72cn (3.31 t; m = 0.27). The (b2) curve, the numerical integration with initial conditions X( t = 0) = 1.72, X( t = 0) = 0.00. In Figure 9 are shown the cases (Oc.2.I.b) and (Oc.Z.II.a), for A/B = -0.952, and z, = 1.2, z2 = -3.11, z3 = 1; from Figure 3, m = 0.59 for a+, and m = 0.70 for a-, in the cnoidal oscillator, and m = O-965 for a- in the dnoidal oscillator. The (a,) curve is the analytical solution X = 1.76 cn (2.30t; m = 0.59). The (az) curve is x= 1.29 cn (1*54?; m = 0.70). The (aj) curve is X = 0.95 dn (0.95t; m = O-965). The (b,) curve is the numerical integration with initial conditions X( t = 0) = 1.76, X(t = 0) = 0.00, the (b2) curve is that with X( t = 0) = 1.29, X( t = 0) = 0.00, and the (b,) curve is that with X( t = 0) = 0.95, X( t = 0) = 0.00. All the analytical results seem to be very good for E s 0.1.
464
J. GARCIA-MARGALLO
AND
J. D. BEJARANO
Figure 8. The (a,) curve is the analytical solution X = 1.47 cn (3.061; m =0.23). The (b,) curve is the numerical integration with X( t = 0) = 1.47, X( f = 0) = 0.00. The (a2) curve is the analytjcal solution X = 1.72 cn (3.3 1 t; m = 0.27). The (b,) curve is the numerical integration with X( r = 0) = 1.72, X( r = 0) = 0.00.
Figure 9. The (a,) curve is the analytical solution X = 1.76 cn (2.30r; m = 0.59). The (a,) curve is X = 1.29 cn (1.54r; m =0.70). The (a,) curve is X = 0.95 dn (0.95r; m = 0.96). The (b,) curve is the numerical integration with X(r=O)= 1.76, k(r=O)=O.OO. The (b2) curve is with X(r=O)= 1.29, X(r=O)=O.OO. The (b,) curve is with X( r = 0) = 0.95, X( r = 0) = 0.00.
GENERALIZED
VAN
DER POL OSCILLATORS
465
6. CONCLUSIONS We have described how to generalize the method of harmonic balance to obtain the first order approximations to the periodic solutions of differential equations of the type (1.4), using elliptic functions, and we have also studied the stability of the fixed points, and of the limit cycles and bifurcations, from a qualitative viewpoint. We have applied the method to the study of the limit cycles of equation (1.5). The presence of zero, one, two, three or four limit cycles depends on the value of the parameters z,, as shown in Figures l-5. We have compared the results with numerical integration results in Figures 6-9 and this comparison shows the analytic approximations to be very good.
ACKNOWLEDGMENTS Thanks are due to the work Dr R. Chatwin has contributed to improvements of this and previous papers [6,7] on non-linear mechanics. We gratefully acknowledge the Comisi6n Asesora de Investigacibn Cientifica y Ttcnica (CAYCIT), Spain, project no. PB87-0007, for financial support and J.D.B. would like to acknowledge the Plan Movilizador, and the E.P. Division at CERN for its hospitality.
REFERENCES 1825 Traits’des functions elliptiques et des intkgrales eule’riennes. Paris: 1. A. M. LEGENDRE Huzard-Courcie. 2. L. EULER 1763 Nov. Corn. Petropotan X pg 4. Id toms VI, VII, 1761. 3. J. D. BEJARANO and A. M. SANCHEZ 1988 Journal of Mathematical Physics 29, 1847-1853. Generalized exponential, circular and hyperbolic functions for non-linear quantum mechanics. 4. J. D. BEJARANO and A. M. SANCHEZ 1988 Journal of Mathematical Physics 30, 1871-1876. Generalised Fourier transforms for nonlinear systems. 5. S. B. YUSTE and J.D. BEJARANO 1986Journal of Sound and Vibration 110,347-350. Construction of approximate analytical solutions to a new class of non-linear equations. and J. D. BEJARANO 1987JournalofSoundand Vibration 116,591-595. 6. J. GARCIA-MARGALLO A generalization of the method of harmonic balance. 7. J. GARCIA-MARGALLO,J.D.BEJARANO~~~S.B.YUSTE~~~~ Journalofsoundand Vibration 125, 13-21. Generalized Fourier series for the study of limit cycles. 8. N. MINORSKY 1962 Nonlinear Oscillations. Malabar: Litton Educational Publishing, Inc., U.S.A. 9. P. HOLMES and D. RAND 1980 Journal of Non-Linear Mechanics 15,449-458. Phase portraits and bifurcations of the non-linear oscillator: jl + (a + yx*),? + /3x + Sx3 = 0. 10. J. GUCKENHEIMER and P. HOLMES 1986 Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. New York: Springer-Verlag. 11. E. KNOBLOCH and M. R. E. PROCTOR 1981 JournalofFluid Mechanics 108,291-316. Nonlinear periodic convection in double-diffusive systems. 12. D. W. JORDAN and P. SMITH 1987 Nonlinear Ordinary Differential Equations. New York: Oxford University Press. 13. R. MICKENS 1987Journal of Sound and Vibration 111,515-518. A generalization of the method of harmonic balance. 14. M. ABRAMOWITZ and J. A. STEGUN (editors) 1972 Handbook of Mathematical Functions. New York: Dover. 15. P. D. BYRD and M. D. FRIEDMAN 1971 Handbook of Elliptic Integrals for Engineers and Scientists. Berlin: Springer-Verlag. 16. S. B. YUSTE. Submittedfor publication. 17. A. M. SANCHEZ and J.D. BEJARANO 1986 Journal of Physics A: Mathematical and General 19, 887. 18. J. D. BEJARANO, A. M. SANCHEZ 85, 5128.
and C. M. RODRIGUEZ
1986 JournalofChemical
Physics
J.
466
GARCIA-MARGALLO
AND
J. D. BEJARANO
APPENDIX I: OSCILLATOR TYPE I-DERIVATION k = 1,2,3
J
237
B,(m)=i
OF THE FUNCTIONS B,(m), 4K
cos4
z
sin z dn u sin (p3dp, = -!-
lr 0 lr I 0 = (4/105~rm~){(9mm,(2m - l)+[lOmm,+4(2m -l)-[lOmm,+4(2m
-(9mm:(2m
cn4 u sn2 u dn* u du - 1)‘](4m - 1)) E(m)
-1)'](2m,-3mm,))K(m)},
271
B,(m)=l
cos'
z
sin z dn u sin p3 dq, = -
cn* u sn2 u dn* u du
r I0 =(4/15~m2)[2(m2+m,)E(m)+m,(m-2)K(m)], sn2 u dn2 u du
B,(m)=:
= (4/3rrm)[(2m
-l)E(m)+
m,K(m)].
Here K(m) and E(m) are the complete elliptic integrals of the first and second kind respectively, and ml = k’* = 1 - m is the complementary parameter. In our treatment in terms of elliptic functions, the argument of sin rp, and cos p3 in the Fourier expansion is, for example, the usual amplitude function cp3 = am (of; m), so that cos 'p3 = cn (wt; m) = cn u and sin cp3 = sn (of; m) = sn u; i.e., cp3 = I dn u du. Observe that in the first order approximation of this paper our method is similar to the principle of energy balance as used in reference [ 1l] and to all orders of approximation to the “cubication” of reference [ 161. The foregoing papers being limited only to the case z, = 0, they need only B2 and s3 with different normalizations and/or notations and factors. The values of these integrals if B = 0, i.e. for m = 0, are given by Minorsky [8], and with our notation and normalizations B, = 5/8, B2 = l/4 and B3 = 1. If A = 0, B # 0 it is well known that for all energies E one has m = l/2 and consequently at all energies of the limit cycles the integrals are also constant: for oscillator type I, B, = 0.0550, B, = 0.2157, B3 = 0.7869. Observe that B3 is also the Q.M. “action integral” 4 p dq for the ASO. In consequence this integral is also found, except for normalization factors, in papers on the Q.M. AS0 oscillators [ 171. B2 is the corresponding integral for the AA0 [ 181. APPENDIX II: OSCILLATOR II-DERIVATION k = 1,2,3
OF THE FUNCTIONS a,(m),
dn4 u sn2 u cn2 u du
-[2m,(3-3m+2m2)K(m)]}, 4K
277
932(m)=L
r
J 0
dn’usnucnusinQ2dyz=~
J
dn’ u sn2 u cn2 u du
0
sn2 u cn2 u du =(4/3nm)[(2-m)E(m)-2m,K(m)].
Fourier expansion, in our treatment in terms of sin Q and cos Q is now calculated with the Q appropriate for this case, i.e., (p2 = k j cn u du, which gives cos cp2 = dn u and sin (p2 = k sn u. Except for the factor m, B2 is the same as for oscillator I.
The