Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Yves Dubief, 2015. NSF for support via NSF-CBET award #1258697.

Crash Course in Python

Python is a powerful programming language that bears many similarities with matlab, C and fortran. The following crash course will take you through some of the key features we will use in ME 249.

1. Simple Mathematical Operation

The python software installed on your computer contains many libraries, the main ones used in ME 249 are:

matplotlib : http://matplotlib.org for examples of plots you can make in python.

math : https://docs.python.org/3/library/math.html for mathematical constants and functions.

numpy : http://docs.scipy.org/doc/numpy/user/index.html Library for operations on matrices and vectors.

scipy : https://www.scipy.org/ Library for fundamental mathematical operations (e.g. ODE solvers).

</ul> Loading a libray in python is done by the command import. The best practice is to take the habit to use

import [library] as [library_nickname]

1.1 Printing $\pi$

The constant $\pi$ is available from all three libraries (math, numpy, scipy). Here I don't give a shorthand name to math, to illustrate how import works with or without as.

The print command here takes a string "pi = %1.15f" where %1.15f signals the print command that a floating number of the form $$ \underbrace{x}_{1}.\underbrace{abcdefghijklmno}_{15} $$ is expected. This number is passed onto the print command via %math.pi.



In [1]:

    
import math

print("pi = %1.15f" %math.pi)
print("pi = %1.15e" %math.pi)

import numpy as np
import scipy as scp
print(" numpy pi = %1.15f" %np.pi)
print(" scipy pi = %1.15f" %scp.pi)









    



pi = 3.141592653589793
pi = 3.141592653589793e+00
 numpy pi = 3.141592653589793
 scipy pi = 3.141592653589793

1.2 Basic operators +,-,\*,/,//,**

A few things you should know about python 3:

2 is an integer, 2.0 is a floating point number, (2.0+3.0j) is a complex.
The division operator returns a floating point number. To perform the division of two integers and expect an integer number as the result, use // or int()
The power operation $x^y$ is x**y or pow(x,y).
a += 1 is a shorthand for a = a + 1, same for subtraction, multiplication and division.



In [2]:

    
print("1+2 = ", 1+2)
print("1.0+2 = ", 1.0+2)
print("1.0+2.0 = ", 1.0+2.0)
print("4/2 = ", 4/2)
print("4//2 = ", 4//2)
print("8//3 = ", 8//3)
print("8/3 = ", 8/3)
print("Integer of 8/3 = ", int(8/3))
print("(2.5,3.0) =", 2.5+3.0j)
print("2^4 = ", 2**4)
print("2.0^4 = ", 2.**4)
print("2.0^4 = ",pow(2.0,4))
a = 10.
a *= 9.
print("a = 10. then a *= 9.=",a)









    



1+2 =  3
1.0+2 =  3.0
1.0+2.0 =  3.0
4/2 =  2.0
4//2 =  2
8//3 =  2
8/3 =  2.6666666666666665
Integer of 8/3 =  2
(2.5,3.0) = (2.5+3j)
2^4 =  16
2.0^4 =  16.0
2.0^4 =  16.0
a = 10. then a *= 9.= 90.0

1.3 Mathematical functions

You may find all mathematical functions here </FONT>: https://docs.python.org/3.5/library/math.html</FONT>.



In [3]:

    
print("e^1=",math.exp(1.))
print("e^(-infinity)=",math.exp(-math.inf))
print("log(e)=", math.log(math.e))
print("log_10(1000.)=",math.log10(1000.))
a = math.cos(math.pi/4.)
print("a = cos(pi/4) = %1.8f" %a)
print("arcos(a)*4./pi = ",math.acos(a)*4./math.pi)









    



e^1= 2.718281828459045
e^(-infinity)= 0.0
log(e)= 1.0
log_10(1000.)= 3.0
a = cos(pi/4) = 0.70710678
arcos(a)*4./pi =  1.0

1.4 Booleans

Logical operations result in True or False outcomes. They mostly drive if and while statements.



In [4]:

    
print("1 equal 2:",1==2)
print("1 non-equal 2:",1!=2)
print("1 > 2:",1>2)
print("1 < 2:",1<2)
print("10.==sqrt(100.):",10.==math.sqrt(100.))









    



1 equal 2: False
1 non-equal 2: True
1 > 2: False
1 < 2: True
10.==sqrt(100.): True



In [5]:

    
a = 0
if a > 0:
    a -= 1
elif a == 0:
    a = 1
print(a)

Let's take a pause here to discuss if, while, for statements. The indentation at the line following if a>0: is critical. It means that all the lines with such indentation following if a>0: are part of the actions needed if a>0 = True. Try the following cells using a = 1 and if a=4. Note also the role of the boolean operators and and or. Note also that each if, while, for statements need to end with : to signal that the following indented lines are part of the statement.



In [6]:

    
a = 4
b = 0
if a > 0 and a < 4:
    a += 1
    b = 3
else:
    a = 0
    b = 1
print(a,b)

0 1



In [7]:

    
a = 4
b = 0
if a > 0 or a < 4:
    a += 1
    b = 3
else:
    a = 0
    b = 1
print(a,b)

5 3

1.5 Loops

A key element of scientific computing is the ability to loop or iterate. You have two options: for and while:

for i in range(lower, upper, step): "range(lower, upper, step)" returns an iterable list of numbers from the lower number to the upper number (-1, if step>0, +1 if step<0), while incrementing by step. If step is not indicated, the default value is 1.

while condition: This loop will run while condition is true.

</ul> Here are some examples:



In [8]:

    
for i in range(0,5,1):
    print(i)



In [9]:

    
for i in range(5):
    print(i)



In [10]:

    
for i in range(3,5):
    print(i)

3
4



In [11]:

    
for i in range(5,2,-1):
    print(i)



In [12]:

    
a = 0
while a <=5 :
    print(a)
    a += 1
print("final a=",a)









    



0
1
2
3
4
5
final a= 6

2 Matrices

Matrices are central to CFD as we will see in future problems. The library or module numpy contains all the definitions and operators necessary for matrix operations.

2.1 Initialization of a matrix

You may set a matrix/vector to zero, unity or the values you may choose using numpy functions zeros, ones, and array. Please note:

Indices always start at 0. More on that later.
For one dimension array (vector) the syntax to initialize an array is np.zeros(n), however for multiple-dimension arrays, the argument to be passed to zeros or ones is a tuple (p,q,r) like this: np.zeros((p,q,r))
dtype allows for the specification of the variable format. There are many ways to pass the format: 'int', 'float32', 'float64', 'float', 'complex64', 'complex128', 'complex'. By default 'float' and 'complex' are 64 bits and 128 bits, respectively, or the same as 'float64' and 'complex128'. Should you be short on memory you may revert to 'float32' and 'complex64' with a loss of accuracy. In the second cell below, hange float64 to 32.



In [13]:

    
import numpy as np

n = 10
a = np.zeros(n)
print("a = ",a)
m = 3
b = np.ones((m,n))
print("b = ",b)
c = np.array([1, 2, 3, 4])
print("c = ",c)
d = np.array([[1, 2], [3, 4]], dtype = np.complex)
print("d = ",d)









    



a =  [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
b =  [[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]
c =  [1 2 3 4]
d =  [[ 1.+0.j  2.+0.j]
 [ 3.+0.j  4.+0.j]]



In [14]:

    
a = np.ones((2,2), dtype = 'float64')
for j in range(2):
    for i in range(2):
        a[i,j] *= np.float64(math.pi/(i+j+1))
print(a)
print("a[0,0] =",a[0,0])









    



[[ 3.14159265  1.57079633]
 [ 1.57079633  1.04719755]]
a[0,0] = 3.14159265359

2.3 The problem with arrays: =

Any programming language has a quirk or two, that may baffle/annoy you. May be it is only me, but I find the equal sign in python dangerous. Here is why. First we define an array $$ a_i = 2^i,\;i=0,\ldots,9 $$



In [15]:

    
n = 10
a = np.zeros(n, dtype = np.float)
for i in range(n):
    a[i] = 2.**i
print("a = ",a)









    



a =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]

Now, I want to modify the last element a[n-1] but in a different array to keep my initial array intact (something quite common in CFD).



In [16]:

    
b = a
b[n-1] -= 500.
print("a = ",a)
print("b = ",b)









    



a =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.   12.]
b =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.   12.]

Both my arrays have been changed!!! The same is not true for simple variables:



In [17]:

    
a = 1.
b = a
b += 1.
print("a=",a)
print("b=",b)









    



a= 1.0
b= 2.0

The fix for arrays is to use np.copy:



In [18]:

    
n = 10
a = np.zeros(n)
for i in range(n):
    a[i] = 2.**i
b = np.copy(a)
b[n-1] -= 500.
print("a = ",a)
print("b = ",b)









    



a =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]
b =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.   12.]

Note that when you perform an operation on all the elements of your array, the problem of linking your old array and new array disappears:



In [19]:

    
n = 10
a = np.zeros(n)
for i in range(n):
    a[i] = 2.**i
b = a * 2.
print("a = ",a)
print("b = ",b)









    



a =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]
b =  [    2.     4.     8.    16.    32.    64.   128.   256.   512.  1024.]

which is equivalent to:



In [20]:

    
n = 10
a = np.zeros(n)
for i in range(n):
    a[i] = 2.**i
b = np.copy(a)
b *= 2.
print("a = ",a)
print("b = ",b)









    



a =  [   1.    2.    4.    8.   16.   32.   64.  128.  256.  512.]
b =  [    2.     4.     8.    16.    32.    64.   128.   256.   512.  1024.]

but this is bad:



In [21]:

    
n = 10
a = np.zeros(n)
for i in range(n):
    a[i] = 2.**i
b = np.zeros(n)
b = a
b *= 2.
print("a = ",a)
print("b = ",b)









    



a =  [    2.     4.     8.    16.    32.    64.   128.   256.   512.  1024.]
b =  [    2.     4.     8.    16.    32.    64.   128.   256.   512.  1024.]

2.4 Array operations

To perform operations an all elements, numpy has many functions. In the following, we

Add 1 to all elements
Multiply two matrices, element by elements
Perform the dot product between two vectors
Perform a mathematical operation



In [22]:

    
n = 4
a = np.zeros(n, dtype = np.float64)
for i in range(n):
    a[i] = 2.**i
print("a=",a)    
b = np.copy(a)
b += 1.
print("b=a+1=",b)
b = np.zeros(n, dtype = np.float64)
for i in range(n):
    b[i] = np.float64(i)
print("a=",a)
print("b=",b)
print("a*b=",a*b)
print("a.b=",np.dot(a,b))
print("exp(-a)=",np.exp(-a))









    



a= [ 1.  2.  4.  8.]
b=a+1= [ 2.  3.  5.  9.]
a= [ 1.  2.  4.  8.]
b= [ 0.  1.  2.  3.]
a*b= [  0.   2.   8.  24.]
a.b= 34.0
exp(-a)= [  3.67879441e-01   1.35335283e-01   1.83156389e-02   3.35462628e-04]

Note that you need to use the exp of numpy and not from math, otherwise you get this:



In [23]:

    
print("exp(-a)=",math.exp(-a))









    



---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-43959dcef9a1> in <module>()
----> 1 print("exp(-a)=",math.exp(-a))

TypeError: only length-1 arrays can be converted to Python scalars

Getting errors is part of code development. Python will provide sometimes short, sometimes way too long error messages. Before you throw your hands in the air, give up and immediately email the TA or the instructor, browse the error message for things that could make sense. Typically, the last line contains useful information. Here it says that math.exp is expecting a scalar but we gave a length-1 array.

2.5 Indexing and Slicing

A very attractive feature of python is the nomenclature for indexing and slicing that enables to maximize the efficiency of vector operations. Here is a summary of the indexing and slicing nomenclature for a vector:



In [32]:

    
%matplotlib inline 
# plots graphs within the notebook
%config InlineBackend.figure_format='svg' # not sure what this does, may be default images to svg format

from IPython.display import Image
from IPython.display import clear_output

