Coupled, First-Order 
Differential Equations

In this laboratory we explore a method that allows us to develop a complete solution to projectile motion (the Lieutenant Chisov problem) among other examples of second order differential equations.

1. An Example.

As we've mentioned before oscillatory systems are found throughout physics (bridges, stars, atomic nuclei). A classic example is the mass on a spring that is described with the following equation derived from Newton's Second Law.

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

The quantity x is the displacement from the equilibrium point of the spring and k is the spring constant. The equation above can be written as a set of coupled equations

                                              m(dv/dt) = -kx
                                              
                                                     v = (dx/dt)
  
where v is the velocity. Applying the centered derivative formula to both expressions generates two finite difference equations that can be solved iteratively.

                                            (Subscript[v, n+1]) = (Subscript[v, n-1]) - ((2hk)/m)(Subscript[x, n])
                                            
                                            (Subscript[x, n+1]) = (Subscript[x, n-1]) + 2(Subscript[hv, n])
                                              
The quantity h is the stepsize. The codes in the following sections present an algorithm for solving these finite difference equations. You might want to experiment with different initial conditions and stepsizes to explore the behavior of the algorithm.

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 15.0 m, the initial speed is -13.6 m/s, the mass is 0.33kg, and the spring constant k is 400 N/m.

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

(* set up the first two points. *)

(* step size *)
step = 0.01;
t1 = t0 + step;
x1 = x0 + v0*step;
v1 = v0 - ( step*kspring*x0/mass);
		
(* 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 = 1.0;

1.2 Getting the velocity.

Load the first two values of the velocity and the position(for t=t0 and t=t0+h) and fill a table consisting of time and velocity (t,v). Plot the results.

(* A centered 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 *)
vminus = v0;   (* initial value of v *)
xmid     = x1;   (* second value of x *)
vmid     = v1;    (* second value of v *)

(* create a table tvel of ordered (t,v). for 
		 each component the next value is calculated
	  first and then the variables are
	  incremented in preparation for the next
	  interation. *)
tvel = Table[
	{t, 		
		vplus = vminus - (2*step*kspring*xmid/mass);
		xplus =  xminus + (2*step*vmid);			
		xminus = xmid;
		xmid = xplus ;
		vminus = vmid;
		vmid  = vplus
	},
	{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 *)
tvel = Prepend[tvel,{t1,v1}];
tvel = Prepend[tvel,{t0,v0}];

ListPlot[tvel,PlotRange->Automatic,Joined->True,Frame->True,FrameLabel->{"t(s)","v(m/s)","Velocity of Mass on a Spring",""},PlotStyle->RGBColor[1,0,0]]

1.3 Getting the position as a function of time.

Reload the first two values of the velocity (for t=t0 and t=t0+h) and the 
position and fill a table consisting of time and velocity (t,x). Note the 
difference with the code from part 1.2. Plot the results.

(* A centered 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 *)
vminus = v0;   (* initial value of v *)
xmid     = x1;
vmid     = v1;

(* create a table tvel of ordered (t,v). for 
		 each component the next value is calculated
	  first and then the variables are incremented 
		 in preparation for the next interation. *)
tpos = Table[
	{t, 		
		vplus = vminus - (2*step*kspring/mass)*xmid;
		xplus =  xminus + (2*step*vmid);			
		vminus = vmid;
		vmid  = vplus;
		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 *)
tpos = Prepend[tpos,{t1,x1}];
tpos = Prepend[tpos,{t0,x0}];

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

2. Testing the Algorithm

In this section we will test the validity of the algorithm in the code presented in Part 1.

a. Write down the differential equations that describe the position and velocity of the mass on a spring as a function of time.

b. Given the following set of initial conditions for a particle of mass m=0.33 kg and a spring with k=400 N/m, what do expect for the behavior of the displacement and velocity as a function? In other words, what are the general and particular solutions of the differential equations and what are the values of the phase angle and the amplitude? Work this out with pencil and paper.

                                            v(t=0) = -13.6 m/s      x(t=0) = 15.0 m
                                            
What will plots of x(t) and v(t) look like? What are the amplitudes and phase angles for the position and velocity curves?

c. Use the three-point formulas to obtain algorithms for calculating the position and velocity as functions of the time using the coupled equations approach and the initial conditions from part 2.b. Execute the commands in Part 1 and print out the table containing the ordered pairs of (t,x). You can do this by entering the name of the table in a new cell and hitting 'Enter'. Using your algorithm calculate by hand (use a hand-held calculator) the first four entries in  the table of time and displacement and make your own (paper) table with the following headings.

                           n             time(s)          (Subscript[x, n-1])(m)             (Subscript[x, n]) (m)                 (Subscript[x, n+1])(m)
                              
Do your results agree with the output of the code in Part 1? If not consult your instructor.
                              
d. Print out the table containing the ordered pairs of (t,v). Calculate by hand (use a calculator again) the first four entries in  the table of time and velocity and make your own (paper) table with the following headings.

                           n             time(s)          (Subscript[v, n-1])(m)             (Subscript[v, n]) (m)                 (Subscript[v, n+1])(m)
                              
Do your results agree with the output of the code in Part 1? If not consult your instructor.

e. Modify the code in Part 1 to plot the difference between the numerical solution and the analytical one. Is there a significant difference between the accuracy of the calculation of x(t) and v(t)?

3. The Full Story of Lieutenant Chisov.

Recall the travails of Lieutenant Chisov.  At an altitude of 21,980 feet he fell from his plane after being hit by German gunfire in 1942. He had no parachute. He landed on the slopes of a snow-covered ravine, slid to the bottom, suffered a fractured pelvis and severe spinal damage, but lived. 

Our model of the drag force acting on a body as it falls through a resistive medium (like air) is given by the following expression.

                                       (Subscript[F, f]) = (1/2) ρ C A (v^2)                .
                                       
Use the discussion in Physics for Scientists and Engineers by Serway and Jewett to choose parameters that might be appropriate for Lieutenant Chisov.                                    

a. What is the pair of coupled, differential equations that describe the velocity and position of the lieutenant as a function of time?

b. Generate an algorithm for solving the equations from Part 3a using the centered derivative method.

c. Use the results of Part 3b to solve the set of coupled equations with a stepsize of  h=0.10 s. Make a plot of the lieutenant's speed and position as a function of time. Assume he is dropped from rest so v(t=0) = 0. 

d. How long did he fall?

e. How fast was he going when he hit?

f. How happy was he during his fall?