Phase Plane Demo

Load the necessary graphics package for plotting vector fields.

<<"VectorFieldPlots`"

Make a plot of the vector field corresponding the right hand side of the 1st order system x' = v, v' = -x . Note the use of options to display axes with labels. (See Help on Options.)

dirPlot=VectorFieldPlot[{v,-x},{x,-2,2},{v,-2,2},Axes->True,AxesLabel->{"x","v"}]

Knowing that the analytic solution is x[t] = A*Cos[t] + B*Sin[t], choose initial conditions, e.g. x[0]=1, x'[0]=0 to obtain a particular solution x[t] = Cos[t] with corresponding velocity v[t] = -Sin[t]. Now plot the phase plane trajectory corresponding to this particular solution.

trajPlot1 = ParametricPlot[{Cos[t], -Sin[t]}, {t,0, 2*Pi}, AspectRatio -> 1, PlotStyle ->AbsoluteThickness[3]]

Now choose another set of initial conditions, e.g. x[0] = 1.5, x'[0] = 0 to obtain another particular solution, x[t]= 1.5*Cos[t] with corresponding velocity v[t] = -1.5*Sin[t] and plot another phase plane trajectory.

trajPlot2 = ParametricPlot[1.5*{Cos[t], -Sin[t]}, {t,0, 2*Pi}, AspectRatio -> 1, PlotStyle ->AbsoluteThickness[3]]

Show the trajectories together with the direction field.

Show[dirPlot, trajPlot1, trajPlot2]

Numerical Version

Now use the numerical differential equation solver NDSolve. Note that we need to work with the 1st order system if we want to have access to the velocity (so we can plot it).

solution=NDSolve[{x'[t]==v[t], v'[t]==-x[t], x[0]==1, v[0] == 0},{x,v},{t,0,2*Pi}]

Let's check if this solution looks reasonable by plotting with the exact solution Cos[t].

Plot[{x[t] /. solution, Cos[t]}, {t,0,2*Pi}]

The plots seem to lie visibly on top of one another. To further quantify the accuracy, plot the difference between the numeric and exact solutions.

Plot[(x[t] /. solution) - Cos[t], {t,0,2*Pi}]

Note that the error is noisy but on the scale of 10^-6 which is probably perfectly acceptable for our purposes.

Plot[(x[t] /. solution) - Cos[t], {t,0,2*Pi},
PlotRange -> {-1,1}, PlotStyle -> AbsoluteThickness[5] ]

Now, plot the phase plane trajectory corresponding to the numerical solution.

numericTrajPlot1 = ParametricPlot[{x[t], v[t]} /. solution, {t,0,2*Pi}, AspectRatio -> 1,
PlotStyle -> AbsoluteThickness[3] ]

Now show this numerical trajectory together with the direction field.

Show[dirPlot, numericTrajPlot1]

Note that the nonlinear version is more interesting, so be sure to plot trajectories with initial amplitudes greater than and less than 2.