Lesson 1

In this lesson, you will write your first python code.

Numbers

>>> 2 + 2
4
>>> 50 - 5*6
20
>>> (50 - 5.0*6) / 4
5.0
>>> 8 / 5.0
1.6
>>> 17 / 3  # int / int -> int
5
>>> 17 / 3.0  # int / float -> float
5.666666666666667
>>> 17 // 3.0  # explicit floor division discards the fractional part
5.0
>>> 17 % 3  # the % operator returns the remainder of the division
2
>>> 5 * 3 + 2  # result * divisor + remainder
17

The equal sign (=) is used to assign a value to a variable. No result is displayed

>>> width = 20
>>> height = 5 * 9
>>> width * height
900

strings

>>> 'spam eggs'  # single quotes
'spam eggs'
>>> 'doesn\'t'  # use \' to escape the single quote...
"doesn't"
>>> "doesn't"  # ...or use double quotes instead
"doesn't"
>>> '"Yes," he said.'
'"Yes," he said.'
>>> "\"Yes,\" he said."
'"Yes," he said.'
>>> '"Isn\'t," she said.'
'"Isn\'t," she said.'

>>> '"Isn\'t," she said.'
'"Isn\'t," she said.'
>>> print '"Isn\'t," she said.'
"Isn't," she said.
>>> s = 'First line.\nSecond line.'  # \n means newline
>>> s  # without print(), \n is included in the output
'First line.\nSecond line.'
>>> print s  # with print, \n produces a new line
First line.
Second line.

(+) operator is available for strings

>>> s = "Hello world."
>>> t = " Said Python"
>>> s+t
'Hello world. Said Python'

In fact, a string is a list object.

lists

Rather than using arrays of a given dimension and size, with list object, you can insert new elements or remove them, iterate, ...

>>> L = [2,3,4,-5,9,10]
>>> L[0]
2
>>> L[:]
[2, 3, 4, -5, 9, 10]
>>> L[2:-2]
[4, -5]

Lists also supports operations like concatenation:

>>> A = [21, 22, 23, 24]
>>> L+A
[2, 3, 4, -5, 9, 10, 21, 22, 23, 24]

Insertion can be done using the append method

>>> L.append(100)
>>> print L
[2, 3, 4, -5, 9, 10, 100]

Getting the index of an element can be done using the index method

>>> L.index(10)
5

Removing an element can be done using the pop or remove methods

>>> L.pop(L.index(10))
>>> L
[2, 3, 4, -5, 9, 100]
>>> L.remove(4)
>>> L
[2, 3, -5, 9, 100]

The list’s length can be accessed using

>>> len(L)
5

numpy

The first thing to do is to import the numpy package:

import numpy as np

If you wan to see what are the functions and modules available in numpy, write np. and press on TAB. The result is:

>>> np.
np.ALLOW_THREADS           np.iscomplex
np.BUFSIZE                 np.iscomplexobj
np.CLIP                    np.isfinite
np.ComplexWarning          np.isfortran
np.DataSource              np.isinf
np.ERR_CALL                np.isnan
np.ERR_DEFAULT             np.isneginf
np.ERR_DEFAULT2            np.isposinf
np.ERR_IGNORE              np.isreal
np.ERR_LOG                 np.isrealobj
np.ERR_PRINT               np.isscalar
np.ERR_RAISE               np.issctype
np.ERR_WARN                np.issubclass_
np.FLOATING_POINT_SUPPORT  np.issubdtype
np.FPE_DIVIDEBYZERO        np.issubsctype
np.FPE_INVALID             np.iterable
np.FPE_OVERFLOW            np.ix_
np.FPE_UNDERFLOW           np.kaiser
np.False_                  np.kron
np.Inf                     np.ldexp
np.Infinity                np.left_shift
np.MAXDIMS                 np.less
np.MachAr                  np.less_equal
np.NAN                     np.lexsort
np.NINF                    np.lib
np.NZERO                   np.linalg
np.NaN                     np.linspace
np.PINF                    np.little_endian
np.PZERO                   np.load
np.PackageLoader           np.loads
np.RAISE                   np.loadtxt
np.RankWarning             np.log
np.SHIFT_DIVIDEBYZERO      np.log10
np.SHIFT_INVALID           np.log1p
np.SHIFT_OVERFLOW          np.log2
np.SHIFT_UNDERFLOW         np.logaddexp
np.ScalarType              np.logaddexp2
np.Tester                  np.logical_and
np.True_                   np.logical_not
np.UFUNC_BUFSIZE_DEFAULT   np.logical_or
np.UFUNC_PYVALS_NAME       np.logical_xor
np.WRAP                    np.logspace
np.abs                     np.long
np.absolute                np.longcomplex
np.add                     np.longdouble
np.add_docstring           np.longfloat
np.add_newdoc              np.longlong
np.add_newdocs             np.lookfor
np.alen                    np.ma
np.all                     np.mafromtxt
np.allclose                np.mask_indices
np.alltrue                 np.mat
np.alterdot                np.math
np.amax                    np.matrix
np.amin                    np.matrixlib
np.angle                   np.max
np.any                     np.maximum
np.append                  np.maximum_sctype
np.apply_along_axis        np.may_share_memory
np.apply_over_axes         np.mean
np.arange                  np.median
np.arccos                  np.memmap
np.arccosh                 np.meshgrid
np.arcsin                  np.mgrid
np.arcsinh                 np.min
np.arctan                  np.min_scalar_type
np.arctan2                 np.minimum
np.arctanh                 np.mintypecode
np.argmax                  np.mirr
np.argmin                  np.mod
np.argsort                 np.modf
np.argwhere                np.msort
np.around                  np.multiply
np.array                   np.nan
np.array2string            np.nan_to_num
np.array_equal             np.nanargmax
np.array_equiv             np.nanargmin
np.array_repr              np.nanmax
np.array_split             np.nanmin
np.array_str               np.nansum
np.asanyarray              np.nbytes
np.asarray                 np.ndarray
np.asarray_chkfinite       np.ndenumerate
np.ascontiguousarray       np.ndfromtxt
np.asfarray                np.ndim
np.asfortranarray          np.ndindex
np.asmatrix                np.nditer
np.asscalar                np.negative
np.atleast_1d              np.nested_iters
np.atleast_2d              np.newaxis
np.atleast_3d              np.newbuffer
np.average                 np.nextafter
np.bartlett                np.nonzero
np.base_repr               np.not_equal
np.bench                   np.nper
np.binary_repr             np.npv
np.bincount                np.number
np.bitwise_and             np.obj2sctype
np.bitwise_not             np.object
np.bitwise_or              np.object0
np.bitwise_xor             np.object_
np.blackman                np.ogrid
np.bmat                    np.ones
np.bool                    np.ones_like
np.bool8                   np.outer
np.bool_                   np.packbits
np.broadcast               np.percentile
np.broadcast_arrays        np.pi
np.byte                    np.piecewise
np.byte_bounds             np.pkgload
np.bytes_                  np.place
np.c_                      np.pmt
np.can_cast                np.poly
np.cast                    np.poly1d
np.cdouble                 np.polyadd
np.ceil                    np.polyder
np.cfloat                  np.polydiv
np.char                    np.polyfit
np.character               np.polyint
np.chararray               np.polymul
np.choose                  np.polynomial
np.clip                    np.polysub
np.clongdouble             np.polyval
np.clongfloat              np.power
np.column_stack            np.ppmt
np.common_type             np.prod
np.compare_chararrays      np.product
np.compat                  np.promote_types
np.complex                 np.ptp
np.complex128              np.put
np.complex256              np.putmask
np.complex64               np.pv
np.complex_                np.r_
np.complexfloating         np.rad2deg
np.compress                np.radians
np.concatenate             np.random
np.conj                    np.rank
np.conjugate               np.rate
np.convolve                np.ravel
np.copy                    np.ravel_multi_index
np.copysign                np.real
np.core                    np.real_if_close
np.corrcoef                np.rec
np.correlate               np.recarray
np.cos                     np.recfromcsv
np.cosh                    np.recfromtxt
np.count_nonzero           np.reciprocal
np.cov                     np.record
np.cross                   np.remainder
np.csingle                 np.repeat
np.ctypeslib               np.require
np.cumprod                 np.reshape
np.cumproduct              np.resize
np.cumsum                  np.restoredot
np.datetime64              np.result_type
np.datetime_               np.right_shift
np.datetime_data           np.rint
np.deg2rad                 np.roll
np.degrees                 np.rollaxis
np.delete                  np.roots
np.deprecate               np.rot90
np.deprecate_with_doc      np.round
np.diag                    np.round_
np.diag_indices            np.row_stack
np.diag_indices_from       np.s_
np.diagflat                np.safe_eval
np.diagonal                np.save
np.diff                    np.savetxt
np.digitize                np.savez
np.disp                    np.savez_compressed
np.divide                  np.sctype2char
np.dot                     np.sctypeDict
np.double                  np.sctypeNA
np.dsplit                  np.sctypes
np.dstack                  np.searchsorted
np.dtype                   np.select
np.e                       np.set_numeric_ops
np.ediff1d                 np.set_printoptions
np.einsum                  np.set_string_function
np.emath                   np.setbufsize
np.empty                   np.setdiff1d
np.empty_like              np.seterr
np.equal                   np.seterrcall
np.errstate                np.seterrobj
np.exp                     np.setxor1d
np.exp2                    np.shape
np.expand_dims             np.short
np.expm1                   np.show_config
np.extract                 np.sign
np.eye                     np.signbit
np.fabs                    np.signedinteger
np.fastCopyAndTranspose    np.sin
np.fft                     np.sinc
np.fill_diagonal           np.single
np.find_common_type        np.singlecomplex
np.finfo                   np.sinh
np.fix                     np.size
np.flatiter                np.sometrue
np.flatnonzero             np.sort
np.flexible                np.sort_complex
np.fliplr                  np.source
np.flipud                  np.spacing
np.float                   np.split
np.float128                np.sqrt
np.float16                 np.square
np.float32                 np.squeeze
np.float64                 np.std
np.float_                  np.str
np.floating                np.str_
np.floor                   np.string0
np.floor_divide            np.string_
np.fmax                    np.subtract
np.fmin                    np.sum
np.fmod                    np.swapaxes
np.format_parser           np.take
np.frexp                   np.tan
np.frombuffer              np.tanh
np.fromfile                np.tensordot
np.fromfunction            np.test
np.fromiter                np.testing
np.frompyfunc              np.tile
np.fromregex               np.timedelta64
np.fromstring              np.timedelta_
np.fv                      np.timeinteger
np.generic                 np.trace
np.genfromtxt              np.transpose
np.get_array_wrap          np.trapz
np.get_include             np.tri
np.get_numarray_include    np.tril
np.get_printoptions        np.tril_indices
np.getbuffer               np.tril_indices_from
np.getbufsize              np.trim_zeros
np.geterr                  np.triu
np.geterrcall              np.triu_indices
np.geterrobj               np.triu_indices_from
np.gradient                np.true_divide
np.greater                 np.trunc
np.greater_equal           np.typeDict
np.half                    np.typeNA
np.hamming                 np.typecodes
np.hanning                 np.typename
np.histogram               np.ubyte
np.histogram2d             np.ufunc
np.histogramdd             np.uint
np.hsplit                  np.uint0
np.hstack                  np.uint16
np.hypot                   np.uint32
np.i0                      np.uint64
np.identity                np.uint8
np.iinfo                   np.uintc
np.imag                    np.uintp
np.in1d                    np.ulonglong
np.index_exp               np.unicode
np.indices                 np.unicode0
np.inexact                 np.unicode_
np.inf                     np.union1d
np.info                    np.unique
np.infty                   np.unpackbits
np.inner                   np.unravel_index
np.insert                  np.unsignedinteger
np.int                     np.unwrap
np.int0                    np.ushort
np.int16                   np.vander
np.int32                   np.var
np.int64                   np.vdot
np.int8                    np.vectorize
np.int_                    np.version
np.int_asbuffer            np.void
np.intc                    np.void0
np.integer                 np.vsplit
np.interp                  np.vstack
np.intersect1d             np.where
np.intp                    np.who
np.invert                  np.zeros
np.ipmt                    np.zeros_like
np.irr

Let’s take the last list and convert it to a numpy array:

>>> x = np.asarray(L)
>>> print x
array([  2,   3,  -5,   9, 100])

You notice that the result is no longer a list, but an array

Arrays can be reshaped using:

>>> L = range(0,12) # create a list or integers from 0 to 11
>>> x = np.asarray(L)
>>> print x
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
>>> print x.shape # access the shape of an array
(12,)
>>> y = x.reshape((3,4)) # reshape
>>> print y
array([[ 0,  1,  2,  3],
[ 4,  5,  6,  7],
[ 8,  9, 10, 11]])
>>> print y.shape
(3, 4)
>>> x.ndim, y.ndim
(1,2)
>>> x.dtype.name
'int64'
>>> x.itemsize, y.itemsize
(8, 8)
>>> x.size, y.size
(12,12)
>>> type(x)
numpy.ndarray

>>> a = np.array(1,2,3,4)    # WRONG
>>> a = np.array([1,2,3,4])  # RIGHT

You can specify the type at the creation, or convert your array:

>>> x = np.array( [ [1,2], [3,4] ], dtype=np.float )
>>> c = np.array( x, dtype=complex )
array([[ 1.+0.j,  2.+0.j],
[ 3.+0.j,  4.+0.j]])

Loop over an array (or list) can be done in different ways:

>>> L = range(0,12) # create a list or integers from 0 to 11
>>> X = np.asarray(L)
>>> # method 1: simple loop
>>> n = X.shape[0]
>>> for i in range(0,n):
>>>     x = X[i]
>>>     print x
>>> # method 2: using enumerate
>>> for i,x in enumerate(X):
>>>     x = X[i]
>>>     print x
>>> # method 3: using intrinsic iterator
>>> for x in X:
>>>     print x

Sometimes you need to loop over two lists (or more), this can be done using zip function:

>>> L = range(0,12) # create a list or integers from 0 to 11
>>> X = np.asarray(L)
>>> Y = X**2 # taking the square of X
>>> for (x,y) in zip(X,Y):
>>>     print x**2-y

Another way of computing the square of X is to use a list:

>>> Z = np.asarray([x**2 for x in X if np.mod(x,2)==0]) # compute squares of even values

From Matlab/Octave to numpy

Arithmetic operators

Description Matlab/Octave Python
Assignment; defining a number a=1; b=2; a=1; b=1
Addition a+b; a+b or add(a,b)
Substruction a-b; a-b or substract(a,b)
Multiplication a*b; a*b or multiply(a,b)
Division a/b; a/b or devide(a,b)
Power a.^b; a**b or pow(a,b)
Remainder rem(a,b); a%b or remainder(a,b)
In place operation a+=1; a+=b or add(a,b,a)
Factorial n! factorial(a);  

More informations can be found at http://mathesaurus.sourceforge.net/matlab-numpy.html

Matplotlib

Let us consider the function

f(x) = sin(2 \pi x), \forall x \in [0,1]

f can be defined in python using:

>>> from numpy import pi, sin, cos
>>> # method 1
>>> def f(x):
>>>     return sin(2*pi*x)
>>> # method 2
>>> f = lambda x: sin(2*pi*x)

The next example shows how to plot f

>>> x = np.linspace(0.,1.,100)
>>> y = f(x)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x,y)
>>> plt.show()
_images/ex1.png

Let’s now take a 2D function

f(x,y) = sin(2 \pi x) sin(4 \pi y), \forall x,y \in [0,1]

f can be defined in python using:

>>> # method 1
>>> def f(x,y):
>>>     return sin(2*pi*x) * sin(4*pi*y)
>>> # method 2
>>> f = lambda x,y: sin(2*pi*x) * sin(4*pi*y)

The next example shows how to plot f

>>> x = np.linspace(0.,1.,100)
>>> y = np.linspace(0.,1.,100)
>>> X,Y = np.meshgrid(x,y)
>>> Z = f(X,Y)
>>> plt.contourf(X,Y,Z)
>>> plt.colorbar()
>>> plt.show()
_images/ex2.png

The following example shows how to plot 3D surface

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
#!/usr/bin/python

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')
X = np.arange(-5, 5, 0.25)
Y = np.arange(-5, 5, 0.25)
X, Y = np.meshgrid(X, Y)
R = np.sqrt(X**2 + Y**2)
Z = np.sin(R)
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False)
ax.set_zlim(-1.01, 1.01)

ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()
_images/ex3.png

scipy

SciPy (pronounced “Sigh Pie”) is open-source software for mathematics, science, and engineering.

Here is a list of some scipy modules

  • Discrete Fourier transforms (scipy.fftpack)
  • Integration and ODEs (scipy.integrate)
  • Interpolation (scipy.interpolate)
  • Input and output (scipy.io)
  • Linear algebra (scipy.linalg)
  • Optimization and root finding (scipy.optimize)
  • Sparse matrices (scipy.sparse)
  • Sparse linear algebra (scipy.sparse.linalg)
  • Compressed Sparse Graph Routines (scipy.sparse.csgraph)
  • Spatial algorithms and data structures (scipy.spatial)
  • Special functions (scipy.special)

We will focus on the use of sparse matrices and their linear solvers.

scipy.sparse

SciPy 2-D sparse matrix package for numeric data.

Available constructors are

  • bsr_matrix(arg1[, shape, dtype, copy, blocksize]) Block Sparse Row matrix
  • coo_matrix(arg1[, shape, dtype, copy]) A sparse matrix in COOrdinate format.
  • csc_matrix(arg1[, shape, dtype, copy]) Compressed Sparse Column matrix
  • csr_matrix(arg1[, shape, dtype, copy]) Compressed Sparse Row matrix
  • dia_matrix(arg1[, shape, dtype, copy]) Sparse matrix with DIAgonal storage
  • dok_matrix(arg1[, shape, dtype, copy]) Dictionary Of Keys based sparse matrix.
  • lil_matrix(arg1[, shape, dtype, copy]) Row-based linked list sparse matrix
>>> import numpy as np
>>> from scipy.sparse import csr_matrix
>>> A = csr_matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])
>>> v = np.array([1, 0, -1])
>>> A.dot(v)
array([ 1, -3, -1], dtype=int64)

Exercise

Let us consider the following advection problem,

\partial_t u + \mathbf{div} (au) = 0, ~ ~ x \in \Omega=(0,2 \pi)
\\
u(t=0,x) = u_0(x), ~ ~ ~ \forall x \in \Omega

where a is a constant or any real function.

Let us denote A^n the resulting finite difference matrix at time n. The spatial grid is defined by

x_i = i  \Delta x, ~ ~ ~ \forall i \in \{ 0, ... , N  \}

Upwind finite difference matrix

We consider an upwind finite difference method for the spatial discretization.

\frac{d}{dt}\mathbf{u} = A(t) \mathbf{u}, ~ ~ ~ \forall t > 0

The \theta-time scheme can be written as

\frac{\mathbf{u}^{n+1} - \mathbf{u}^{n}}{\Delta t} = \theta A^{n+1} \mathbf{u}^{n+1} + (1-\theta) A^{n} \mathbf{u}^{n}, ~ ~ ~ \forall t > 0

Reordering the last equation, leads to

\{ 1 - \theta \Delta t A^{n+1} \} \mathbf{u}^{n+1} = \{ 1 + (1-\theta) \Delta t A^{n} \} \mathbf{u}^{n}

In the case of dirichlet boundary, the vector \mathbf{u}^{n} is

\mathbf{u}^{n} = \{ u_1, u_2, ..., u_{N-1} \}

The finite difference matrix is therefor

A^{n}_{i,i}   =   \frac{a_i^{n}}{\Delta x}
\\
A^{n}_{i,i-1} = - \frac{a_{i-1}^{n}}{\Delta x}
\\
A^{n}_{i,j}   = 0, ~ ~ ~ \mbox{if} ~ ~ j \notin \{ i-1, i \}

In the case of periodic boundary, the vector \mathbf{u}^{n} is

\mathbf{u}^{n} = \{ u_1, u_2, ..., u_{N-1}, u_{N} \}

and the finite difference matrix is

A^{n}_{i,i}   =   \frac{a_i^{n}}{\Delta x}
\\
A^{n}_{i,i-1} = - \frac{a_{i-1}^{n}}{\Delta x}
\\
A^{n}_{1,0}   := A^{n}_{1,N}
\\
A^{n}_{i,j}   = 0, ~ ~ ~ \mbox{if} ~ ~ j \notin \{ i-1, i \}

  1. write the python script for these problems
  2. study the stability of the \theta -scheme
  3. study the consistence of the \theta -scheme
  4. check the expected analytical CFL number, for the periodic case

Dissipation and dispersion

We consider the periodic boundary condition and we assume that a is a constant. We also consider the Euler centered-explicit scheme, which can be written as

u_{j}^{n+1} = u_{j}^{n} - \frac{\lambda a}{2} (u_{j+1}^{n} - u_{j-1}^{n})

Note

any finite difference scheme can be written in the form u_{j}^{n+1} = u_{j}^{n} - \lambda(h_{j+\frac{1}{2}}^{n} - h_{j-\frac{1}{2}}^{n}). Where h_{j+\frac{1}{2}}^{n} = h(u_j^n, u_{j+1}^n), h(\cdot, \cdot) is the numerical flux.

In the case of an Euler center-explicit scheme, the numerical flux is given by h_{j+\frac{1}{2}}^{n} = \frac{a}{2} (u_{j+1}^n+u_j^n)

Let us consider the exact solution

u(t^n,x) = u_0(x - a n \Delta t), ~ \forall n \geq 0, ~ \forall x \in \Omega

We have

u(t^n,x_j) = \sum_{k} \alpha_k e^{ikj \Delta x} g_k^n, ~~~ \mbox{where} ~~ g_k = e^{-iak\Delta t}

Let us define \phi_k=k\Delta x and \lambda = \frac{\Delta t}{\Delta x}

\phi_k is the k-harmonic phase.

  1. Show that \forall n \geq 0, we have

u_j^n = \sum_{k} = \alpha_k e^{ikjh} \gamma_k^n

where \gamma_k = 1-\frac{a \Delta t}{\Delta x} i sin(k\Delta x) is the amplification coefficient for the k -harmonic

Note

The coefficient \gamma_k is the analogue of g_k for the numerical scheme.

Notice that |g_k|=1 while |\gamma_k| \leq 1.

The amplification error for the k -harmonic is defined by \epsilon_a(k) = \frac{|\gamma_k|}{|g_k|}

  1. plot the amplification error for the Euler centered explicit scheme

Let us define thee propagation velocity \frac{\omega}{k} as

\gamma_k = | \gamma_k | e^{-i\omega \Delta t} = | \gamma_k | e^{-i \frac{\omega}{k} \lambda \phi_k}

The dispersion error is defined by

\epsilon_d(k) = \frac{\omega}{ka} = \frac{\omega \Delta x}{\phi_k a}

  1. Plot the dispersion error for the Euler centered explicit scheme
  2. Do the same study for the following schemes

Lax-Friedrichs

u_{j}^{n+1} = \frac{1}{2} (u_{j+1}^n + u_{j-1}^{n}) - \frac{\lambda}{2}(u_{j+1}^n - u_{j-1}^n)

can also be defined using the numerical flux

h_{j+\frac{1}{2}}^{n} = \frac{1}{2} \{ a(u_{j+1}^{n} + u_j^n) - \lambda^{-1} (u_{j+1}^n - u_j^n) \}

Lax-Wendroff

h_{j+\frac{1}{2}}^{n} = \frac{1}{2} \{ a(u_{j+1}^{n} + u_j^n) - \lambda a^2 (u_{j+1}^n - u_j^n) \}

Decentered Explicit Euler

h_{j+\frac{1}{2}}^{n} = \frac{1}{2} \{ a(u_{j+1}^{n} + u_j^n) - |a| (u_{j+1}^n - u_j^n) \}