The phugoid equations

= = - sin - *Rv*^{2}

run into trouble when

> R:=0.3: DEplot([diff(theta(t),t)=(v(t)^2 - cos(theta(t)))/v(t), diff(v(t),t) = -sin(theta(t)) - R*v(t)^2], [theta,v], t=0..10, [seq([theta(0)=0,v(0)=2+i/10],i=0..10)], theta=-1.5..3, v=-0.05..3, title="R=0.3",linecolor=black,stepsize=0.1);

Something looks very wrong with the solution which starts at
= 0, *v* = 2.5. Notice how the solution makes a sharp turn, heading
in a completely different direction from the vector field. Why is
this happening?

Just to verify that it isn't a fluke, let's examine several solutions with
initial conditions
(0) = 0 and
2.4 < *v*(0) < 2.6, looking in
the problem area. We'll use obsrange=false to tell Maple we
want it to continue computing the solutions, even when they go outside
of the viewport.

> DEplot( [diff(theta(t),t) = (v(t)^2 - cos(theta(t)))/v(t), diff(v(t),t) = -sin(theta(t)) - R*v(t)^2 ], [theta,v], t=0..10, [seq([theta(0)=0,v(0)=2.4+.02*i],i=0..10)], theta=1..2, v=-0.05..0.5, obsrange=false, title="R=0.3",linecolor=black,stepsize=0.1);

Even worse! Not only does the solution for *v*(0) = 2.5 go the wrong
way, the one with *v*(0) = 2.52 also turns the wrong way,
crossing over another solution.
Since our equations are autonomous, this is *never* supposed to
happen.

We can remedy this by decreasing the stepsize significantly. Decreasing the stepsize by a factor of 10 (to 0.01) in the above example gives solutions which do the right thing. But a significant price in computation must be paid, and this doesn't really solve the problem. There are other initial conditions nearby which lead to the same trouble.

In order to decide how to fix the problem, we must first understand
what is going wrong. The first problem we noticed happened for a
solution with
(*t*_{i}) 1.44, *v*(*t*_{i}) 0.09, using a
stepsize *h* = 0.1. Since Maple is using a Runge-Kutta method as
described in §4.2, let's calculate what is
happening.

Recall that in Runge-Kutta, we evaluate the vector field at four
points, and average them together to get the next point. We denoted
these vectors as
, , , and
. Since we have to calculate the vector field at four
points, we'll write a little function to do this for us, and use the
fact that Maple lets us do aritmetic on lists as though they are
actual vectors; that is, we can add them and multiply by
scalars.^{3.7}

> h:=0.1: VF:=p->[(p[2]^2 - cos(p[1]))/p[2], -sin(p[1])-R*p[2]^2]:p:=[1.44, 0.9]; k[l]:=VF(p); k[c]:=VF(p+h/2*k[l]); k[m]:=VF(p+h/2*k[c]); k[r]:=VF(p+h*k[m]); step:=h/6*(k[l] + 2*k[c] + 2*k[m] + k[r]);

The problem here is the vector , which has a huge
component, and points in a different direction from the
others. This happens because *k*_{r} is the value of the vector field
calculated at
1.35,
*v* - 0.00315. Not only
is the value of *v* very small here (giving a huge
component), it is outside the range of allowed values, and the vector
field points the wrong way. Similar problems will occur for
any solution which comes too close to the *v* = 0 axis with *v*(*t*)
decreasing.

Now the question is how to avoid this? One not very satisfying answer
is just to use a very small step size. But the problem with this idea
is that we then must do a huge amount of computation, even when the
vector field is behaving nicely.

A better solution is to use an *adaptive stepsize*: when the
vector field is ``nice'', we can use a large step, and when the vector
field is being troublesome, we use a small step. Maple has such
methods built in to it, such as the Runge-Kutta-Fehlberg method, which
we can use by appending method=rkf45 to the list of options to
DEplot.

Sometimes, it is not always possible to change the numerical method.
In the case of the phugoid, we can accomplish the same goal (in fact,
gain something extra) by adjusting the length of the vectors that make
up the vector field.

If we think of our solution curves as being parametric curves with
tangent vectors which coincide with the vector field, then changing
the lengths of the vectors (but not their direction) will not change
the solution *curves*, only the parameterization. That is, we can
rescale time so that the vectors no longer blow up as solutions
approach the *v* = 0 line.

Multiplying each vector by *v*(*t*) cancels out the
problems in
, and gives a system which has the same
solution curves as long as *v* > 0. Thus, we get the new system

= *v*^{2} - cos = - *v* sin - *Rv*^{3} = *v*^{2}cos = *v*^{2}sin

This new system of equations is defined and continuous for all values
of
Furthermore, if we want to keep track of the original time variable
*t*, we can add a new variable *T*. In the original, non-rescaled
equations, the time could be described as *T*(*t*) = *t*, or
equivalently,
= 1. After rescaling the vector field by the
factor of *v*, this gives
= *v*.

> R:='R'; vphug:= [ diff(theta(t),t) = v(t)^2 - cos(theta(t)), diff(v(t), t) = -v(t)*sin(theta(t)) - R*v(t)^3, diff(x(t), t) = v(t)^2 *cos(theta(t)), diff(y(t), t) = v(t)^2 *sin(theta(t)), diff(T(t), t) = v(t) ]:R:=0.3:

DEplot(vphug, [theta, v, x, y, T], t=0..20, [seq([theta(0)=0, v(0)=2+i/10, x(0)=0, y(0)=5, T(0)=0],i=0..10), [theta(0)=0, v(0)=2.51, x(0)=0, y(0)=5, T(0)=0], [theta(0)=0, v(0)=2.52, x(0)=0, y(0)=5, T(0)=0], seq([theta(0)=0, v(0)=2.51+i/1000, x(0)=0, y(0)=5, T(0)=0],i=2..5) ], theta=-1.5..3, v=-0.05..3, x=0..10, y=0..5, T=0..20, obsrange=false, linecolor=black, scene=[theta,v], stepsize=0.1);

> zoom(

As you can see, using the desingularized equations solves the problems
for small values of *v*.

In addition, we can see additional structure not apparent in the
original equations. There are new fixed points at *v* = 0,
= ±/2. These fixed points are saddles, as we can see by
looking at the Jacobian. In the desingularized equations, the
Jacobian is

While it is possible to deduce this behavior from the original equations, it is much more obvious one the singularity is filled in. Similar techniques are used in other branches of mathematics (blowing up a singularity in algebraic geometry, construction of a collision manifold in celestial mechanics) with great success.

- ...
scalars.
^{3.7} - This ability is not present in releases prior to Maple 6, making this computation a little more cumbersome.

2002-08-29