The goal is to have working python with the following packages and bindings:
numpy (commonly used data types and algorithms in numerical computation)scipy (commonly used algorithms in scientific computation)matplotlibipython notebook or jupyter notebookIn the following we will explain how to get a native installation on your system. Beside that it is also possible to use the prepared VirtualBox image, follow the Windows/macOS section for this.
Install the packages that are provided by your linux distribution, for Ubuntu up to 16.10:
sudo apt install python-numpy python-scipy python-matplotlib python-pip cython ipython-notebook build-essential libblas-dev liblapack-devwhile for Ubuntu 17.04:
sudo apt install python-numpy python-scipy python-matplotlib python-pip cython jupyter-notebook build-essential libblas-dev liblapack-devThis gives you numpy, scipy, matplotlib and the ipython notebook respective jupyter notebook.
We recommend using Anaconda which bundles all required packages. Installation instructions are provided on the Anaconda website.
We suggest to start an ipython/jupyter notebook and do some interactive experiments. Execute on the command line (Ubuntu <= 16.10, macOS):
ipython notebookrespective (Ubuntu 17.04)
jupyter notebookand click on New Notebook.
It is good practice to first include all python packages that are going to be used.
We will likely always use numpy that provides some numerical datatypes and standard numerical algorithms as well as the print function from Python 3 replacing the print statement and the / function from Python 3 that automatically type converts integer to floats to avoid counterintuitive behaviour.
After you have entered commands into a cell, press Shift + Enter to execute it.
In [1]:
from __future__ import print_function, division
import numpy as np # numpy will be used a lot, thus it is convenient to address it with np instead of numpy
In [2]:
a = 1 # new integer variable a that is assigned the value 1
b = 2 # new integer variable b that is assigned the value 2
c = 2.0 # new floating point variable c that is assigned the value 2.0
d = True # new boolean variable d that is assignes the value True
In [3]:
print(b*c) # automatic type conversion to float
In [4]:
print(a/b) # what you expect with automatic type conversion, beware: this is different in standard python 2!
In [5]:
print(a//b) # integer division
Note that all these datatypes are immutable, this will be explained in the functions paragraph.
Furthermore there are the composed datatypes list and dict (short for dictionary).
A list is a ordered collection of arbitrary python objects, a dictionary a mapping between arbitrary python objects (the object on the source side have to be hashable).
In contrast to the previously encountered datatypes, list and dict are mutable, explained in the function section.
In [6]:
a = [1, 2, 14.0, 3] # declaration of a list, note that members of the collection may have different types
print(a[0], a[2], a[-1]) # indexing is zero based: this gives the first, third and last element of the list
print(len(a)) # len gives the length of a list
In [7]:
b = list(range(10)) # range is a function that constructs a generator giving a list with integers
c = list(range(3,7))
d = [x*x for x in b if x % 2 == 0] # list comprehension is a method to construct new lists out of a given one,
# here all even squares < 100
print(b, c, d)
In [8]:
cd = c + d # concatenate lists c and d
print(cd)
In [9]:
print(c)
c.append(5) # extend the list c with another element 5
print(c)
In [10]:
e = 'hello world' # declares a string which is a list of characters
print(e, e[3])
In [11]:
f = { # declare a dictionary with the keys 'algorithm', 'tolerance', 'suspicious iterates'
'algorithm': 'downhill',
'tolerance': 1e-6,
'suspicious iterates': [3,7,11]
}
print(f['algorithm'], f.keys()) # access mapped elements of dictionary by their key, the function keys gives a list of keys
In [12]:
for elem in a:
print(elem, type(elem))
Note that the second line that constitutes the inner block of the for loop has been indented by 4 characters.
This is one syntax rule in python, which you will also encounter in if-else and while statements:
Blocks are announced by a : and have to be indented. End of indentation indicates end of block. In Java/C++ you would use { to indicate the begin of a block and } to indicate its end.
In [13]:
# less elegant, but same functionality as in previous example
for ii in range(len(a)):
print(a[ii], type(a[ii]))
In [14]:
for key in f:
print(key, ':', f[key])
In [15]:
# sum all cubes of integers from 0,...,99
cubesum = 0
for ii in range(100):
cubesum += ii*ii*ii
print(cubesum)
In [16]:
print(3 == 7)
print(7 < 7)
print(7 <= 7)
print(3 in range(0,9))
print(4 in ['apple', 'pear'])
Booleans can be composed with and and or and negated with not:
In [17]:
print(cubesum % 100 == 0 and cubesum < 0)
print(cubesum % 100 == 0 or cubesum < 0)
print(not cubesum in ['apple', 'pear'])
In [18]:
a = 2.0
x = 2.0
# Newton iteration to compute x = sqrt(a) with a residual <= 1e-8
while( abs(x*x-a) > 1e-8 ):
x = .5*(x+a/x)
print(x, x*x-a)
Customary with loops are the break and continue statements that exit the loop execution if a certain condition is satisfied or jump to the next loop execution:
In [19]:
a = 2.0
x = 2.0
count = 0
# Newton iteration to compute x = sqrt(a) with a residual that cannot be achieved in double precision
while True:
if abs(x*x-a) <= 1e-30: # exit loop because precision is reached
print("Very precise square root computed")
break
if count > 100: # exit loop because iteration limit is reached
print("Iteration limit exceeded")
break
x = .5*(x+a/x)
count += 1
print(count, x, x*x-a)
In [20]:
# sum all primes from 10,...,99 (this is inefficient!)
primesum = 0
for ii in range(10,100):
# test if ii is prime
prime = True
for jj in range(2,ii):
if ii % jj == 0:
prime = False
# continue with next iteration if ii is not a prime number
if not prime:
continue
primesum += ii
print(primesum)
In [21]:
def newton_squareroot(a, x0): # declare a function newton_squareroot that takes two arguments a and x0
"""
Compute square root of a by Newton iteration, with initial guess x0
"""
x = 1.0
it = 0
itmax = 100
tol = 1e-12
while (abs(x*x-a) > tol and it < itmax):
x = .5*(x+a/x)
it += 1
return x # return the computed square root x
In [22]:
print(newton_squareroot(2.0, 1.0), newton_squareroot(9.0, 1.0))
In [23]:
entropy = lambda x: -x*np.log(x) # declare a function entropy with one argument by a lambda expression
print(entropy(0.2))
In [24]:
# declare a function density with three arguments by a lambda expression
density = lambda x, m, s: np.exp(-(x-m)*(x-m)/(2.0*s*s))/(np.sqrt(2.0*np.pi)*s)
print(density(0.3,0.1,1.2))
In functions declared with def it is also possible to have optional arguments.
They have to appear after nonoptional arguments:
In [25]:
def newton_squareroot(a, x0 = 1.0, tol = 1e-12, itmax = 100):
"""
Compute square root of a by Newton iteration, with initial guess x0
"""
x = x0
it = 0
while (abs(x*x-a) > tol and it < itmax):
x = .5*(x+a/x)
it += 1
return x, it # return the computed square root x and the number of iterations
In [26]:
print(newton_squareroot(3.4))
Now the difference between mutable and immutable can be explained: arguments passed to a function that are mutable can be changed during function execution, while immutable arguments can not:
In [27]:
def double(x):
x = x + x
def appendzero(x):
x.append(0)
In [28]:
# for an immutable argument
x = 1.0
double(x) # inside the function x = x+x, so x = 2.0, but as x is immutable the value does not change outside
print(x)
In [29]:
# for a mutable argument
x = [1]
appendzero(x) # inside the function x will be appended with 0, so x = [1, 0]
print(x)
In [30]:
a = np.array([3.0, 4.2, -1.0]) # converting an python array of floats to a numpy vector
b = np.zeros([4]) # a numpy vector containing 4 zeros
c = 3.0*np.ones([13]) # a numpy vector containing 13 times the entry 3.0
print(a, b, c)
In [31]:
d = np.zeros_like(a) # a numpy vector containing as many zeros as a has entries
e = np.ones_like(b) # a numpy vector containing as many ones as b has entries
print(d, e)
In [32]:
f = np.linspace(3.0, 7.0, 101) # a numpy vector with 101 elements, spaced equidistantly over [3.0, 7.0]
print(f)
Matrices can be constructed in a similar way:
In [33]:
A = np.array([ # convert a python array of arrays of floats to a numpy matrix
[3.0, 0.5],
[3.7, 1.2],
[5.0, -3.0]
])
print(A)
In [34]:
B = np.diag(a) # create a diagonal matrix with diagonal given by the array a
print(B)
some operations with matrices and vectors
In [35]:
# get dimensions
print(a.shape)
print(A.shape)
In [36]:
# transposition
print(A.transpose())
print(A.transpose().shape)
In [37]:
# matrix vector product
print(B.dot(a))
In [38]:
# access first row of matrix A:
print(A[0,:])
In [39]:
# access second column of matrix A:
print(A[:,1])
The previous examples are a special case of so called slicing in numpy to obtain subvectors and submatrices, further examples:
In [40]:
print(f[:4]) # the first four elements of f, equivalent to f[0:4]
print(f[12:18]) # elements f[12], ..., f[17]
print(f[::20]) # every twentieth element of f
print(f[1::20]) # every twentieth element of f, but starting at the second position
print(f[[0,5,8]]) # the first, sixth, ninth element of f
print(f[85:-4]) # all elements of f starting with the 86st and omitting the four last elements
It is sometimes necessary to "vectorize" matrices, forming the vector with all matrix rows appended.
In [41]:
vecA = A.reshape(-1) # vectorize A
print(vecA)
reshapedA = vecA.reshape([3,2])
print(reshapedA) # change shape of A to 3 x 2 instead 3 x 2, order of elements is preserved
In [42]:
# compute trace of matrix
print(B.trace())
In [43]:
# compute sum and norms of vector
print(u'\u03a3\u2096a\u2096 =', a.sum())
print(u'\u2016a\u2016\u2082 =', np.linalg.norm(a))
print(u'\u2016a\u2016\u2081 =', np.linalg.norm(a, 1))
print(u'\u2016a\u2016\u221e =', np.linalg.norm(a, np.inf) )
In [44]:
# solve linear system (backslash in matlab)
np.linalg.solve(B, np.ones([3])) # find x with B*x = [1,1,1]
Out[44]:
In [45]:
# compute inverse of A (you should always prefer solve or getting a decomposition)
print(np.linalg.inv(B))
It is often necessary to construct vectors and matrices from vectors and matrix of smaller dimensions by stacking or forming block matrices. This can be done with the concatenate and bmat in numpy.
In [46]:
ab = np.concatenate([a, b]) # vector obtained by stacking a and b
print(ab)
K = np.bmat([[np.eye(3), A], [A.transpose(), np.zeros([2,2])]]) # block matrix, eye(n): identity in R^n
print(K)
Many of the standard functions are available in numpy and vectorized:
They can be applied to a numpy vector and return a vector with the function applied to the components.
This is very efficiently implemented:
In [47]:
print(np.sin(a))
print(np.array([np.sin(x) for x in a])) # same result, but slower
To verify that the first variant is faster, we use %timeit, that determines the average execution time of a python statement. A command starting with % is a IPython magic command that only works in ipython and ipython notebooks, but not in python scripts.
In [48]:
%timeit test = np.sin(np.linspace(0.0, 1.0, 1000))
In [49]:
%timeit test = np.array([np.sin(ii*1e-3) for ii in range(1000)])
To get the execution time of a code block in python, you can use the time.time() of the time module that gives you the current system time:
In [50]:
def fib(n):
if n in [0, 1]:
return 1
return fib(n-1) + fib(n-2)
import time
zerotime = time.time()
print('30th fibonacci number:', fib(30))
executiontime = time.time() - zerotime
print('took {:.2f} s to compute using recursion'.format(executiontime))
You can also create a custom vectorized function out of a function that is defined for scalar values:
In [51]:
newton_squareroot_vectorized = np.vectorize(newton_squareroot)
In [52]:
print(newton_squareroot_vectorized(np.linspace(1.0, 10.0, 5)))
In [53]:
%timeit test = newton_squareroot_vectorized(np.linspace(1.0, 10.0, 100))
In [54]:
%timeit test = [newton_squareroot_vectorized(x) for x in np.linspace(1.0, 10.0, 100)]
In [55]:
import matplotlib.pyplot as plt
%matplotlib notebook
In [56]:
xs = np.linspace(-5.0, 10.0, 1000)
ys = np.sin(xs)
xs2 = np.linspace(-2.0, 8.0, 25)
ys2 = np.cos(xs2)
plt.figure()
plt.title('A sample plot')
plt.plot(xs, ys)
plt.plot(xs2, ys2, 'o')
plt.show()
In [57]:
# Compute iterates in Newton iteration for x^2 = 2
iterates = [1.0]
a = 2.0
while abs(iterates[-1]*iterates[-1] - a) > 1e-15:
iterates.append(.5*(iterates[-1]+a/iterates[-1]))
residuals = [abs(it*it-a) for it in iterates]
plt.figure()
plt.title(r'Convergence Newton Iteration $ x^2 = 2 $')
plt.xlabel(r'iteration $k$')
plt.ylabel(r'residual $\vert x_k^2 - 2 \vert$')
plt.semilogy([0, len(iterates)-1], [1e-15, 1e-15], '--')
plt.semilogy(range(len(iterates)), residuals, '-o')
plt.show()
Initial value problems as considered in the lecture are of the form $$ \dot x(t) = f(t, x(t), p), \quad x(t_0) = x_0 $$ where $ p $ is a parameter.
Scipy provides with ODE a class to solve ordinary differential equations.
The following function provides a convenience interface to solve an initial value problem as specified above:
In [58]:
import scipy.integrate
class IVPResult:
pass
def solve_ivp(f, ts, x0, p=None, integrator='dopri5', store_trajectory=False):
"""
Solve initial value problem
d/dt x = f(t, x, p); x(t0) = x0.
Evaluate Solution at time points specified in ts, with ts[0] = t0.
"""
ivp = scipy.integrate.ode(f)
ivp.set_integrator(integrator)
if store_trajectory:
times = []
points = []
def solout(t, x):
if len(times) == 0 or t != times[-1]:
times.append(t)
points.append(np.copy(x[:x0.shape[0]]))
ivp.set_solout(solout)
ivp.set_initial_value(x0, ts[0])
ivp.set_f_params(p)
result = IVPResult()
result.ts = ts
result.xs = np.zeros([ts.shape[0], x0.shape[0]])
result.success = True
result.xs[0,:] = x0
for ii in range(1,ts.shape[0]):
ivp.integrate(ts[ii])
result.xs[ii,:] = ivp.y[:x0.shape[0]]
if not ivp.successful():
result.success = False
break
if store_trajectory:
result.trajectory_t = np.array(times)
result.trajectory_x = np.array(points)
return result
As an example, we consider a harmonic oscillator that is driven by a periodic external force with frequency $ \omega $. The system is described by the differential equation $$ \ddot x + k^2 x = F \sin(\omega t) $$ with $ k = \sqrt{D/m} $ where $ D $ is the spring constant and $ m $ mass, furthermore we set $ F = 1 $.
Reformulation as a system of first order is given by $$ \begin{pmatrix} \dot x_0 \\ \dot x_1 \end{pmatrix} = \begin{pmatrix} x_1 \\ -k^2 x_0 + \sin(\omega t) \end{pmatrix}. $$
We have to define the function $ f $ for the right hand side in python where we set $ p = \begin{pmatrix} k \\ \omega \end{pmatrix} $.
In [59]:
def harm_osc(t, x, p):
k = p[0]
omega = p[1]
return np.array([x[1], -k*k*x[0]+np.sin(omega*t)])
First we use $ k = 1, \omega = 2 $ and solve the initial value problem with $ x(0) = 1, \dot x(0) = 0 $ as initial value on the time horizon $ [0, 25] $ and store the trajectory as computed by the initial value solver:
In [60]:
result = solve_ivp(harm_osc, np.array([0., 25.]), np.array([1., 0.]), p=np.array([1., 2.]), store_trajectory=True)
print(result.xs)
fig = plt.figure()
plt.plot(result.trajectory_t, result.trajectory_x)
plt.show()
Now the case of resonance disaster: we set $ k = \omega = 1 $:
In [61]:
result = solve_ivp(harm_osc, np.array([0., 25.]), np.array([1., 0.]), p=np.array([1., 1.]), store_trajectory=True)
print(result.xs)
fig = plt.figure()
plt.plot(result.trajectory_t, result.trajectory_x)
plt.show()
SLSQP provided by scipyWe use the problem HS71 from the Hock-Schittkowski collection of optimization problems. This is given by $$ \begin{array}{ll} \min_{x \in \mathbb R^4} & x_0 x_3 (x_0+x_1+x_2) \\ \text{s.t.} & x_0 x_1 x_2 x_3 -25 \ge 0 \\ & \Vert x \Vert_2^2 - 40 = 0 \\ & 1 \le x_i \le 5 \quad (i=0,\ldots,3) \end{array}. $$
Now we define functions
ffcn that computes $ f(x) := x_0 x_3 (x_0+x_1+x_2) + x_2 $,dfcn that computes $ \nabla f(x) = \begin{pmatrix} x_3 (x_0+x_1+x_2) + x_0 x_3 \\ x_0 x_3 \\ x_0 x_3 + 1 \\ x_0 (x_0+x_1+x_2) + x_0 x_3 \end{pmatrix} $cfcn_ieq that computes $ c_I(x) := x_0 x_1 x_2 x_3 - 25 $cfcn_eq that computes $ c_E(x) := \Vert x \Vert_2^2 - 40 $jfcn_ieq that computes $ Jc_I(x) = \nabla c_I(x)^T = \begin{pmatrix} x_1 x_2 x_3 & x_0 x_2 x_3 & x_0 x_1 x_3 & x_0 x_1 x_2 \end{pmatrix}$jfcn_eq that computes $ Jc_E(x) = \nabla c_E(x)^T = \begin{pmatrix} 2 x_0 & 2 x_1 & 2 x_2 & 2 x_3 \end{pmatrix}$
In [62]:
ffcn = lambda x: x[0]*x[3]*np.sum(x[:3]) + x[2]
dfcn = lambda x: np.array([x[3]*np.sum(x[:3])+x[0]*x[3], x[0]*x[3], x[0]*x[3]+1, x[0]*np.sum(x[:3])])
cfcn_ieq = lambda x: np.array([np.prod(x)-25.0])
cfcn_eq = lambda x: np.array([np.inner(x,x)-40])
def jfcn_ieq(x):
return np.array([
[np.prod(x[1:]), np.prod(x[[0,2,3]]), np.prod(x[[0,1,3]]), np.prod(x[:3])]
])
jfcn_eq = lambda x: np.array([2*x])
For sanity reasons we test first if our implementation has no obvious bugs by comparing with finite differences at a random point.
In [63]:
x0 = np.random.randn(4) # x0 is an array with 4 entries chosen from a normal distribution
h = 1e-8 # perturbation in finite differences
# get unit vectors
unit = np.eye(4)
# compute finite difference approximation to gradient in x0
fx0 = ffcn(x0) # store f(x0) to compute it only once
dfx0 = np.array([ffcn(x0+h*unit[ii,:])-fx0 for ii in range(4)])/h # approximation to gradient
print(np.linalg.norm(dfx0 - dfcn(x0)), np.inf)
# compute finite difference approximation to inequality jacobian in x0
cix0 = cfcn_ieq(x0)
jcix0 = (np.array([cfcn_ieq(x0+h*unit[ii,:])-cix0 for ii in range(4)])/h).transpose()
print(np.linalg.norm(jcix0 - jfcn_ieq(x0)), np.inf)
# compute finite difference approximation to equality jacobian in x0
cex0 = cfcn_eq(x0)
jcex0 = (np.array([cfcn_eq(x0+h*unit[ii,:])-cex0 for ii in range(4)])/h).transpose()
print(np.linalg.norm(jcex0 - jfcn_eq(x0)), np.inf)
The resulting deviations from the finite difference approximations give some confidence that we have done a correct implementation of the gradient and jacobian. Let us now solve the optimization problem with initial guess $ x^0 = \begin{pmatrix} 1 \\ 5 \\ 5 \\ 1 \end{pmatrix} $. Using the specified options shows a one line summary per SQP iteration.
In [64]:
import scipy.optimize
n = 4 # number of variables
m = 2 # number of constraints
l = np.ones([n]) # lower bound variables
u = 5.0*np.ones([n]) # upper bound variables
x0 = np.array([1., 5., 5., 1.]) # initial guess
result = scipy.optimize.minimize(ffcn, x0, method='SLSQP', jac=dfcn,
bounds=[(l[i], u[i]) for i in range(l.shape[0])],
constraints=[
{'type': 'eq', 'fun': cfcn_eq, 'jac': jfcn_eq},
{'type': 'ineq', 'fun': cfcn_ieq, 'jac': jfcn_ieq}
],
options={'disp': True, 'iprint': 2})
result is now a dictionary containing the solution vector and some status information:
In [65]:
print(result)
First and foremost it is useful to consult the documentation of the packages that are used.
There is a plethora of manuals, tutorials, introductions on the web for python and its usage in scientific computing. For example the recent notebooks by Robert Johannson. Using google is usually a good idea.
If you are fluent in german and looking for another basic introduction, consider the tutorial of the introductory numerical analysis lecture by Andreas Potschka.
In [ ]: