CoCalc Public FilesFull_Phugoid-01.ipynb
Author: Ed Johnson
Views : 64
Description: Jupyter notebook Full_Phugoid-01.ipynb
Compute Environment: Ubuntu 18.04 (Deprecated)
In [2]:
from math import sin, cos, log, ceil
import numpy
from matplotlib import pyplot
%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16

In [3]:
# model parameters:
g=9.8   # gravity in m s^{-2}
v_t = 30.0 # trim velocity in m s^{-1}
C_D = 1/40   # drag coefficient --- or D/L if C_L=1
C_L = 1      # for convenience, use C_L = 1

### set initial conditions ###
v0 = v_t      # start at the trim velocity (or add a delta)
theta0 = 0    # initial angle of trajectory
x0 = 0      # horizontal position is arbitary
y0 = 1000   # initial altitude

In [4]:
def f(u):
"""Returns the righ-hand side of the phugoid system of equations.
Parameters
----------
u : array of float
array containing the solution at time n.

Returns
-------
dudt : array of float
array containing the RHS given u.
"""

v = u[0]
theta = u[1]
x = u[2]
y = u[3]
return numpy.array([-g*sin(theta) - C_D/C_L*g/v_t**2*v**2,
-g*cos(theta)/v + g/v_t**2*v,
v*cos(theta),
v*sin(theta)])

In [5]:
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_plus_1 : array of float
approximate solution at the next time step.
"""

return u + dt * f(u)

In [6]:
T = 100         # final time
dt = 0.1        # time increment
N = int(T/dt) + 1        # number of time-steps

# initialise the array containing the solution for each time-step
u = numpy.empty((N, 4))
u[0] = numpy.array([v0, theta0, x0, y0])  # fill 1st element with initial values

# time loop - Euler method
for n in range(N-1):

u[n + 1] = euler_step(u[n], f, dt)

In [7]:
# get the glider's position with respect to the time
x = u[:,2]
y = u[:,3]

In [8]:
# viisualisation of the path
pyplot.figure(figsize=(8,6))
pyplot.grid(True)
pyplot.xlabel(r'x', fontsize=18)
pyplot.ylabel(r'y', fontsize=18)
pyplot.title('Glider trajectory, flight time = %.2f' % T, fontsize=18)
pyplot.plot(x, y, 'k-', lw=2);

In [9]:
dt_values = numpy.array([0.1, 0.05, 0.01, 0.005, 0.001])

u_values = numpy.empty_like(dt_values, dtype=numpy.ndarray)

for i, dt in enumerate(dt_values):

N = int(T/dt) + 1    # number of time-steps

### discretize the time t ###
t = numpy.linspace(0.0, T, N)

# initialize the array containing the solution for each time-step
u = numpy.empty((N, 4))
u[0] = numpy.array([v0, theta0, x0, y0])

# time loop
for n in range(N-1):

u[n+1] = euler_step(u[n], f, dt)   ### call euler_step() ###

# store the value of u related to one grid
u_values[i] = u

In [10]:
def get_diffgrid(u_current, u_fine, dt):
"""Returns the difference between one grid and the fine one using L-1 norm.

Parameters
----------
u_current : array of float
solution on the current grid.
u_finest : array of float
solution on the fine grid.
dt : float
time-increment on the current grid.

Returns
-------
diffgrid : float
difference computed in the L-1 norm.
"""

N_current = len(u_current[:,0])
N_fine = len(u_fine[:,0])

grid_size_ratio = ceil(N_fine/N_current)

diffgrid = dt * numpy.sum( numpy.abs(\
u_current[:,2]- u_fine[::grid_size_ratio,2]))

return diffgrid

In [11]:
# compute difference between one grid solution and the finest one
diffgrid = numpy.empty_like(dt_values)

for i, dt in enumerate(dt_values):
print('dt = {}'.format(dt))

### call the function get_diffgrid() ###
diffgrid[i] = get_diffgrid(u_values[i], u_values[-1], dt)

dt = 0.1 dt = 0.05 dt = 0.01 dt = 0.005 dt = 0.001
In [12]:
# log-log plot of the grid differences
pyplot.figure(figsize=(6,6))
pyplot.grid(True)
pyplot.xlabel('$\Delta t$', fontsize=18)
pyplot.ylabel('$L_1$-norm of the grid differences', fontsize=18)
pyplot.axis('equal')
pyplot.loglog(dt_values[:-1], diffgrid[:-1], color='k', ls='-', lw=2, marker='o')

[<matplotlib.lines.Line2D at 0x7fa9d0317978>]
In [13]:
r = 2
h = 0.001

dt_values2 = numpy.array([h, r*h, r**2*h])

u_values2 = numpy.empty_like(dt_values2, dtype=numpy.ndarray)

diffgrid2 = numpy.empty(2)

for i, dt in enumerate(dt_values2):

N = int(T/dt) + 1   # number of time-steps

### discretize the time t ###
t = numpy.linspace(0.0, T, N)

# initialize the array containing the solution for each time-step
u = numpy.empty((N, 4))
u[0] = numpy.array([v0, theta0, x0, y0])

# time loop
for n in range(N-1):

u[n+1] = euler_step(u[n], f, dt)         ### call euler_step() ###

# store the value of u related to one grid
u_values2[i] = u

#calculate f2 - f1
diffgrid2[0] = get_diffgrid(u_values2[1], u_values2[0], dt_values2[1])

#calculate f3 - f2
diffgrid2[1] = get_diffgrid(u_values2[2], u_values2[1], dt_values2[2])

# calculate the order of convergence
p = (log(diffgrid2[1]) - log(diffgrid2[0])) / log(r)

print('The order of convergence is p = {:.3f}'.format(p));

The order of convergence is p = 1.014
In [14]:
from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, "r").read())

--------------------------------------------------------------------------- FileNotFoundError Traceback (most recent call last) <ipython-input-14-fde599084f9c> in <module>() 1 from IPython.core.display import HTML 2 css_file = '../../styles/numericalmoocstyle.css' ----> 3 HTML(open(css_file, "r").read()) FileNotFoundError: [Errno 2] No such file or directory: '../../styles/numericalmoocstyle.css' 
In [ ]: