Maple Worksheet for the Solution for Bonus Problem No. 1

A Trajectory with Quadratic Drag that is also air-density dependent,

and compared with no dependence on air density as well as with no drag at all

    We begin the problem by setting up the entire system of differential equations that needs to be solved.

It is best to write them all as first-order differential equations;

therefore, we want an equation for

      1) the time-derivative of the x-component of the velocity, which we denote by vx,

      2) the time-derivative of the y-component of the velocity, which we denote by vy,

      3) the time derivative of x, which should equal vx,

      4) the time derivative of y, which should equal vy,

and of course, also the dependence of the quadratic drag coefficient, which depends on y.

We are interested in three different cases, actually,

   a) when the drag coefficient depends on y, i.e., it's some constant multiplied by exp(-y/lambda);

   b) when the drag coefficient is just the constant, that one gets by setting y=0 in the equation above,

   c) and when the drag coefficient is simply zero, as if there were no air resistance.

We begin, then, by setting up that drag coefficient, where c0 is a constant, depending on the resistance

of the air and the mass of the object, and lambda, which I will just call L in the typescript,

 is the atmospheric scale length, given to us as 7.600 kilometers.

In principle the coefficient for the drag should be called c; however, when it is inserted into the actual

differential equations that constitute Newton's second Law, that c will be divided by m, so that we need

some new quantity, c/m, which I will call w, which is then the constant w0 times the exponential decrease.

Since we know that---when measured in MKS units---the quantity c is given by D*D/4, we may re-write

the mass m as the density, rho, multiplied by the volume, which is (Pi/6)*D*D*D, we can write out explicitly

what the constant w0 should be that multiplied the exponential decrease with height, namely

w0 = (D^2/4)/[rho*(Pi/6)*D^3] = (3/2*Pi)/(rho*D).

Since this just multiplies the square of the velocity to give the acceleration, it must, again in MKS units,

have dimensions of 1/meter.  Therefore, if I calculate it that way, and then multiply it by the factor (1000 meter/1 km)

I will have the result in kilometers per second, which is the way I want to measure things.

This gives us  the value for w0.  I want to insert that value as late as possible, so I will first calculate it here,

and call it w00.  I use the given values for the density and diameter, but all measured in MKS units:

>	rho:=7800; Dd:=0.15;

>	w00:=(3/(2*evalf(Pi)))*1000/(rho*Dd);

>	w:=w0*exp(-y(t)/L);

>	wwo:=limit(w,L=infinity);

The next pair of equations are just the x- and y-components of Newton's Second Law for this situation,

with all terms moved to the same side of the equation, so that the right-hand side of this equation

should just be zero, which is what Maple will assume when nothing else is given.

>	xdf:=diff(vx(t),t)+w*sqrt(vx(t)^2+vy(t)^2)*vx(t);

>	ydf:=diff(vy(t),t)+g+w*sqrt(vx(t)^2+vy(t)^2)*vy(t);

Lastly, we have to tell the system how the derivatives of x and y are related to vx and vy

>	ddx:=diff(x(t),t)-vx(t); ddy:=diff(y(t),t)-vy(t);

That constitutes the equations, so we now put them all together:

>	syst:={xdf,ydf,ddx,ddy};

   The next thing to do is to take the limit as L goes to infinity, on these equations, which will give us the "second set"

of equations, where the change in air resistance due to altitude is being ignored.

This will make the air resistance somewhat higher at high altitudes than it is the equations above.

>	consyst:=map(limit,syst,L=infinity);

We might as well also now go ahead and create the equations that constitute the case when one totally ignores the air resistance,

which we can do by setting w0 to zero.  Of course in that case we could have solved the differential equations and presented

their solutions and simply graphed them.  But, we can do it this way as well:

>	nosyst:=subs(w0=0,consyst);

The next command sets up the initial conditions for this system of first-order ode's, using the initial speed and initial angle,

and beginning the motion from the origin:

>	init:=[vx(0)=v0*cos(theta0),vy(0)=v0*sin(theta0),x(0)=0,y(0)=0];

   Next, we put in a choice of parameters for the problem at hand, where I am measuring everything in kilometers, and seconds.

>	g:=9.8/1000; v0:=0.400; theta0:=evalf(Pi)/3; L:=7.600; w0:=w00;

The command just below tells Maple to load a "package" of commands that will allow us to plot

the solutions of our differential equations.

>	with(plots):

