First Order Differential 
Equations

Explore the solutions of a variety of first order differentail equations using the numberical methods developed in class.

1. An Example.

Consider the differential equation shown below which is reminiscent of the nuclear decay law.

                                              (dy/dt) = - y
            
 Applying the right derivative formula to the left side generates a finite difference equation that can be solved iteratively.
 
                                             (Subscript[y, n+1]) = (1-h)(Subscript[y, n])
                                             
The code below presents an algorithm for solving this finite difference equation. You might want to experiment with different initial conditions and stepsizes to explore the behavior of the algorithm.

(* Initial Conditions. *)

y0 = 5.0 ; (* starting value of y *)
t0 = 0.0;  (* starting value of t *)

(* The algorithm
	
	   Use a right difference formula to generate a
	   finite difference equation to solve the nuclear
	   decay law. *)
			
(*	 Set up the starting point. *)
h = 0.1; (* stepsize *)
yn=y0;   (* initial value of y *)

(* the limits of the iterations. since we already
	  have y(t=0) given, then the first value of 
		 y we calculate is for t=h.  *)
tmin = h; (* first calculated point *)
tmax = 10.0; (* maximum t *)
stepsize = h; (* step size *)

(* create a table t1 of ordered (t,y) *)
t1 = Table[
	{t,
		yplus = yn*(1-stepsize); (* get the next value of y *)
		yn = yplus   (*increment yn for the next step *)
	},
	{t,tmin,tmax,h}
];
(* the table contains points that start at t=h so we have to add the starting point to the table *)
t1 = Prepend[t1,{t0,y0}];

p1=ListPlot[t1,PlotRange->{{t0,tmax},{0,1.1` y0}},Joined->True,Frame->True,FrameLabel->{"t","y","",""},PlotStyle->RGBColor[1,0,0]]

2. Problems in Paradise.

Consider the example shown part 1. We now want to explore the accuracy of the algorithm used to solve the nuclear decay law described with the equation below.

                                              (dy/dt) = - y
            
a. What is the exact solution of the differential equation above?

b. Generate a plot of the value of y as a function of time up to t=10 (i.e., just use the stuff in part 1).

c. Modify your code in part 2b above to calculate the difference between the estimated solution you calculate with the algorithm from part 1 and the exact solution. Where is the numerical solution valid? Does the accuracy get better as the time t get larger? Save your plot.

d. Modify your code to calculate the normalized difference or

                                             δ((Subscript[t, n])) = (((Subscript[y, n]) - (Subscript[y, 0])(e^(-(Subscript[t, n]))))/((Subscript[y, 0])(e^(-(Subscript[t, n])))))                .

Where do you think the numerical solution is valid now? Does the accuracy get better as the time t get larger? Can you explain what is happening?

3. Poor Lieutenant Chisov.

In January, 1942 a Soviet Ilyushin 4 flown by Lieutenant I.M.Chisov was badly damaged by German gunfire. At an altitude of 21,980 feet Lieutenant Chisov fell from the plane. Unfortunately, he did not have a parachute on when he fell. He landed on the slopes of a snow-covered ravine and slid to the bottom. He suffered a fractured pelvis and severe spinal damage, but lived. 

When an object like the good lieutenant falls through a resistive medium (like air) its acceleration is retarded by a force proportional to the square of its speed.

                                       (Subscript[F, f]) = (1/2) ρ C A (v^2)                .
                                       
Consider an object falling through the air accelerated by the force of gravity and retarded by air friction. Use the discussion in Principles of Physics by Serway and Jewett to choose parameters that might be appropriate for Lieutenant Chisov.

a. What is the total force on the Lieutenant Chisov?

b. Express your result from part a. as a first-order, linear, ordinary differential equation where the speed v is a function of the time t.

c. Solve the equation from part b using the right derivative or two-point method and a stepsize of  h=0.1 s. Make a plot of the object's speed as a function of time. Assume the object (Lieutenant Chisov) is dropped from rest so v(t=0) = 0. 

d. Objects falling through a resistive medium typically reach a terminal velocity. What is that velocity for skydivers? Does the speed of the object in your model approach a limiting value? What is that value? Is it consistent with the observed terminal velocity? 

e. Find the speed at t=25  s for several values of h in the range (10^(-5 ))s to 1.0 s. Is the solution stable over this range of stepsizes? A semilog plot of your results may be appropriate to adequately view the results.

f. How fast was Lieutenant Chisov moving when he hit the ravine? 

g. Have we done enough to figure out how long he had to contemplate the imminent impact? If not, what is missing?