## Robert C. Hilborn

Print publication date: 2000

Print ISBN-13: 9780198507239

Published to British Academy Scholarship Online: January 2010

DOI: 10.1093/acprof:oso/9780198507239.001.0001

Show Summary Details

# (p.589) Appendix I The van der Pol Oscillator

Source:
Chaos and Nonlinear Dynamics
Publisher:
Oxford University Press

# I.1 The van der Pol model

Limit cycles describe the spontaneous occurrence of periodic time-dependent behavior in some models. Since this behavior may be at odds with our linearly trained intuition, we examine a two-dimensional model in some detail to understand at a more physical, intuitive level how the various parts of the system interact to produce limit cycle behavior. The less-mathematically inclined reader should feel free to skim this appendix.

The model we shall describe has a venerable history in nonlinear dynamics. It was originally developed by van der Pol in the 1920s (VDP26) to describe the dynamics of a triode electronic oscillator. (A triode is an electronic vacuum tube with three elements.) We will not describe the details of van der Pol's derivation. Instead, we will try to make some plausibility arguments and then see how we can understand the appearance of limit cycle behavior. The rest of this appendix will be devoted to working through some of the analytic methods that can give us an approximate description of the van der Pol oscillator.

In the van der Pol model, the electrical charge [denoted by q(t)] passing through the triode tube is assumed to be described by an equation that is similar to that for a linear, damped, simple harmonic oscillator:

(I-1)
where γ is the so-called damping rate (representing a damping or energy loss mechanism) and ω is the frequency with which the charge would oscillate in the absence of damping. Van der Pol's insight was to model the behavior of the triode tube by allowing the damping parameter to depend on the amount of charge q. For small q the tube would tend to increase the amplitude of the oscillation due to the circuit's behavior as an amplifier. For large q, however, the amount of oscillating charge q would be limited by so-called saturation effects in the tube and the associated circuitry. (In rough terms, the tube elements and associated circuitry can supply only so much charge in a given period of time.) This behavior is modeled by making γ depend on q in such a way that for small q, γ is less than 0. A negative value for “dissipation” means that the amplitude of the oscillations grows, rather than decays. For large q, however, γ becomes positive, representing a dissipation of energy from the oscillating charge. Van der Pol chose the simple function (p.590)
(I-2)
We take γo > 0. Then for q < q0, the damping parameter γ is negative, and for q > q0, we have γ positive. If we now use Eq. (I-2) in Eq. (I-1) and introduce the variables
(I-3)
the van der Pol equation becomes
(I-4)
We now put Eq. (I-4) into our standard first-order form by introducing the variable U = Q = dQ/dτ to obtain
(I-5)
(I-6)

Next we find the fixed points for the system by setting f 1(Q, U) = 0 and f 2(Q, U) = 0. It is obvious that the only fixed point is U = 0, Q = 0, which corresponds to the no oscillation condition. Is this point stable or unstable? To answer this question, we evaluate the Jacobian matrix for the system of Eqs. (I-5) and (I-6).

(I-7)
Thus we see that the determinant Δ = 1 at Q = 0, U = 0 and that TrJ = R. Since R is positive by definition, we see that the no-oscillation fixed point is unstable. The characteristic values for this fixed point are
(I-8)
Hence, for R < 2, the fixed point is a spiral repellor. For R > 2, the fixed point is a simple repellor. The behavior for R = 0.3 is shown in Fig. I.1.

(p.591)

Fig. I.1. A state space trajectory starting near U = 0.1, Q = 0.1 for the van der Pol oscillator with R = 0.3. The fixed point at U = 0, Q = 0 is clearly a spiral repellor. Trajectories approach the limit cycle as t → ∞.

We should remind ourselves that the analysis of the fixed point's stability does not tell us what happens to those trajectories repelled by the fixed point. However, our intuition tells us that the trajectories cannot get too far from (0, 0) because eventually the damping term becomes positive and the corresponding dissipation of energy will limit the size of the trajectory.

We gain some insight into what happens by considering the time dependence of the energy associated with the charge oscillations. We can write this energy as a sum of terms that are analogous to the kinetic energy and potential energy for a mechanical oscillator ([Berg, Pomeau and Vidal, 1986], pp. 28-29):

(I-9)
Here L represents the inductance of the oscillator circuit, and C is its capacitance. We can choose units such that L = 1 and C = 1. We then find that the rate of change of this energy is given by
(I-10)
Using Eqs. (I-5) and (I-6) for the time derivatives of U and Q, we obtain (p.592)
(I-11)
We now average this change of energy over one oscillation period
(I-12)
The bar symbol indicates that we are taking the time average of each term. (Strictly speaking, we are averaging the equation over a time of about one oscillation period, which is assumed to be short compared to the damping time 1/ γo). The first term on the right-hand side of the previous equation represents the (time average) rate of energy generation in the circuit. Of course, this energy is provided by the “power supply” in the associated circuitry. The second term with its negative sign represents the rate of energy dissipation. A steady-state is achieved when these terms balance.

If we now assume that the circuit settles into sinusoidal oscillations (this in fact occurs for small values of R as we might guess from the spiral nature of the repellor), we may write

(I-13)
where Q0 is the amplitude of the charge oscillations and ω is the oscillation frequency. We can use this ansatz to determine how the amplitude of the oscillations depend on R. First we evaluate the time-averaged quantities:
(I-14)
(I-15)

We set dW/dt = 0 in Eq. (I-12) for a steady-state oscillation and then use Eqs. (I-14) and (I-15) to find

(I-16)

We see that in the limit of small R, we expect to have sinusoidal oscillations, represented in state space by a circle whose radius is given by . For larger values of R, the oscillations become non-sinusoidal. Typical behavior is shown in (p.593)

Fig. I.2. Behavior of the van der Pol oscillator for R = 3.0. On the left are plotted Q(τ) and U(τ). On the right is the Q–U state space behavior. Notice the repelling nature of the fixed point at the origin.

Fig. I-2 for R = 3.0. These oscillations, which switch rapidly from one extreme value to another, were called “relaxation oscillations” by van der Pol. The state space trajectory is shown on the right-hand side of Fig. I-2.

We now show how to analyze the stability of the limit cycle. That is, we want to know if trajectories near the limit cycle are attracted toward it or are repelled from it. The procedure we shall use is called “the method of slowly varying amplitude and phase.” It finds many applications in the study of nonlinear dynamics (see, for example, [Sanders and Verhulst, 1984]). The method is also called the KBM averaging method after the mathematicians Krylov, Bogoliubov, and Mitropusky, who developed the general formalism.

Let us begin by rewriting the differential equation (I-4):

(I-17)
If the right-hand side of the previous equation were 0, then Q would oscillate sinusoidally in time. This observation leads us to introduce the following expressions for Q and its time derivative (keeping in mind the small R limitation):
(I-18a)
(p.594)
(I-18b)

Here a is a time-varying amplitude, and Φ a time-varying phase. Note that we define a(τ) and Φ(τ) so that the previous equations are true. We do not calculate by taking the derivative of the first expression. However, for these definitions to be consistent with the usual derivative, we must have

(I-19)
or
(I-20)
We now compute the second derivative of Q by taking the derivative of the second equation in (I-18) to obtain
(I-21)
Next, we use Eqs. (I-21) and (I-18) in Eq. (I-17) to find
(I-22)
Now we want to find separate equations for ȧ and . To do that, we first multiply Eq. (I-20) by sin(τ + Φ) and Eq. (I-22) by cos(τ + Φ) and add the resulting two equations. We then multiply Eq. (I-20) by cos(τ + Φ) and Eq. (I-22) by −sin(τ + Φ) and add those. After all these algebraic manipulations we finally arrive at the desired results:
(I-23)

We should point out that Eqs. (I-23) are exactly equivalent to Eq. (I-17); we have just implemented a change of variables. Up to this point, no approximations have been made. Now we want to invoke the following crucial notion: When a trajectory gets near the limit cycle, its amplitude a and its phase Φ vary slowly over the time scale of the period of oscillation. Hence, the time derivatives of these quantities are nearly constant over one period of oscillation. If these arguments indeed apply to trajectories near the limit cycle, then we can get approximate equations for the amplitude and phase by integrating the right-hand sides of Eqs. (I-23) over one period and treating in those integrations the amplitude and phase as constants. In carrying out those integrations, we make use of the following integrals:

(I-24a)
(p.595)
(I-24b)
(I-24c)
After evaluating those integrals, we arrive at the approximate equations
(I-25)

Note that the limit cycle is reached when ȧ = 0; that is, when the same value we found before. The present method, however, allows us to find the rate at which a nearby trajectory approaches the limit cycle. To find this rate of approach, we expand the right-hand side of the first of Eqs. (I-25) in a Taylor series about the limit cycle value :

(I-26)
where the last equality follows because the derivative of f evaluated at a* is equal to −R. Equation (I-26) tells us that the trajectory approaches the limit cycle exponentially. If we let d be the difference between the trajectory amplitude a(t) and the limit cycle amplitude a*, we see that
(I-27)
where d0 is the value of that difference at time τ = 0. Thus, we conclude that the limit cycle is stable because trajectories on either side of the limit cycle approach it as time goes on.

We can use what we have just learned to construct an approximate Poincaré map function for the van der Pol oscillator. If we choose the Poincaré section to be the positive Q axis in state space, then the amplitude a(τ) gives us the location of the point at which the trajectory crosses that section. Since the time between crossings is in our units for time, Eq. (I-27) tells us that the Poincaré map function, expressed in terms of the distance from the limit cycle amplitude, must be

(I-28)
(p.596)

Fig. I.3. Sketch of the Poincaré map function for the van der Pol oscillator. We have assumed that R < 1.

and we see that the limit cycle crossing point, for which dn = 0, is a stable fixed point of the Poincaré map function.