Chapter 5
SecondOrderPlot

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.

5.0  Overview

This chapter discusses the use of SecondOrderPlot to study the solutions to one or more second-order equations. This is really the culmination of the whole package since this function makes use of many of the other comprehensive functions such as PhasePlot and SystemSolutionPlot, and also because this class of equation encompasses so many fascinating examples.
FUNCTIONS DISCUSSED.  ToSystem, PhasePlot, PoincareSection, SecondOrderPlot.

5.1  A Single Second-Order Equation

SecondOrderPlot accepts a second-order equation in the form x''[t]==x[t]+t, or x''[t]-x[t]==t. One important distinction is that the second argument must consist of the dependent function (or list of functions in the case of a system). The simplest usage is as follows. No plotting domain is given and the default is to provide graphs of x and x' in a single image.

Get["VisualDSolve`"]

SecondOrderPlot[x''[t]==-x[t],x[t],{t,0,2 π},InitialValues->{0,1}]

Note that the initial values must always appear in the order (x, x'), or (x, y, …, x', y',…) in the case of a system, and the t-value at which the initial values are set is the first t-value in the t-iterator. This is all a little less flexible than the initial value possibilities for PhasePlot, but the user can always transform the equation(s) to a system and use PhasePlot (an example where this is useful occurs in the helium balloon section that follows).

If separate images are desired, that is accomplished through the PlotVariables option, which is passed through to SystemSolutionPlot (we can use the PlotLabels option to identify the graphs).

GraphicsRow[SecondOrderPlot[x''[t]==-x[t],x[t],{t,0,2 π},
InitialValues->{0,1},PlotVariables->{{x},{x'}},PlotLabels->{"x vs. t","x' vs. t"}]]

If only one graph is wanted, then that is obtained by feeding only a single plotting variable as an argument; it can be either x or x'. In these cases the plotting range would be controlled by the PlotRange option (unlike VisualDSolve, which used an iterator, {x,a,b}, to control the window height).

SecondOrderPlot[x''[t]==-x[t],x[t],{t,0,2 π},
PlotVariables->{{x'}},InitialValues->{0,1},FrameLabel->{"Time","Speed"}]

Note that SecondOrderPlot can interpret, say, x' properly in its role as, for example, a plotting variable. When the plotting variables are given in the form of iterators, as opposed to just variables, then they form the basis of a phase plot. Often for a second-order equation one will want to see x vs. x'. The initial values may be set to any of the forms suitable for PhasePlot.

SecondOrderPlot[x''[t]==-x[t],x[t],{t,0,2 π},{x,-1.5,1.5},{x',-1.5,1.5},InitialValues->{{0,1},{0,0.5,{0,π}}},FrameLabel->{"x","x'"}, RotateLabel->False]

Here are some views of the phase space for a standard pendulum. We begin with the simple case where there is no friction. The undamped pendulum is a Hamiltonian system (see Chapter 10), and that allows us to figure out some initial values on the orbits that separate the closed orbits from the open ones. Thus the image includes orbits forming the separatrix. The equation is in the package as PendulumEquation.

pendulumEqn=PendulumEquation[{8, 1}, x[t]]

SecondOrderPlot[pendulumEqn,x[t],{t,0,4.2},{x,-10,10},{x',-10,10},Nullclines->True,NullclineShading->True, StayInWindow->True,ShowEquilibria->True,EquilibriumPointStyle->{PointSize[0.02],White},InitialValues->Join[{{0,-4 (Sqrt[2])},{0,4 (Sqrt[2])},{2 π,-4 (Sqrt[2])},{2 π,4 (Sqrt[2])},{-2 π,-4 (Sqrt[2])},{-2 π,4 (Sqrt[2])}},Table[{π,y},{y,-6.5,-1.5}],Table[{π,y},{y,1.5,6.5}],Table[{x,0},{x,-2 π,-π,(π/4)}],Table[{x,0},{x,0,π,(π/4)}],Table[{x,0},{x,2 π,3 π,(π/4)}]],PlotStyle->Black, AspectRatio->1/GoldenRatio,FrameTicks->{{Automatic,None},{Range[-3π,3π,π],None}}]

The next example shows the effect of adding friction. Now all the orbits eventually settle into one of the equilibria on the x-axis. Note how useful the shaded nullcline view is: the regions allow us to see at a glance where the equilibria are and the general shape of the orbits. This image also gives a sense of the shape of the basin of attraction of a single equilibrium by looking at the band of orbits that end up there. We use the ForcedDampedPendulumEquation in the package.

dampedPendulumEqn=ForcedDampedPendulumEquation[{1, 0, 10, 1}, x[t]]

SecondOrderPlot[dampedPendulumEqn,x[t],{t,0,15},{x,-16,16},{x',-16,16},Nullclines->True,NullclinePlotPoints->30,
NullclineShading->True,FastPlotting->True,InitialValues->Join[Table[{x,15},{x,-14,-7,0.8}],Table[{x,-15},{x,7.6,14,0.8}]],PlotStyle->Black,AspectRatio->1/GoldenRatio,FrameTicks->{{Automatic,None},{Range[-5π,5π,π],None}}]

Another tool that can be used with a single second order differential equation is PoincareSection. Here is an example that comes from the Duffing equation, which is studied in much more detail in Chapter 15.

pSec=PoincareSection[x''[t]==(1/4) Cos[2 t]-x[t]-((x[t])^3),{t,0,50 π},{x,2,2.15},{x',-1.15,1.15},{2.03284,0},TimeInterval->π,Ticks->{{2.02,2.06,2.1,2.14},Range[-1,1,0.5]},PointStyle->PointSize[0.02],AspectRatio->1/GoldenRatio]

We can examine the orbit, though the view we get gives no indication of where the heart-like shape lies.

orbit=SecondOrderPlot[x''[t]==(1/4) Cos[2 t]-x[t]-((x[t])^3),x[t],{t,0,10 π},{x,-2.2,2.2},{x',-4,4},InitialValues->{2.03284,0},FrameLabel->{"x","x'"}, AspectRatio->1/GoldenRatio,RotateLabel->False]

Combining the two images shows approximately where the points of the Poincaré section lie.

Show[orbit,pSec]

A closer look is called for.

Show[pSec,orbit,Axes->Automatic,AxesOrigin->{1.9,0},Ticks->{{2,2.15},Range[-2,2]},PlotRange->{{1.9,2.15},{-2,2}}]

5.2  The Ups and Downs of a Helium Balloon

Helium balloons provide an interesting example because they will oscillate up and down to an equilibrium. Usually such a balloon, when released, travels high into the sky until it finds its equilibrium point or, much more likely, pops because of the reduced air pressure. But if you attach a moderately long length of string to the balloon, you can see the equilibrium occurring at eye height. Here we will show how to model the balloon-and-string with a second-order differential equation that has a few subtle points.

An interesting question is this: Imagine a helium balloon on a string in an environment with no air resistance. We can't use a vacuum, because there is no air! But we can imagine a very, very thin balloon. There will be an equilibrium height where the weight of the string balances the buoyant force due to the helium. But there is no air resistance. So will the balloon oscillate around the equilibrium and converge to it, or will it oscillate around the equilibrium but never get close to it? We will answer that question here.

Suppose a helium-filled balloon attached to a long string with linear density ρ grams/cm. When the balloon descends because of the weight of the string, the system will lose mass as some string falls to rest on the ground; this changing-mass aspect is what makes the problem interesting. Surprisingly, the problem was first studied by J. Buquoy in 1814 [SP]. We denote gravitational acceleration (980 cm/(sec^2)) by g. We use y(t) for the height of the balloon's low point off the ground at time t and v(t) for the velocity of the balloon (same as y'(t)).

Let M denote the mass of the system, which is the variable mass of the part of the string that is off the ground (ρ·y) plus the mass of the balloon and the helium ((Subscript[M, b])=3 grams). Then the usual force of gravity is: (Subscript[F, g])=-g·M - g·ρ·y(t). Let (Subscript[F, h]) denote the
constant upward force from the buoyancy of the helium in the balloon. Assume the balloon displaces 13 liters. Using the specific gravity of air at room temperature (1.161 grams/liter) and of helium (0.16 grams/liter), we can approximate the
upward force due to the helium’s buoyancy as

(Subscript[F, b])=g·(mass air displaced-mass helium)=(13·1.161-13·0.16)g=13ggrams·cm/(sec^2).

Let (Subscript[F, d]) be the drag, or air resistance, on the balloon, which is a function of the velocity of the balloon. Following [YSCD], we take (Subscript[F, d]) to be -sign(y'(t))(Subscript[c, d])((y'(t))^2) where (Subscript[c, d]) is the coefficient of drag. This force is always opposite to the direction of motion and that is what the negative of the sign function accomplishes. To start we will ignore drag, using (Subscript[c, d])=0. Newton’s second law of motion — the time derivative of momentum is equal to the sum of the forces — will yield a second-order differential equation. But something atypical arises in this problem as there is an important difference between the ascent and descent phases. Newton’s second law says

(*)      total force = F = (Subscript[F, g])+(Subscript[F, b])+(Subscript[F, d])= change in momentum =(d/(dt))(M·v)=M v'+v M'.

For a system in which the mass does not change, M'=0 and we have the familiar F=M a. But for the balloon on a string, the mass is changing on descent as some of the string strikes the ground, going instantaneously from speed v (t) to speed 0. And on ascent the same happens in reverse as some string is pulled up from rest to assume speed v (t). We will assume that these changes are instantaneous. Such an instantaneous change in acceleration causes infinite jerk ("jerk" is the third derivative of position). That cannot happen in reality, but it is a simplifying assumption here that leads to no difficulty. In some situations infinite jerk is a disaster (e.g., if a parachute opens too quicky the jerk will kill the human; see Chapter 7).

The term v M' in (*) needs careful study. On ascent, this gain in mass results from an upward pull by the buoyant force of the helium; because it involves new mass, it is an interaction that will decelerate the system. This can be seen more clearly by isolating y'': y''=(F/M)-((v )/M); here v is positive and F is decreasing (because of the increasing mass) and so the system decelerates. Thus the force law in (*) applies as written. But on descent the loss of mass and momentum due to the bottom of the moving string crashing into the ground is an interaction between the ground and the bottom part of the string. It does not involve the system above the bottom part of the string and so has no effect on the acceleration of the system. Therefore it must be ignored and removed from (*). Another way to achieve the same result is to say that, on descent, there is a new force provided by the ground, one that cancels the mass-changing term: (Subscript[F, ground])=v M'=(ρ v)v=ρ (v^2). This is different than the normal force of the ground, which supports the part of the string that is at rest; that force and that string can be ignored. For more details and some variations in the approach see [SP], [YSCD], [KRR], [DKS].

So we will have a piecewise defined differential equation, but that is no problem. We first define all forces; it is convenient to use Fground and define it by an If-statement to account for the difference in ascent and descent.

g=980;ρ=0.1;Mballoon=3;
Fgravity=-g(Mballoon+ρ y[t]);
Fbuoyant=13g;
Fdrag[cd_]:=-Sign[y'[t]]cd ((y'[t])^2);
Fground=If[y'[t]<0,ρ ((y'[t])^2),0];

We can now set up the differential equation using (*) and the product of mass and velocity, ((Subscript[M, b]) + ρ y)y', to get the left-hand side (change in momentum). We are here using Mathematica to differentiate the abstract expression representing momentum with respect to t.

balloonEqn[cd_]:=D[(Mballoon+ρ y[t]) y'[t],t]==Fbuoyant+Fgravity+Fground+Fdrag[cd];
balloonEqn[0]

SecondOrderPlot[balloonEqn[0],y[t],{t,0,6},{y,0,175},{y',-160,260},InitialValues->{0,0},AspectRatio->1,ShowInitialValues->True,FrameLabel->{"y","y'"}, RotateLabel->False]

We would like to rotate this so that the vertical axis represents the height of the balloon. However, SecondOrderPlot is finicky in that it wants its initial values to have the function values first, followed by the derivative values; and the iterators have to respect this order as well. So to do what we want we resort to transforming the second-order equation to a system and using PhasePlot. This is not too onerous because we can use ToSystem (see section 2.5) to do the transformation for us. And we use the AuxiliaryVariables option to force v to be the auxiliary variable that represents y'.

balloonSystem[cd_]:=ToSystem[balloonEqn[cd],y[t],t,AuxiliaryVariables->v];
balloonSystem[cd]

Now we can feed balloonSystem to PhasePlot with v, the velocity, coming ahead of y at all steps, including the equation and the initial values. We start the balloon at ground level.

PhasePlot[Reverse[balloonSystem[0]],{v[t],y[t]},{t,0,15},{v,-250,250},{y,0,175},InitialValues->{0,0},AspectRatio->1,Nullclines->True,NullclineShading->True,FlowField->True,FieldColor->GrayLevel,FrameLabel->{"y'","y"},FieldLogScale->Automatic,PlotStyle->{Black,Thickness[0.01]},RotateLabel->False]

We can now add air resistance on the model by specifying (Subscript[c, d]) to be positive. It is not hard to set up a demonstration that allows several parameters to be controlled, though here we control only the coefficient of drag and whether or not to include the force of the ground.

Manipulate[PhasePlot[
{(Mballoon+ρ y[t])v'[t]==13g-Sign[v[t]]cd ((v[t])^2)+If[v[t]<0&&groundQ,ρ ((v[t])^2),0]-ρ ((v[t])^2)-g(Mballoon+ρ y[t]),y'[t]==v[t]},{v[t],y[t]},{t,0,15},{v,-250,250},{y,0,175},InitialValues->{0,0},PlotPoints->100,AspectRatio->1,NullclineShading->True,Nullclines->True,FrameLabel->{"y'","y"},PlotStyle->{Black,Thickness[0.01]},RotateLabel->False],{{cd,0.05,"coefficient of drag"},0,0.1},
{{groundQ,True,"include force of ground"},{True,False}},TrackedSymbols:>{cd,groundQ},
Initialization:>(Needs["VisualDSolve`VisualDSolve`"];ρ=0.1;Mballoon=3;g=980)]

It is interesting to look at the energy loss more closely. Energy is lost as the string suddenly decelerates and strikes the ground; the crash landing causes heat to form among the molecules in the string. Similarly, the violent upward jerk from rest to a nonzero speed causes a disturbance that leads to an energy loss as heat. The kinetic energy is just (1/2)m (v^2). For the potential energy we can take the highest point of the string as the base point and use my. But the buoyant force of the helium provides a negative potential energy; that is, this force can be viewed as a gravitational force that acts in the opposite of the usual direction. For the string, the potential energy is computed as if all the string is at its center of mass, which is y/2.

First generate a solution and store it in ySol and ySolDer.

PhasePlot[balloonSystem[0],{y[t],v[t]},{t,0,15},{y,0,175},{v,-250,250},InitialValues->{0,0}, SaveSolution->True];
{ySol,ySolDer}=VDSSolution

Now define the two energy functions and the total energy.

kineticE=(1/2) (ρ ySol+Mballoon) (ySolDer^2);
potentialE=Mballoon g ySol-Fbuoyant ySol+(1/2) ρ g (ySol^2)  ;
totalEnergy = kineticE+potentialE;

For a more useful graph it is useful to make the times when the balloon is at a local maximum or minimum.

extrema=FindAllRoots[ySolDer, {t,0.1, 14}]//Quiet
equilibrium = y[t]/. Solve[(Mballoon+ρ y[t]) ==Fbuoyant+Fgravity/. g-> 980 /. ρ-> .1, y[t]][[1]]

Here is a plot of the balloon’s height.

Plot[ySol, {t,0,8}, Frame->True, GridLines->{Table[{x,{Red,AbsoluteThickness[1]}}, {x,extrema}],{{equilibrium, AbsoluteThickness[1]}}}, FrameLabel->{"time", "y"}, RotateLabel->False]

And here is the total energy.

Plot[totalEnergy,{t,0,8},Frame->True,GridLines->{Table[{x,{Red,AbsoluteThickness[1]}}, {x,extrema}],{}},PlotRange->All, FrameLabel->{"time", "total\nenergy"}, RotateLabel->False]

So we see that energy is lost throughout the oscillation. In [YSCD] it is shown that the rate of energy loss is the same in both the ascent and descent cases: that rate equals to ((1/2)ρ+ (Subscript[C, d]))((Abs[v])^3), where (Subscript[C, d]) is a scalar multiple of (Subscript[c, d]). We can easily verify that result symbolically.

Clear[ρ, g]
Fgravity=-g(Mballoon+ρ y[t]);
Fbuoyant=13g;
Fdrag[cd_]:=-Sign[y'[t]]cd ((y'[t])^2);
Fground=If[y'[t]<0,ρ ((y'[t])^2),0];

Here is the total energy in symbolic form.

totalEnergySymbolic=(1/2) (ρ y[t] ((y'[t])^2)+ Mballoon ((y'[t])^2))+Mballoon g y[t]-Fbuoyant y[t]+(1/2) ρ g ((y[t])^2)

And here we look at the derivative of the energy, eliminate y'' by using the basic differential equation.

solAscent=Solve[(∂_t)((Mballoon+ρ y[t])y'[t])==Fbuoyant+Fgravity+Fdrag[cd],y''[t]][[1]];
solDescent=Solve[(∂_t)((Mballoon+ρ y[t])y'[t])==Fbuoyant+Fgravity+Fdrag[cd]+ρ ((y'[t])^2),y''[t]][[1]];
Collect[Expand[{
Simplify[(∂_t)totalEnergySymbolic/. solAscent,y'[t]>0],
Simplify[(∂_t)totalEnergySymbolic/. solDescent,y'[t]<0]}],y'[t]]

These formulas are identical to formulas 2a and b in [YSCD]. It is surprising that the rate of energy loss is identical in the two types of motion, and is simply proportional to the cube of the speed.

5.3  Second-Order Systems

We illustrate second-order systems by looking at Newton’s laws applied to planetary motion. We can consider a planar cross-section, with an eccentricity of 0.5. Since X''=-X/|X((|)^3) governs the motion, where X represents the position vector (x,y), this is easily handled by SecondOrderPlot. Since we want to see an x-y plot, we use those two variables in the iterators. And, to illustrate Kepler’s law (a planet sweeps out equal area in equal time), we color the orbit in terms of its speed. The Epilog option is used to add the sun as a large point. Note that the initial values (four in this case) are understood to be in the order x, y, x', y'.

SecondOrderPlot can accept a differential equation in vector form; the package identifies the proper dimension from the initial conditions. This leads to elegant code for situations such as Newton's law. Note that if the unknown vector function is VEC then the two coordinate functions, needed for the window specification, are VEC[1] and VEC[2].

e=0.5;
SecondOrderPlot[
X''[t]==((-X[t])/((Norm[X[t]])^3)),X[t],{t,0,2 π},{X[1],-2,1},{X[2],-1.4,1.4},InitialValues->{1-e,0,0,(Sqrt[((1+e)/(1-e))])},ParametricPlotFunction->ColorParametricPlot,PlotStyle->Thickness[0.015],Speed->True,WindowShade->Gray,Epilog->{Yellow,PointSize[0.03],Point[{0,0}]}]

Recalling that red means fast and blue slow, we see that the planet hits maximum speed when it is closest to the sun, as it must do because of Kepler’s law, which states that the planet sweeps out equal areas in equal time.

We can illustrate Kepler’s law more convincingly by drawing lines from the sun at the origin to the planet at equal time intervals. To do this we store the solution in orbit (via the SolutionName option) and then do some graphics programming to get the lines we want. Since SaveSolution (which is automatically turned on by the use of SolutionName) causes only the plotted variables to be saved, substituting a Range for t in orbit will create a pair: some X-data and some Y-data. Transposing gives us the pairs we want, and we can Map the pure function that draws a line from the origin to a point onto this set of points.

orbitImage[e_]:=SecondOrderPlot[X''[t]==-X[t]((Norm[X[t]])^(-3)),
X[t],{t,0,2 π},{X[1],-2,1},{X[2],-1.4,1.4}, InitialValues->{1-e,0,0,(Sqrt[((1+e)/(1-e))])},PlotStyle->Thickness[0.008],SolutionName->"orbit",Epilog->{AbsolutePointSize[6],Point[{0,0}]}];

Show[orbitImage[0.5],Graphics[Table[(Line[{{0,0},tt}]), {tt,Transpose[orbit/.t->Range[0,2 π-(π/15),(π/15)]]}]]]

Another example of a second-order system arises from the flight of a thrown object, such as a ball. In the absence of air resistance, the model is quite simple. In the following example, the object is thrown from a height of 1.5 meters with initial velocity vector (20,20); this is the same as saying the initial velocity is 28.28 meters per second, at an angle of 45°. The result is the familiar parabolic trajectory. We present the example in two forms, the second using the vector approach.

g=9.8;
SecondOrderPlot[{x''[t]==0,y''[t]==-g},{x[t],y[t]},{t,0,10},{x,0,100},{y,0,25},InitialValues->{0,1,20,20}]

g=9.8;
SecondOrderPlot[X''[t]=={0,-g},X[t],{t,0,10},{X[1],0,100},{X[2],0,25},InitialValues->{0,1,20,20}]

We can now add air resistance to the model. A typical assumption is that the resistance force varies with the square of the speed. The velocity vector at time t is given by (x'(t),y'(t)), and the drag force is in the opposite direction (hence the minus signs). We use 0.002 as the resistance coefficient. Again, we give both the coordinate and vector forms.

k=0.002;g=9.8;
SecondOrderPlot[{x''[t]==-k (Sqrt[((x'[t])^2)+((y'[t])^2)]) x'[t],y''[t]==-g-k (Sqrt[((x'[t])^2)+((y'[t])^2)]) y'[t]},
{x[t],y[t]},{t,0,10},{x,0,100},{y,0,22},InitialValues->{0,1,20,20}]

k=0.002;g=9.8;
SecondOrderPlot[X''[t]=={0,-g}-k X'[t] Norm[X'[t]], X[t],{t,0,10},{X[1],0,100},{X[2],0,22},InitialValues->{0,1,20,20}]

A much more complicated  model of projectile motion is discussed in Chapter 13, where the flight of a discus is analyzed.