ON THE APPLICATION OF IHB TECHNIQUE TO THE ANALYSIS OF NON-LINEAR OSCILLATIONS AND BIFURCATIONS

Article No. sv981641

ON THE APPLICATION OF IHB TECHNIQUE TO

THE ANALYSIS OF NON-LINEAR OSCILLATIONS

AND BIFURCATIONS

S. S

Department of Civil Engineering, Eastern Mediterranean University, G. Magosa, Kibris, Mersin 10, Turkey

 K. H†

Department of System Design Engineering, University of Waterloo, Waterloo, Ontario, Canada N2L 3G1

(Received 5 September 1997, and in final form 17 February 1998)

The main objective of this paper is to introduce certain refinements and alternative formulations, which enhance the applicability and availability of the intrinsic harmonic balancing technique. This is achieved by considering certain illustrative examples concerning non-linear oscillations and dynamic bifurcation phenomena. Indeed, the bifurcation behaviour of a harmonically excited non-autonomous system is analyzed conveniently, with reference to the corresponding autonomous system, by applying the IHB technique, which yields the bifurcation equation as an integral part of the perturbation procedure. A symbolic computer language, namely MAPLE, facilitates the analysis as well as the verification of the ordered approximations to the solutions. The methodology lends itself to MAPLE readily, which in turn, enhances the applicability of the IHB technique. 7 1998 Academic Press

1. INTRODUCTION

Non-linear oscillations and bifurcation problems can be analyzed via a variety of methods [1–3] such as averaging techniques, multiple time scaling, harmonic balancing, etc. Averaging method, for example, yields a lowest order approximation conveniently, but higher order calculations become lengthy and complicated. The method of harmonic balancing is conceptually simple and it leads to algebraic equations only; however, the results may be inconsistent [3].

The intrinsic harmonic balancing technique was introduced earlier [4–6] in order to overcome certain observed inconsistencies in the application of the conventional harmonic balancing method. The method has been applied effectively to the analysis of non-linear vibrations and dynamic bifurcation problems systematically. In addition, the basic concepts and the methodology of the IHB technique is generalized and adopted for the analysis of non-linear forced oscillations [7, 8].

Nevertheless, it seems that there are a number of issues regarding the perturbation procedure that require further clarification and refinement. This is true for non-linear

† Also associated with Eastern Mediterranean University, Department of Civil Engineering, G. Magosa, Mersin 10, Turkey.

oscillations as well as bifurcation problems. In the former case, for example, the evaluation of perturbation equations for the non-linear system may be referred to the linearized system or directly to the origin of the system in which case the evaluations of the ordered perturbation equations can be performed rapidly and conveniently. In this paper, alternative formulations, certain clarifications, and variations of the methodolgy are discussed with the aid of illustrative examples. The bifurcation analysis of a non-autonomous system is performed through the IHB technique and alternative formulations are discussed. It is observed that the IHB technique, as applied to this bifurcation problem here, provides a very simple treatment compared to other methods. It is expected that the exposition presented in this paper will enhance the application of the IHB technique to a variety of specific problems in many fields.

2. AUTONOMOUS SYTEMS

Consider an autonomous system generally described by

x¨ + g(x, x˙,o) = 0, (1)

where g is a polynomial function of x, x˙ and o; dots on x indicate differentiation with

respect to time and g(x, x˙, 0) gives the corresponding linear system. The solution of

equation (1) is sought in a parametric form x = x(t,o) where o is a small parameter.

A series of perturbation equations can be obtained by introducing the assumed solution

x = x(t,o) back into equation (1), differentiating with respect to o successively, and

evaluating these equations ato = 0:

o0_{: x¨ + g = 0,} o1_{: x¨' + g} xx' + gx˙x˙' + go= 0, o2: x¨0 + [g xxx' + gxx˙x˙' + 2gxo]x' + gxx0 + [gx˙xx' + gx˙x˙x˙' + 2gx˙o]x˙' + gx˙x˙0 + goo= 0, (2) etc.,

where the primes and subscripts on g denote differentiation with respect too and related

variables, respectively, and all perturbation equations are evaluated ato = 0.

Further, the assumed solution x = x(t,o) can be represented by a Fourier series of the

form [4–6]

x(t,o) = s

M m = 0

[pm(o) cos (mv(o)t) + rm(o) sin (mv(o)t)], (3)

which is substituted back into the perturbation equations sequentially. At each step, balancing the harmonics, one obtains the derivatives of the amplitudes to construct the Taylor’s expansion of the amplitudes to a desired order as

pm(o) = p0m+ p'mo +21pm0o2+ O(o3),

rm(o) = r0m+ r'mo +21rm0o2+ O(o3),

and the solutions are envisaged as

x(t,o) = x0(t, 0) +ox'(t, 0) + . . . , (4)

where x0(t, 0) is the solution of the linearized equation.

this may require a further step in the perturbation process. This can be done by introducing appropriate scales such that the solution is expressed in the form of

x(t,o) = ox'(t, 0) +1

2!o

2_x0(t, 0) + · · · , ₍₅₎

where o = 0 identifies the origin.

Obviously, the non-linear system is now referred to the origin rather than to the solution

of the linearized system since x(t, 0) = 0 and one has pm(0) = rm(0) = 0 in the Fourier

series (3).

As an example, consider a system described by

dx2_/dt2_{+ x = a +}ox2_{, (6)}

subject to initial conditions x(0) = A, xt(0) = 0. Hereo is a small positive parameter, o = 0

giving the linearized equation. System (6) was solved before [4, 6, 9] on the basis of the linearized equation.

In order to obtain the solution in the form of equation (5), the perturbation procedure may be facilitated by introducing certain scaling as

a =ob, A =oB. (7)

Further, one may introduce the time scaling t = v(o)t, in which case equation (6) takes

the form

v2_x

tt+ x =ob + ox2, (8)

where the subscriptt indicates differentiation with respect to t. The periodic solution may

be expressed in a parametric form in terms of o,

x = x(t; o), v = v(o),

and the assumed solution will be in the form of equation (3) with

pm(0) = rm(0) = 0, [m,

since x(t, 0) = 0.

The solution is now 2p-periodic in t, with v(0) = vc= 1. A sequence of perturbation

equations can now be generated by differentiating equation (8) with respect to o and

evaluating the derivatives at o = 0. Thus, one obtains

x'tt+ x' = b (first order), (9)

4v'x'tt+ x0tt+ x0 = 0 (second order), (10)

6v0x'tt+ 6v'x0tt+ 6(v')2x'tt+ x1tt+ x1 = 6(x')2 (third order), (11)

etc.

Substituting the assumed solution (3) into equation (9), one obtains sM

m = 0

[(1 − m2_)p'

mcos (mt) + (1 − m2)r'msin (mt)],

and balancing the harmonics yields the derivatives p'0= b and p'm= r'm= 0, me 2, which

give

Substituting initial conditions x'(0, 0) = B and x't(0, 0) = 0 yields p'1= B − b and r'1= 0. Using the above information one has

x'(t, 0) = b + (B − b) cos t

and the first order solution can be written as

x(t, o) = ob + o(B − b) cos t + O(o2_).

Similarly, substituting the assumed solution and the first order solution x'(t, 0) into the

second order perturbation equation (10) and balancing all the harmonics, one obtainsv'

and all the coefficients except p01 and r01 as zero. Further, using the initial conditions

x0(0, 0) = 0 and x0t(0, 0) = 0 yields p01= r01 = 0, which results in x0(t, 0) = 0.

Keeping in mind the scaling introduced to a and A, for the first order approximation one should go to the third perturbation. Substituting the assumed solution, first and second order solutions (x'(t, 0) and x0(t, 0)) once more into the third order perturbation equation (11) and balancing all the harmonics, one obtains the non-zero coefficients:

p01= 6b2+ 3(B − b)2, p21= −(B − b)2, v0 = −2b,

and using the third derivative of amplitudes yields

x1(t, 0) = 6b2_{+ 3(B − b)}2+ p1

1 cos (t) + (B − b)2cos (2t) + r11 sin (t),

where p11 and r11 are obtained from the initial conditions x1(0, 0) = 0 and x1t (0, 0) = 0

as

p11= −[6b2+ 2(B − b)2] and r11 = 0, respectively.

The solution, which is expressed as

x(t, o) = ox'(t; 0) +_2!1 o2_x0(t; 0) +1 3!o

3_x1(t; 0) + O(o4_{), (12)}

now takes the form

x(t, o) = a + (A − a) cos t + o

$

0

1

%

after rescaling (i.e., returning to the original variables). Moreover, approximation to the

v(o) can be written as

v(o) = 1 − ao + O(o2_).

Note that equation (12) is the third order solution of the scaled system and after rescaling, the corresponding solution becomes equation (13).

system or to the origin. If scalings similar to equations (7) are introduced, the non-linear system will be referred to the origin. For example consider the system

dx2_/dt2_{+ x = a + cx}2_{, (14)}

subject to initial conditions: x(0) = A, xt(0) = 0. Here, c is not small and cannot be treated

as a perturbation parameter. A small unidentified parameter should be introduced in order to facilitate the perturbation procedure and similar scalings as in equations (7) are

necessary for the application of the IHB technique. Let the unidentified small parameter

bem, and introduce the scaling a = mb, A = mB, together with the time scaling t = v(m)t. Then equation (14) becomes

v2_x

tt+ x =mb + cx2. (15)

Upon writing equation (14) in the first order form x˙ = y, y˙ = a − x + cx2_{, equilibrium}

points can be obtained as y = 0 and the roots of cx2_{− x + a = 0. Periodic solutions may}

exist if the roots of cx2_{− x + a = 0 are real. In addition, attention here is focused on small} vibrations in the vicinity of an equilibrium point (centre), so that for an asymptotic solution A and a can be assumed to be small. Then, the procedure of the IHB technique, as described above, is followed and the periodic solution of system (14) is obtained as

x(t, m) = [b + (B − b) cos (t)]m + [cb2₊1 2c(B − b) 2_{− (cb}2₊1 3c(B − b) 2_{) cos (t)} − (1 6c(B − b)2) cos (2t)]m2+ [2c2b3+ c2b(B − b)2− c2b2(B − b) − 1 3c2(B − b)3 +1_6(−12c2_b3_{− 4c}2_{b(B − b)}2_{+ 4c}2_b2_{(B − b) +}29 24c2(B − b)3) cos (t) +1_6(−2c2_{b(B − b)}2_{+ 2c}2_b2_{(B − b) +}2 3c2(B − b)3) cos (2t) + 1 48c2(B − b)3cos (3t)]m3. (16)

Also, the amplitude–frequency relation can be constructed as

v(m) = 1 − cbm −1_2(3c2_b2₊5 6c

2_{(B − b)}2₎m2_{. (17)}

Solution (16) may also be obtained by referring the system to the linearized equation by introducing the additional scaling x:mz, so that system (15) may be written as

v2_z

tt+ z = b +mcz2, (18)

and the initial conditions become z(0) = B, zt(0) = 0. One obtains the solution of system

(18) by substituting the assumed solution (3) sequentially into the zero order, first order,

etc., perturbation equations, and an approximate solution z(t, m) consistent to a desired

order can be constructed. The solution of the original system can be written after back

scaling (returning to the original variables),mz:x, which leads to the same solution (16).

The solution obtained by referring the system to the linearized equation reduces the number of perturbations and it seems an advantage. However, the former solution procedure, where the perturbation equations are evaluated at the origin is more convenient since the solution is referred to the origin as mentioned before.

2.1.   

The solutions obtained above are given in an ordered form of approximations. To verify the solution, one direct way is to substitute this solution back into the original differential equation [10]. In general, if the Nth order approximate solution is substituted back into the original equation and it yields a result in the order O(eN + 1_{), i.e., x = · · · + O(e}N + 1_),

words, the solution is a consistent approximation. If substitution, which contains lower order terms than N + 1, the approximation is inconsistent.

The third order asymptotic solution (12) is verified by substituting the solution into equation (8) with the aid of MAPLE, which yields a result of O(o4_{). Similarly, the first order} solution (13), obtained after returning to the original variables (rescaling) is substituted into the original equation (6) and a result of O(o2_{) is obtained. The approximate solution}

(16) of the system (15) is verified by substituting the solution (16) and v(m) back into

equation (15) and a result of O(m4_{) is obtained, confirming the consistency of the}

approximation (16).

3. NON-AUTONOMOUS SYSTEMS

The IHB technique has been adapted for the analysis of non-linear forced oscillations [7, 8, 11]. Here, the technique will be applied to a bifurcation problem under external harmonic excitation. The analysis will be carried out without the aid of multiple time scaling [11], and a number of important aspects of the procedure will be pointed out. Consider a specific, harmonically excited (non-autonomous) bifurcation problem given by

x¨ +v2

cx −o(h − x2)x˙ = F sin (Vt), (19a)

where o q 0 is a small parameter. It is assumed that natural frequency vc and external

frequency V satisfy the non-resonance relationship l1vc+ l2V $ 0 where l1, l2 are any

positive or negative integers. This system exhibits bifurcations from a periodic solution to a two-frequency quasi-periodic solution. The system is weakly non-linear and the solutions can be obtained without introducing scaling, by referring the non-linear system to the linearized one. If one uses the scale F:oC, the solution will be referred to the origin. However, the application here will be carried out by referring the system to the linearized equation without introducing any scaling.

Obviously, F = 0 gives the corresponding autonomous system and x = 0 is the

equilibrium position. System (19a) with F = 0 exhibits Hopf bifurcation as h passes

through zero. In order to demonstrate this one can write the above second order system as the first order one

x˙1= x2, x˙2= −v2cx1+o(h − x21)x2+ F sin (Vt). (19b)

The state defined by the origin, x1= x2= 0, is now identified as an equilibrium state of

the F = 0 system. The Jacobian is evaluated at the origin, and its eigenvalues are given by A =

$

%

Clearly,h Q 0 gives complex conjugate eigenvalues with negative real part and when h

is greater than zero real part of eigenvalues become positive. At h = hc= 0 the pair of

complex conjugate eigenvalues becomes an imaginary pair and Hopf bifurcation occurs. The family of periodic or quasi-periodic solutions is expressed in the parametric form

where t1=Vt and t2=v(o)t are introduced so that equation (19a) becomes

V2_x

11+ 2Vvx12+v2x22+vc2x −o(h − x2)(vx2+Vx1) = F sint1. (21)

Subscripts 1 and 2 on x indicate differentiation with respect to t1 and t2, respectively.

Further, one may assume a solution in the form of a generalized Fourier series (with two frequencies), x(t1,t2;o) = s M m = 0 m1+ m2= m pm1,m2(o) cos (m1t1+ m2t2) + rm1,m2(o) sin (m1t1+ m2t2), (22)

where m1and m2may be chosen positive or negative; and M is an arbitrary positive integer.

It is noted that equation (22) reduces to the ordinary Fourier series in the case of m10 0

or m20 0. m10 0 describes periodic solutions of the associated autonomous system while

m20 0 denotes periodic solutions which are purely excited by the external force F sin (Vt). As a matter of fact, assumption (22) embraces equilibrium, and periodic solutions as well as quasi-periodic motions, thus enabling one to identify bifurcations from one solution to the other.

A series of perturbation equations is obtained by substituting equation (20) into equation (21), differentiating with respect too and evaluating at o = 0. It is noted that o = 0 is the corresponding linear system so that one obtains zeroth, first, second, etc., order perturbation equations as V2_x 11+ 2Vvcx12+v2cx22+vc2x = F sin (t1) (0th order), (23) V2_x' 11+ 2Vv'x12+ 2Vvcx'12+ 2vcv'x22+v2cx'22+vc2x' − (h − x2₎₍v cx2+Vx1) = 0 (1st order), (24) etc.,

where primes indicate differentiation with respect too evaluated at o = 0.

Substituting the series solution (22) into equation (23) and balancing the harmonics yields the non-zero coefficients as

p0 01= p01(0)$ 0, r010 = r01(0)$ 0, r010= r10(0) = F/(v2c−V2), which yield x0₍t_1,t_{2; 0) = p}0 01cost2+ r001sint2+ F (v2 c−V2) sint1. Here, it is understood that the solution is envisaged as

x(t1,t2,o) = x0₍t_1,t_{2, 0) +}ox'(t_1,t_{2, 0) + O(}o2_).

Similarly, substituting the assumed solution and the zeroth order solution x0₍t_1,t_{2, 0)} into the first order perturbation equation (24) and balancing the harmonics one obtains the following non-zero coefficients.

sin (−t1+ 2t2)c r'−12= 2Fp0 01r 0 01(2vc−V) 4(v2 c−V2)(4vcV − V2− 3v2c) , cos (t1)c p'10= hVF (v2 c−V2)2 − F 3V 4(v2 c−V2)4 −[(p 0 01)2_{+ (r}0 01)2]VF 2(v2 c−V2)2 , cos (2t1+t2)c p'21= F2_r01(20 _{V + v} c) 4(v2 c−V2)2(−4V2− 4Vvc), sin (2t1+t2)c r'21= −F2_p0_{01(2V + v} c) 4(v2 c−V2)2(−4V2− 4Vvc), cos (t1+ 2t2)c p'12= F[−(p0 01)2+ (r01)0 2](2vc+V) 4(v2 c−V2)(−4vcV − V2− 3v2c) , sin (t1+ 2t2)c r'12= −2Fp0 01r01(20 vc+V) 4(v2 c−V2)(−4vcV − V2− 3v2c) , cos (3t1)c p'30= V 4(v2 c− 9V2)

0

1

Note that, amplitudes p01(o) and r01(o) are envisaged in the form of Taylor’s expansions

given by p01(o) = p0 01+ p'01o +1 2!p010o 2_{+ O(}o3_{), r01(}o) = r0 01+ r'01o +1 2!r010o 2_{+ O(}o3_).

Depending on the ordered form of the results, p01(o) and as r01(o) may be represented by

the first, second, etc., order terms in the above expansions. The important point here is to keep the consistency of approximations with regard to the ordered form of the solutions.

Thus, upon using the above information, the first perturbation gives

x'(t1,t2; 0) = p'2 − 1cos (2t1−t2) + r'2 − 1sin (2t1−t2) + p'−12cos (−t1+ 2t2) + r'−12sin (−t1+ 2t2) + p'10cos (t1) + p'21cos (2t1+t2) + r'21sin (2t1+t2) + p'12cos (t1+ 2t2) + r'12sin (t1+ 2t2) + p'30cos (3t1) + p'03cos (3t2) + r'03sin (3t2),

and the first order approximation for the solution of equation (19a) can be expressed as

x(t1,t2;o) = p01cos (t2) + r01sin (t2) + r010sin (t1)

+o{p'2 − 1cos (2t1−t2) + r'2 − 1sin (2t1−t2) + p'−12cos (−t1+ 2t2) + r'−12sin (−t1+ 2t2) + p'10cos (t1) + p'21cos (2t1+t2) + r'21sin (2t1+t2) + p'12cos (t1+ 2t2) + r'12sin (t1+ 2t2) + p'30cos (3t1) + p'03cos (3t2)

On the other hand, comparing the coefficients of cos (t2) and sin (t2) leads to cos (t2)c1 2 F2_r0 01 (v2 c−V2)2 −hr0 01+ 1 4(p 0 01)2_r0 01+ 1 4(r 0 01)3_{− 2}v'p0 01= 0, (26) sin (t2)c−1 2 F2_p0 01 (v2 c−V2)2 +hp0 01− 1 4p 0 01(r0₀₁₎2₋1 4(p 0 01)3− 2v'r 0 01= 0, (27)

which can be used to construct 1 2z(p₀₁₎0 2_{+ (r01)}0 2*{−p 0 01*[equation (26)] − r0_{01*[equation (27)]} = 0,} 1 2z(p₀₁₎0 2_{+ (r01)}0 2*{−r 0 01*[equation (26)] + p001*[equation (27)]} = 0. The first equation above yields

av' = 0, (28)

and the second equation gives the important result

a 2

$

0

1

%

01 and r001 represent p01(o) and r01(o), respectively, keeping the consistency of approximations. Equation (29) yields two distinct, steady state solutions:

Solution (I), a = 0 and Solution (II), h =1₄a2₊1 2

0

1

A careful inspection of equation (25), with the aid of the derivatives of the amplitudes

listed above, reveals that a0 0 represents a periodic solution (with frequency V) and

Solution (II) represents quasi-periodic motions on an invariant torus (with frequenciesV

and v). Indeed, Solution (I) is described by

x(t1;o) = r0

10sin (t1) +o{p¯'10cos (t1) + p'30cos (3t1)},

where p¯'10is p'10with p010 = r010 = 0(a = 0) and Solution (II) by equation (25) with a$ 0. On

a plot of a versush (see Figure 1), one clearly observes a bifurcation phenomenon, Solution

(II) bifurcating from Solution (I) at the critical point

hc=1 2

0

1

On the other hand, equation (28) yieldsv' = 0 for a $ 0, and one has to carry out higher

order perturbations ifv0, v1, etc., are needed in constructing

v(o) = vc+v'o +12v0o

Quasi-periodic solution, , a c Periodic solution,

Figure 1. Plot of a versush: hc= 1/2 (F/(vc2−V2))2.

It is noted that equation (29) is a bifurcation equation, providing valuable information about the behaviour of the system, and it has been obtained through the application of the IHB technique only. A formal stability analysis is not among the objectives of this

paper, but it can be shown (see the Appendix) that Solution (I) is stable for h Q hcr and

unstable for h q hcr. On the other hand, Solution (II), bifurcating from Solution (I) at

h = hcr, represents a stable family of quasi-periodic motions.

It was demonstrated earlier that the corresponding autonomous system, described by

equation (19a) with F = 0, exhibits a Hopf bifurcation ash passes through zero. It is now

observed that, the introduction of the external harmonic excitation results in a shift of

bifurcation point along theh-axis to h = hcdefined by equation (30), and the bifurcation

takes place from a family of periodic motions to a family of quasi-periodic motions, as compared to Hopf bifurcation from an equilibrium path to a family of periodic solutions, associated with the corresponding autonomous system.

Finally, it is noted that the consistency of solutions (I) and (II) as well as that of the general solution (25), up to first order approximations, has been verified by substituting these results into equation (21) and following the procedure described earlier with the aid of MAPLE.

4. CONCLUSIONS

motions bifurcates from a family of periodic motions at a critical value of the parameter. It is also observed that the external harmonic excitation results in a shift of the bifurcation point along the parameter axis, compared to the Hopf bifurcation associated with the corresponding autonomous system. A formal stability analysis is not the objective of this paper; however, a brief discussion concerning the stability properties of the solution is presented in the Appendix for completeness.

A verification scheme is outlined and applied to ascertain the validity and consistency of the ordered approximations. A symbolic computer language, namely MAPLE, is used extensively to obtain and verify various solutions. It is observed that the method lends itself conveniently to this process.

REFERENCES

1. C. H. H 1964 Nonlinear Oscillations in Physical Systems. New York: McGraw Hill. 2. D. W. J and P. S 1987 Nonlinear Ordinary Differential Equations. Oxford: Oxford

University Press; second edition.

3. A. H. N and D. T. M 1979 Nonlinear Oscillations. New York: John Wiley. 4. A. S. A and K. H 1982 Journal of Sound and Vibration 85, 129–131. A note on

‘‘Uniformly valid asymptotic solution for d2_y/dt2_{+ y = a +oy}2_’’.

5. A. S. A and K. H 1982 Journal of Sound and Vibration 95, 525–530. An intrinsic method for harmonic analysis for nonlinear oscillations.

6. K. H 1986 Multiple Parameter Stability Theory and Its Applications. Oxford: Oxford University Press.

7. K. H and S. W 1991 Journal of Sound and Vibration 148, 361–363. An extension of the Intrinsic Harmonic Balancing technique to the analysis of forced vibrations.

8. S. W and K. H 1992 Mathematical and Computer Modeling 16, 49–57. ‘‘MAPLE’’ analysis of nonlinear oscillations.

9. R. E. M 1981 Journal of Sound and Vibration 76, 150–152. A uniformly valid asymptotic solution for d2_y/dt2_{+ y = a +}oy2_.

10. K. H and S. W 1992 Bulleting of the Technical University of Istanbul 45, 87–98. A perturbation analysis of nonlinear forced oscillations.

11. K. H and S. W 1992 World Congress of Nonlinear Analysts, Tampa, FL. On the analysis of dynamic bifurcations and stability of solutions.

12. K. H and P. Y 1988 Applied Mathematical Modelling 12, 189–201. On bifurcations into nonresonant quasi-periodic motions.

APPENDIX

It can be shown that the local dynamics in the vicinity of the critical point is governed by the first order differential equations

da dt= oa 2

$

0

1

%

In order to prove the above relations, one considers perturbing the general solution (25) (x(t1,t2;o)) in the vicinity of the critical point. All amplitudes in general solution (25) are in terms of P01and r01, as indicated on the list in the text. To this end, consider the first order system (19b), where the solution of the first order system can be written as

and x2=V dx(t1,t2;o) dt1 +vc dx(t1,t2;o) dt2 .

Here, for simplicity, r01can be assumed to be zero, as in reference [11]. Suppose now that

time-invariant constant p01is a function of time denoted by a. Further,t2is also assumed

to be u(t). After these transformations, the solution takes the form

xi= fi(a,u, Vt). (A1)

By substituting equation (A1) into equation (19b) and solving for da/dt and du/dt, the

rate equations are obtained after truncating the higher order terms. This procedure has been described in a number of earlier papers for autonomous systems (see for example, Appendix D in reference [12]) and now it is applied to the non-autonomous system considered in this paper.

The stability of the steady states can be examined by considering the Jacobian of da/dt,

J = d da

0

1

$

0

1

%

Thus, evaluating the Jacobian for Solution (I) (a = 0) it is concluded that Solution (I) is stable (unstable) forh Q hc (h q hc). Similarly, evaluating the Jacobian for Solution (II)