I'm using Mathematica to prepare some teaching materials based on the discussion of the damped, driven pendulum in Chapter 12 of Taylor's Classical Mechanics. This involves using NDSolveValue to solve the non-linear ODE $$ \ddot{\phi} + 2 \beta \dot{\phi} + \omega_0^2 \sin \phi = \gamma \omega_0^2 \cos \omega t. $$ We then want to create a phase plot, which is a parametric curve of the form $(\phi(t), \phi'(t))$ as a function of $t$. This is easy enough to implement:

DDPeq = D[ϕ[t], {t, 2}] + 2 β D[ϕ[t], t] + Subscript[ω, 0]^2 Sin[ϕ[t]] == γ Subscript[ω, 0]^2 Cos[ω t]
paramvals = {ω -> 2 π, Subscript[ω, 0] -> 3 π, β -> 3 π/8, γ -> 1.5}
tmin = 0;
tmax = 200;
ϕsol1[t_] = NDSolveValue[{DDPeq /. paramvals, ϕ[0] == -π/2, ϕ'[0] == 0}, ϕ[t], {t, tmin, tmax}];
ParametricPlot[{ϕsol1[t], ϕsol1'[t]}, {t, 0, tmax/4}, AspectRatio -> 1/GoldenRatio, AxesLabel -> {ϕ, \!\(\*OverscriptBox[\(ϕ\), \(.\)]\)}, ImageSize -> Full, PlotPoints -> 100]

enter image description here

However, the angle of the pendulum is really only defined modulo $2 \pi$; ideally, what we'd like to do is have Mathematica plot with the horizontal axis "wrapped around". We can use Mod to calculate this, and we can use Exclusions to ensure that the opposite edges of the graph aren't connected.

ParametricPlot[{Mod[ϕsol1[t], 2 π, -π], ϕsol1'[t]}, {t, tmax/10, tmax}, AspectRatio -> 1/GoldenRatio,  AxesLabel -> {ϕ, \!\(\*OverscriptBox[\(ϕ\), \(.\)]\)}, ImageSize -> Full,  
  PlotPoints -> 100, Exclusions -> Mod[ϕsol1[t], 2 π] == π, PlotStyle -> Thin]

enter image description here

The problem: Mathematica also seems to be "breaking" the plots when $\phi = 0$ (modulo $2 \pi$); you can see that the lines are broken whenever the curve crosses the vertical axis. Increasing PlotPoints makes this problem less noticeable, by narrowing the span of the curve that is omitted (at the edges as well as at the vertical axis); however, the gaps are still present.

Why is Mathematica excluding points that I didn't want it to exclude? Also, is there a better way to create this plot? (I'm still using Mathematica 10.4.1, though my college now has a site license for 11, so if there's a nifty feature in 11 that fixes this problem, I'm all ears.)

  • $\begingroup$ Apologies for the nasty-looking code. In a recent upgrade I lost my browser plugin that re-inserts special characters, etc. into code samples on StackExchange, and now I can't find it again. If anyone wants to prettify the code (or point me towards I can download that browser extension so I can do it myself), feel free. $\endgroup$ – Michael Seifert Apr 10 '18 at 16:33

Why not plot the phase space with the true topology?

ParametricPlot3D[{50 Cos[ϕsol1[t]], 50 Sin[ϕsol1[t]], ϕsol1'[t]}, {t, tmax/10, tmax},
                 PlotPoints -> 1000, PlotStyle -> Thin]



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.