Chapter 15
The Duffing Equation

This notebook is from VisualDSolve, Visualizing Differential Equations with Mathematica, by Antonín Slavík, Stan Wagon, and Dan Schwalbe, copyright 2022. Unauthorized reproduction is not permitted.

The Duffing equation is a much studied differential equation arising from a forced system (typically an oscillator subject to a periodic force; for physical details see [BCB] or [Str] and for more on the fascinating mathematics of this equation see [GH, chap. 2]). Duffing orbits exhibit a fascinating variety of behaviors typical of chaos. In this chapter we show some interesting orbits, and discuss ways to recognize and deal with sensitive situations. The techniques that can help us recognize an inaccurate solution are quite general, and should be used whenever a software user seeks evidence regarding the accuracy of the output of a numerical algorithm.

15.1 The Stable Case

The Duffing equation is x''-a x+b(x^3)+c x'=dcos(et). The equation is included in the VisualDSolve package.

Get["VisualDSolve`"]

DuffingEquation[{a,b,c,d,e}, x[t]]

We restrict our analysis to the case when the parameters a, b, c, d, e are nonnegative. The simplest situation occurs for c=d=0 (unforced and undamped case).

SecondOrderPlot[DuffingEquation[{0.5,0.5,0,0,1}, x[t]],x[t],{t,0,30},{x,-1.7,1.7},{x',-0.8,0.8},InitialValues->Table[{x0,0},{x0,{-(Sqrt[2]),-1.3,-1.2,-1.1,1.1,1.2,1.3,(Sqrt[2]),1.5,1.6}}],FrameLabel->{"x","(x^(′))"},NullclineShading->True,Nullclines->True,ShowEquilibria->True,RotateLabel->False]

From the corresponding pair of first-order equations, we immediately see that the system has three equlibria: (0,0) and (±(Sqrt[(a/b)]),0).

duffingSystem=ToSystem[DuffingEquation[{a,b,0,0,1}, x[t]],x[t],t,AuxiliaryVariables->{y}]

The system possesses a Hamiltonian function:

hamiltonian=HamiltonianFunction[duffingSystem,{x[t],y[t]}]

Hence, the orbits are indeed closed and coincide with the contours of the Hamiltonian function.

The separatrix passes through the origin and corresponds to the zero value of the Hamiltonian function. It intersects the x-axis at  (±(Sqrt[((2a)/b)]),0).

Solve[hamiltonian==0/. y->0,x]

15.1.1  Detailed Analysis

Solutions of the preceding unforced and undamped Duffing equation can be expressed analytically in terms of the Jacobi elliptic functions. Because the problem is symmetric with respect to the y-axis and each orbit intersects the x-axis, it is enough to consider initial-value problems with x(0)=(Subscript[x, 0]), x'0)=0, where (Subscript[x, 0])>0. We follow the approach presented in [Mey].
For t=0, the value of the Hamiltonian is -((a (Subscript[x, 0]^2))/2)+((b (Subscript[x, 0]^4))/4). Since the value is preserved along the orbit, we get -((a ((x(t))^2))/2)+((b ((x(t))^4))/4)+((x'(((t))^2))/2)=-((a (Subscript[x, 0]^2))/2)+((b (Subscript[x, 0]^4))/4), or equivalently ((x'(t))^2)=-a (Subscript[x, 0]^2)+((b (Subscript[x, 0]^4))/2)+a ((x(t))^2)-((b ((x(t))^4))/2).

Let us start with the separatrix, where (Subscript[x, 0])=(Sqrt[((2a)/b)]). The previous equation reduces to ((x'(t))^2)=a ((x(t))^2)-((b ((x(t))^4))/2), and Mathematica finds the solution easily.

FullSimplify[DSolve[{((x'[t])^2)==a ((x[t])^2)-(1/2) b ((x[t])^4),x[0]==(Sqrt[((2 a)/b)])},x[t],t]]//Quiet

Here is an illustration:

With[{a=0.5,b=0.5},Plot[(Sqrt[2]) (Sqrt[(a/b)]) Sech[(Sqrt[a]) t],{t,-8,8}]]

Next, consider the orbits inside the right loop of the separatrix. Each orbit intersects the x-axis at a point ((Subscript[x, 0]),0), where (Sqrt[(a/b)])≤(Subscript[x, 0])<(Sqrt[((2a)/b)]). Mathematica expresses the solution in terms of the Jacobi cd function.

FullSimplify[DSolve[{((x'[t])^2)==-a (x0^2)+((b (x0^4))/2)+a ((x[t])^2)-(1/2) b ((x[t])^4),x[0]==x0},x[t],t]]//Quiet

Note that the term under the square root is negative, i.e., the first argument of the cd function is imaginary. Since we prefer to stay in the real domain, let us derive the solution in another form by hand calculation. In the upper half of the orbit, we have x'(t)≥0, and consequently x'(t)=(Sqrt[-a (Subscript[x, 0]^2)+((b (Subscript[x, 0]^4))/2)+a ((x(t))^2)-((b ((x(t))^4))/2)]). This is a first-order separable differential equation. Integration gives

({{({{t=(∫_0^x)((du)/(Sqrt[-a (Subscript[x, 0]^2)+((b (Subscript[x, 0]^4))/2)+a (u^2)-((b (u^4))/2)]))=(∫_0^(x/(Subscript[x, 0])))(((Subscript[x, 0])dv)/(Sqrt[-a (Subscript[x, 0]^2)+((b (Subscript[x, 0]^4))/2)+a (Subscript[x, 0]^2)(v^2)-((b (Subscript[x, 0]^4)(v^4))/2)]))=(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[-a +((b (Subscript[x, 0]^2))/2)+a (v^2)-((b (Subscript[x, 0]^2)(v^4))/2)]))=(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[-a(1-(v^2)) +((b (Subscript[x, 0]^2))/2)(1-(v^4))]))f=(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[(1-(v^2)) (-a+((b (Subscript[x, 0]^2))/2)(1+(v^2)))]))=(1/((Sqrt[(b/2)])(Subscript[x, 0])))(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[(1-(v^2)) ((v^2)-((a-((b (Subscript[x, 0]^2))/2))/((b (Subscript[x, 0]^2))/2)))])).}})}})

The inverse function to the elliptic integral x↦(∫_0^x)((dv)/(Sqrt[(1-(v^2)) ((v^2)-(k^2))])) is the Jacobi elliptic function dn(x|1-(k^2)). From our previous calculation, it follows that (x/(Subscript[x, 0]))=dn((Sqrt[(b/2)])(Subscript[x, 0])t|1-((a-((b (Subscript[x, 0]^2))/2))/((b (Subscript[x, 0]^2))/2))), i.e., x=(Subscript[x, 0])dn((Sqrt[(b/2)])(Subscript[x, 0])t|2-((2a)/(b (Subscript[x, 0]^2)))). 
The period of the function u↦dn(u|m) is 2K(m), where K is the complete elliptic integral of the first kind. Also, dn is even about the point K(m). Hence, our formula for x is correct not only in the upper half-plane, but remains valid for all times.

We now use the explicit formula for x to plot graphs of several solutions inside the separatrix. Note how their periods depend on the choice of the initial value.

With[{a=0.5,b=0.5},Plot[Evaluate[Table[x0 JacobiDN[(Sqrt[(b/2)]) x0 t,2-((2 a)/(b (x0^2)))],{x0,{1,1.1,1.2,1.3}}]],{t,0,10}]]

The period of the solution starting at (Subscript[x, 0]) is 2K(2-((2a)/(b (Subscript[x, 0]^2))))/((Sqrt[(b/2)])(Subscript[x, 0])). The next plot shows the dependence of the period on (Subscript[x, 0]) for all orbits inside the separatrix.

With[{a=0.5,b=0.5},
Plot[2EllipticK[2-((2 a)/(b (x0^2)))]/(Sqrt[(b/2)]),{x0,(Sqrt[(a/b)]),(Sqrt[((2 a)/b)])},AxesLabel->{"(Subscript[x, 0])","period"}]]

It seems that the period goes to infinity as we move towards the separatrix, but has a positive lower bound when we approach the equilibrium. The next calculation confirms these hypotheses.

Assuming[{a>0,b>0},Limit[2EllipticK[2-((2 a)/(b (x0^2)))]/(Sqrt[(b/2)]),x0->(Sqrt[((2a)/b)])]]

Assuming[{a>0,b>0},Limit[2EllipticK[2-((2 a)/(b (x0^2)))]/(Sqrt[(b/2)]),x0->(Sqrt[(a/b)])]]

In a similar way, we can analyze orbits outside the separatrix. We now look for the solution starting at ((Subscript[x, 0]),0), where (Subscript[x, 0])>(Sqrt[((2a)/b)]). Proceeding as before, we get

({{({{t=(∫_0^x)((du)/(Sqrt[-a (Subscript[x, 0]^2)+((b (Subscript[x, 0]^4))/2)+a (u^2)-((b (u^4))/2)]))=(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[(1-(v^2)) (-a+((b (Subscript[x, 0]^2))/2)(1+(v^2)))]))=(1/(Sqrt[b (Subscript[x, 0]^2)-a]))(∫_0^(x/(Subscript[x, 0])))((dv)/(Sqrt[(1-(v^2)) (((-a+((b (Subscript[x, 0]^2))/2))/(b (Subscript[x, 0]^2)-a))+(((b (Subscript[x, 0]^2))/2)/(b (Subscript[x, 0]^2)-a))(v^2))])).}})}})

This time, we use the fact that the inverse function to the elliptic integral x↦(∫_0^x)((dv)/(Sqrt[(1-(v^2)) (1-(k^2)+(k^2)(v^2))])) is the Jacobi elliptic function cn(x|(k^2)). This leads to the result x=(Subscript[x, 0])cn((Sqrt[b (Subscript[x, 0]^2)-a])t|((b (Subscript[x, 0]^2))/(2b (Subscript[x, 0]^2)-2a))). Since the period of the function u↦cn(u|m) is 4K(m), the period of x is 4K(((b (Subscript[x, 0]^2))/(2b (Subscript[x, 0]^2)-2a)))/(Sqrt[b (Subscript[x, 0]^2)-a]). Here is a plot of the period as a function of (Subscript[x, 0]):

With[{a=0.5,b=0.5},
Plot[((4EllipticK[((b (x0^2))/(2 b (x0^2)-2 a))])/(Sqrt[b (x0^2)-a])),{x0,(Sqrt[((2 a)/b)]),8},AxesLabel->{(Subscript[x, 0]),"period"},AxesOrigin->{(Sqrt[((2 a)/b)]),0}]]

The next calculation shows that the period approaches infinity as we move towards the separatrix, and decreases to zero when (Subscript[x, 0]) goes to infinity.

Assuming[{a>0,b>0},Limit[((4 EllipticK[((b (x0^2))/(2 b (x0^2)-2 a))])/(Sqrt[b (x0^2)-a])),x0->(Sqrt[((2a)/b)])]]

Assuming[{a>0,b>0},Limit[((4 EllipticK[((b (x0^2))/(2 b (x0^2)-2 a))])/(Sqrt[b (x0^2)-a])),x0->∞]]

We next consider a case that has an additional damping term (c>0), but no forcing (d=0).

duffingStable=DuffingEquation[{0.5,0.5, 0.015, 0, 1}, x[t]]

SecondOrderPlot[duffingStable,x[t],{t,0,50},{x,-1.7,1.7},{x',-0.65,0.65},InitialValues->{{1,0.54},{1,0.56}},FrameLabel->{"x","(x^(′))"},NullclineShading->True,Nullclines->True, PlotStyle->{{},{Red}},ShowEquilibria->True,RotateLabel->False]

Note how the nearby initial values lead to orbits that spiral toward different equilibria. We can use ProjectionPlot3D to get a 3-dimensional view of one of these orbits, together with its x, y, and x-and-y projections.

SecondOrderPlot[duffingStable,x[t],{t,0,60},{x,-2,2},{x',-1,1},InitialValues->{1,0.54},ParametricPlotFunction->ProjectionPlot3D,Axes->True,AxesLabel->{"x","(x^(′))","t"}]

There is no indication of instability in this orbit, but we may as well discuss now a method that can be used to try to detect bad solutions. Sophisticated numerical algorithms try to estimate the error at each step; then they do additional work until this error is below a prescribed tolerance. A typical goal, and this is the default for Mathematica’s NDSolve, is to attempt to get this local error less than (10^(-8)) (in a relative sense, of course). This means that, after the very first step, one can expect an error of (10^(-8)) to be introduced. We can anticipate this error by perturbing the initial conditions by (10^(-8)). In short, a good general rule is to perturb the initial conditions by an amount on the order of the local truncation error. Let’s do it.

ϵ=(10^(-8));
SecondOrderPlot[duffingStable,x[t],{t,0,100},{x,-1.7,1.7},{x',-0.65,0.65},InitialValues->{{1,0.54},{1+ϵ,0.54+ϵ}},SaveLastPoint->True, PlotStyle->{{Blue, Thickness[0.015]}, {Yellow, Thickness[0.005]}}]

NumberForm[VDSLastPoint, 10]

The results differ by an amount roughly similar to the perturbation, which indicates the absence of a large error-propagation effect. And the plot shows that the two trajectories agree. A large discrepancy, on the other hand, would indicate that the error propagation was large, and so the results cannot be trusted.

It is natural to wonder which initial conditions lead to which equilibrium points. Such computations of “basins of attraction” can be very time-consuming when the basins are complicated, but they do lead to very striking images, since the basins of attraction often have fractal boundaries (see section 15.4). In the present case, the system is autonomous and one can begin to investigate the basins by simply drawing lots of orbits and checking where they go.

The code that follows uses the auxiliary function GetPts from the VisualDSolve package, which turns an interpolating function into the list of interpolating points. We can look at the last point to see which of the two equilibria has attracted the orbit, and therefore which color to use. The curve function that we now define does all this; that is, it takes an InterpolatingFunction and returns an appropriately colored Line object.

curve[pair_]:=(temp=Transpose[GetPts/@pair];{{If[temp[[-1,1]]>0,Blue,Yellow],Line[temp]}})

We now set some parameters so that we get five orbits. These parameters were chosen after some trial and error.

ymin=0.555;ymax=0.6175;
s=((ymax-ymin)/4);
{a,b}={1.6,0.73};
SecondOrderPlot[duffingStable,x[t],{t,0,400},{x,-a,a},{x',-b,b},InitialValues->Table[{1,y},{y,ymin,ymax,s}],SolutionName->"DuffingOrbits"];

Finally, we apply curve to the functions stored in DuffingOrbits to get an informative image.

Graphics[{{GrayLevel[0.3],Rectangle[{-a,-b},{a,b}]},
Thickness[0.0015], curve/@DuffingOrbits},Frame->True,PlotRange->{{-a,a},{-b,b}}]

This image shows that the two basins of attraction are intertwined, with smooth boundaries. But let’s not forget about the separatrices: orbits that go exactly to the saddle equilibrium at the origin. In fact, we can generate the two separatrices by considering two initial values, one just above and one just below the origin. Then a little graphics programming using the points on the separatrices will yield two polygon objects that capture the basins perfectly.

First we generate the separatrices. Note how the initial conditions specify that the initial value is at t=0, but the plotting intervals are taken to be negative. Thus we are reversing time to get these orbits.

SecondOrderPlot[duffingStable,x[t],{t,-66.7,0},{x,-1.7,1.7},{x',-0.83,0.83},InitialValues->{{0,{0,-0.001}},{0,{0,0.001},{-61.9,0}}},PlotStyle->{{Blue},{Yellow}},FastPlotting->True,WindowShade->GrayLevel[0.4],SolutionName->"DuffingSeparatrices"]

The programming that follows is a little intricate, since we must get the exact points to match up in order to get both Polygon objects. In particular, we need to locate which point in the blue interpolating function (in the preceding image) is near (-0.64,0.69). It turns out to be the 139th point.

{a,b}={1.66,0.85};
highData=Transpose[GetPts/@DuffingSeparatrices[[1]]];
lowData=Transpose[GetPts/@DuffingSeparatrices[[2]]];
Graphics[{{GrayLevel[0.4],Rectangle[{-a,-b},{a,b}]},
{Blue,Polygon[Join[Reverse[highData],lowData]]},{Yellow,Polygon[Join[Reverse[lowData],Drop[highData,138]]]},{PointSize[0.02],Red,Point/@{{-1,0},{1,0},{0,0}}}},Frame->True,FrameTicks->{{Automatic,None},{Automatic,None}},PlotRange->{{-a,a},{-b,b}}]

This image shows the basins quite nicely. Of course, generating such an image for a general autonomous system might well involve some graphics tricks peculiar to the situation. But the ideas above should work in many other autonomous cases as well. Nonautonomous cases are a different story!

15.2 The General Duffing Equation

We now study the impact of a forcing term. And we continually use the perturbation idea just discussed as a way of monitoring the likelihood that the results are accurate. The forced system — the forcing term is 0.3 cos t — is not autonomous, so we cannot look at nullclines or equilibrium points.

duffingForced = DuffingEquation[{1,1,0.15,0.3,1}, x[t]]

ϵ=(1/(10^8));
SecondOrderPlot[duffingForced,x[t],{t,0,6},{x,-2.2,2.2},{x',-2.2,2.2},InitialValues->{{0.6,1.3},{0.6+ϵ,1.3+ϵ}},FastPlotting->True,AspectRatio->(1/GoldenRatio)]

SecondOrderPlot[duffingForced,x[t],{t,0,40},{x,-2.2,2.2},{x',-2.2,2.2},InitialValues->{{0.6,1.3},{0.6+ϵ,1.3+ϵ}},FastPlotting->True,SaveLastPoint->True,AspectRatio->(1/GoldenRatio)]

In this case the solution seems very close to being periodic. To examine the result of the perturbation, we first define a distance function.

distance[{p_,q_}]:=Norm[p-q]
VDSLastPoint

distance[VDSLastPoint]

Again, the results are good in that the distance between the last points is roughly the same size as the perturbation. Of course, this is not a proof, but it is good evidence. What would happen if we perturbed by (10^(-3)), or (10^(-12))? We will leave such computations as exercises. But because the precision goal is only 8 digits, an 8-digit error will arise at the first step anyway, so a (10^(-12)) perturbation will yield a discrepancy similar to that for a (10^(-8)) perturbation. In section 5 we will discuss how and when to change the precision goal.

Now we take the same Duffing equation but move the initial value to (-1,1).

{x0,y0}={-1,1};
SecondOrderPlot[duffingForced,x[t],{t,0,100},{x,-1.7,1.7},{x',-1.7,1.7},InitialValues->{x0,y0},ShowInitialValues->True,FastPlotting->True,AspectRatio->(1/GoldenRatio)]

Here we see that the orbit has become quite disorganized. It looks like there are two regions that the orbit has an inclination to visit, and that it jumps from one to the other somewhat randomly. It will help to have a physical model to think about. Imagine a bead sliding on a wire in the shape illustrated. This is called a two-well model because, if a damping term is added, the bead will settle into one of the wells, the bump being an unstable equilibrium.

Plot[1-2 (x^2)+(x^4),{x,-1.5,1.5},AspectRatio->0.3,Axes->None,PlotRange->{-0.2,1.6},Epilog->{Thickness[0.01],Line[{{1,-0.08},{1,0.08}}],Line[{{-1,-0.08},{-1,0.08}}],Line[{{0,0.92},{0,1.08}}]}]

But if the whole apparatus is rigidly shaken from left to right with a periodic forcing function, then the bead can slide chaotically from well to well. The equations for this forced apparatus are quite similar (though not identical) to the Duffing equation (see [Str, §12.5]). Thus the images of section 15.2 can be thought of as showing the two wells, with the bead moving around in one of the wells for a while and then jumping to the other one.

Of course, when the motion is this complicated, we must wonder whether the solution is really correct. As usual, we try a perturbation of (10^(-8)).

ϵ=(10^(-8));
SecondOrderPlot[duffingForced,x[t],{t,0,100},{x,-1.7,1.7},{x',-1.5,1.5},InitialValues->{{x0,y0},{x0,y0}+ϵ},ShowInitialValues->True,SolutionName->"DuffingAndPerturbed",PlotStyle->{{},{Red}},SaveLastPoint->True, AspectRatio->1/GoldenRatio]

It is clear that the orbits have diverged somewhere near the end of the time interval.

distance[VDSLastPoint]

There are two approaches to getting more legible and more informative images. One is to plot only x against t. We could do this using SecondOrderPlot in a way that calls SystemSolutionPlot, but then the PlotStyle option only works to set styles for different plotting variables, rather than different initial values. So instead we use the fact that we have just stored the solution in DuffingAndPerturbed.

Plot[{DuffingAndPerturbed[[1,1]],DuffingAndPerturbed[[2,1]]},{t,0,100},PlotStyle->{{Blue,Thickness[0.008]},{Red}},Ticks->{Range[10,100,10],Automatic}]

This shows that visual divergence occurs near t=90, but that is a little misleading. Plotting the difference between the two solutions is better.

Plot[DuffingAndPerturbed[[1,1]]-DuffingAndPerturbed[[2,1]],{t,0,100}, FrameTicks->{{Automatic,{}},{Range[10,100,10],{62}}}, PlotRange->All, Frame->True, Axes->False, GridLines->{{62},{0}}]

Even clearer is to use logarithms. It is natural to use Log10[Abs[data]], but this is exactly what the RealExponent function does.

Plot[RealExponent[DuffingAndPerturbed[[1,1]]-DuffingAndPerturbed[[2,1]]],{t,0,100}, FrameTicks->{{Automatic,{}},{Range[10,100,10],{62}}}, PlotRange->All, Frame->True, Axes->False, GridLines->{{62},{-2}}]

This view is much clearer and shows that the small initial error gets slowly magnified to visible size after 62 seconds. We can also deduce from this the time limit for dependability. The evidence is that the original solution is valid (meaning: accumulated error is less than 1/100)  up to time 62. And one concludes that neither solution is accurate beyond t=62. This approach is not a proof, but it does provide evidence that no serious error has contaminated the basic solution up to time 60. We can check this with a high-precision plot; for that it is important that the initial values have infinite precision: (-1,1) is fine but (-1.,1.) would not be. We set the working precision of 30; this causes the precision goal to be set to 15.

duffingForcedRational = DuffingEquation[{1,1,15/100,3/10,1}, x[t]];

SecondOrderPlot[duffingForcedRational,x[t],{t,0,100},{x,-1.7,1.7},{x',-1.5,1.5},InitialValues->{x0,y0},SolutionName->"DuffingTrue",WorkingPrecision->30];

A comparison of the machine precision solution and the high-precision solution shows visible error at around 65.

Plot[{DuffingAndPerturbed[[1,1]],DuffingTrue[[1]]},{t,0,100},PlotStyle->{{Blue},{Red}}, PlotRange->All, Frame->True,FrameTicks->{{Automatic,{}},{Range[10,100,10],{}}}]

And a plot of the difference shows the real story: the machine precision solution has definite error after time 62, as predicted by the perturbation method.

Plot[{Abs[DuffingAndPerturbed[[1,1]]-DuffingTrue[[1]]]},{t,0,100},PlotStyle->{{},{Red}}, PlotRange->All, Frame->True,FrameTicks->{{Automatic,{}},{Range[10,100,10],{}}}, GridLines->{{62},{}}]

The error plot here shows that the true error from the high-precision solution mimics the error seen using only machine precision with a (10^(-8)) perturbation.

One can also look at the phase plane image, but restrict the plotting domain. Now that we know that t=60 is critical, we can look at the orbits from, say, 40 to 100. It is even more easy to read if we do it in steps of 2.5 seconds. The efficient way to do this is to first generate the solution out to 100 (which we have done already), and then use it to generate a few plots.

Several books have published orbits of the Duffing equation with t running up to 10000 or more. But the above computations show that, if machine precision is used (as is the case with almost all differential equations software), then the results are not correct.

15.3 Flying Duffing Circles

Returning to the more modest range of t between 0 and 10, we can gain some understanding of why the extreme sensitivity exists by making some movies. The nullcline shading idea does not apply to a nonautonomous system, but we can apply it at a given time instant. In other words, the underlying flow field changes as a function of t (because of the cos t term), but by fixing t we can look at the instantaneous nullclines. Similarly, there are no equilibrium points, but we can look at the “instantaneous equilibrium points” (which correspond to the two wells and the bump).

So here’s the plan: consider a small circle (radius (1/10)) of initial values around (-1,1), and make a movie that shows where these initial values go as t increases. In each frame we can include the instantaneous equilibria (the wells and bump) and nullcline shading.

An Insensitive Case

In the insensitive case — initial values (0.6,1.3) — we see that the thumping term does not greatly distract the flying circle. Using the two-well analogy, the ball seems to have enough energy to fly around the surface in a way that avoids the wells and the bump between, and it is unaffected by the forcing term.

We first compute several solutions and save them.

SecondOrderPlot[duffingForced,x[t],{t,0,4 π},{x,-2,2},{x',-2,2},InitialValues->Table[{0.6,1.3}+0.1 {Cos[s],Sin[s]},{s,0,((15 π)/8),(π/8)}],SolutionName->"SixteenDuffings1"];

Define a function to get all 16 points for a given time.

data1[tval_]:=SixteenDuffings1/. t->N[tval];

Define a frame of the desired movie.

frame1[time_]:=frame1[time]=SecondOrderPlot[x''[t]==-0.15 x'[t]+x[t]-((x[t])^3)+0.3 Cos[time],x[t],{x,-2,2},{x',-2,2},NullclineShading->True,Nullclines->True,ShowEquilibria->True,
EquilibriumPointStyle->{PointSize[0.03],White},Epilog->{PointSize[0.015],Point[data1[time]]}]

Compute all the frames in advance and generate the movie.

ListAnimate[Table[frame1[u],{u,0,4 π,(π/8)}],AnimationRepetitions->1]

A Sensitive Case

However, if we start at (-1,1), then the nearby initial conditions do interact with the bump. The animation shows how the bump breaks up the circle, causing some orbits to visit one well while the rest visit another. And keep in mind that this behavior is occurring over the modest time range of [0,4π].

SecondOrderPlot[duffingForced,x[t],{t,0,4 π},{x,-2,2},{x',-2,2},InitialValues->Table[{-1,1}+0.1 {Cos[s],Sin[s]},{s,0,((15 π)/8),(π/8)}],SolutionName->"SixteenDuffings2"];

data2[tval_]:=SixteenDuffings2/. t->N[tval];

frame2[time_]:=frame2[time]=SecondOrderPlot[x''[t]==-0.15 x'[t]+x[t]-((x[t])^3)+0.3 Cos[time],x[t],{x,-2,2},{x',-2,2},NullclineShading->True,Nullclines->True,ShowEquilibria->True,
EquilibriumPointStyle->{PointSize[0.03],White},Epilog->{PointSize[0.015],Point[data2[time]]}]

ListAnimate[Table[frame2[u],{u,0,4 π,(π/8)}],AnimationRepetitions->1]

The shape of the image of the initial circle in the sensitive case is intriguing in the way that it folds back upon itself. In fact, this fractal-like folded-over shape occurs again and again in orbits of the Duffing equation, and is related to something called the Smale horseshoe map (see [Str, chap. 12]). It is strongly suspected (though not proved) that a fractal-like strange attractor exists for the Duffing equation. Briefly, this refers to a very complicated set, with a fractal nature, that attracts even chaotic orbits.

Another way of gaining appreciation of the evolution of Duffing orbits is to look at a single image that shows what happens to a circle of initial conditions. In the code that follows we look at orbits out to t=7 and connect the resulting points (32 of them, from a circle of 32 initial points) with lines so as to see what happens to the initial circle. Note the folding that takes place. Note also that the area within the evolving circle seems to shrink. These behaviors are connected to the Lyapunov exponents, which can be estimated [KW].

SecondOrderPlot[DuffingEquation[{1,1,0.15,0.3,1},x[t]],x[t],{t,0,7},{x,-2,2},{x',-2,2},InitialValues->Prepend[Table[{-0.6,1.3}+0.1 {Cos[s],Sin[s]},{s,0,2 π-(π/32),(π/32)}],{-0.6,1.3}],SolutionName->"DuffingOrbits1"];

image1=Show[ParametricPlot[Evaluate[First[DuffingOrbits1]],{t,0,7},PlotStyle->{{Thickness[0.004],Red}}],
Graphics[Table[{EdgeForm[Thickness[0.006]],FaceForm[],Polygon[Rest[DuffingOrbits1]]},{t,0,7}]],PlotRange->{{-2,2},{-1.4,1.4}},Frame->True, Axes->None];

SecondOrderPlot[DuffingEquation[{1,1,0.15,0.3,1},x[t]],x[t],{t,0,7},{x,-2,2},{x',-2,2},InitialValues->Prepend[Table[{-1,1}+0.1 {Cos[s],Sin[s]},{s,0,2 π-(π/32),(π/32)}],{-1,1}],SolutionName->"DuffingOrbits2"];

image2=Show[ParametricPlot[Evaluate[First[DuffingOrbits2]],{t,0,7},PlotStyle->{{Thickness[0.004],Red}}],
Graphics[Table[{EdgeForm[Thickness[0.006]],FaceForm[],Polygon[Rest[DuffingOrbits2]]},{t,0,7}]],PlotRange->{{-2,2},{-1.4,1.4}},Frame->True, Axes->None];

GraphicsRow[{image1, image2}, ImageSize->500]

15.4 Forced Attraction

If we vary the parameters a little (reduce the forcing amplitude to 0.25 and increase the damping coefficient to 0.25), then the Duffing orbits are much more well behaved: they end up in one of the gravity wells near (-1,0) or (1,0). The following example is typical: after some oscillations, it falls into the left well and stays there.

duffingForced2= DuffingEquation[{1,1,0.25,0.25,1}, x[t]];

SecondOrderPlot[duffingForced2,x[t],{t,0,60},{x,-1.7,1.7},{x',-1.7,1.7},InitialValues->{1,1.1},FastPlotting->True, AspectRatio->1/GoldenRatio]

A plot of x vs. t shows the behavior a little more clearly.

SecondOrderPlot[duffingForced2,x[t],{t,0,150},PlotVariables->{{x}},InitialValues->{1,1.1},FastPlotting->True,PlotRange->All]

It is not an easy task to determine which well will attract a given initial value. Unlike the autonomous case, the phase plane is not simply partitioned by the orbits. So there is nothing for it but to look at thousands of initial conditions and simply keep track of where they end up. We proceed by defining AttractingWell to return ±1 depending on where the orbit gets to. Then we make an array of such values and use it to make a colored plot of the two basins. While some cases are stable for t out to 100 or more, other cases are not. Here we restrict ourselves to t≤8 π. Thus we must be aware that some orbits will not have settled down yet.

The straightforward way to do this is as follows. The order of the iterators is reversed from one would first do it, but this is how to get ListDensityPlot to interpret the values properly. One could double-check this by storing triples {a,b,sign}, but that uses more memory.

AttractingWell[a_,b_]:=Sign[NDSolveValue[{duffingForced2,x[0]==a,x'[0]==b},x[t],{t,0,8 π}]/. t->8. π];
res=40;

In the following we show how to use Monitor to monitor a computation. This one is very fast, but the technique is useful.

Monitor[data=Table[AttractingWell[a,b],{b,-2.5,2.5,5/res},{a,-2.5,2.5, 5/res}];,StringForm["``% done",Round[100 (b+2.5)/5]]]

And we show the points as follows, tweaking the FrameTicks setting to get the true values.

ListDensityPlot[data,ColorFunction->(If[#==1,Yellow,Blue]&),FrameTicks->{{({(1/5) (#+2.5) res,#}&)/@Range[-2,2],{}},{({(1/5) (#+2.5) res,#}&)/@Range[-2,2],{}}},InterpolationOrder->0]

The preceding approach will work with a resolution of 400, though it takes a while. A useful feature of Mathematica is the ability to set up a parallel computation, provided the platform supports multiple kernels. On a MacBook Pro with Intel 2.66 GHz CPU processor this works fine as follows. First launch four kernels.

LaunchKernels[4]

$KernelCount

Now we need to distribute the key definitions to the parallel kernels. And we set the resolution to 400.

res=400;
DistributeDefinitions[AttractingWell,duffingForced2,res];

And the ParallelTable command then gets the result in almost half the time of Table. We set the resolution to 400; the computation takes takes about two and a half minutes in absolute time. Monitor appears not to work in the parallel environment.

AbsoluteTiming[pData=ParallelTable[AttractingWell[a,b],{b,-2.5,2.5, 5/res},{a,-2.5,2.5,5/res}]]//First

A quick check that we have what we want.

{Length[Flatten[pData]],(401^2)}

And we can now look at the nicely resolved basin image. This too takes some time to generate.

ListDensityPlot[pData,ColorFunction->(If[#==1,Yellow,Blue]&),FrameTicks->Table[{({ (((#+2.5))/5) res,#}&)/@Range[-2,2],{}},{2}], ImageSize->400]

15.5 What To Do?

The sensitivity we have seen is not simply an anomaly of the numerical methods. It is in fact true that very small differences in the Duffing initial conditions can lead to very large differences in the orbits 100 seconds later. This is the often mentioned “butterfly effect”, which refers to the fact that a complicated system of differential equations, such as those that occur in weather prediction, might (note: “might”, not “will”) be sensitive to changes on the order of (10^(-10)) or smaller. There are several different ways to deal with this problem.

Physical Reality

If it is not the Duffing equation we are so interested in, but the underlying physical phenomenon that it models, then we can take the point of view that, because we cannot measure initial conditions to more than one part in (10^10), it is somewhat meaningless to ask for the true orbit.

Change the Question

The subsequent comments, which are related to those in the previous paragraph, were pointed out to us by Robert Corless, who has written on the usefulness of residuals [Cor].

Look at the residual (the ResidualPlot function, discussed in Chapter 2, can be used to generate its graph), denoted r(t). Then the computed solution (an interpolating function) is a completely accurate solution to the perturbed differential equation, in which the d cos t term is replaced by d cos t+r(t). If r(t) is less than, say, (10^(-10)), then we can’t be certain that it is not present as noise in the physically real situation. This gives us a way of changing the question so that we are absolutely certain that the computed answer is the true answer to the new question.

The Shadow Knows!

In some sensitive cases it has been proven that computed orbits are actually true orbits, but for slightly perturbed initial conditions. While this phenomenon, called shadowing, has not been strictly proved for the Duffing equation, it is suspected to be true. From this point of view, one is seeing a true orbit, but not for the precise initial values presented. One can in fact find explicit shadows using high-precision optimization; see [Wag, chap. 13].

Don’t Give Up

All this might seem hopeless to the idealistic mathematician who wants the true solution (for some infinite precision initial conditions such as (-1,1)). And in a sense it is hopeless! But in fact, Mathematica’s high-precision capabilities can give us some additional useful information, as in the example in section 15.2. We must emphasize that our point here is not so much the solving of the Duffing equation, but the discussion of general techniques that can be used on any equation to help determine if the numerical algorithm is doing the right thing.

Our basic suggestion is to increase the precision goals. If we work with, say, 20-digit working precision throughout and sharpen the local truncation error to, say, (10^(-10)), then we will maintain more precision at each step and get past the t=60 barrier that stops machine precision cold. True, we won’t get very much farther. But the idea is instructive to follow through. Instead of using VisualDSolve, it is simpler to use the built-in NDSolve command.

Here are two important points regarding increased precision:

When using higher precision it is crucial that no floating-point reals be used in the differential equation. That is why we will use rational coefficients such as 3/10 and 15/100.

When the working precision is set to a value (such as 24), then the precision goal, if it is set to Automatic, is automatically set to half that value (12).

Do[wpsol[wp]=NDSolveValue[{x''[t]-x[t]+((x[t])^3)+((15 x'[t])/100)==((3 Cos[t])/10),x[0]==-1,x'[0]==1},x[t],{t,0,100},WorkingPrecision->wp],{wp,16,32,4}];

We can now check the number of steps used. It is comforting to see this number grow at a reasonable rate.

Table[Length[GetPts[wpsol[i]]],{i,16,32,4}]

Now we look at a single perturbed solution, using working precision 20, for which the precision goal is 10.

HighPrecPerturbed=NDSolveValue[{x''[t]-x[t]+((x[t])^3)+(15/100)x'[t]==((3 )/10)Cos[t],x[0]==-1+(10^(-10)),x'[0]==1+(10^(-10))},x[t],{t,0,100},WorkingPrecision->20];

When we plot the two solutions for x(t) we see that they are in agreement out to about t=90, well past the t=60 barrier for machine precision.

Plot[{wpsol[20],HighPrecPerturbed},{t,0,100},PlotStyle->{{},{Red}}]

The error exists, but it does not become visually apparent until about t=90.

Plot[{wpsol[20]-HighPrecPerturbed},{t,0,100},PlotStyle->{{},{Red}}, PlotRange->All, Frame->True, Axes->False]

Using these ideas it is possible to get somewhat past t=100. However, some authors have computed Poincaré sections of the Duffing equation to t=10000 and beyond. Rob Knapp and Stan Wagon [KW] have used Lyapunov exponents to estimate that 800 digits of precision would be required if the computation were to pass the present tests for accuracy. This is not in itself an insurmountable obstacle. But they also estimated that (10^66) steps would be needed!

We conclude by presenting the orbit out to t=100, generated using 40 digits of working precision and a precision goal of 20. We believe that this orbit, unlike previously published ones, is an accurate representation to within 1 part in (10^(-6)) of the solution for the initial conditions (-1,1). We use a VisualDSolve approach, since working precision can be passed to its functions.

SecondOrderPlot[duffingForcedRational,x[t], {t,0,100},{x,-1.7,1.7},{x',-1.7,1.7},
 InitialValues->{-1,1},WorkingPrecision->40,PrecisionGoal->20, MaxSteps->20000]

15.6 Strange Attraction

Numerical solutions to the Duffing equation show the apparent existence of a strange attractor. One sees this by looking at a Poincaré section in steps of 2π. An important point is that it has not been proven that such an attractor exists for the ideal equation. It is possible (though unlikely) that the attractor is an artifact of the numerical work. Discussing this point Guckenheimer and Holmes [GH, p. 91] state that, “There is a substantial theoretical question as to whether this ‘strange attractor’ of the Duffing equation is an artifact of the noise and is absent in the ideal deterministic system.”

Here is a Poincaré section out to time 10000 for the starting values (-1,1).

PoincareSection[DuffingEquation[{1,1,0.15,0.3,1}, x[t]],{t,0,50000},{x,-1.6,1.6},{x',-0.85,1.3},{-1,1},TimeInterval->2 π,Axes->False,Frame->True,FrameTicks->{{Range[-1,1,1],{}},{Range[-1.5,1.5,.5],{}}},PointStyle->{Yellow,PointSize[0.01]},
WindowShade->Gray,MaxSteps->∞]

What is remarkable is that these points appear to be nearing a shape that does not change even as the initial position changes. The following demonstration allows one to move the locator defining the start position and watch the effect on the Poincaré section. For more on strange attractors, such as the precise definition, see [GH].

Manipulate[PoincareSection[DuffingEquation[{1,1,0.15,0.3,1},x[t]],{t,0,2 π n},{x,-1.6,1.6},{x',-0.85,1.3},initVal,TimeInterval->2 π,Axes->False,Frame->True,FrameTicks->{{Range[-1,1,1],{}},{Range[-1.5,1.5,0.5],{}}},PointStyle->{Yellow,PointSize[size]},WindowShade->Gray,MaxSteps->∞],{{initVal,{-1,1}},Locator},{{n,500,"number of points in Poincaré section"},10,1000},{{size,0.015,"point size"},0.001,0.03},TrackedSymbols->{size,n,initVal},Initialization:>Needs["VisualDSolve`VisualDSolve`"]]

The next demo shows how the Poincaré section is formed from the solutions, which are plotted in 3-space so as to include time, but reduced mod 2π. One can change the initial conditions by the sliders and see how the points at multiples of 2π are attracted to the attractor.

Manipulate[sol=NDSolveValue[{DuffingEquation[{1,1,0.15,0.3,1},x[t]],x[0]==a,x'[0]==b},x[t],{t,0,tmaxForDE},MaxSteps->∞];
Show[
ParametricPlot3D[Evaluate[{Mod[t,2π],sol,(∂_t)sol}],{t,0,tmaxForDE},PlotStyle->{Purple,Thickness[0.004]}, Exclusions-> Table[2 π i, {i, 0, 100}]],
Graphics3D[{{PointSize[0.002],Red,Point[ps]},
{Blue,PointSize[0.02],Point[{0, a,b}],Point[Table[Evaluate[{ 2π,sol,(∂_t)sol}/.t->2πi],{i,(tmaxForDE/(2 π))}]]},
{Gray,Line[Table[{{0,j,-1.2}, {0,j,1.2}}, {j,-1.5,1.5,.5}]],
Line[Table[{{0,-1.4,j}, {0,1.7,j}}, {j,-1,1,.25}]]}}],ImageSize->400,PlotRange->{{0.001, 6.28},{-1.4,1.7},{-0.7,1.2}},BoxRatios->{6.4, 3,2.6},
Ticks->{Range[0,2π,π], Range[-1,1], Range[-.5,1,.5]},AxesLabel->{"t (mod 2 π)", "x", "(x^(′))"},ViewPoint->3{-1.6, -2, .4}, AxesEdge->{{-1,-1},{-1,-1},{-1,1}}],
{{a,0, "x(0)"}, -1,1},{{b,0, "x'(0)"},-1,1},{{tmaxForDE,15,"max time for solution"},0.001,150},TrackedSymbols->{a,b,tmaxForDE},TrackedSymbols->Manipulate, Initialization :> (Needs["VisualDSolve`VisualDSolve`"];PoincareSection[DuffingEquation[{1,1,0.15,0.3,1},x[t]],{t,0,2 π 3000},{x,-1.5,1.5},{x',-1.3,1.3},{-1,1},TimeInterval->2 π,SavePoints->True,MaxSteps->∞];ps =Prepend[#, 2π]& /@  VDSPoints)]