ESE - Numerical Methods I: Exercises



In [1]:

    
%pylab inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as scip
import scipy.optimize as scop
matplotlib.rcParams['figure.figsize'] = 14,7
matplotlib.rcParams['font.size'] = 14
matplotlib.rcParams['font.family'] = 'serif'









    



Populating the interactive namespace from numpy and matplotlib

Interpolation - 1D continued



In [2]:

    
#define/load data
x=np.array([0,1,2,3])
y=np.array([-5,-6,-1,16])

#create lin space for plotting
x_int = np.linspace(np.min(x), np.max(x), 100)

#one way to do it:
#define interpolating functions
polynomial_interpolation = scip.BarycentricInterpolator(x,y)
linear_interp1d = scip.interp1d(x, y, kind='linear')
quadratic_interp1d = scip.interp1d(x, y, kind='quadratic')
cubic_interp1d = scip.interp1d(x, y, kind='cubic')
#evaluate
y_poly = polynomial_interpolation(x_int)
y_linear = linear_interp1d(x_int)
y_quadratic = quadratic_interp1d(x_int)
y_cubic = cubic_interp1d(x_int)

#another way to do it - use 'convenience functions'
y_poly2 = scip.barycentric_interpolate(x, y, x_int)
y_pchip = scip.pchip_interpolate(x,y, x_int)

#generate the lagrangian polynomial
print 'Lagrangian polynomial'
print scip.lagrange(x,y)

# plot
plt.scatter(x, y, marker='x', s=45, label='data')
plt.plot(x_int, y_poly,'-', color='green', label='poly')
plt.plot(x_int, y_pchip,'--', color='red', label='pchip')
plt.plot(x_int, y_linear,'--', color='blue', label='linear')
plt.plot(x_int, y_quadratic,'-', color='black', label='quadratic')
plt.plot(x_int, y_cubic,'--', color='gray', label='cubic')
plt.legend(loc=2)
plt.show()









    



Lagrangian polynomial
   3
1 x - 2 x - 5

Problems

Based on data published by Kerry Emanuel in Nature, try to predict Hurricane energy for a mean sea temperature of 30 degrees Celsius! - download data (click raw)
Compare the linearly extrapolated sea level in 2050 between Aberdeen and Newlyn graphically! - download data (click raw)
Based on well production data, create a map to find regions for infill wells. - download data (click raw)



In [3]:

    
# AD 1)

# load data
data = np.loadtxt('hurricane_energy_emanual_2003.txt', skiprows=1)
x = data[:,2]
y = data[:,1]

# plot data
plt.scatter(x, y, marker='x',s=45, label='data')
plt.xlabel('Temperature [deg C]')
plt.ylabel('Hurricane Energy [-]')

# function to be fitted
def exponential(x, a, b):
    return a + b**x

# get parameters for best fit
po, pc = scop.curve_fit(exponential, x, y)

# 100 points to evaluate fitted function at (x-values)
x_fit = np.linspace(np.min(x), 30., 100)

# evaluate fit (y-values)
y_fit = exponential(x_fit, *po)

# plot fit and precition
plt.plot(x_fit, y_fit, lw=2, ls='--', label='fit')
plt.scatter(30, exponential(30, *po), marker='x', s=200, label='prediction')
plt.legend(loc=4, scatterpoints=1)
plt.show()



In [4]:

    
# AD 2)

# load data
data_a = np.loadtxt('aberdeen_msl.txt', skiprows=1)
msl_a = data_a[:,2]
year_a = data_a[:,0]
data_n = np.loadtxt('newlyn_msl.txt', skiprows=1)
msl_n = data_n[:,2]
year_n = data_n[:,0]

# plot data
plt.scatter(year_a, msl_a, marker='1', s=45, label='Aberdeen data', alpha=0.5)
plt.scatter(year_n, msl_n, marker='x', color='red', s=45, label='Newlyn data', alpha=0.5)
plt.xlabel('Year')
plt.ylabel('Mean Sea Level [m]')

# function to be fitted
def linear(x, a, b):
    return a + b*x

# get parameters for best fit
po_a, pc_a = scop.curve_fit(linear, year_a, msl_a)
po_n, pc_n = scop.curve_fit(linear, year_n, msl_n)

# 2 points to evaluate fitted function at (x-values)
year_pred = 2050
x_fit = np.array([np.min(year_a), year_pred])

# plot prediction
plt.plot(x_fit, linear(x_fit, *po_a), lw=3, ls='--', label='Aberdeen fit')
plt.scatter(year_pred, linear(year_pred, *po_a), marker='1', s=200, label='Aberdeen prediction')
plt.plot(x_fit, linear(x_fit, *po_n), color='red', lw=3, ls='--', label='Newlyn fit')
plt.scatter(year_pred, linear(year_pred, *po_n), marker='x', s=200, color='red', label='Newlyn prediction')

plt.legend(loc=4, scatterpoints=1)
plt.show()



In [5]:

    
# AD 3)

# data points
prod = np.loadtxt('prod_map.txt', skiprows=1)
x = prod[:,0]
y = prod[:,1]
z = prod[:,2]

# 1. points to evaluate interpolation at
interpolated_points = 200
xspace = np.linspace(np.min(x), np.max(x), interpolated_points)
yspace = np.linspace(np.min(y), np.max(y), interpolated_points)
X, Y = np.meshgrid(xspace, yspace)

# 2. cubic interpolation in 2D
Z = scip.griddata((x,y), z, (X,Y), method='cubic')
    
# plot
matplotlib.rcParams['figure.figsize'] = 14,7
# isodepth colors
img = plt.contourf(X, Y, Z, 7, zorder=0)
# data point positions
plt.plot(x, y, 'x', mew=2, ms=10, color='black')
plt.ylabel('Northing')
plt.xlabel('Easting')
# contour lines
ctour = plt.contour(X, Y, Z, 7, colors='black')
plt.show()



In [5]: