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

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

# 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:

*γ*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)

*γ*> 0. Then for

_{o}*q < q*, the damping parameter

_{0}*γ*is negative, and for

*q > q*, we have

_{0}*γ*positive. If we now use Eq. (I-2) in Eq. (I-1) and introduce the variables

*U = Q = dQ/dτ*to obtain

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).

*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

*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.

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):

*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

*U*and

*Q*, we obtain (p.592)

*γ*). 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.

_{o}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

*Q*is the amplitude of the charge oscillations and

_{0}*ω*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:

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

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)

*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):

*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):

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

*Q*by taking the derivative of the second equation in (I-18) to obtain

*ȧ*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:

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:

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 :

*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

*d*is the value of that difference at time

_{0}*τ*= 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 *2π* 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

*d*= 0, is a stable fixed point of the Poincaré map function.

_{n}# I.2 Further Reading

B. van der Pol, “On Relaxation Oscillations,” *Phil. Mag.* (7) **2**, 978–92 (1926). This paper describes the original van der Pol oscillator.

# I.3 Computer Exercises

CEI-1. Use *Chaos Demonstrations* to study the van der Pol equation limit cycles in state space. Vary the parameter *h* (equivalent to the parameter *R* used in the text) to see how the oscillations change from simple harmonic (for small values) to relaxation oscillations for larger values.

(p.597)
CEI-2. Use *Chaotic Dynamics Workbench* to study the Shaw-Van der Pol Oscillator with the force term set to 0 (to make the state space two-dimensional). Observe the time dependence of the dynamical variables and the state space diagrams as the coefficient *A* (corresponding to *R* in the text) increases.