We'll look at the same phugoid model as in the first lesson, but instead of representing the velocity as speed & angle, we'll write everything in $x$ and $y$ coordinates. The results ought to be the same in the end.
In Figure 1, $\vec L$ is the lift vector, $\vec W$ is the weight vector, $\vec v$ is the velocity vector of the plane. We'll write $L, W$, and $v$ for the lengths of these vectors, and $L_x$, $L_y$, etc., for their Cartesian $x$ and $y$ components.
The total force on the plane is $\vec F = \vec L + \vec W$. Once we know the total force, Newton's second law, ``$F=ma$,'' tells us how the velocities change with time, $$ m\frac{dv_x}{dt} = F_x \qquad m\frac{dv_y}{dt} = F_y, $$ while the $(x, y)$ coordinates of the plane change according to $$ \frac{dx}{dt} = v_x \qquad \frac{dy}{dt} = v_y. $$
The weight vector $\vec W$ is easy: its horizontal component $W_x$ is zero, its vertical component is $W_y=-mg$. The magnitude of the weight is $W=mg$.
To find the lift vector we first find its magnitude, and then its direction. Recall, from the first lesson, that the strength of the lift is proportional to the squared velocity $v^2$, and that the trim velocity $v_t$ is defined to be the velocity at which lift and weight are equal (in size, not necessarily direction). This led us to the following formula for the strength of the lift: $$ L = mg\frac{v^2}{v_t^2}.\tag{1} $$ To find the direction of the lift look at figure 1: the two yellow triangles are similar, so we have $$ L_x : L_y = -v_y :v_x. $$ Thus $$ \frac{L_x}{-v_y} = \frac{L_y}{v_x}. $$ If we call this ratio $a$ then we can write the components of the lift as $L_x = -av_y$, $L_y= av_x$. The length of the lift is given by (1), so we get $$ L = \sqrt{L_x^2+L_y^2} = \frac{mg}{v_t^2} \bigl(v_x^2+v_y^2\bigr). $$ Substituting $L_x = -av_y$, $L_y= av_x$, allows us to solve for $a$. The upshot is that $$ L_x = -\frac{mg}{v_t^2} vv_y, \qquad L_y = \frac{mg}{v_t^2} vv_x \tag{2} $$ (where $v=\surd(v_x^2+v_y^2)$).
The combined forces on the paper plane are $$ F_x = -\frac{mg}{v_t^2} vv_y, \qquad F_y = -mg + \frac{mg}{v_t^2} vv_x.\tag{3} $$
By combining Newton's law and the total force (3) we get the following four equations that describe the motino of the paper plane: $$ \left. \begin{aligned} \frac{dv_x}{dt} & = -\frac{g}{v_t^2}\, vv_y\\ \frac{dv_y}{dt} & = -g + \frac{g}{v_t^2}\, vv_x\\ \frac{dx}{dt} & = v_x\\ \frac{dy}{dt} & = v_y. \end{aligned} \right\}\tag{4} $$ Here, as before, $v$ is an abbreviation for $\surd(v_x^2+v_y^2)$.
Adapt the code from the first few Phugoid notebooks to solve the differential equations (4), using as parameters $g=9.8$ m/sec, and as trim velocity $v_t=4.9$ m/sec over a time interval of about 8 seconds.
Initial conditions: we choose everything zero, i.e. suppose you start with $v_x=v_y=0$, what solution do you get?
Use both Euler's method and the Midpoint method. Plot approximate solutions in one figure.
Begin by importing numpy, matplotlib, etc...
In [3]:
from math import sqrt
import numpy
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
Set the physical parameters. Adjust $g$ to your planet and $v_t$ to your glider.
In [4]:
# Physical model parameters:
g = 9.8 # gravity in m s^{-2}
vt = 4.9 # trim velocity in m s^{-1}, value from the paper plane model
Define the function f() from the right hand side in the diffeqs (4), and define the two functions that perform one Euler step, and one midpoint step.
In [5]:
def f(u):
"""Returns the right-hand side of the phugoid system of equations.
input: u, an array containing the solution at time n.
output: dudt=f(u), an array containing the derivatives
"""
vx,vy = u[0],u[1]
x,y = u[2],u[3]
v = sqrt(vx*vx+vy*vy)
return numpy.array([-g*v*vy/(vt*vt),
-g+g*v*vx/(vt*vt),
vx,
vy ])
def euler_step(u, f, dt):
"""Returns the solution at the next time-step using Euler's method.
Parameters
----------
u : array of float
solution at the previous time-step.
f : function
function to compute the right hand-side of the system of equation.
dt : float
time-increment.
Returns
-------
u_{n+1} : array of float
approximate solution at the next time step.
"""
return u + dt * f(u)
def midpoint_step(u, f, dt):
"""Returns the solution at the next time-step using the midpoint method.
Parameters
u : array of float
solution at the previous time-step.
f : function
function to compute the right hand-side of the system of equation.
dt : float
time-increment.
Returns
u_{n+1} : array of float
approximate solution at the next time step.
"""
umid = u+dt*f(u)/2
return u+dt*f(umid)
Choose simulation parameters: time interval, number of time steps. Then choose the initial conditions.
It turns out that we will need about 1000 time steps to keep Euler's method from running into trouble.
In [6]:
### Choose time interval and approximation parameters
T = 8.0 # final time
N = 1000 # number of time-steps
dt = T/N # time increment
t = numpy.linspace(0.0, T, N+1) # time discretization
### set initial conditions ###
vx0 = 0.0 # start at the trim velocity (or add a delta)
vy0 = 0.0 # initial vertical velocity of trajectory
x0 = 0.0
y0 = 0.0
# initialize the array containing the solution for each time-step
u = numpy.empty((N+1, 4))
v = numpy.empty((N+1, 4))
u[0] = numpy.array([vx0, vy0, x0, y0]) # fill 1st element with initial values
v[0] = numpy.array([vx0, vy0, x0, y0]) # u&v have same initial values
# time loop - Euler method
for n in range(N):
u[n+1] = euler_step(u[n], f, dt)
v[n+1] = midpoint_step(v[n], f, dt)
In order to plot the path of the glider, we need the location (x, y) with respect to time. That information is already contained in our NumPy array containing the solution; we just need to pluck it out.
Make sure you understand the indices to u, below, and the use of the colon notation. If any of it is confusing, read the Python documentation on Indexing.
In [7]:
# get the glider's position with respect to the time
xEuler = u[:,2]
yEuler = u[:,3]
xMidp = v[:,2]
yMidp = v[:,3]
Time to plot the path of the glider and get the distance travelled! Which of the two methods does a better job at preserving energy?
In [11]:
# visualization of the path
plt.figure(figsize=(12,5))
plt.grid(True)
plt.xlabel(r'x', fontsize=18)
plt.ylabel(r'y', fontsize=18)
plt.title('GLIDER TRAJECTORY: flight time T = %.2f, Euler & midpoint time step dt= %4f' %(T,dt), fontsize=16)
plt.plot(xEuler, yEuler, 'r', label="Euler")
plt.plot(xMidp, yMidp, 'g', label="Midpoint")
plt.ylim(-4,1.0)
plt.legend();
In [1]:
from IPython.core.display import HTML
css_file = '../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())
Out[1]: