ESE - Numerical Methods I

Overview

The objective of this course is to understand what a numerical method is and what it can be used for. We will study together a variety of numerical methods that are used to solve equations and analyse experimental data. We will approach each topic from a practical perspective, covering for each motivation, method, and application.

The course will make extensive use of Python, specifically, we will learn how to use SciPy, Matplotlib, and Pandas. The course will be driven using iPython Notebooks which you will be able to run and edit locally on your computer (no worries, you are an expert on iPython Notebooks as you will have taken G. Gorman's excellent 'Programming for Geoscientists' class!).

We will have eight sessions, every Monday from 2:00 to 4:45pm. Every lecture we will take a coffee break (we will try to break after 2pm and before 3pm). Each session will blend theory and practice. I will introduce concepts, then we will use python to understand how the methods work and how and when we can apply them. We will complete all the work during class, and there will be only one (!) final assignment.

All course material will be posted online through ESESIS. If you have any questions after class, please contact me by e-mail (apaluszn@imperial.ac.uk). Office hours will be Tuesday morning 10:00 - 11:00am (at RSM 2.48).

What is a numerical method?

Numerical methods are used to analyse numerical data. In many cases, methods are used to find solutions that cannot be computed "analytically" (by using ready-made equations), or to analyse "discrete" data such as the output of a lab experiment or the observations over time of some property. Numerical methods are techniques used to approximate solutions, so they rarely provide an exact solution - but a "good" estimate of the solution. In many cases, it is the only way of obtaining a solution!

In dealing with numerical methods, the main concerns are

  • performance (how long solutions take to compute) and
  • accuracy of the solution (what is the error in our approximation).

Specific numerical methods can be used to interpolate data, compute approximate derivatives, integrate functions, among others. There are special libraries dedicated to implement these "solutions", in Python the library is called SciPy. In this class you will learn

  • what numerical methods are available and what are their advantages and disadvantages
  • how to use SciPy to solve a variety of problems numerically
  • how to evaluate the numerical solutions that you obtain

Numerical Methods will be very useful in your career in almost any path you take! May it be when analysing lab data, making sense of observations in the field, or writing complex simulators; chances are you will end up using them on a daily basis.

Outline

  1. Loading and Plotting data. Numerical Errors. Examples.
  2. Roots of nonlinear equations. Methods: Incremental Search, Bisection, False Position, Secant, and Newton-Raphson. (optimize)
  3. Curve Fitting. Interpolation. (optimize, interpolate). General Least Squares Fitting. Fitting functions: polynomial, exponential, non-linear. (leastsq)
  4. Solving simultaneous linear equations. Gaussian elimination. Pivoting. Tri-diagonal systems. Norms and Condition Number. (linalg)
  5. Numerical Integration - Quadrature rules: rectangle, trapezoid, Simpson’s Rule. (integrate) - with Practical 1 and Practical 2

Tools: (Python) SciPy and MatPlotLib

We will learn how numerical methods work, and what they are used for. We will use the Python SciPy library to execute implementations of the numerical methods that we discuss (in other words: SciPy is a large collection of programmes that execute numerical methods, like a super-powerful calculator!). We will use the Python library MatPlotLib to plot the data we analyse.

We will use NumPy as an underlying library to define some of the objects we use, such as polynomials (poly1d) and arrays.

From the SciPy Manual: "SciPy is a collection of mathematical algorithms and convenience functions built on the Numpy extension of Python. It adds significant power to the interactive Python session by providing the user with high-level commands and classes for manipulating and visualizing data. With SciPy an interactive Python session becomes a data-processing and system-prototyping environment rivaling sytems such as MATLAB, IDL, Octave, R-Lab, and SciLab." An important source of information on SciPy is the Reference Manual and its Tutorials, see: SciPy Tutorial.

From the MatPlotLib website: "matplotlib is a python 2D plotting library which produces publication quality figures in a variety of hardcopy formats and interactive ". In particular, we will be using PyPlot, a collection of specific commands that make plotting simpler (see PyPlot Tutorial).

In this course, we assume that you already have some basic Python skills (see: How to install Python). However, for most exercises example code will be provided. The two exercise sessions on Week 4 and 8 will be dedicated to practising what we learn about Numerical Methods, and learning how to apply them to analyse data.

Setting up NumPy, SciPy, and MatPlotLib


In [1]:
%pylab inline
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt


Populating the interactive namespace from numpy and matplotlib

SciPy packages

SciPy is organised into packages. We will only review the following:

- io
- interpolate
- optimize
- linalg

Arrays (numpy)

Arrays can be used to store data. For example:


In [2]:
a=np.array([100,3,-5,0.15])
print a


[ 100.      3.     -5.      0.15]

In [3]:
b=np.array([40,13,-5,7,9])
print b


[40 13 -5  7  9]

In [4]:
plt.plot(a)
plt.ylabel('some numbers')
plt.show()


The command 'plot' can be used to plot data onto the screen. Its documentation can be found here: http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.plot. In general, plot takes either one argument (data points directly), or two arguments (two arrays of datapoints - such as x,y), and it takes successively parameters that define the colors and symbols to be used in the plot. So, adding 'r' will make the data red, adding 'b' will make it blue, adding '-' will make sure its a line, adding 'o' will make it plot as circles, and so on. And if you want to plot multiple datasets on one plot, you just continue to add information onto the command line.


In [5]:
plt.plot(a, 'ro', b, 'b-')
plt.show()



In [6]:
c=np.arange(10)
print c


[0 1 2 3 4 5 6 7 8 9]

In [7]:
d=np.arange(0,10,0.5)
print d


[ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5  7.
  7.5  8.   8.5  9.   9.5]

In [8]:
e=np.random.rand(20)
print e


[ 0.27177748  0.30166386  0.9577788   0.74645191  0.25017223  0.43295796
  0.56341258  0.96661882  0.31703513  0.87176957  0.88314891  0.16863722
  0.91600137  0.07045681  0.77285472  0.62513014  0.10897923  0.61134761
  0.23990097  0.82380901]

In [9]:
plt.plot(d, e, '-', np.random.randint(0,10,100), np.random.rand(100), 'ro')
plt.ylabel('y')
plt.xlabel('x')
plt.show()


Polynomials (numpy.poly1d)

Data and functions can also be respresented using polynomials. A simple way of expressing a polynomial using Python is to use the poly1d class (from Numpy). This class constructs a polynomial of one variable, the poly1d polynomials can be initialised using a list of the polynomial coefficients. The polynomial can then be added, substracted, multiplied, differentiated, integrated, evaluated, etc. using in-built functions. Poly1D objects can even be printed to screen like a polynomial:

Poly1D Documentation can be found here: poly1d

A polynomial can be defined using its coefficients:


In [10]:
p = poly1d([10,-2,3])
print p


    2
10 x - 2 x + 3

In [11]:
q = poly1d([3,4,5],True)
print q


   3      2
1 x - 12 x + 47 x - 60

Polynomials can be 'operated', they can be summed, substracted, multiplied with each other.


In [12]:
print p*q


    5       4       3       2
10 x - 122 x + 497 x - 730 x + 261 x - 180

In [13]:
print p+q


   3     2
1 x - 2 x + 45 x - 57

In [14]:
print p.deriv()


 
20 x - 2

In [15]:
print p([5])


[243]

In [16]:
print p([5,6,7,8,9,10])


[243 351 479 627 795 983]

In [17]:
p.r


Out[17]:
array([ 0.1+0.53851648j,  0.1-0.53851648j])

In [18]:
p.c


Out[18]:
array([10, -2,  3])

In [19]:
x=np.arange(0,10,0.5)
print x


[ 0.   0.5  1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5  7.
  7.5  8.   8.5  9.   9.5]

In [20]:
p(x)


Out[20]:
array([   3. ,    4.5,   11. ,   22.5,   39. ,   60.5,   87. ,  118.5,
        155. ,  196.5,  243. ,  294.5,  351. ,  412.5,  479. ,  550.5,
        627. ,  708.5,  795. ,  886.5])

Furthermore, we can plot our polynomial evaluated over a certain array, using the MatPlotLib library (that we imported earlier on) as:


In [21]:
plt.plot(x, p(x))
plt.ylabel('y')
plt.xlabel('x')
plt.xlim(0,10)
plt.ylim(0,1000)
plt.show()



In [21]: