Duffing Equation

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

The Duffing equation (Georg Duffing, 1918)

(((∂^2) x)/(∂(t^2))) + 2γ ((∂x)/(∂t)) + k x +α (x^3) ==f Cos[ω t]

can be used to describe a driven anharmonic oscillator with linear damping.  The anharmonicity parameter is α, the damping parameter is γ, the driving amplitude is f, and the driving frequency is ω.   In the absence of damping or anharmonicity the natural oscillation frequency would be (Subscript[ω, 0]) = (Sqrt[k]).  The solutions to this equation display a dazzling variety of patterns depending upon the values of these parameters.  In this notebook we present relatively elegant functions which solve this equation and produce phase diagrams.  We also introduce the Poincaré map, which provides a tool for analyzing structures that are sometimes hidden within chaotic trajectories.  The basic idea is strobe the phase trajectory at the frequency of the driving force.  Thus, we examine the set of points

P  = { {x[(Subscript[t, k])], (x^′)[(Subscript[t, k])]},  
        k = 1, ∞}
(Subscript[t, k]) = 
      (Subscript[t, 0])  +  ((2  π k)/ω)

where (Subscript[t, 0]) is an arbitrary offset.  Although the points within such sets usually depend strongly upon both (Subscript[t, 0]) and the initial conditions, the qualitative structure or topology of the Poincaré map depends primarily upon the dynamical parameters, here {k, α, γ, f, ω}, and often provides considerable insight into the behavior of the system.  We will use a few representative examples to illustrate some of the possible behaviors but will not attempt an exhaustive study of this very rich system.

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1]

$DefaultFont={"Times", 12}; 
$TextStyle={FontFamily->"Times", FontSize->12, FontSlant->"Italic"};

Needs["Graphics`"];
Needs["Utilities`Notation`"];
Needs["Utilities`FilterOptions`"]

Symbolize[(Subscript[ω, 0])]; Symbolize[(Subscript[x, 0])]; Symbolize[(Subscript[v, 0])]

Implementation

In this section we construct functions which define and solve the Duffing equation given its parameters and initial conditions.

\(\(duffingEq::usage\  = \ \*""duffingEq[k,α,γ,f,ω] defines the Duffing equation for x[t], \n (x^″) + 2γ (x^′) + k x + α (x^3)) == f Cos[ω t], \n given elasticity parameter k, nonlinearity parameter α, damping rate γ, driving amplitude f, and driving frequency ω."";\)
\)
duffingEq[k_, α_, γ_, 
        f_, ω_]  := (x^″)[t]  +  
          2  γ (x^′)[t]  +  
          k x[t] + α x[t]^3 == 
        f Cos[ω t]

DSolve[duffingEq[k, 0, γ, f, ω], x, t]

It is useful to provide a place for optional parameters when constructing a function which invokes a Mathematica command to perform a task which may require some attention by the user.  For example, in this case it may be necessary to increase MaxSteps if the range of integration is large enough to exceed the default value for NDSolve.  Placing opts___Rule last in the argument list with the triple-underscore pattern indicates that zero or more replacement rules may be included there.

Clear[duffingSol];
\(\(duffingSol::usage\  = \*""duffingSol[k,α,γ,f,ω,(Subscript[x, 0]),(Subscript[v, 0]),tmax,opts] uses NDSolve to produce a numerical solution to the Duffing equation given initial conditions x[0]==(Subscript[x, 0]) and (x^′)[0]==(Subscript[v, 0]) for {t,0,tmax}.  See duffingEq for other parameter definitions.  The optional parameters opts are passed to NDSolve."";\)
		\)
duffingSol[k_, α_, γ_, 
          f_, ω_, (Subscript[x, 0]) : _, (Subscript[v, 0]) : _, tmax_, opts___Rule][
        t_Symbol]  := {x[t], (x^′)[
            t]} /. NDSolve[{duffingEq[k, α, γ, f, ω], 
              x[0] == (Subscript[x, 0]), (x^′)[0] == (Subscript[v, 0])}, 
            x, {t, 0, tmax}, 
            opts][[1]]

Since duffingPlot calls duffingSol to produce the solution it seeks to plot, and both may benefit from options, it is convenient to employ FilterOptions, found in the package  Utilities`FilterOptions`, to choose those options appropriate to NDSolve and those appropriate to ParametricPlot.  The With statement makes local assignments based upon the output of FilterOptions and then passes those option lists to the appropriate commands.  However, those lists are not defined outside this context.  (Check for yourself that optPlot and optND remain undefined after executing duffingPlot.)  Thus, this method allows options for both commands to be included at the end of the argument list.  Note that the options may occur in any order and need not be grouped by command.  Also note that options are applied in order but only once.  Therefore, by placing the user-defined options first in ParametricPlot defaults for labels can be included at the end without overriding the user's selections, if present.

Clear[duffingPlot];
\(\(duffingPlot::usage\  = \ \*""duffingPlot[k,α,γ,f,ω,(Subscript[x, 0]),(Subscript[v, 0]),tmin,tmax,opts] produces a phase diagram for the Duffing equation.  The equation is solved for {t,0,tmax}, but the plot may be restricted to {t,tmin,tmax}.  The optional parameters opts can regulate either NDSolve or ParametricPlot as appropriate.  See duffingEq and duffingSol for parameter definitions."";\)
\)
duffingPlot[k_, α_, γ_, 
        f_, ω_, (Subscript[x, 0]) : _, (Subscript[v, 0]) : _, tmin_, tmax_, 
        opts___Rule]  :=  
      With[{optPlot  =  FilterOptions[ParametricPlot, opts] , 
          optND  = FilterOptions[NDSolve, opts]},  
        ParametricPlot[
          Evaluate[duffingSol[k, α, γ, f, ω, (Subscript[x, 0]), 
                (Subscript[v, 0]), tmax, optND][t]], {t, tmin, tmax},  optPlot, 
          AspectRatio -> 1,  Frame -> True, 
          FrameLabel -> {"position",  "velocity"},  
          PlotLabel -> "Duffing Equation"]]

Finally, the following function constructs a Poincaré set from a phase trajectory given as a pair of functions.  The user must choose the proper strobe frequency.

poincare[{x_[t_], v_[t_]}, ω_, tmin_, tmax_] := 
      Table[{x[t], v[t]}, {t, tmin, tmax, ((2  π)/ω)}]

Examples

Conservative systems

In the absence of damping or driving forces, the Duffing equation describes the motion of a particle in a potential well which conserves energy.  The behavior can then be understood using energy diagrams.  The shape of the potential is obtained by integrating the restoring force.

ν[k_, α_][
        x_]  = (∫_0^x)(k x + α (x^3)) \[DifferentialD]x   // Expand

Below we plot three cases for k > 0.  The red curve with α = 0 represents the harmonic potential, the blue curve with α = +1 is more strongly confining for large displacements, while the green curve with α = -1 is less strongly confining for moderate displacements and becomes repulsive for large displacements.  The motion for the strictly confining potentials will be periodic between the turning points determined by the condition that the total energy is potential, whereas character of the motion for potentials with α < 0 depends upon both the initial position and total energy.  If the system begins with total energy and position within the confined region the motion will be periodic, but if the energy is too large or the position is outside the well the particle escapes to infinity.

Plot[
        Evaluate[
          Table[ν[1, α][x], {α, -1, 
              1}]], {x, -2, 2}, 
        PlotRange -> {Automatic, {-1, 1}}, 
        PlotStyle -> {RGBColor[0, 1, 0], RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}];

The phase diagram for simple harmonic motion in a parabolic well is circular as both position and velocity are represented by complementary sinusoidal functions (circular functions), but for nonparabolic wells the phase diagrams are slightly deformed and the periods differ from 2  π/(Subscript[ω, 0]).