Warning, the name changecoords has been redefined

And now we can begin the solution of the three different systems of differential equations:

Maple has various different ways to do this.  One of the best is just "dsolve" and "odeplot" which I am using below.

I note that the command "op" when acting on a list---in brackets---or a set---in (curly) braces---simply removes those brackets

or braces, leaving the result as a sequence, so that it can more easily be combined with some other sequence.

So the use of dsolve creates a Maple procedure to do a numerical solution of the specific differential equation with initial

conditions.  Then odeplot uses that procedure to create numerical values over the range that it needs for the plot desired.

We first do that for the complete problem, with altitude-dependent air resistance:

>	nutri:=dsolve({op(syst),op(init)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwww:=odeplot(nutri,[x(t),y(t)],0..41,refine=2,color=blue,thickness=2):

Next we do the same thing when the air resistance does not decrease with altitude:

>	conn:=dsolve({op(consyst),op(init)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwcon:=odeplot(conn,[x(t),y(t)],0..40,color=red,thickness=1):

And then yet again the same thing, but now with NO air resistance

>	nonn:=dsolve({op(nosyst),op(init)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwnon:=odeplot(nonn,[x(t),y(t)],0..71,color=black):

Having saved all three of those graphs, we now plot them together.

>	display([wwww,wwcon,wwnon],title=`at 60 degrees, initial angle`);

Now, let's change the initial angle to 30 degrees, and look at the results.

Recall that for no air resistance the range is the same for 60 and 30 degrees, since they are the same

distance on both sides from 45 degrees, although, of course, the heights are very different.

>	theta0:=evalf(Pi)/6;

>	init2:=[vx(0)=v0*cos(theta0),vy(0)=v0*sin(theta0),x(0)=0,y(0)=0];

>	nutri2:=dsolve({op(syst),op(init2)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwww2:=odeplot(nutri2,[x(t),y(t)],0..26,refine=2,color=blue,thickness=2):

>	conn2:=dsolve({op(consyst),op(init2)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwcon2:=odeplot(conn2,[x(t),y(t)],0..26,color=red,thickness=1):

>	nonn2:=dsolve({op(nosyst),op(init2)},[x(t),y(t),vx(t),vy(t)],type=numeric,range=0..100):

>	wwnon2:=odeplot(nonn2,[x(t),y(t)],0..41,color=black):

>

>	display([wwww2,wwcon2,wwnon2],title=`at 30 degrees, initial angle`);

The two sets of graphs appear "almost" the same; there is a very slight difference in the two graphs for the disparity

between the red and the blue curves.

However, this similarity is misleading because Maple has changed the vertical scale on the two graphs.

Let's plot all 6 on the same figure:

>	display([wwww2,wwcon2,wwnon2,wwww,wwcon,wwnon]);

For amusement we can also use the calculations we have made to plot, for instance the horizontal component of

velocity as a function of time, for the 60 degree launch.

We see that it slows down, in a fairly nonlinear way, toward zero.

It really comes eventually to zero; however, looking at a graph somewhat further down in this presentation,

we see that it strikes the ground about t=40, so that the remainder of the graph does not actually happen unless

there's a high cliff on which all of this is happening.

>	odeplot(nutri,[t,vx(t)],0..100,refine=2,color=green);

As a last comment, related to the fact that the horizontal velocity went to zero, is the graph below,

which shows that x(t) started out nicely, and then is levelling off, never getting beyond 3.5 kilometers.

Compare the blue curve above, where we also see that effect.

However, again, it hits the ground again at about t=40.

>	odeplot(nutri,[t,x(t)],0..100,refine=2,color=red);

Deciding that we should also wonder about the vertical direction, the plot below is

the vertical velocity, which of course starts out positive, decreases to zero about t=19,

and then continues to decrease, overshoots, and starts coming back toward a final,

terminal velocity in that direction.

Actually in this graph it's not completely clear that it's not headed back to zero;

therefore, let's look at a much larger time, in the next figure, even though it doesn't fit with the physics

of this problem, where it hits the ground at about t=40.

>	odeplot(nutri,[t,vy(t)],0..100,refine=2,color=blue);

>	odeplot(nutri,[t,vy(t)],0..1400,refine=2,color=blue);

Lastly, let's look at a graph of the vertical location itself, where we can, finally,

see that it hits the ground again at about t=41.

>	odeplot(nutri,[t,y(t)],0..50,refine=3,color=blue);

>