In this section, we analyze differential equation modeling Duffing oscillator,
named after Georg Duffing (1861--1944). It provides an example of a periodically forced
oscillator with a nonlinear elasticity. The Duffing equation is an example of
a dynamical system that exhibits chaotic behavior. For this type of system,
there are frequencies at which the vibration suddenly jumps-up or down, when
it is excited harmonically with slowly changing frequency. The frequencies at
which these jumps occur depend upon whether the frequency is increasing or
decreasing and whether the nonlinearity is hardening or softening. Between
these frequencies, multiple solutions exist for a given frequency of
excitation, and the initial conditions determine which of these solutions
represents the response of the system.
Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part III of the course APMA0340
Introduction to Linear Algebra with Mathematica
In order to analyze chaos in a simple system, we will consider the
basic equation, called after the German engineer and freelance
researcher Georg Wilhelm Christian Caspar Duffing
(1861--1944). Duffing had a heart valve defect that freed him from military service.
The Duffing equation was originally derived as a model for describing the forced vibrations of inductrial machinery. It provides a very good approximation of the motion of a damped, driven inverted pendilum with a torsional restoring force. The unforced pendulum equation is
where θ is the angle of inclinantion of the bob of mass m from the downward vertical position, k is the torsional spring constant, and γ is the damping coefficient. Here dot denotes the derivative with respect to time t (Newton's notation). Upon approximation the sine function by its two term Maclausin polynomial, we get
where we keep the same notations for coefficients k and γ. Linear terms can be united.
Nowadays, the term "Duffing equation" is used for any equation that
describes an oscillator that has a cubic stiffness term, regardless of
the type of damping or excitation.
This, however, was not the case in Dufﬁng’s original work. The above
equation can display chaotic behavior. For ω_{0}²>0, the
Duffing oscillator can be interpreted as a forced oscillator with a spring
whose restoring force is written as F = - ω_{0}²x - βx^{3}. When β>0, this equation
represents a "hard spring," and for β<0, it represents a "soft
spring." If β<0, the phase portrait curves are closed.
We analyze the unforced system first.
For ω_{0}²<0 (so ω_{0} is pure imaginary),
the Duffing oscillator describes the dynamics of a point mass in a double well
potential, and it can be regarded as a model of a periodically forced steel
beam which is deflected toward the two magnets. It is known that chaotic
motions can be observed in this case.
When there is no resistance proportional to the velocity (γ = 0), the
Duffing equation can be integrated upon multiplication by the velocity:
The function in parenthesis \( H= \frac{1}{2}\, \dot{x}^2 +
\frac{1}{2}\,\omega_0^2 x^2 + \frac{1}{4}\, \beta\,x^4 \) is called
the Hamiltonian for the Duffing equation. Then
For positive coefficients ω²_{0} and β, the solution
is bounded: \( |x| \le \sqrt{2H/\omega_0^2} \) and
\( |\dot{x}| \le \sqrt{2H} . \)
When γ≥0, the function E(t) satisfies
therefore, E(t) is a Lyapunov function, and every trajectory moves on
the surface of E(t) toward the equilibrium position---the origin. When
the Duffing equation has nonzero coefficients, there exist stationary
solutions that are obtained upon solving the cubic equation
\[
\omega_0^2 x + \beta\,x^3 =0 \qquad \mbox{or} \qquad x \left( \omega_0^2 +
\beta\,x^2 \right) =0 .
\]
So we get two other equilibrium solutions \( x =
\pm \sqrt{- \omega_0^2 / \beta} , \) when \(
\omega_0^2 \beta < 0 . \) To analyze their stability, we apply
the linearization procedure, so we calculate the Jacobian matrix
and it is found that this equilibrium is stable for \(
\omega_0^2 \ge 0 . \) On the other hand, the eigenvalues of the
equilibria \( x = \pm \sqrt{- \omega_0^2 /\beta} \)
are
Elliptic functions are functions of two variables. The first variable
might be given in terms of the amplitude φ. The second variable might
be given in terms of the parameter m, or as the elliptic
modulus k, where k² = m, or in terms of the
modular angle α, where m = sin² α. Out of twelve Jacobi
elliptic functions, there are two other elliptic functions closely
related to the previous one: the elliptic sine
sn and the delta amplitude dn. They are
defined by
Mathematica has build-in corresponding
functions: JacobiSN, JacobiCN, and JacobiDN.
The number m = k² (0 < m < 1)
is called elliptic modulus and the number φ is
called Jacobi amplitude and it is denoted by am(t, m).
The following identities are hold:
This solution is periodic and bounded if \( \omega_0^2 +
\beta\,x_0^2 > 0 \) for any β. In case when
\( \omega_0^2 +
\beta\,x_0^2 < 0 , \) solutions are bounded when β > 0;
other wise, they are unbounded.
When the initial displacement is zero, x_{0} = 0 but
the velocity is not x_{1} ≠ 0, the solution of the Duffing
equation is expressed via elliptic sine function because it is a solution of
the second order differential equation:
This solution is bounded and periodic with period 2K(1/2) ≈
3.3715, where K is the complete elliptic integral
of the first kind: \( \displaystyle K(m) = \int_0^{\pi /2}
\frac{{\text d}\phi}{\sqrt{1 - m\, \sin^2 \phi}} = \int_0^1 \left( 1 - s^2 \right){-1/2} \left( 1 - m \,s^2 \right){-1/2} {\text d}s . \)
This solution is unbounded and periodic with period 2K(3^{½}/2) ≈ 4.313.
Here Jacobi elliptic function nc has default notation in Mathematica
as JacobiNC.
We choose the units of length so that the minima are at \( x = \pm 1 , \) and the
units of energy so that the depth of each well is at -1/4. The potential is given by \( V(x) = \frac{x^4}{4} - \frac{x^2}{2} . \) The force is given by \( F(x) = - {\text d}V/{\text d}x = x- x^3 . \) As usual we solve the second order differential equation F = ma by
expressing it as two first order differential equations
Because of energy conservation, one can clearly never get chaos from the
motion of a single degree of freedom. We therefore add both a driving force
and damping, in order to remove energy conservation. The equations of motion
become
\[
\frac{{\text d}x}{{\text d}t} = v , \qquad \frac{{\text d}v}{{\text d}t} = x -
x^3 - \gamma \,v + A\, \cos \left( \omega t \right) ,
\]
where γ is the friction coefficient, and A is the strength of the
driving force which oscillates at a frequency ω. We will
see that a transition to chaos now occurs as the strength of the driving force
is increased. We will fix coefficients as \( \gamma = 0.05, \ \omega = 1.1 , \) and initially set
A = 0.1 (which will be in the non-chaotic regime).
We now start the particle off at rest at the origin and integrate the
equations of motion. We will go up to t = 200. In order to have
sufficient accuracy we need to set the option
MaxSteps of NDSolve to a larger value than
the default (otherwise Mathematica will grumble).
This looks complicated, but in fact, most of the plot shows the initial period
of time during which the motion is approaching its final behavior which is
much simpler. The early behavior is called an "initial transient". One now
sees a simple periodic orbit around the minimum of the potential at
x = 1. We now investigate how the phase space plots change when the
strength of the driving force is changed. It is convenient to define a
function, which we call "solution[A, tmax]", which will integrate the
equations for a given value of the forcing, A, up to a time tmax. We
use a delayed assignment because we want the evaluation only to be done when
we call the function later with particular values of A and tmax.
Note that the particle has moved through both of the wells. However, again, most of this complexity is due to an initial
transient. If we look at the behavior at later times, we take from 740 to 802, we see a much simpler curve:
graph[740, 802]
Next let's increase the driving force to 0.338.
sol4 = solution[0.4, 1000]
graph[700, 800]
Finally, we present animation of periodic change the velocity versus
displacement of the chaotic attractor of the Duffing oscillator for
β = 1, ω_{0}² = 1, γ =
0.2, A = 0.3, and ω = 1.
A useful way of analyzing chaotic motion is to look at what is called the Poincaré section. Rather than considering the
phase space trajectory for all times, which gives a continuous curve, the Poincaré section is just the discrete set of
phase space points of the particle at every period of the driving force, i.e. at \( t= 2\pi /\omega , \ 2\pi /\omega , 3\pi /\omega , \ldots , \) etc. Clearly for a periodic orbit the Poincaré section is a single point, when the period has doubled it consists of two points, and so on.
We define a function, poincare, which produces a Poincaré section for given values of A, γ, and ω, in which the first ndrop periods are assumed to be initial transient and so are not plotted, while the subsequent ndrop eriods are plotted. The point size is given by the parameter psize. Note that the function g[{xold, vold}] maps a point in phase space {xold, vold} at time t to the point in phase space {x , v} one period T later
This strange diagram is the strange attractor.
It is the limiting set of points to which the trajectory tends (after the
initial transient) every period of the driving force.
We now want to study the structure of the strange attractor in more detail. To do this we will blow up the region in the
shaded box in the above figure. To do so, we define a function
poincare2, in which we specify the range of x and v to
be included in the plot:
We zoom in to the region marked by the box in the figure above (this is
quite slow because I had to increase the time of integration to 10^{5} to get enough points in this region)
Zooming in to the region in the box in the above figure, shows many of the same features as in the figure. (This takes a very long time, because I had to increase the length of the integration to 10^{6} to get a significant density of points in this small region. You may want to skip the execution of the rest of this notebook, and just look at it.)
Integrating a differential equation, as we have done here, is much more time consuming than iterating a map, such as the logistic map. People have therefore investigated maps which have similar behavior to that of driven, damped
differential equations like the Duffing equation. One popular choice is the Hénon map:
in which two variables, x and y, are iterated. The parameters a and b can be adjusted to get a transition to chaos. In the chaotic regime, the points converge to a strange attractor similar to the one found for the Duffing equation. Note, in particular, the way it folds back on itself. A discussion of using
Mathematica to display the Hénon map is given in
Zimmerman and Olness, Mathematica for Physics, Pearson, 2nd edition, 2002. ISBN-13: 978-0805387001
Going back to the Duffing equation, you can try different values of the parameters gamma and omega and see where the period doubling transition to chaos occurs. The function bifurcation below, is similar to poincare except that it scans a range of values of A, and gives the x-values on the attractor for each A.
The following command scans 200 values of A in the range from 0.3 to 0.35, with γ =0.1 and ω =1.4. It starts with v = 0 and x = 1, i.e.at rest in the x = 1 minimum. The equations are integrated for 1500 periods of the driving force with the first 1000 being ignored and the value of x plotted at the end of each of the remaining 500 periods.
One clearly sees period doubling leading into chaos. With the parameters chosen, in the region of limit cycles the
system is either in the well at positive x or in the well at negative x, depending on the precise value of A, but does not hop between wells. For example, for A < 0.309 there are actually two separate fixed points rather than a period two limit cycle. (This accounts for the “spotty” nature of the lines, when one is blank the other is colored because the particle has converged on to one fixed point or the other.) However, in the chaotic region, the system goes between different wells.
bifurcation2, below, is similar to bifurcation, except that it also specifies the range of values of x (the vertical axis) that will be plotted.
Byrd, P.F. and Friedman, M.D., Handbook of Elliptic Integrals for
Engineers and Scientists, Springer, 1971.
Duffing, G. (1918), Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung [Forced oscillations with variable natural frequency and their technical relevance; PhD Thesis] (in German), Heft 41/42, Braunschweig: Vieweg.
J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer-Verlag, 1983.
Holmes, P.J. and Rand, D.A., The bifurcations of Duffing's equation: An application of catastrophe theory, Journal of Sound and Vibration, 1976, 44, pp. 237--253.
Holmes, P., A nonlinear oscillator with a strange attractor, Philosophical Transactions of the Royal Society A, 292, 419-448, 1979.
Holmes, C. and Holmes, P., Second order averaging and bifurcations to subharmonics in Duffing's equation, Journal of Sound and Vibration, 1981, 78, pp. 161--174.
Holmes P. and Whitley, D., On the attracting set for Duffing's equation, II. A geometrical model for moderate force and damping, Physica, 1983, Vol. 7D, pp. 111--123.
Kovacic, I. and Brennan, M.J., The Duffing Equation: Nonlinear
Oscillators and their Behaviour, Wiley, 2011, ISBN: 978-0-470-97783-5
J.M.T. Thompson and H.B. Stewart, Nonlinear Dynamics and Chaos (2nd edition), John Wiley & Sons, 2002.
Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions