Second Order 
Differential Equations

In this laboratory we explore a method that allows us to directly solve second-order differential equations. We examine first our dear old friend the simple harmonic oscillator and then consider the modifications necessary to investigate the nonlinear physical pendulum.

1. More on Simple Harmonic Motion.

We are becoming disgustingly familiar with the simple harmonic oscillator exemplified by the omnipresent mass on a spring. The differential equation that describes the system is 

                                            F = m (((d^2)x)/(dt^2)) = -kx       

where k is the spring constant, m is the mass, and x is the displacement from equilibrium. Apply the centered formula for the second derivative and show that

                                      (Subscript[x, n+1]) = 2(Subscript[x, n]) - (Subscript[x, n-1]) - (k/m)(h^2)(Subscript[x, n])          .
                                     
After generating the algorithm shown above read the remaining portions of Section 1 and execute the commands in each cell.

1.1 Setting things up.

The first code presented here defines the initial parameters and the starting point for the algorithm describing the motion of a mass on a spring. The initial position is  x(0) = 0.0 m, the initial speed is v(0) = 2.0 m/s, the mass is m=0.33 kg, and the spring constant is k=0.33 N/m.

(* Initial conditions and parameters
	
   initial position in meters *)
x0 = 0.10;
(* initial velocity in m/s *)
v0 = 2.0;
(* initial time in seconds *)
t0 = 0.0;
(* the mass in kg *)
mass = 0.33;
(* spring constant in N/m *)
kspring = 0.33;

(* set up the first two points. *)

(* step size *)
step = 0.1;
t1 = t0 + step;
x1 = x0 + v0*step;
		
(* limits of the iterations. since we already
				have y(t=0) and we have calculated y(t=step),
				then the first value in the table will be 
				for t=2*step. *)
tmin = 2*step;
tmax = 50.0;

1.2 Solving it.

What do you expect the plot of displacement versus time to look like? Load the first two values of the position(for t= (Subscript[t, 0])  and t = (Subscript[t, 0]) + h) and fill a table consisting of time and displacement (t,x). Plot the results. Does the plot fulfill your expectations.

(* A centered second derivative formula is used to
	  generate a iterative solution for the mass 
	  on a spring.
			
   first load the starting point values for the 
	  algorithm. *)
xminus = x0;   (* initial value of x *)
xmid      = x1;   (* second value of x *)

(* create a table tdispl of ordered (t,x). for 
		 each component the next value is calculated
	  first and then the variables are
	  incremented in preparation for the next
	  interation. *)
tdispl = Table[
	{t, 		
		xplus =  2*xmid - 
			                  xminus -  
			                 (kspring*step*step*xmid/mass);			
		xminus = xmid;
		xmid = xplus 
	},
	{t,tmin,tmax,step}
];

(* the table contains points that start 
		 at t=h and t=2h so we have add the 
			starting	point to the	table *)
tdispl = Prepend[tdispl,{t1,x1}];
tdispl = Prepend[tdispl,{t0,x0}];

Now make a plot of your results.

ListPlot[tdispl,PlotRange->Automatic,Joined->True,Frame->True,FrameLabel->{"t","x","Displacement of Mass on a Spring",""},PlotStyle->RGBColor[1,0,0]]

1.3 Phase Space, the Final Frontier.

What do you expect a plot of the velocity versus displacement to look like? Clear the table from the previous section (to save memory) and create a new one that includes (x,v) at each time step.

(* Initial conditions and parameters are the
	  same as before.
	
   initial position in meters *)
x0 = 0.10;
(* initial velocity in m/s *)
v0 = 2.0;
(* initial time in seconds *)
t0 = 0.0;
(* the mass in kg *)
mass = 0.33;
(* spring constant in N/m *)
kspring = 0.33;

(* set up the first two points. *)

(* step size *)
step = 0.1;
t1 = t0 + step;
x1 = x0 + v0*step;
x2 = 2*x1 - x0 - (kspring*step*step*x1/mass);
v1 = (x2 - x0)/(2*step);
		
(* limits of the iterations. since we already
				have y(t=0) and we have calculated y(t=step),
				then the first value in the table will be 
				for t=2*step. *)
tmin = 2*step;
tmax = 50.0;

(* get rid of the previous results for the table 
	  and proceed *)

Clear[tdispl];

(* A centered second derivative formula is used to
	  generate a iterative solution for the mass 
	  on a spring.
			
   first load the starting point values for the 
	  algorithm. *)
xmid      = x1;   (* second value of x *)
xplus   = x2;   (* third value of x *)

(* create a table of ordered (x,v). for 
		 each component the next value is calculated
	  first and then the variables are
	  incremented in preparation for the next
	  interation. *)
tdispl = Table[
	{xminus = xmid;
		xmid      = xplus;
		xplus   =  2*xmid - 
			                    xminus -  
			                   (kspring*step*step*xmid/mass);			
		xmid,
		vmid = (xplus - xminus)/(2*step)
	},
	{t,tmin,tmax,step}
];

(* the table contains points that start 
		 at t=h and t=2h so we have add the 
			starting	point to the	table *)
tdispl = Prepend[tdispl,{x1,v1}];
tdispl = Prepend[tdispl,{x0,v0}];

Now make a plot of your results.

ListPlot[tdispl,PlotRange->Automatic,Joined->True,Frame->True,FrameLabel->{"x","v","Phase Space of the Mass on a Spring",""},PlotStyle->RGBColor[1,0,0]]

1.4 Seeing the time dependence

If we now add the time component to the precious table of points what do you expect this three-dimensional trajectory in phase space to look like? As above, clear the table and view the phase space trajectory with the time included in 3-D. The plotting of a three-dimensional phase trajectory is not as simple in Mathematica as one might like so carefully look over the plotting commands at the end of this section.

(* Initial conditions and parameters are the
	  same as before.
	
   initial position in meters *)
x0 = 0.10;
(* initial velocity in m/s *)
v0 = 2.0;
(* initial time in seconds *)
t0 = 0.0;
(* the mass in kg *)
mass = 0.33;
(* spring constant in N/m *)
kspring = 0.33;

(* set up the first two points. *)

(* step size *)
step = 0.2;
t1 = t0 + step;
x1 = x0 + v0*step;
x2 = 2*x1 - x0 - (kspring*step*step*x1/mass);
v1 = (x2 - x0)/(2*step);
		
(* limits of the iterations. since we already
				have y(t=0) and we have calculated y(t=step),
				then the first value in the table will be 
				for t=2*step. *)
tmin = 2*step;
tmax = 20.0;

(* get rid of the previous results for the table 
	  and proceed *)

Clear[tdispl];

(* A centered second derivative formula is used to
	  generate a iterative solution for the mass 
	  on a spring.
			
   first load the starting point values for the 
	  algorithm. *)
xmid      = x1;   (* second value of x *)
xplus   = x2;   (* third value of x *)

(* create a table tdispl of ordered (t,x,v). for 
		 each component the next value is calculated
	  first and then the variables are
	  incremented in preparation for the next
	  interation. *)
tdispl = Table[
	{t,
		xminus = xmid;
		xmid      = xplus;
		xplus   =  2*xmid - 
			                    xminus -  
			                   (kspring*step*step*xmid/mass);			
		xmid,
		vmid = (xplus - xminus)/(2*step)
	},
	{t,tmin,tmax,step}
];

(* the table contains points that start 
		 at t=h and t=2h so we have add the 
			starting	point to the	table *)
tdispl = Prepend[tdispl,{t1, x1,v1}];
tdispl = Prepend[tdispl,{t0, x0,v0}];

Now make a plot of your results. There are a couple of ways to make 3-dimensional plots in Mathematica. We display two of them here. The simplest command is the ListPointPlot3D command, but it will not join adjacent points to make a pretty curve so the second command is presented if you want to do that.

ListPointPlot3D[tdispl,
	           BoxRatios->{1.5,1,1},
	           Axes->True,
	           AxesLabel->{"Time","x(m)","v(m/s)"},
	           PlotLabel->"HO Phase Trajectory",
	           ViewPoint->{2.002, -2.494, 1.104}
          ]

Show[ Graphics3D[ { RGBColor[1,0,0],
			                                       Thickness[0.01], 
			                                        Line[tdispl ]
                                         }],
	           BoxRatios->{1.5,1,1},
	           Axes->True,
	           AxesLabel->{"Time","x(m)","v(m/s)"},
	           PlotLabel->"HO Phase Trajectory",
	           ViewPoint->{2.002, -2.494, 1.104}
          ]

2. Beyond General Physics, the Non-linear Physical Pendulum.

In our study of second order differential equations we have encountered the physical pendulum. The 'model' that describes the motion of the physical pendulum is

                                                       I(((d^2)θ)/(dt^2)) =  -mgh sinθ

where I is the moment of inertia, θ is the angle with respect to the vertical, m is the mass, and h is the distance from the pivot point to the center of mass (not the time step size).

a. Use the centered, second-derivative formula to generate an algorithm for solving the differential equation above.

b. Consider a physical pendulum consisting of a thin rod of length L=0.6 m and pivoted about one end. It is held at an angle of (10^o) and then released from rest. Using the examples in Section 1 as a guide, generate a plot of angular displacement  θ versus time. Use a stepsize of about 0.01 s.

c. Compare your plot of angular displacement versus time with what you expect for the linear oscillator (in other words, find the particular solution for this set of initial conditions for the simple harmonic oscillator). How do the plots differ? How are they similar? What happens if you increase the initial angle to (12^o)? What happens if you increase the inital angle to (75^o).

d. Again using the examples above generate a phase-space plot of velocity versus angular displacement for an inital angle of (10^o) and  (75^o). Compare your result with what you expect for the linear oscillator. How do the plots differ? How are they similar? What happens when you change the initial angle of release by a small amount? Is this system chaotic? Explain.