plot1  =  
        duffingPlot[1, 0, 0.0, 0, 0, 0.5, 0, 0, 10,  
          PlotLabel -> "Free Conservative Oscillations", 
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
plot2  =  
        duffingPlot[1, 1, 0.0, 0, 0, 0.5, 0, 0, 10, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
plot3  =  
        duffingPlot[1, -1, 0.0, 0, 0, 0.5, 0, 0, 10, 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
Show[{plot1, plot2, plot3}, 
        DisplayFunction -> $DisplayFunction];

It is also of interest to consider the case k < 0, α > 0, which produces double-well potentials.  The motion is still confined and periodic, but the phase trajectory is more interesting.

Plot[
        Evaluate[
          Table[ν[-1, α][x], {α, (1/3), 1, 
              (1/3)}]], {x, -3, 3}, 
        PlotRange -> {Automatic, {-1, 1}}, 
        PlotStyle -> {RGBColor[0, 1, 0], RGBColor[1, 0, 0], 
            RGBColor[0, 0, 1]}];

plot1  =  
        duffingPlot[-1, (1/3), 0.0, 0, 0, 0, 1, 0, 10,  
          PlotLabel -> "Free Conservative Oscillations", 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
plot2  =  
        duffingPlot[-1, (2/3), 0.0, 0, 0, 0, 1, 0, 10, 
          PlotStyle -> RGBColor[1, 0, 0], 
          DisplayFunction -> Identity];
plot3  =  
        duffingPlot[-1, 1, 0.0, 0, 0, 0, 1, 0, 10, 
          PlotStyle -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
Show[{plot1, plot2, plot3}, 
        DisplayFunction -> $DisplayFunction];

Driven undamped oscillator

Harmonic oscillator

If the natural and driving frequencies are commensurate (in other words, their ratio is rational), we find closed curves in phase space with multiple periodicities.

duffingPlot[1, 0, 0, 1, 0.25, 0, 0, 0, 50, MaxSteps -> 2000, 
        PlotPoints -> 100];

duffingPlot[1, 0, 0, 1, 3, 0, 0, 0, 50, MaxSteps -> 2000, 
        PlotPoints -> 100];

If the natural and driving frequencies noncommensurate, the motion appears to be space filling — as tmax increases, the phase trajectory tends to completely fill a circle in phase space whose radius depends upon the total energy specified by the initial conditions.  (Verify by make plots for several values of tmax!)

duffingPlot[1, 0, 0, 1, π, 0, 0, 0, 200, MaxSteps -> 5000, 
        PlotPoints -> 100];

Nevertheless, we see below that the Poincaré map is just a circle.

Clear[sol];
sol[
          t_] = duffingSol[1, 0, 0, 1, π, 0, 0, 1000, 
            MaxSteps -> (10^5)][t];
pts = poincare[sol[t], π, 0, 1000];
\(\(ListPlot[pts, AspectRatio -> Automatic, Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k]^′)""}];\)\)

However, the sequence of points is fairly complicated — the phase points hops from one place to another rather than filling the circle in an orderly fashion.  Below we display sequence of positions and velocities at each (Subscript[t, k]); to minimize clutter, we show just the first 100 points.

\(\(ListPlot[Take[First /@ pts, 100], Frame -> True, 
        FrameLabel -> {"k", \*""\!\(x\_k\)""}];\)\)
\(\(ListPlot[Take[Last /@ pts, 100], Frame -> True, 
        FrameLabel -> {"k", \*""(Subscript[x, k]^′)""}];\)\)

Both quantities have show 3 interleaved curves, with the phase point hopping from one curve to the next.  To demonstrate this, we plot every third point.  Similar figures are obtained starting with either the second or third points.

\(\(ListPlot[First /@ Take[pts, {1, Length[pts], 3}], Frame -> True, 
        FrameLabel -> {"k", \*""\!\(x\_\(3  k\)\)""}];\)\)

Although each of these curves appears to be sinusoidal, the interleaving of three such curves makes the mapping from (Subscript[x, k]) to (Subscript[x, k + 1]) very complicated — we know that the pairs {(Subscript[x, k]), (Subscript[x, k + 1])} fall on a simple curve but they hop around that curve almost unpredictably.  If we take enough samples, this ellipse will eventually fill in.  We leave the interesting question of whether this Poincaré map is truly continuous in the limit t -> ∞ to mathematicians.  This question cannot be answered by numerical simulation and requires more sophisticated mathematical analysis.

\(\(ListPlot[Partition[First /@ pts, 2], AspectRatio -> Automatic, 
        Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k + 1])""}];\)\)

Anharmonic oscillator

Next, consider a driven but undamped anharmonic oscillator.  The phase trajectories are not closed and fill regions of the phase space that depend upon the dynamical parameters.  A frequency scan is shown below for initial conditions {(Subscript[x, 0]), (Subscript[v, 0])} = {0, 1}.  You should try other scans as well.

Table[
        duffingPlot[1, 1, 0, 1, ω, 1, 0, 1, 100, 
          MaxSteps -> (10^4)], {ω, 0.5, 1.5, 0.5}];

Once again, we use the Poincaré map to analyze the phase-space filling pattern.  The figures below show that the Poincaré map for {k, α, γ, f, ω} = {1, 1, 0, 1, 
          0.5} is again a relatively simple closed curve but it is populated by 9 interleaved qausisinusoidal recurrence relations (Subscript[x, k]) -> (Subscript[x, k + 9]).

Clear[sol];
sol[
          t_] = duffingSol[1, 1, 0, 1, 0.5, 0, 1, 10000, 
            MaxSteps -> 5×(10^5)][
          t];
pts = poincare[sol[t], 0.5, 0, 10000];
\(\(ListPlot[pts, Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k]^′)""}];\)\)

\(\(ListPlot[First /@ pts, Frame -> True, 
        FrameLabel -> {"k", \*""\!\(x\_k\)""}];\)\)
\(\(ListPlot[Last /@ pts, Frame -> True, 
        FrameLabel -> {"k", \*""(Subscript[x, k]^′)""}];\)\)

\(\(ListPlot[First /@ Take[pts, {1, Length[pts], 9}], Frame -> True, 
        FrameLabel -> {"k", \*""\!\(x\_\(6  k\)\)""}];\)\)

Finally, consider a double-well potential with modest driving force.  The phase trajectory appears to be rather complicated, with peculiar kinks and loops, and tends to fill an area of phase space.

duffingPlot[-0.25, 1, 0, 0.5, 2, 1, 0, 0, (10^2), 
        MaxSteps -> (10^4)];

In this case the Poincaré map consists of three disconnected closed curves between which the phase point hops in the sequence indicated.  Even though the phase trajectory is ultimately space filling, there is underlying order.

Clear[sol];
sol[
          t_] = duffingSol[-0.25, 1, 0, 0.5, 2, 1, 0, (10^3), 
            MaxSteps -> (10^5)][t];
pts = poincare[sol[t], 2, 100, (10^3)];
\(\(ListPlot[pts, Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k]^′)""}, 
        Epilog -> {Text["1", {\(-0.8\), 0.25}], 
            Text["2", {0.85, 0.1}], 
            Text["3", {\(-0.25\), \(-0.8\)}]}];\)\)

\(\(Table[
        ListPlot[Take[pts, {ℓ, Length[pts], 3}], Frame -> True, 
          FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k]^′)""}], {ℓ, 1, 3}];\)\)

Damping

When γ > 0, damping removes energy from the system such that without a driving force oscillations gradually decay away and the phase diagram becomes an inward spiral.

duffingPlot[1, 0, 0.1, 0, 0, 1, 0, 0, 30];

Driven, damped, anharmonic oscillator

A complete catalog of the complex motions that become possible when both driving and damping forces are applied to an anharmonic oscillator is beyond the scope of this notebook, but it is fun to  present a few interesting cases.

Modest driving and damping forces

We begin at the origin, for simplicity, and allow the driving force to seek a periodic motion.  The initial transients are suppressed by the damping forces and the system eventually settles upon a periodic orbit through phase space.  The shape of the limiting orbit depends upon α, but for a wide range of initial conditions is independent of those initial conditions; hence, the limiting orbit is called an attractor and the domain of initial conditions which lead to that attractor is called the basin of attraction.  The basin of attraction is clearly limited for negative α, but for positive α is the entire phase plane.

duffingPlot[1, 1, 0.05, 0.1, 0.5, 0, 0, 0, 100, 
        MaxSteps -> 3000, PlotPoints -> 100, 
        AspectRatio -> Automatic];

Below we generate a table of plots that show the limiting orbits for an array of initial conditions, with the early transient behavior suppressed.  Plotting this array together shows that the limiting orbit for specified {α, γ, f}is independent of initial conditions.

stuff = 
        Table[duffingPlot[1, 1, 0.05, 0.1, 0.5, x0, v0, 150, 300, 
            MaxSteps -> 3000, PlotPoints -> 100, 
            DisplayFunction -> Identity], {x0, 0, 1, 0.2}, {v0, 0, 1, 
            0.2}];

Show[stuff, DisplayFunction -> $DisplayFunction, 
        AspectRatio -> Automatic];

The period is determined by the driving force, so the Poincaré map quickly settles upon a single point as transients decay.

Clear[sol];
sol[
          t_] = duffingSol[1, 1, 0.05, 0.1, 0.5, 0, 1, 3000, 
            MaxSteps -> (10^5)][t];
pts = poincare[sol[t], 0.1, 0, 3000]

Strong driving, modest damping

In the example below, the driving frequency is just above the natural frequency, but its amplitude is quite large.  Nevertheless, after some rather chaotic transients, the system appears to settle upon a stable attractor.

duffingPlot[1, 1, 0.05, 10, 1.25, 0, 0, 0, 100, 
        MaxSteps -> 4000, PlotPoints -> 100];

The shape of the attractor can be clarified by using tmin to suppress transients.

duffingPlot[1, 1, 0.05, 10, 1.25, 0, 0, 100, 200, 
        MaxSteps -> 8000, PlotPoints -> 100];

The figure below illustrates the dependence of the limiting orbit upon α.

stuff = 
        Table[{duffingPlot[1, α, 0.05, 10, 1.25, 0, 0, 100, 200, 
              MaxSteps -> 8000, PlotPoints -> 100, 
              PlotStyle -> Hue[α], 
              PlotRange -> {{-22, 22}, {-22, 22}}, 
              DisplayFunction -> Identity], 
            Graphics[
              Text[ToString[α], {10, 6 - 12  α}, 
                TextStyle -> {FontColor -> 
                      Hue[α]}]]}, {α, 0, 1, 0.2}];

Show[stuff, DisplayFunction -> $DisplayFunction, 
        PlotRange -> {{-22, 22}, {-22, 22}}, 
        AspectRatio -> Automatic];

For low-frequency, large-amplitude driving forces the attractor has a rather complicated shape in which the system jumps from one side to the other for no apparent reason (other than symmetry, perhaps) after executing some rather intricate loops.

duffingPlot[1, 1, 0.05, 10, 0.1, 0, 0, 0, 100, MaxSteps -> 6000, 
        PlotRange -> All, PlotPoints -> 100];

Although the limiting orbit remains quite complicated even after suppressing initial transients, graphical experimentation suggests that it is independent of initial conditions nonetheless.  However, formal proof of this hypothesis is beyond the scope of this presentation.

stuff = 
        Table[duffingPlot[1, 1, 0.05, 10, 0.1, x0, v0, 200, 300, 
            MaxSteps -> 10000, PlotRange -> All, 
            PlotPoints -> 100, DisplayFunction -> Identity, 
            Frame -> False, PlotLabel -> " ", 
            Axes -> False], {x0, -1, 1, 1}, {v0, -1, 1, 
            1}];

Show[GraphicsArray[stuff]];

Double-well potential with strong driving force and weak damping

Next, consider a double-well potential including both damping and driving forces.

duffingPlot[-1, 0.5, 0.1, 1, 0.5, 0. , 0.2, 0, 500, 
        MaxSteps -> (10^6), PlotRange -> All, 
        PlotPoints -> 100];

This system does not appear to approach a periodic trajectory.  Again, we use the Poincaré map to search for underlying periodicities but this time we obtain a fractal structure, described as a strange attractor, instead of a closed curve.  The phase point jumps around unpredictably on this fractal from one period of the driving force to the next.  Thus, the behavior of this system is truly chaotic.  We suppressed the first (Subscript[t, min]) = (10^3) time steps to permit transients to decay.  Although the detailed shape of the strange attractor depends upon the value of (Subscript[t, min]), the observation that it is not a simple closed curve does not.

Clear[sol];
sol[
          t_] = duffingSol[-1, 0.5, 0.1, 1, 0.5, 0, 0.2, (10^5), 
            MaxSteps -> (10^7)][t];
pts = 
        poincare[sol[t], 0.5, 1000, (10^5)];
\(\(ListPlot[pts, Frame -> True, 
        FrameLabel -> {\*""(Subscript[x, k])"", \*""(Subscript[x, k]^′)""}];\)\)

\(\(ListPlot[First /@ pts, Frame -> True, 
        FrameLabel -> {"k", \*""\!\(x\_k\)""}];\)\)
\(\(ListPlot[Last /@ pts, Frame -> True, 
        FrameLabel -> {"k", \*""(Subscript[x, k]^′)""}];\)\)

(((∂^2) x)/(∂(t^2)))+γ ((∂x)/(∂t))- x (1 - α (x^2))==f Cos[ω t]

(((∂^2) x)/(∂(t^2)))+ ((Subscript[ω, 0])^2) x +((Subscript[γ, 1])-((((Subscript[γ, 2])(((∂x)/(∂t)))))^2)) ((∂x)/(∂t)) == 0