topic:  vector fields, flow lines, and Euler's method
description:  discusses Euler's method and the approximation of flow lines in vector fields.  Contains applications to free fall and to pendulum motion.
file name:  vectorfield.nb
directions:  Read the text cells for background information.  Execute any existing input cells as you come to them.  Answer questions or work exercises as you come to them.  Modify the existing input cells as needed to work exercises requiring the use of Mathematica.

Computer Assignment #6
Euler's Method and Vector Field Flow

Special thanks to Calculus&Mathematica (Addison-Wesley, 1994) for their support of these assignments.  Each is very much in the style of C&M, and some problems and commands are adaptations of C&M problems and commands.

Why should I have to do this computer assignment anyway?

Vector fields and vector field flow are the geometric tools which we use in the study of differential equations.  Differential equations, in turn, provide us with models for a large number of the observable phenomena in the "real-world."
	In this assignment, you will learn about flow lines (which are the solution curves to differential equations) and how they help us understand and visualize the world around us.  You will also gain some experience with a numerical method for "solving" differential equations known as Euler's Method.

Note:  Since we will be discussing real-world applications in this assignment, it will be important to note the units involved.  Be certain to include the appropriate units in your answers.

Part 1 -- Euler's Method

Euler's method is one of the most commonly-used numerical methods for APPROXIMATING flow lines through a vector field.  In English, the method uses the following algorithm:

EULER'S METHOD:
	Look at your current position in the vector field.  Associated with that position is a corresponding field vector.  That corresponding field vector gives you the velocity (direction and speed) at which you should begin traveling.  Travel for some small amount of time with that velocity.  This takes you to a new position in the vector field.  Repeat the process by looking at the field vector associated with your new position and traveling with that as your new velocity.  Continue in this fashion.

	Euler's method is analogous to finding your way through a forest with a compass.  Rather than traveling toward distant landmarks, which you are apt to lose sight of during your travel, you should periodically sight on nearby landmarks, taking new compass readings once each landmark is reached.

Example 4 on Page 335, Section 17.2

The following cells show how to code Euler's method for Example 4 on Page 335 of the textbook:

First, we define and plot the field.

field[x_,y_] = {y^2,2x^2};
Needs["Graphics`PlotField`"];
fieldplot = 
	PlotVectorField[field[x,y],{x,0,2},{y,1,3},
		Axes->True,
		AxesLabel->{"x","y"}];

Next, we apply Euler's method.  The example asks us to approximate the flow line through the point (1,2).  The book chooses 9 time steps, each of size 0.02 time units.

Clear[next,x,y,point,k];
start = {1,2};
timestep = 0.02;
steps = 9;

next[{x_,y_}] = {x,y} + timestep * field[x,y];

point[0] = start;
point[k_] := point[k] = N[next[point[k-1]]];

path = Graphics[
     {Thickness[0.01],
      Line[Table[point[k],{k,0,steps}]]}];

Show[path,fieldplot,
	Axes->Automatic,
	AxesLabel->{"x","y"},
	AspectRatio->Automatic,
	PlotRange->{{0,2},{1,3}}];

That flow line seems to "follow" the vector field!

Finally, here is a recreation of Table 17.2 in the book:

Table[point[k],{k,0,steps}]

a.  Use Mathematica to work Problem 8 on Page 337 of the textbook.  Show a vector field plot as well as a table of values.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

b.  How does your approximate flow line (for 0 < x < 2, 1 < y < 3) compare with the approximate flow line of Example 4 on Page 335?  Is it a better or a worse approximation of the true flow line than that of the example?  How do you know?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

c.  The flow line which passes through the point (0,2) also passes through the point (1.86,____).  What value (approximately) should go in the blank to make the statement true?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

Part 2 -- An Application to Free Fall

When an object falls under the influence of gravity, its height after t seconds is described by

Clear[h,t,v0,h0,g];
h[t_] = g/2 * t^2 + v0 * t + h0

where v0 is the initial velocity of the object, h0 is the initial height of the object, and g is the acceleration due to gravity.  (You have undoubtedly seen this function numerous times, both in physics and in calculus!)  Let's use vector fields to gain an understanding of how objects behave in free fall.

The velocity is given by

Clear[v,t];
v[t_] = h'[t]

and the acceleration is given by

Clear[a,t];
a[t_] = v'[t]

Suppose we were interested just in the relationship between a falling object's height and its velocity.  Whenever we're interested in the relationship between two variables, a fairly natural thing to do is to sketch the relationship in a plane.  So, we could represent the velocity vs. the height in a two-dimensional coordinate system, called a phase plane.  Every possible height and velocity combination for a falling object corresponds to exactly one point in the phase plane.
	For example, suppose that at time t = 0, we throw a football from ground level (h0 = 0) with an initial upward velocity of v0 = 15 m/sec.  We could represent this information as follows:

ListPlot[{{0,15}},
	Axes->True,
	AxesLabel->{"h","v"},
	PlotRange->{{-25,25},{-25,25}},
	AspectRatio->Automatic,
	PlotStyle->{PointSize[.05]}];

Now, what will happen to the football as time progresses?  Obviously, it is going to change both its height and its velocity!  But HOW?!
	Vectors are the answer!  Our location in the phase plane at time t is given by the vector {h(t),v(t)}.  Our change in location in the phase plane will be given by the vector derivative:  {h'(t),v'(t)}.
	BUT, h'(t) = v(t), by definition, and v'(t) = g.  So, we have a vector field which tells us how our location in the phase plane changes.
	Let's work under the influence of earth's gravity, so that g = -9.80 m/sec^2:

Clear[fallfield,h,v];
fallfield[h_,v_] = {v,-9.80};
Needs["Graphics`PlotField`"];
fieldplot = 
	PlotVectorField[fallfield[h,v],{h,-25,25},{v,-25,25},
		Axes->True,
		AxesLabel->{"h","v"}];

Okay, so WHAT'S THE POINT?!
	Well, you've probably been asked at some point in your life a question such as the following:

QUESTION:  A football is thrown from ground level with an upward velocity of 15 m/sec.
	a.	Find the maximum height that the football reaches.
	b.	Find the speed of the football when it is 5 meters off the ground.

Now you have a visual tool for answering these two questions, quickly and easily!  The answer depends on the flow line through the phase plane!  Let's use Euler's method:

Clear[next,x,y,point,k];
start = {0,15};
timestep = .1;
steps = 30;

next[{x_,y_}] = {x,y} + timestep fallfield[x,y];

point[0] = start;
point[k_] := point[k] = N[next[point[k-1]]];

path = Graphics[
     {Thickness[0.01],
      Line[Table[point[k],{k,0,steps}]]}];

Show[path,fieldplot,
	Axes->Automatic,
	AxesLabel->{"h","v"},
	AspectRatio->Automatic,
	PlotRange->{{-25,25},{-25,25}}];

Now you should be able to see that the maximum height which the football reaches is roughly 12 meters.

a.  What is the (approximate) speed of the football when it is 5 meters off the ground?  (Note that there are TWO times when it is 5 meters off the ground!)  Remember to include your units.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

b.  Suppose we drop an object from 20 meters off the ground.  (Note that "drop" means the initial velocity is 0.)  Plot the flow line corresponding to the object's change in height and velocity.  What is the object's approximate velocity on impact with the ground?  Remember to include your units.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

Part 3 -- An Application to Pendulum (Periodic) Motion

In Physics 215, (you will learn/you have learned) the equation of motion for a simple pendulum.  Here it is:  If a pendulum consisting of a point mass of m kilograms is suspended by a light string of length L meters, and theta(t) represents the angle in radians which the string makes with the vertical at time t seconds, then, for small values of the angle theta(t), the motion of the pendulum is governed by the differential equation theta''(t) = -g/L theta(t).  (We can define negative theta values as representing angles to the left of vertical as we look at the pendulum.)
	Suppose we have a string of length 4.9 meters.  Then, we have the simple differential equation theta''(t) = -2 theta(t).
	Let's study the relationship between the angle theta(t) and the time rate of change of the angle, given by theta'(t).  Our phase plane is the plane of coordinates {theta(t),theta'(t)}.
	As in Part 2, our vector field is given by the time derivative of the phase plane coordinates:
{theta'(t),theta''(t)}.  So, we define our vector field as follows:

Clear[pendulumfield,theta,thetaprime];
pendulumfield[theta_,thetaprime_] =
			{thetaprime,-2 theta};
Needs["Graphics`PlotField`"];
fieldplot = 
	PlotVectorField[pendulumfield[theta,thetaprime],
		{theta,-0.2,0.2},{thetaprime,-0.2,0.2},
		Axes->True,
		AxesLabel->{"theta[t]","theta'[t]"}];

Let's release (initial angular speed 0) the pendulum from 0.1 radians from vertical and see what happens:

start = {0.1,0};

timestep = .02;
steps = 60;

Clear[next,point,k];
next[{x_,y_}] = {x,y} + timestep pendulumfield[x,y];

point[0] = start;
point[k_] := point[k] = N[next[point[k-1]]];

path = Graphics[
     {Thickness[0.01],
      Line[Table[point[k],{k,0,steps}]]}];

Show[fieldplot,path];

a.  At what rate is the angle changing when the pendulum is hanging straight up-and-down vertical?  Remember to include your units.

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

Change the line steps = 60 in the code above to steps = 120.

b.  What type of common curve does the flow line appear to follow, based on the above plot?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

Change the line steps = 120 to steps = 240.

c.  The flow line seems to be spiraling outward.  What does an outward-spiraling flow line say about the motion of the pendulum?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

d.  The answer to Part c should be disturbing to you.  Explain why the flow line above seems to spiral outward.  How would you "fix" it, or at least make it better?  (Hint:  Review Part 1, Question b.)

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

In reality, pendulums will not swing forever.  They experience friction.  One type of friction which they experience is air resistance.  The air resistance which a pendulum experiences is proportional to how fast it is traveling through the air.  This is reflected by the refined differential equation
theta''(t) = -b theta'(t)  - g/L theta(t).
	Let's introduce a drag coefficient of b = 0.5 into our previous example.  Then our vector field changes as follows:

Clear[pendulumfield,theta,thetaprime];
dragfield[theta_,thetaprime_] = 
			{thetaprime,-0.5 thetaprime - 2 theta};
Needs["Graphics`PlotField`"];
fieldplot = 
	PlotVectorField[dragfield[theta,thetaprime],
		{theta,-0.2,0.2},{thetaprime,-0.2,0.2},
		Axes->True,
		AxesLabel->{"theta[t]","theta'[t]"}];

e.  Describe the flow line passing through the point (0.1,0).  What does this say about the motion of the pendulum?

<<ANSWER QUESTION IN YOUR COPY OF THE FILE COMP6ANS.MA>>

The differential equations above are the governing principles behind periodic motion, both undamped (without friction) and damped (with friction).  The universe around us is quite fond of periodic motion.  The motion of pendulums and springs are just two examples from mechanics.  The orbits of planets are another.  We also find abundant periodic behavior in electromagnetic physics.  RLC circuits are governed by the exact same form of differential equation as above.
	You will do more with these types of differential equations if/when you take Math 245.