from IPython.core.display import HTML
def header(text):
    raw_html = '<h4>' + str(text) + '</h4>'
    return raw_html

def box(text):
    raw_html = '<div style="border:1px dotted black;padding:2em;">'+str(text)+'</div>'
    return HTML(raw_html)

def nobox(text):
    raw_html = '<p>'+str(text)+'</p>'
    return HTML(raw_html)

def addContent(raw_html):
    global htmlContent
    htmlContent += raw_html
    
class PDF(object):
  def __init__(self, pdf, size=(200,200)):
    self.pdf = pdf
    self.size = size

  def _repr_html_(self):
    return '<iframe src={0} width={1[0]} height={1[1]}></iframe>'.format(self.pdf, self.size)

  def _repr_latex_(self):
    return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)


PDF('figures/indexing.pdf',size=(800,600))









    Out[32]:



In [46]:

    
n = 6
a = np.zeros(6,dtype=np.float)

for i in range(n):
    a[i] = 2**i
print("a[0:n] = ",a[0:n])
print("a[:] = ",a[:])
print("a[1:]  = ",a[1:])
print("a[1:n] = ",a[1:n])









    



a[0:n] =  [  1.   2.   4.   8.  16.  32.]
a[:] =  [  1.   2.   4.   8.  16.  32.]
a[1:]  =  [  2.   4.   8.  16.  32.]
a[1:n] =  [  2.   4.   8.  16.  32.]



In [47]:

    
print("a[0:n]   = ",a[0:n])
print("a[0:n-1] = ",a[0:n-1])
print("a[0:-1]  = ",a[0:-1])
print("a[:-1]   = ",a[:-1])









    



a[0:n]   =  [  1.   2.   4.   8.  16.  32.]
a[0:n-1] =  [  1.   2.   4.   8.  16.]
a[0:-1]  =  [  1.   2.   4.   8.  16.]
a[:-1]   =  [  1.   2.   4.   8.  16.]



In [48]:

    
print("a[2:] = ",a[2:])
print("a[:-3]= ",a[:-3])









    



a[2:] =  [  4.   8.  16.  32.]
a[:-3]=  [ 1.  2.  4.]



In [52]:

    
b = np.zeros(n,dtype=np.float)
b[:-1] = a[1:] - a[:-1]
print("a[1:]  =",a[1:])
print("-")
print("a[:-1] =",a[:-1])
print("       =",b[:-1])









    



a[1:]  = [  2.   4.   8.  16.  32.]
-
a[:-1] = [  1.   2.   4.   8.  16.]
       = [  1.   2.   4.   8.  16.]



In [ ]:

    
n = 100000
mu = 0.
sigma = 1.
a = np.random.normal(loc = mu, scale = sigma, size = n)

print("a[:-(n-10)] = ", a[:-(n-10)])
print("a[0:10] = ", a[0:10])



In [ ]:

    
import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(a, 50, normed=True)
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * \
         np.exp( - (bins - mu)**2 / (2 * sigma**2) ),\
         linewidth=2, color='r')
plt.show()



In [ ]:

    
n = 100
mu = 0.
sigma = 0.2
x = np.linspace(0,1.,n)
a = np.sin(2.*math.pi*x)
plt.plot(x,a, lw = 2)
plt.show()
a_noisy = np.copy(a)
n_noise = 4
k_noise = 50
noise_amplitude = np.random.normal(loc = mu, scale = sigma, size = n_noise)
noise_phase = 2*math.pi*np.random.normal(loc = mu, scale = 1.0, size = n_noise)
for i in range(n_noise):
    a_noisy += noise_amplitude[i]*np.sin(2*(k_noise+i)*math.pi*x + noise_phase[i])
plt.plot(x,a_noisy, lw = 2)
plt.show() 
print(noise_amplitude)
print(noise_phase)



In [ ]:

    
a_filter = np.zeros(n)
a_filter[1:-1] = (a_noisy[0:-2] + 2. * a_noisy[1:-1] + a_noisy[2:])/4.
a_filter[0] = (a_noisy[n-2] + 2. * a_noisy[0] + a_noisy[1])/4.
a_filter[n-1] = a_filter[0]
plt.plot(x,a_noisy, linestyle="solid", lw = 1, label = r"noisy signal")
plt.plot(x,a_filter, color = "red",linestyle="solid", lw = 2, label = r"filtered signal")
plt.plot(x,a, color = "green", linestyle="None", marker = "o", markeredgecolor = "green", markeredgewidth = 1,\
         markerfacecolor = "None", label = r"original signal")
plt.legend(loc="lower left", bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.xlabel(r"$x$", fontdict = fontlabel)
plt.ylabel(r"$a(x)$", fontdict = fontlabel)
plt.show()



In [ ]:

    
plt.semilogy(x,np.abs(a_noisy - a_filter), lw = 2, label=r"\vert a_{noisy}-a_{filter}\vert")
plt.semilogy(x,np.abs(a - a_filter), color = "red", lw = 2, label=r"\vert a-a_{filter}\vert")
plt.legend(loc="lower left", bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.xlabel(r"$x$", fontdict = fontlabel)
plt.ylabel(r"difference", fontdict = fontlabel)



In [ ]:

    
def explicit_filter(a):
    n = len(a)
    a_filter = np.zeros(n)
    a_filter[1:-1] = (a[0:-2] + 2.*a[1:-1] + a[2:])/4.
    a_filter[0] = (a[n-2] + 2.*a[0] + a[1])/4.
    a_filter[n-1] = a_filter[0]
    return a_filter



In [ ]:

    
a_filter = explicit_filter(a_noisy)
plt.semilogy(x,np.abs(a_noisy - a_filter), lw = 2, label=r"$\vert a_{noisy}-a_{filter}\vert$")
plt.semilogy(x,np.abs(a - a_filter), color = "red", lw = 2, label=r"$\vert a-a_{filter}\vert$")
plt.legend(loc="lower left", bbox_to_anchor=[0, 1],
           ncol=2, shadow=True, fancybox=True)
plt.xlabel(r"$x$", fontdict = fontlabel)
plt.ylabel(r"difference", fontdict = fontlabel)



In [ ]:



In [ ]:

    
%matplotlib inline 
# plots graphs within the notebook
%config InlineBackend.figure_format='svg' # not sure what this does, may be default images to svg format

from IPython.display import display,Image, Latex
from __future__ import division
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')

import time

from IPython.display import display,Image, Latex

from IPython.display import clear_output

#import SchemDraw as schem
#import SchemDraw.elements as e

import matplotlib.pyplot as plt
import numpy as np
import math
import scipy.constants as sc

import sympy as sym

from IPython.core.display import HTML
def header(text):
    raw_html = '<h4>' + str(text) + '</h4>'
    return raw_html

def box(text):
    raw_html = '<div style="border:1px dotted black;padding:2em;">'+str(text)+'</div>'
    return HTML(raw_html)

def nobox(text):
    raw_html = '<p>'+str(text)+'</p>'
    return HTML(raw_html)

def addContent(raw_html):
    global htmlContent
    htmlContent += raw_html
    
class PDF(object):
  def __init__(self, pdf, size=(200,200)):
    self.pdf = pdf
    self.size = size

  def _repr_html_(self):
    return '<iframe src={0} width={1[0]} height={1[1]}></iframe>'.format(self.pdf, self.size)

  def _repr_latex_(self):
    return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)

class ListTable(list):
    """ Overridden list class which takes a 2-dimensional list of 
        the form [[1,2,3],[4,5,6]], and renders an HTML Table in 
        IPython Notebook. """
    
    def _repr_html_(self):
        html = ["<table>"]
        for row in self:
            html.append("<tr>")
            
            for col in row:
                html.append("<td>{0}</td>".format(col))
            
            html.append("</tr>")
        html.append("</table>")
        return ''.join(html)
    
font = {'family' : 'serif',
        #'color'  : 'black',
        'weight' : 'normal',
        'size'   : 12,
        }
fontlabel = {'family' : 'serif',
        #'color'  : 'black',
        'weight' : 'normal',
        'size'   : 16,
        }

from matplotlib.ticker import FormatStrFormatter
plt.rc('font', **font)