1 Projectile Motion
2 The equations of motion with air resistance
3 Numerical solutions of differential equations in Mathematica
4 Baseball home run with air resistance
5 Other exercises

How to' s 		in Lesson 04

How to 		solve coupled nonlinear differential equations
			NDSolve [ eqs , vars , { x , x1 , x2 } ] 
				
How to 		make parametric plots
			ParametricPlot [ { x[t] , y[t] }, { t , t1 , t2},
				PlotRange -> { {x1 , x2} , {y1 ,y2} } ]

Lesson 04
Week of Monday January 28
Projectile Motion -- Part 2

Outline
1- Projectile Motion
2- The Equations of Motion
3- Numerical Solutions of Differential Equations in Mathematica
4- A baseball home run with air resistance
5- Other exercises

EXERCISE 0 (review exercise)
Use ParametricPlot to make a graph of y versus x;
	x(t) = cos(80t) – (cos^3)( t) and y(t) = sin(t) – (sin^3)(80*t)
	(t runs from 0 to 7)
SKETCH THE GRAPH ON THE ANSWER SHEET (a ROUGH SKETCH IS SUFFICIENT)

E0

ParametricPlot[{Cos[80*t]-Cos[t]^3,Sin[t]-Sin[80*t]^3},
{t,0,100},PlotPoints->10000,PlotStyle->AbsoluteThickness[1]]
Rasterize[%,RasterSize->400];

/1/ Projectile Motion

A mass moving in Earth's gravity near the surface of the Earth is called a projectile.

The motion of a projectile is affected by three forces:
	the force of Earth's gravity  =    (Subscript[(OverVector[F]), g])  = – m (Overscript[g, ->]) ;
	buoyancy (usually very small)  = (Subscript[(Overscript[F, ->]), b]) = + V (Subscript[ρ, air]) (Overscript[g, ->]) ;
	the aerodynamic force exerted by the air =    (Subscript[(OverVector[F]), air])   .

In general the motion is three-dimensional, so we could set up a 3D Cartesian coordinate system, (x,y,z).  However, if there is no sideward aerodynamic force then the particle will only move on a vertical plane, which we’ll take to be the xy-plane; x = horizontal, y = vertical.

In today’s calculations we’ll assume that the only aerodynamic force is aerodynamic friction, also called air resistance.

The force of air resistance has these properties ...

...  the direction is opposite to the direction of the velocity;

...  the magnitude depends on the speed, and on the size and shape and surface roughness of the projectile.

A pretty good approximation for the air resistance force acting on a sphere moving with velocity (OverVector[v])  is
		(Subscript[(OverVector[F]), air]) = – 0.25  ρ  A |(OverVector[v])| (OverVector[v]) 
	where ρ = the density of air = 1.225  kg/(m^3) at sea level;
	A = π (R^2) = cross sectional area of the sphere and R = radius of the sphere.
The magnitude of the force is  (Subscript[F, air]) = 0.25 ρ A (v^2) ;
the direction is parallel to   –(OverVector[v]) , i.e., opposite to the velocity.
We’ll use this model for our calculations today.

/2/  The Equations of Motion

The equation of motion for a sphere moving in the atmosphere near the surface of the Earth is, in vector form,
		m (Overscript[r, ••]) = (Subscript[(OverVector[F]), g]) + (Subscript[(Overscript[F, ->]), b]) + (Subscript[(OverVector[F]), air])

Set up a 2D Cartesian coordinate system , x and y , where y is vertically upward and x is horizontal. The projectile will move in the xy-plane. The equations of motion, in coordinate form, are
		m (Overscript[x, ••]) = – 0.25 ρ A (Sqrt[(Subscript[v, x]^2)+(Subscript[v, y]^2)]) (Subscript[v, x])
		m (Overscript[y, ••]) = – (Subscript[m, e]) g – 0.25 ρ A (Sqrt[(Subscript[v, x]^2)+(Subscript[v, y]^2)]) (Subscript[v, y])  where (Subscript[m, e]) = m – Vρ ≈ m

These amount to two coupled second-order differential equations.

To make it simpler, we’ll set (Subscript[m, e]) = m, which should be a good approximation.

It will be useful to convert the two 2nd-order equations into four 1st-order equations.
The equations are

(Overscript[x, •])(t)  =  v(t)    and    (Overscript[v, •])(t) = – 0.25 (ρA/m)  (Sqrt[(v^2)+(w^2)]) v
	(Overscript[y, •])(t)  =  w(t)    and    (Overscript[w, •])(t) = – g – 0.25 (ρA/m) (Sqrt[(v^2)+(w^2)]) w

Here the velocity vector is  (v,w) = ( (Subscript[v, x]) , (Subscript[v, y]) ) = ( (Overscript[x, •]) , (Overscript[y, •]) ) .

In addition to these differential equations, we must also have some boundary conditions; for example, we might know the coordinates and velocity components at some initial time, t = 0.

The problem is, solving these nonlinear coupled equations analytically, is very difficult.

/3/  Numerical Solutions of Differential Equations in Mathematica

The Mathematica commands for solving differential equations numerically are NDSolve (most general, but a little tricky to use) and NDSolveValue (less general, not always applicable, but easier to use).

Get some help on NDSolve. Type ?NDSolve (Shift-Enter).

?NDSolve

Do you get it?

To use NDSolve, you will type a command like this,

	NDSolve[ eqs, u, { x , xmin , xmax } ]

where eqs is a list of equations (both diff. eqs. and boundary conditions); u is a list of functions (the functions that the equations are supposed to determine); x is the independent variable (i.e., the derivatives in the diff. eqs. are derivatives with respect to x); and xmin ⟶ xmax is the  range of the independent variable.

The output from the command NDSolve will not be a function! The output will be a list of “replacement rules” from which you can extract one or more solution functions. You need to learn how to extract the functions. , using slash-dot. The extracted function will not be an analytic function, but rather an “interpolation function”. That is because NDSolve solves the diff. eqs. numerically.

NDSolveValue is a little easier to use because the output from the command will be the interpolation function directly—no need for the extraction from a replacement rule. However, it can only give one solution function.

Example 1 = Exponential Decay

Let’s start with a really simple differential equation,

   (dP/dt) = – 0.3 P.
	
We are looking for the function P(t). We have a first-order differential equation, so we need one boundary condition. Let’s say P(0) = 500.  Then what is P(t)? I believe that you can figure out the answer yourself—without Mathematica—analytically.

But for practice, let’s get the answer using NDSolve. Type these commands in an Input Cell:

eqs={P'[t] == - 0.3*P[t], P[0] == 500}        (* eqs is the list of equations *)
NDSolve[ eqs , P , {t,0,10}]        (* and now hit Shift-Enter *)

What did Mathematica give us?
The first output cell is just repeating what we typed for eqs.
The second output cell is the result from NDSolve.

The result from NDSolve is a list with only one element. That element is called a “replacement rule”. The replacement rule, used with slash-dot,  will replace the variable P by an “Interpolating Function”.
The replacement is not an analytic function! Instead, because Mathematica solved the equation numerically, the replacement rule is an Interpolating Function, valid for t in the requested range (0,10).
Now let's make a plot of the result.

First run this ...

Remove["Global`*"]
(* now type the equations *)
eqs={P'[t] == - 0.3*P[t], P[0] == 500} 
(* next execute NDSolve and name the result "reprule" *)
reprule=NDSolve[ eqs , P , {t,0,10}] 
(* next extract the solution function using "slash-dot", and give it the name "soln" *)
soln=P /. reprule[[1]]
(* Evaluate some points to make sure it's working *)
{soln[0.0],soln[1.0],soln[2.0],soln[10.0]}

Now make the plot ...

Plot[ soln[t] , {t,0,10}, GridLines->Automatic ]

(* As you probably know,
the analytic solution is 500*Exp[-0.3*t] ; 
let's verify it. *)
Plot[   500*Exp[-0.3*t]  ,{t,0,10},GridLines->Automatic]

EXERCISE 1
WRITE YOUR ANSWERS ON THE ANSWER SHEET; ACCURCY COUNTS!
What is the analytic solution to Example 1? 
Calculate (Subscript[P, exact])(1) and  (Subscript[P, exact])(2) to 6 place accuracy.
Now check that NDSolve gave the same results.

E1

500*Exp[-0.3*t]  /. {t->1}
500*Exp[-0.3*t]  /. {t->2}
{soln[1],soln[2]}

Example 2 = numerical solution of an unknown diff. eq.

Now let’s try a more difficult example,

   (((d^2)f)/(du^2)) + [ 1+u f(u) ]  (df/du) + 3 (f^2) = 0.			(☺)
	
We are looking for the function f(u). We have a second-order differential equation, so we need two boundary conditions. Let’s say f(0) = 10 and f ’(0)= 0.  Then what is f(u) for u from 0 to 0.7?

I’m pretty sure you can’t figure out the answer yourself , i.e., analytically.

So let’s get the answer using NDSolve. This time, we’ll use the simpler form NDSolveValue.

?NDSolveValue

First run this ...

Remove["Global`*"]
eqs={ f''[u]+(1+u*f[u])*f'[u]+3*f[u]^2 == 0 ,
f[0] == 10, f'[0]==0}
f2=NDSolveValue[ eqs , f , {u,0,0.75}] 
{f2[0.0],f2[0.1],f2[0.2]}

Now make the plot ...

Plot[ f2[u], {u,0,0.75}, GridLines->Automatic,PlotStyle->Red]

How do we know that the answer is correct?

We cannot compare the Interpolating Function  f2[x] to the analytic solution, because we don’t know the analytic solution ! However, we can verify that the Interpolating Function obeys the differential equation (☺):

Plot[f2[u],{u,0.0,0.75},GridLines->Automatic]    
(* this is the solution from NDSolveValue *)

Plot[{f2''[u],(1+u*f2[u])*f2'[u],3*f2[u]^2}, 
{u,0.0,0.75},GridLines->Automatic,
PlotStyle->{Pink,Pink,Pink}]  
(* these are the three parts of the diff. eq. *)

Plot[f2''[u]+(1+u*f2[u])*f2'[u]+3*f2[u]^2, 
{u,0.0,0.75},
PlotRange->{ {0.0,0.75},{-0.002,0.002} },
PlotStyle->Blue]
(* this is the left-hand side of the diff. eq. (the sum of the three parts); it should be 0 *)

(* Now hit Shift-Enter *)

The function plotted in blue should be 0. It is not exactly 0, but it is very small. So, the numerical solution is not perfect, but it is close. The error is less than 0.002 out of 10, so that is pretty close; error < 2 × (10^(-4)).

EXERCISE 2
WRITE YOUR ANSWERS ON THE ANSWER SHEET; ACCURACY COUNTS.
What is the solution to Example 2?  Show the answer as a sketch of the graph of f(u).  How do you know that NDSolveValue gave the correct result?

/4/  A baseball home run with air resistance

WRITE AND DRAW YOUR ANSWERS TO THE QUESTIONS (A) — (F) ON THE ANSWER SHEET.

EXERCISE 3 : A baseball home run with air resistance
A baseball batter hits a home run.
As the ball leaves the bat, the position coordinates are ((Subscript[x, 0]) , (Subscript[y, 0]) ) and the components of the velocity vector are ( (Subscript[v, 0x]) , (Subscript[v, 0y]) ). The initial direction of the velocity vector is at angle (Subscript[θ, 0]) above the horizontal.
(A) DRAW A PICTURE THAT SHOWS THE INITIAL STATE OF THE BASEBALL.
Label these quantities on the picture:  (Subscript[y, 0]) , (Subscript[v, ox]) , (Subscript[v, oy]) , (Subscript[θ, 0]) .

(B) WRITE THE EQUATIONS THAT YOU NEED FOR THE CALCULATION.
(Hint: Refer back to Section 4.2; don’t forget the initial conditions.)
Have someone verify that you have the right equations.

Parameters of a regulation baseball:
mass = 0.145 kg
radius = 3.28 cm = 0.0328 m

(C) Assume these initial values:
	(Subscript[x, 0]) = 0 ;  (Subscript[y, 0]) = 3 ft ; (Subscript[v, 0]) = 140 ft/sec ; (Subscript[θ, 0]) = 45 degrees = π/4 radians.
Use NDSolve to determine x(t) and y(t).
Then plot the trajectory of the baseball on the xy-plane with ParametricPlot.
[  Hint: Before you start, convert the distances to meters. 1 m = 3.28 feet.
Then do the calculation in SI units. But finally convert back to feet for the unit of length.
The plot of the trajectory should be vertical axis = y (in feet) versus horizotal axis = x (in feet). ]

(D) Calculate the time (Subscript[t, grd]) when the ball will hit the ground.
(E) Calculate the horizontal distance (Subscript[x, grd]) ( in feet )  where the ball will hit the ground.
(F) Is the effect of air resistance negligible? Explain clearly and accurately. For example, compare the range with and without air resistance.

/5/ Other Exercises

EXERCISE 4: A home run in Denver
The density of air at sea level is (Subscript[ρ, 0]) = 1.225 kg/(m^3).
The density of air in Denver is (Subscript[ρ, 1]) = 0.98 kg/(m^3).   
Consider the same home run initial values as in Exercise 3.
Calculate the horizontal distance where the ball will hit the ground for a home run in Denver, and compare the result to sea level.
WRITE THE ANSWER ON THE ANSWER SHEET; ACCURACY COUNTS.

EXERCISE 5: A ballistic missile test; Part A
A country has developed a ballistic missile.
To test the missile they fire it straight vertically upward.
The initial velocity is
	(Subscript[(Overscript[v, →]), 0]) = (Subscript[v, 0])  (OverHat[y])  	where 	(Subscript[v, 0]) =  6264.2 m/s ,
at ground level ( i.e., at (Subscript[y, 0]) = 0 and t = 0 ).
(A) First approximate the force of gravity by (Subscript[(Overscript[F, →]), g]) = – m g  (OverHat[y])  .
( g = 9.81 m/(s^2)  ;  we can ignore air resistance here because most of the time the missile is outside Earth’s atmosphere.)
Solve the equations using NDSolve or NDSolveValue.
Plot the trajectory (ParametricPlot).
By graphical analysis calculate the maximum altitude ( (Subscript[y, max]) ) that the rocket will reach;
and also the time ( (Subscript[t, 1]) ) when it will return to the ground.
WRITE THE ANSWERS ( (Subscript[y, max]) and (Subscript[t, 1]) ) ON THE ANSWER SHEET; ACCURACY COUNTS!

EXERCISE 5: A ballistic missile test; PART B
(B) Now use the correct Newtonian formula for the force of gravity,
	(Subscript[(Overscript[F, →]), g])  = – GMm /(((R+y))^2)    (OverHat[y])  . 
( GM = g (R^2)   and  R = the radius of the Earth = 6400.0 km. )
Determine the maximum altitude ( (Subscript[y, max]) ) that the rocket will reach;
and also the time ( (Subscript[t, grd]) ) when it will return to the ground.
WRITE THE ANSWERS ( (Subscript[y, max]) and (Subscript[t, grd]) )  ON THE ANSWER SHEET; ACCURACY COUNTS!

ANSWER SHEET FOR LESSON 04

NAME _____________________________

HAND IN DATE _______________________

EXERCISE 0
 Sketch the parametric plot

EXERCISE 1

analytic solution:

(Subscript[P, exact])(1) and (Subscript[P, exact])(2) = 

(Subscript[P, NDSolve])(1) and (Subscript[P, NDSolve])(2) =

EXERCISE 2
Sketch the plot of f(u)

EXERCISE 3A
Diagram of the problem

EXERCISE 3B
The differential equations

EXERCISE 3C
Sketch the parametric plot; label the axes: x [in feet] and y [in feet]

EXERCISE 3D

(Subscript[t, grd]) =

EXERCISE 3E

(Subscript[x, grd]) =

EXERCISE 3F
Is the air resistance important?

EXERCISE 4
In Denver, (Subscript[x, grd]) =

EXERCISE 5A
Neglecting air resistance ...

Estimate (Subscript[y, max]) ≈

Estimate (Subscript[t, grd]) ≈

EXERCISE 5B
Including air resistance ...

Estimate (Subscript[y, max]) ≈

Estimate (Subscript[t, grd]) ≈

My Calculations

Exercise 0 = review exercise

Remove["Global`*"]
x[t_]=Cos[80t] - Cos[t]^3
y[t_]=Sin[t] - Sin[80 t]^3
ParametricPlot[{x[t],y[t]},{t,0,7}]

Exercise 1 = exponential decay

Exercise 2 = an unknown diff. eq.

Remove["Global`*"]
eqs={ f''[u]+(1+u*f[u])*f'[u]+3*f[u]^2 == 0 ,
f[0] == 10, f'[0]==0}
f2=NDSolveValue[ eqs , f , {u,0,0.75}] 
{f2[0.0],f2[0.1],f2[0.2]}

Plot[ f2[u], {u,0,0.75},
Frame->True,FrameLabel->{"u","f(u)"},
 GridLines->Automatic,PlotStyle->Red]

Exercise 3 (baseball with air resistance)

(* FIRST convert everything to SI units *)
Remove["Global`*"]; 
{m, g} = {0.23, 9.81};
{y0, v0, θ0} = {0.99 , 30.0, Pi/4};
{ρ, R} = {1.225, 0.05};
const=ρ*Pi*R^2;

bbeq={
x'[t]==v[t],v'[t]==(1/m)*(-0.25*const*Sqrt[v[t]^2+w[t]^2]*v[t]),
y'[t]==w[t],w'[t]==(1/m)*(-m*g-0.25*const*Sqrt[v[t]^2+w[t]^2]*w[t]),
x[0]==0,v[0]==v0*Cos[θ0],y[0]==y0,w[0]==v0*Sin[θ0]}
replacementrules=NDSolve[bbeq,{x,y},{t,0,5}]

xCoord=x /. replacementrules[[1]]
yCoord=y /. replacementrules[[1]]
fig3=ParametricPlot[{3*xCoord[t],3*yCoord[t]},{t,0,5},
Frame->True,FrameLabel -> {"x [ft]","y [ft]"},
GridLines->Automatic]

Exercise 4 (baseball home run in Denver)

(* FIRST convert everything to SI units *)
(* Remove["Global`*"];   *)
{m, g} = {0.23, 9.81};
{y0, v0, θ0} = {0.99 , 30.0, Pi/4};
{ρ, R} = {0.98, 0.05};
const=ρ*Pi*R^2;

bbeq={
x'[t]==v[t],v'[t]==(1/m)*(-0.25*const*Sqrt[v[t]^2+w[t]^2]*v[t]),
y'[t]==w[t],w'[t]==(1/m)*(-m*g-0.25*const*Sqrt[v[t]^2+w[t]^2]*w[t]),
x[0]==0,v[0]==v0*Cos[θ0],y[0]==y0,w[0]==v0*Sin[θ0]}
replacementrules=NDSolve[bbeq,{x,y},{t,0,5}]

xCoord=x /. replacementrules[[1]]
yCoord=y /. replacementrules[[1]]
fig4=ParametricPlot[{3*xCoord[t],3*yCoord[t]},{t,0,5},
Frame->True,FrameLabel -> {"x [ft]","y [ft]"},
GridLines->Automatic]

Show[fig3,fig4]

Exercise 5 (the ballistic missile test)

EXERCISE 5 Part A
A country has developed a ballistic missile. To test the missile they fire it straight vertically upward.
The initial velocity is    (Subscript[(Overscript[v, →]), 0]) = (Subscript[v, 0])  (OverHat[y])  where   (Subscript[v, 0]) =  6264.2 m/s , at ground level ( i.e., at (Subscript[y, 0]) = 0 and t = 0 ).
(A) First approximate the force of gravity by (Subscript[(Overscript[F, →]), g]) = – m g  (OverHat[y])  .
( g = 9.81 m/(s^2)  ;  we can ignore air resistance here because most of the time the missile is outside Earth’s atmosphere.) Solve the equations using NDSolve or NDSolveValue. Plot the trajectory (ParametricPlot). By graphical analysis calculate the maximum altitude ( (Subscript[y, max]) ) that the rocket will reach; and also the time ( (Subscript[t, 1]) ) when it will return to the ground. WRITE THE ANSWERS ( (Subscript[y, max]) and (Subscript[t, 1]) ) ON THE ANSWER SHEET; ACCURACY COUNTS!

v0=6264.2;g=9.81;
meqs={y''[t]==-g, y'[0]==v0, y[0]==0}
yA=NDSolveValue[meqs,y,{t,0,1400}]
partA=Plot[yA[t]/1000,{t,0,1400},
Frame->True,FrameLabel->{"t [sec]","y [km]"}]

t1=t /. FindRoot[yA[t]==0,{t,1200}]
ta=t /. FindRoot[yA'[t]==0,{t,800}]
ymax=yA[ta]/1000

EXERCISE 5 PART B
(B) Now use the correct Newtonian formula for the force of gravity,
(Subscript[(Overscript[F, →]), g])  = – GMm /(((R+y))^2)    (OverHat[y])  . 
( GM = g (R^2)   and  R = the radius of the Earth = 6400.0 km. )
Determine the maximum altitude ( (Subscript[y, max]) ) that the rocket will reach; and also the time ( (Subscript[t, 1]) ) when it will return to the ground. WRITE THE ANSWERS ( (Subscript[y, max]) and (Subscript[t, 1]) )  ON THE ANSWER SHEET; ACCURACY COUNTS!

v0=6264.2;g=9.81;R=6400*^3;
GM=g*R^2;
meqsB={y''[t]==-GM/(R+y[t])^2, y'[0]==v0, y[0]==0}
yB=NDSolveValue[meqsB,y,{t,0,2200}]
partB=Plot[yB[t]/1000,{t,0,2200},
Frame->True,FrameLabel->{"t [sec]","y [km]"}]

t1=t /. FindRoot[yB[t]==0,{t,2000}]
ta=t /. FindRoot[yB'[t]==0,{t,1000}]
ymax=yB[ta]/1000

Show[partB,partA]

ANSWER SHEET FOR LESSON 04

NAME ___________________________      DATE HANDED IN___________________________

EXERCISE 1

EXERCISE 2

EXERCISE 3

EXERCISE 1

EXERCISE 1

EXERCISE 1

EXERCISE 1

EXERCISE 1