In §2, we exhibited a number of solutions to
the phugoid model with various initial conditions and values for the drag,
as plots of the parametric curves
[(t), v(t)]. From these, we
could deduce an approximate path of the glider through space via reasoning
like ``
(t) is monotonically increasing, so the glider must be doing
loop-the-loops''. We will now discuss how to solve for the position
[x(t), y(t)] directly.
First, notice that if we know the value of v and at some time
t, we also have
and
at that t, since
Aside from the fact that there is an easier way, there is a serious
practical problem with this approach. When we get a numerical solution
[(t), v(t)], this is a collection of points
(
, v0, t0),(
, v1, t1), ... ,(
, vn, tn)
along a curve. This means we only know the vector field for (x, y) at
those points. It is highly unlikely that these will be the points we need
to use Runge-Kutta. We could use these to approximate our (x, y)-solution
using Euler's method, although even then we would have to program it
ourselves rather than use Maple's built-in method.
If you think for just a moment, you should realize that there is no reason
we must find (t) and v(t) before attempting to find x(t) and
y(t). Instead, we can augment the original system, adding in the new
equations for
and
. Then, we can use numerical methods
to solve for
[
(t), v(t), x(t), y(t)] all in one step. The methods
we discussed in §4.2 apply equally well to vector
fields in any number of variables.
Consequently, we rewrite our phugoid system as
Asking Maple to solve this new system is almost the same as before, but we
must specify the additional variables and use the scene option to
DEplot to tell Maple which projection we want to see. Note that we
reset R to ensure that it appears in the equations as an arbitrary
constant, rather than with any previously assigned value.
>
R:=`R`:
xphug :=[ diff(theta(t),t) = (v(t)^2 - cos(theta(t)))/v(t),
diff(v(t),t) = -sin(theta(t)) - R*v(t)^2 ,
diff(x(t),t) = v(t)*cos(theta(t)),
diff(y(t),t) = v(t)*sin(theta(t))];
Now we ask Maple to plot two solution curves in both , v coordinates
and x, y. One is started with an initial angle of
= 0 and velocity
v = 1.5 at a height y(0) = 1; the other starts with
= 0, v = 2.5, and
initial height y = 2. We also use display with an array to
show the two scenes side by side.
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[x,y]),
])));
>
R:=0;
plots[display](array([
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[theta,v]),
If we change the value of R to 0.2 and repeat the command, we can see
how things change when there is friction. Notice how spiraling towards the
fixed point in the -v coordinates corresponds to a plane which is
diving with an oscillatory motion.
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[x,y]),
])));
>
R:=0.2;
plots[display](array([
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[theta,v]),
Finally, we redo the pictures with a large amount of friction (R = 3). Note how the glider acts more like a rock than a glider, and there is no oscillation at all.
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[x,y]),
])));
>
R:=3;
plots[display](array([
DEplot(xphug, [theta,v,x,y], t=0..15,
[[theta(0)=0,v(0)=1.5,x(0)=0,y(0)=1],[theta(0)=0,v(0)=2.5,x(0)=0,y(0)=2]],
theta=-Pi..5*Pi/2, v=0.1..3, x=-1..10, y=0..4,
linecolor=black, stepsize=0.1, scene=[theta,v]),