Numerical Techniques
Numerical Techniques
Chemical Engineers
Solving chemistry and engineering problems using Python
Publisher
First published in November 2020 by F.C. Grozema
Contents
Contents iii
6 Curve Fitting 36
6.1 Linear regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
6.2 Polynomial regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
6.3 Fitting to any function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
6.4 Some remarks on fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.5 Exercise: fitting to a Gaussian . . . . . . . . . . . . . . . . . . . . . . . . . . 43
11 Numerical differentiation 86
11.1 First derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
11.2 Second derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
11.3 Partial derivatives in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
11.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
11.5 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
Second derivative . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
12 Non-linear equations 94
12.1 The bisection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
Recursion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
Recursive bisection method . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
12.2 The Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Two coupled equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
12.3 Reduced Newton Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
Multiple equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
12.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
The bisection method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
The Newton method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
Newton for two equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
12.5 Extra exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Gradient descent in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
Index 157
Basic Introduction to Programming in
Python
Introduction to Python 1
Computers are exceptionally good at performing repetitive task at a 1.1 Getting Started with Python . 2
very high speed, and with virtually zero probability of mistakes. The Installing Python . . . . . . . . 3
word computer already implies that the main task of computers is to Documentation . . . . . . . . . 3
1.2 Basic calculations in Python . 3
do computations or calculations, which is the kind of task that occurs
Variables . . . . . . . . . . . . . 4
very often in engineering disciplines but also in basic sciences, including
Comments . . . . . . . . . . . . 5
chemistry and physics. Using computers, a wide variety of mathematical
Printing and formatting text . 5
problems can be solved, generally not by deriving analytical solution 1.3 Some remarks on syntax . . . 7
to a certain mathematical problem but by coming up with a numerical Equality vs. assignments . . . 7
solution. While analytical solutions are sometimes preferred, in a majority Precedence of operators . . . . 8
of cases it is not even possible to get to such a solution without making 1.4 Importing modules . . . . . . . 8
very severe approximations. In order to use a computer to help us solve 1.5 Exercise: Height of ball . . . . 9
out problems in science and engineering, we need a way of instructing
the computer what to do; i.e. we need a programming language. In
this course we use a language called Python, which is an easily usable
high-level programming language that is quickly becoming a standard
in scientific computing. One of the big advantages is that it is open
source (free) and can be used on virtually all platforms (Windows, Mac,
Linux). The language is easy to learn and programs written in Python are
generally easy to read and often much shorter than comparable programs
written in languages such as C or Fortran. A major difference between
Python and these languages is that C or Fortran codes have to be compiled
into an executable format before they can be run, while in Python the
code is interpreted while it is running. The disadvantage of this is that
computationally intensive programs run a lot slower than pre-compiled
ones written in C or Fortran. However, Python comes packed with an
enormous range of standard modules that contain many precompiled
functions and algorithms for many common (numerical) tasks. Examples
include opening and reading files, complicated mathematical functions,
plotting data in graphs and advanced numerical methods. The version
of Python that is currently used, is version 3.7 and that is what we
will use in this course (however most of the examples will also run
in Python 2.7). In this chapter we will give a general overview of the
programming language Python, covering most of the basic features that
all programming languages share and that will be used in the chapters
that follow.
Installing Python
Documentation
The first way in which Python can be used is as a basic calculator, where
you can just print the result of a simple equation, for instance to calculate
the height of a ball from the surface of the earth with an initial speed
𝑣0 :
1 2
𝑦(𝑡) = 𝑣0 𝑡 − 𝑔𝑡 (1.1)
2
Here 𝑦(𝑡) is the height of the ball, 𝑔 the gravitational constant, and 𝑡 the
time elapsed after the ball has been launched. We can simply compute
the height of the ball at 𝑡 = 0.6 in a Python program as:
1 print(5*0.6 - 0.5*9.81*0.6**2)
The line of code in the grey area also illustrates how pieces of example
code in this text can be recognised. The words in red are the so-called
reserved words in Python. They indicate a specific action that you want
to happen in the code, in this case printing the answer to the calculation
on the screen. These reserved words can not be used for anything else in
your code, for instance as the name of a variable.
Variables
While the piece of code above does something useful, it is not easy to
recognise the different parameters, let alone vary them easily. For instance
selecting a different time, t, or a different initial velocity requires us to
know the equation and what each number stands for. It is much easier
to change the different parameters if the equation would be written in
terms of the variables. In programming languages, and Python is no
exception, it is possible to define 𝑣 0 , 𝑔 , 𝑡 and 𝑦 as variables and combine
them into to the right hand side of the equation for our ball. The result
can be written as:
1 v0 = 5
2 g = 9.81
3 t = 0.6
4 y = v0*t - 0.5*g*t**2
5 print(y)
Comments
So far we have only the discussed the actual program statements that
actually ’do’ something. For a readable code it is advisable to add
comments that do not actually do something but they describe the code.
This makes it a lot easier to follow what is going on in the code and will
make it easier to understand and modify a code at a later stage. In Python,
comments start with the hashtag character. All text after this character is
a comment and will not be considered when the code is executed. Our
code with comments looks like this:
1 # Program to calculate the height of a ball that moves in the vertical
direction
2 v0 = 5 # Initial velocity
3 g = 9.81 # acceleration of gravity
4 t = 0.6 # time
5 y = v0*t - 0.5*g*t**2 # calculation of the position of the ball
6 print(y)
This program again does exactly the same thing as before, but it is much
easier to understand for other people, they may actually figure out what
is going on in the code when they read it without explanation by the
programmer. Any code that consist of more than a few lines will greatly
benefit from the inclusion of comments and a smart choice of the names
of variables.
In the programs we have looked at so far the output was just a single
number, which basically does the job we would like but it is in many cases
nice to have a somewhat more informative line of output. Ideally, we
would also like to have some control over the formatting of the number
in the output, for instance the number of decimal spaces. An example of
such output is given below:
At t=0.6 s, the height of the ball is 1.23 m.
Such kind of output can be achieved using the print statement that we
have seen before but the way it is used is a bit more complicated. The
formatting that is used here is the so-called ’printf’ formatting that has
its origin in the programming language C. A line of code that supplies
such output is the following:
1 print(’At t=%g s, the height of the ball is %.2f m.’ % (t, y))
1 Introduction to Python 6
This new syntax uses the method format() of the strings. Each slot where
we want to insert the value of a variable is now given by a {what:how}
where ’what’ is what we want to print and ’how’ specifies the format.
In the example above the whats are 0 and 1 and refers to the first and
second argument of the format() method. Alternatively we could have
written:
1 " At t = {a:1.3f} sec., the height of the ball is {b:8f}".format(a=t, b=y)
Clearly, it is not the best idea to collect all statements on a single line;
this line is almost impossible to read! In this line of code, the spaces, for
instance around ’=’ are omitted, which is also true for other mathematical
operations. A general convention is to use one blank space around
=, - and + and no spaces around *, / and **. It should be noted that
this is just a matter of readability, omitting the spaces is in principle
correct. For guidelines about syntax conventions you can refer to the
PEP8 documents at [Link]
This series of conventions are widely adopted in the Python community
and significantly improve the readability of codes.
When we first wrote the line of code that contained the equation for
calculating the height of a ball at time 𝑡 the result looked very much
like a mathematical equation. In mathematics the ’=’ sign means that
the expressions on both sides of this sign are in fact equal. In Python
(and most other programming languages) the ’=’ indicates what we call
an assignment. An assignment basically means the following: calculate
everything on the right-hand side of the ’=’ and store the result in the
variable on the left-hand side. Although this may seem very similar to
the mathematical meaning it is in fact very different. Consider the line
below:
1 y = y + 3
This program will result in the output of three numbers: 3, 7 and 49.
1 Introduction to Python 8
Precedence of operators
Formulas in Python are evaluated in the same way as they would normally
be evaluated in mathematics. This means that power operations, 𝑎 4 coded
as a**4 have precedence over multiplications and division, which in turn
have precedence over addition and subtraction. We can use parentheses
in our code to change the way an expression is evaluated.
The basic types of operations that are present in Python are very limited,
but it is already possible to perform many interesting calculations. As
we have seen, the basic arithmetic operations are present, but if we
need to calculate a square root for instance this is not included in the
standard operations. The good thing about Python is that there is a wide
range of modules that can be included that contain useful functions. For
instance the module ’math’ contains a square root function (sqrt) and
many other mathematical functions, including sin, cos, exp and log. In
order to use these functions it is necessary to ’import’ these functions or
a whole module into you Python program, which is done by the import
statement:
1 import math
2 v0 = 9.0
3 answer = [Link](v0)
4 print(answer)
As can be seen in this code, if the module is imported in this way, the
function name requires a prefix which is the name of the module where
it comes from. This prefix can be a bit annoying, especially if the function
is used very often. An alternative import syntax makes it possible to skip
the prefix, and only specific functions can be imported:
1 from math import sqrt
2 v0 = 9.0
3 answer = sqrt(v0)
4 print(answer)
After this, the function ’sqrt’ can be used without the prefix. It is also
possible to import more than one function at once:
1 from math import sqrt, exp, sin
In the last statement all functions in the math module are imported and
can be used in your code without the prefix. In general it is advised
to import only the functions are actually used in a program because
importing all results in a lot of names in the program that are not used,
but can also not be used for variables anymore. However it is sometimes
quite convenient, especially in the case of the math module, to just import
all functions.
Sometimes it can be convenient to import modules and functions and
give them new names in one go:
1 Introduction to Python 9
1 import math as m
2 #now the math module is called m in the program
3 v = [Link]([Link]) #pi is a constant in the math module that is now called
’m’
Again, all the red words are reserved words in Python. You can also
change the name of a given function that you are importing from a
module. This is done in a very similar way than for the changing the
name of the module itself
1 from math import sin as theSinus, cos as theCosinus, log as theLog
2 x = 5.0
3 v = theLog(x)
4 v = theSinus(x)*theCosinus(x) + theLog(x)
However this is generally not a very good idea as you can loose track of
where a given function is defined. Besides writing [Link](x) is generally
much faster much more readable as theSinus(x).
The math module is only one of the many modules that can be used in
Python, the virtually endless list contains modules related to designing
web interfaces, doing operations on all sorts of graphical formats (pic-
tures), interacting with the operating system of a computer, etc. In this
course we are dealing mostly with the use of Python in a scientific/engi-
neering environment and we focus on numerical mathematics techniques.
For the purpose of this course there are two additional modules that
we will encounter in the chapters that follow. One is a collection of a
wide variety of pre-programmed numerical methods to solve all sorts of
problems, including integration, matrix operations, solving differential
equations, etc. The library that contains these functions is called numpy.
The second module that is useful is called matplotlib, which contains
functions to make a graphical representation of data on screen or in a
picture. These modules will be introduced in later chapters but they are
imported in the same way as the math module here.
For the first exercise you will reproduce the example presented in
section 1.2. You will therefore write a small program named ’[Link]’
that automatically computes the height of a ball that falls vertically
following:
1 2
𝑦 = 𝑣0 𝑡 − 𝑔𝑡 (1.2)
2
and print the answer. You’ll define the variable as shown in page 8. We
now ask the question : How long does it take the ball to reach a certain
height 𝑦 𝑐 ? Posing 𝑦 = 𝑦 𝑐 in the equation above and rearranging the
terms we find
1 2
𝑔𝑡 − 𝑣0 𝑡 + 𝑦 𝑐 = 0 (1.3)
2
1 Introduction to Python 10
You can solve this equation using the well known formula for a quadratic
equation and find the expression of the two roots
q
𝑡± = 𝑣 0 ± 𝑣02 − 2 𝑔 𝑦 𝑐 /𝑔 (1.4)
There are two times because the ball can reach 𝑦 𝑐 on its way up or down.
Update your code so that it also computes the values of 𝑡 𝑝𝑚 for 𝑦 𝑐 = 0.2.
What happens if you take 𝑦 𝑐 = 2.0 ? Change your code to intercept the
error message.
Note: the square root function can be loaded from the 𝑚𝑎𝑡 ℎ module.
1 import math as m
2 [Link](4.)
Loops and lists 2
As mentioned in Chapter 1, computers are exceptionally well-suited 2.1 While-loops . . . . . . . . . . . 11
to perform (boring) repetitive tasks in a automatic way. In computer 2.2 Boolean expressions . . . . . 13
programs such tasks can conveniently be performed in loops of which 2.3 Lists . . . . . . . . . . . . . . . 13
2.4 For-loops . . . . . . . . . . . . 14
different types exist. Additionally, what is often useful for large collections
2.5 The range construction . . . 15
of data (a row of numbers) to store them in a smart way, for instance
2.6 Branching: if-else statements 16
in a vectors or an array. In Python, such storage can be done in a ’list’.
2.7 Exercise: Money in the bank . 17
Loops and lists, together with functions and if-statements form the basic 2.8 Exercise: Cooking a perfect egg18
core of programming operations that are used in this course and a
thorough understanding of these concepts is essential for understanding
the following chapters.
2.1 While-loops
This does the job perfectly, twelve almost identical lines of code (with
three statements on each line) that give the required result. Now imagine
wanting a table in steps 0.01 degrees? It is clear who is doing all the hard
work here, in terms of typing the lines of code. While it works perfectly
well, it isn’t using the capabilities of computers to the fullest. The main
issue is: all lines do the same thing, only the value in the variable C
changes. All programming languages are equipped with a construct to
deal with such problems in a very easy way, so-called loops. In Python
there are two kinds of loops, while-loops and for-loops.
A while loop is used to repeat a set of statements as long as a certain
condition is true. We will use this type of loop here as an example to
generate the list temperatures for our conversion table. A sequence of
operations, or an algorithm for this is listed below in text form.
2 Loops and lists 12
This small code contains one of the most important features of program-
ming in Python: the block of statements that is to be executed in the while
loop has to be indented by the same amount. All operations under point
three in the algorithm above start in the code on the same position of the
line. In this way, the Python interpreter knows that these three statement
should be ’inside the loop’. At the top of the loop, line 4, Python checks
whether the condition is true (C <= 40), which is certainly true at the
beginning since the initial value of C is -20. Subsequently it executes
the three lines inside the loop (including printing the table entry) and
goes back to the top of the loop, line 4. Inside the loop the value of C has
changed so the check has to be performed again. This cycle of events is
repeated until, at some point, the condition is not true anymore: C has
the value of 45 after twelve cycles. At this point, the program does not go
into the loop anymore but skips to the first statement after the loop, line
8. The colon : behind the while statement on line 4 is essential because it
marks the start of the block of code on the next lines to be executed. It is
interesting to experiment with this code, for instance giving the last line
the same indentation as those inside the loop. In this case the lines in the
table will be separated by a line of dashes.
It is easy to see from this code that this way of solving this problem is
much more convenient and powerful than typing a separate line of code
for each temperature. On top of that, it is very easy now to make a much
longer list, incremented every 1 degree for instance by just changing the
values of dC. Note the use in the piece of code above of the printing
syntax:
1 print(’-’*20)
An important aspect of while loops is that at the top of the loop there
is an expression of which it is evaluated whether it is true. This type
of expression is called a boolean expression that can either be true or
false. Other comparisons that can be useful are listed in the following
fragment:
1 C == 40 # C equals 40
2 C != 40 # C does not equal 40
3 C >= 40 # C is larger than or equal to 40
4 C > 40 # C is larger than 40
5 C < 40 # C is smaller than 40
Boolean expressions can also be combined with the keywords ’and’ and
’or’:
1 while x > 0 and y <= 1:
2 print(x, y)
If both separate conditions are true the overall result is ’True’, and in all
other cases the result will be ’False’.
2.3 Lists
Until now all the calculations we have made using variables contained
a single number in an isolated variable. In many cases, especially in
mathematics, it is more natural if these number are grouped together
in a list, or an array as it is called in many text about programming. An
illustration is the list of temperatures in the examples above. In Python a
list can be made just assigning a variable name to a list of numbers that
is contained in square brackets, separated by commas:
1 C = [-20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40]
The variable C now refers to a list object that contains 13 elements and
all these elements are integer numbers. Each element of the list can be
assessed separately using an index that refers to its position in the list.
The index runs from 0 to 12 in this case. To indicate the third element in
the list we can write C[2], which refers to an object of the type int, with
the value -10.
Some basic operations can be performed on lists for instance appending
elements at the end:
1 C = [ -10, -5, 0, 5, 10, 15, 20, 25, 30]
2 [Link](35) #add a new element at the end
1 C = [ -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40, 45]
2 [Link](0, 15) # insert 15 a new element on index 0
There are many more operations that can be performed on lists that may
be handy in certain situations but the ones we will use in this course are
quite limited. The combination of a while (or for) loop and lists can be
very powerful since the elements can be addressed one by one in a loop
in a very short piece of code. For instance in the code below, a list is filled
with temperatures starting from -50 up to 200 degrees in steps of 2.5
degrees:
1 C =[]
2 C_value = -50
3 C_max = 200
4 while C_value <= C_max:
5 [Link](C_value)
6 C_value += 2.5
In the last line we have used the operator +=, which has the equivalent
effect of
1 C_value = C_value + 2.5
2.4 For-loops
This for loop runs over all elements in the list and in every pass the
variable C takes the value of one of the elements of the list. Again, the for
specification ends with a colon : and the lines that follow below with the
same indentation are executed in each cycle of the loop. Using the for
loop we can now take a list of temperature in Celcius and convert them
to Fahrenheit:
1 Cdegrees = [-20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40]
2 for C in Cdegrees:
3 F = (9.0/5)*C + 32
4 print(C, F)
Executing this piece of code gives a rather ugly looking table, but we
have seen before that the way numbers are printed can be influenced in
the print statement very easily:
2 Loops and lists 15
1 Cdegrees = [-20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40]
2 for C in Cdegrees:
3 F = (9.0/5)*C + 32
4 print(’%5d %5.1f’ % (C, F))
Finally let’s mention that if the sole purpose of the for loop is to create
the Fdegrees list, this can be achieved in Python using only one line !
1 Cdegrees = [-20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40]
2 Fdegrees = [ (9.0/5)*C + 32 for C in Cdegrees ]
This is very Pythonic ! In the second line a new list is create where each
element depends is built up from the list Cdegrees.
We have now encountered two different kinds of loops, for-loops and
while-loops. These loops differ in the way they are implemented but in
principle, either of the two can be used to produce any loop. In fact, in
many programming languages only a single type of loop construct is
available, and this works perfectly fine. The for-loop just above can be
programmed using a while-loop with exactly the same result:
1 Cdegrees = [-20, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35, 40]
2 index = 0
3 while index < len(Cdegrees):
4 C = Cdegrees[index]
5 F = (9.0/5)*C + 32
6 print(’%5d %5.1f’ % (C, F))
7 index += 1
This code is slightly longer any some would argue that it is less clear than
the implementation with the for loop, but both give exactly the same
result. Therefore, it is up to the programmer to decide which construct
to use and it is advised to take the version that is the most intuitive for
the specific application.
Using the range construct a for loop over a range of integers can be
written as:
1 for i in range(start, stop, step):
2 ...
This code will will print the value of z, for every combination of x and y,
where x and y range from 1 to 10.
The two print statement on lines 2 and 3 are only executed if C < 273.15 is
True, otherwise it is skipped and the code on lines 5 and 6 is executed. The
final printing statement on line 7 is always executed since the indentation
tells Python that it is not inside the block of statements in the if-else
statement. The else can also be skipped :
1 if C < -273.15:
2 print(’%g degrees Celcius is non-physical!’ % C)
3 F = 9.0/5*C + 32
4 print(F)
5 print(’end of program’)
In this case the conversion in Fahrenheit is not in the if-else block and is
therefore always executed, regardless of whether it makes sense or not.
2 Loops and lists 17
These statements have only given the choice between two options but it
is in fact possible using the elif keyword (short for else if) to ’branch’ the
program into as many different flows are needed:
1 if condition1:
2 <block of statements>
3 elif condition2:
4 <block of statements>
5 elif condition3:
6 <block of statements>
7 else:
8 <block of statements>
9 <next statement>
There are many other ways of using the if statement and some of these
will appear in the code that is used in further chapter in this course. In
most cases it is easy to understand how these constructs work because
as with a lot of statements in Python: they actually make some sense
in normal language too. It is also a good exercise to find your own
background information about the if-else statement since the online
documentation is virtually endless, including many examples. This is
true for all other Python constructs discussed in this chapter.
Let 𝑝 be a bank interest rate in percent per year. An initial amount 𝐴 has
then grow to
𝑝 𝑛
𝐴𝑛 = 𝐴0 (1 + ) (2.1)
100
after 𝑛 years. Make a program called ’[Link]’ that computes and print
on the screen how much money 𝐴0 = 1000 euros have grown to after
𝑛 = 1 , 2, 3, 4..., 𝑁 years with a 5% interest rate. To do that, make a for
loop over the value of 𝑛 and store the successive values of in a list
𝐴𝑛 = [1000 , ......]. Take 𝑁 = 50. If you don’t know how to start, try
something like the following:
1 A0 = ....
2 N = ...
3 p = ...
4 A = []
5 ....
6 for i in range(N):
7 [Link]( .... )
8 print(’year = %d \t money = %1.6f’ % (i, A[i]))
1 n = 0 # counter
2 N_MAX = 1000 # maximum number of iterations
3 An = 0. # initial value
4
As an egg cooks, the proteins first denature and then coagulate. When the
temperature exceeds a certain critical point, reactions begin and proceed
faster as the temperature increases. In the egg-white the proteins start to
coagulate at temperatures above 63 °C, while in the yolk the proteins start
to coagulate for temperatures above 70 °C. For a soft boiled egg, the white
needs to have been heated long enough to coagulate at a temperature
above 63 °C, but the yolk should not be heated above 70 °C. For a hard
boiled egg, the center of the yolk should be allowed to reach 70 °C. The
following formula expresses the time 𝑡 it takes (in seconds) for the center
of the yolk to reach the temperature 𝑇𝑦 (in Celsius degrees):
𝑀 2/3 𝑐𝜌1/3 𝑇𝑜 − 𝑇𝑤
𝑡= ln 0.76 (2.2)
𝐾𝜋2 (4𝜋/3)2/3 𝑇𝑦 − 𝑇𝑤
As an example consider the code below that calculates and prints the
factorial of the range of integer numbers from 1 to 10:
1 for n in range(1, 11):
2 result = 1
3 for m in range (1, n+1):
4 result = result*m
5 print(n, result)
In the previous subsection we have seen how to use one simple argument
as input of a function. Of course a function can take as many arguments as
you want. For example we can define a function that takes two numbers
as argument and return their product. This is simply done by the small
snippet of code below:
1 def myproduct(n, m):
2 return(n*m)
One very useful feature of Python is that some arguments can be defined
as optional. To do that you simply have to give a value to the optional
argument in the definition of the function.
1 def myfunc(w, a=0, b=2*[Link], n=101):
2 # the parameters a,b and n are here optional
3 # default values are used if these options are not specified
4 dx = (b-a)/(n-1)
5 for i in range(n):
6 x = a + i*dx
7 [Link]([Link](w*x))
8 return S
9
10 w = 0.1
3 Structured programming and functions 21
11 x = myfunc(0.1)
12 y = myfunc(0.1, b=4*[Link])
As you can see in the definition of the function the parameters 𝑎 ,𝑏 and 𝑛
are given a default value. Therefore we can omit them during the call of
the function as in line 8. In that case the default values of these parameter
are assumed. If on the other hand we want to use different values than
the default ones we just have to specify these values when calling the
function as on line 9. Note however that to specify a value we need to
explicitly write the name of the argument (here 𝑏 ) in the function call.
Such argument is generally called a keyword argument in Python
1 2
𝑦(𝑡) = 𝑣0 𝑡 − 𝑔𝑡 (3.1)
2
𝑦 0(𝑡) = 𝑣0 − 𝑔𝑡 (3.2)
Note that you also have to specify two values in the call of the function.
Note as well the use of the optional argument 𝑔 whose default values is
9.81 (like on Earth) but that we change to 3.711 (like on Mars) when we
call the function.
5
𝑇𝐶 = (𝑇𝐹 − 32) (3.3)
9
9
𝑇𝐹 = 𝑇𝐶 + 32 (3.4)
5
3 Structured programming and functions 22
The functions in the previous exercise take a single value as their argument
and return a single converted temperature. It is also possible to pass a list
of values to the functions, which can then also return a list of converted
temperatures. In this exercise you can rewrite one or both of the functions
written above so that they can convert a whole list of temperatures at
once. To achieve this the function should:
I establish the number of elements in the list that was passed to the
function: use the ’len()’ function.
I use a loop to convert all the elements in the list to other units and
store them all in a (second) list.
I return the list with converted values.
Mathematics and data analysis using
Python
Mathematics using Numpy 4
In the last chapters we have seen how to use computer programs to solve 4.1 Basics of Numpy . . . . . . . 24
simple problems. In this chapter we will turn to more dedicated methods 4.2 Numpy arrays . . . . . . . . . 24
that are specialized for scientific computing. Python is very well suited for 4.3 Arrays indexing and slicing 26
4.4 Vectorization . . . . . . . . . . 27
scientific use because specific modules have been written called 𝑛𝑢𝑚𝑝 𝑦
4.5 The copy issue . . . . . . . . . 29
and 𝑠𝑐𝑖𝑝 𝑦 . These modules contain a large number of functions and
4.6 Exercise: Slicing and Copying 29
methods to model and solve complex engineering problems. Together
with the visualization module 𝑚𝑎𝑡𝑝𝑙𝑜𝑡𝑙𝑖𝑏 these modules offer a very
powerful environment for scientific use.
The module 𝑛𝑢𝑚𝑝 𝑦 provides Python with very data structures that
can be used in mathematical programming. In particular, it allows the
definition of vectors and matrices in an easy way. On top of that it contains
a multitude of linear algebra routines. In this chapter we present the
main features of the numpy module. Additional functionality will be
introduced throughout the next chapters when needed. You can also
find all the information related to the numpy module in the online
documentation [Link]
As all modules, 𝑛𝑢𝑚𝑝 𝑦 needs to be imported before it can be used.
Therefore to use the numpy module in your code the following line must
be present before any call to any numpy routine
1 import numpy as np
Note that in the code listing above we specifically specify the type of
numbers, i.e. 64-bit floating point numbers, even though the numbers in
4 Mathematics using Numpy 25
16 print(c)
17 print(C)
18
19 C *=10
20
21 print (C)
In this code two Python lists are defined and these are subsequently used
to create two numpy-arrays by passing them to the [Link] function.
The difference between the arrays and lists becomes clear if they are
added up. The sum of two lists stored in the variable 𝑐 leads to a list
containing six numbers, 𝑐 = [1 , 2 , 3 , 4 , 5 , 6], i.e, the concatenation and not
the mathematical sum of the two list. In general, when doing mathematics,
we prefer our arrays of numbers to behave as vectors and this is not what
happens if we use lists. This is solved in numpy arrays that have many of
the characteristics of vectors (and matrices for two-dimensional arrays).
In this case the sum becomes 𝐶 = [5 , 7 , 9] as expected mathematically.
This example nicely illustrates what numpy does: it provides numerical
and mathematical operations for Python.
A large number of methods can be used to create 1D arrays. A few of
them are shown in the snippet of code below
1 v = [Link](5) # [0 0 0 0 0]
2 v = [Link](5) # [1 1 1 1 1]
3 v = [Link](5) # [0 1 2 3 4]
4
Note the extra set of brackets: the list contains three lists each contained
in their own brackets. This line of code creates the matrix
1 2 3
mat = 4 5 6® (4.2)
© ª
«7 8 9¬
where start/end define the first and last index in the range and increment
the distance between two elements of the range. For example
1 Y[0:10:2]
In the last line no indexes are give, just the ’:’, which means that all
elements are addressed. Finally, we can also use negative indexes to start
counting from the end of the array, for instance the last, and second-last
elements are accessed by:
1 Y[-1] # last element
2 Y[-2] # second to last
3 ...
4.4 Vectorization
In this code, the complete array is passed to the function [Link] and the
function returns a numpy-array with all the results.
While this leads to the exact same results as in the implementation using
a for-loop, this vector implementation is much faster. To illustrate this,
the time required to fill up the the array 𝑌 for different number of points
is listed in Table 4.1.
It is obvious from this table that it becomes quite attractive to use this
vectorized approaches, especially when the number of elements in an
array is very large.
4 Mathematics using Numpy 29
This code creates an array 𝑋 and second one 𝑌 that contains the 10 first
elements of 𝑋 . Note here the use of the method 𝑎𝑠𝑡 𝑦𝑝𝑒() to force the
numbers to be floating-point numbers. The instruction + = is encountered
for the first time here and simply results in adding a certain number
(here 1) to all the elements of 𝑌 . The results of this code is that we have
indeed added one to all the elements of the array 𝑌 , however, this way of
creating 𝑌 does not result in a physical copy in the computer memory.
Instead 𝑌 simply points to the same physical positions in memory as the
first ten elements of 𝑋 . Therefore, if we now print 𝑋 we will see that its 10
first elements have also been incremented by 1! If we want to circumvent
this issue we have to make sure that we explicitly reserve new memory
positions to store an actual copy of the numbers. This can be achieved in
multiple ways as shown in the code below:
1 # these are example of safe copies
2
3 # create an array
4 A = [Link](20).astype(’float64’)
5
We have barely scratched the surface of what 𝑛𝑢𝑚𝑝 𝑦 is capable of. Many
more possibilities will be during in the next chapters to complete the
presentation of numpy. A companion module of 𝑛𝑢𝑚𝑝 𝑦 , called 𝑠𝑐𝑖𝑝 𝑦
will also be introduced. Scipy provides even more sophisticated methods
for scientific computing, as for example the manipulation of sparse
matrices, and is intensively used.
Use a for-loop to print a table containing the first 10 values in the arrays
𝑇 , 𝑋 and 𝑌 on the screen. If you inspect the numbers in the arrays X and
Y you will find something is wrong. How can you fix this problem?
Plotting data with Matplotlib 5
5.1 Plotting data with Matplotlib 31
5.1 Plotting data with Matplotlib 5.2 Simple Curves . . . . . . . . . . 31
5.3 Contour plots . . . . . . . . . 32
The plotting of numerical data is a very important aspect of science 5.4 Exercise: Plotting a function 33
general and in scientific computing in particular. Plotting data in graphs 5.5 Extra exercises . . . . . . . . . 34
is a very quick way to check and analyze the results of experiments or Fourier series . . . . . . . . . . 34
Taylor expansion . . . . . . . 34
simulations. The module Matplotlib is one of many libraries that adds
plotting capabilities to Python. The most general set of plotting functions
for numerical data is found in the submodule pyplot, which will we use
here. Like any other (sub-)module, pyplot library has to be imported in
the python code before we can use its functions:
1 import [Link] as plt
Matplotlib is a very large library, so large that entire books have been
written to document the different styles and capabilities it offers. A
lot of the information about the module can be found in its online
documentation [Link] The documentation
contains many coding examples and tutorials to help you mastering
Matplotlib. In this chapter a few basic examples of of making plots are
discussed. In the following chapters a few more of the pyplot capabilities
will become clear, but we will use only a very small subset of the many
possibilities that are available.
4 X = [Link](0, 2, 100)
5
11 [Link](’Time’, fontsize=20)
12 [Link](’F(t)’, fontsize=20)
13
14 [Link](’My plot’)
15 [Link](loc=2)
16
17 [Link]()
The simple plot shown above is the simplest kind of plot that can be made
with Matplotlib and it will be sufficient for the most common situations.
However, Matplotlib contains a wealth of other plotting options. An
example is a so-called contour-plot that allows to plot a 3D function that
is a function of two arguments. An example of such a function is given
by:.
𝑥
𝑓 (𝑥, 𝑦) = (1 − + 𝑥 2 + 𝑦 3 ) exp(−𝑥 2 − 𝑦 2 ) (5.1)
2
9 n = 256
10 x = [Link](-2, 4, n)
11 y = [Link](-2, 4, n)
12
13 X, Y = [Link](x, y)
14
(𝑥 − 𝑏)2
𝑦(𝑥) = 𝑎 · exp − (5.2)
𝑐2
Fourier series
1 0 ≤ 𝑡 < 𝑇/2
𝑓 (𝑡) = 0 𝑡 = 𝑇/2 (5.3)
−1 𝑇/2 < 𝑡 ≤ 𝑇
𝑛
2(2 𝑖 − 1)𝜋𝑡
4X 1
𝑆(𝑡, 𝑛) = sin (5.4)
𝜋 𝑖=1 2 𝑖 − 1 𝑇
1. A function that calculates 𝑓 (𝑡). You will take 𝑇 = 2𝜋 and use 251
points between 𝑡 = 0 and 𝑡 = 𝑇 . This is the true step-function.
2. A function that calculates the fourier series 𝑆(𝑡, 𝑛) fan any value of
𝑛.
3. The calculation of 𝑆(𝑡, 𝑛) for n=2,5,10,30 and 100. Plot these different
approximations and compare to the true 𝑓 (𝑡) function in the same
graph.
Taylor expansion
𝑛
𝑥 2 𝑗+1
(−1) 𝑗
X
sin(𝑥) ≈ 𝑆(𝑥 ; 𝑛) = (5.5)
𝑗=0 (2 𝑗 + 1)!
5 Plotting data with Matplotlib 35
between the linear fit and the data points. In other words, we would like
to minimise the least squares error, 𝐸 𝑙𝑠 defined as:
𝑛
X
𝐸 𝑙𝑠 = (𝑦 𝑖 − 𝑝(𝑥 𝑖 ))2 (6.1)
𝑖=1
The function 𝑝(𝑦 𝑖 ) can be a polynomial of any order but as for the
purpose of this illustration we will restrict ourselves to a linear function.
We can therefore replace the function 𝑝(𝑦 𝑖 ) by the linear equation in
terms of 𝑐 0 and 𝑐 1 . The goal of linear regression is to minimise 𝐸 𝑙𝑠 , i.e. to
determine the values of the coefficients for which 𝐸 𝑙𝑠 is smallest. Since
eq. 6.1 is a quadratic equation we know that the minimum occurs where
the derivative of this function is zero. Therefore, we take the derivative
of eq. 6.1 with respect to the two unknown parameters:
𝑛
𝜕𝐸 𝑙𝑠 X
= −2 (𝑦 𝑖 − 𝑐 1 𝑥 𝑖 − 𝑐 0 ) = 0 , (6.2)
𝜕𝑐 0 𝑖=1
𝑛
𝜕𝐸 𝑙𝑠 X
= −2 ((𝑦 𝑖 − 𝑐 1 𝑥 𝑖 − 𝑐 0 ) · 𝑥 𝑖 ) = 0. (6.3)
𝜕𝑐 1 𝑖=1
𝑛
X 𝑛
X
(𝑦 𝑖 ) = 𝑐 1 (𝑥 𝑖 ) + 𝑛 · 𝑐0 , (6.4)
𝑖=1 𝑖=1
𝑛
X 𝑛
X 𝑛
X
(𝑦 𝑖 𝑥 𝑖 ) = 𝑐1 (𝑥 2𝑖 ) + 𝑐 0 (𝑥 𝑖 ). (6.5)
𝑖=1 𝑖=1 𝑖=1
6 Curve Fitting 38
This results in two linear equations with two unknown variables, 𝑐 0 and
𝑐 1 that can be solved easily by hand, giving:
P𝑛
(𝑥 𝑖 ) 𝑛𝑖=1 (𝑦 𝑖 )
P
P𝑛
𝑖=1 (𝑦 𝑖 𝑥 𝑖 )
𝑖=1
− 𝑛
𝑐1 = (6.6)
( 𝑛𝑖=1 (𝑥 𝑖 ))2
P
P𝑛 2
𝑖=1 (𝑥 𝑖 ) − 𝑛
P𝑛 P𝑛
𝑖=1 (𝑦 𝑖 − 𝑐1 ) 𝑖=1 (𝑥 𝑖 )
𝑐0 = (6.7)
𝑛
Although these equations look rather fearsome, in a computer program
it is very easy to calculate all the summations that are involved. In the
code below the summations are evaluated using the numpy routine
[Link], while summations over the product of the elements in
the arrays are evaluated using [Link]. The sum() function simply
returns the sum of all the elements in an array, while dot() calculates
the dot-product of two vectors, which is the sum of the products of the
elements of two vectors. After sorting out the sums, the evaluation of
the two coefficients and hence the best possible linear fit is calculated in
two lines of code (lines 16 and 17). After calculation the best possible
linear fit the code prints the two coefficients to the screen and plots both
the experimental data (using [Link]) and the linear fit. The
output of the program is shown in Figure ?? (right) and it is clear that it
yields a nice linear fit. The resulting fit is the best possible in the sense
that it minimises the least square difference between the measured data
and the fitted line.
1 import numpy as np
2 import [Link] as plt
3
10 Xsum = [Link](temp)
11 Ysum = [Link](yld)
12 XYsum = [Link](temp, yld)
13 X2sum = [Link](temp, temp)
14
19 print(c0, c1)
20
21 regression = c0 + c1*temp
22
This piece of code uses two internal numpy functions: [Link](X) and
[Link]. The first one computes the sum of all the elements contained in
the array X. The second one, i.e. [Link](X,Y) computes the dot product
of the two vectors passed in argument. In terms of equations we have
X
[Link](X) = 𝑋𝑛 (6.8)
𝑛
X
[Link](X,Y) = 𝑋𝑛 × 𝑌𝑛 (6.9)
𝑛
P𝑛 P𝑛
𝑛 𝑐0
𝑖=1
(𝑥 𝑖 ) 𝑖=1
(𝑦 𝑖 )
P𝑛 P𝑛 · = P𝑛 (6.10)
𝑖=1 (𝑥 𝑖 )
2
𝑖=1 (𝑥 𝑖 ) 𝑐1 𝑖=1 (𝑦 𝑖 𝑥 𝑖 )
21 print(c)
22
While it is very instructive to write the code for fitting a function to exper-
imental data, this is usually not necessary. In many software packages
(Origin, Igor, MATLAB) there are standard routines available for fitting
functions to any kind of data. Similarly, in the Python modules Numpy
and Scipy there is a variety of routines available for curve fitting. In the
section above we have written our own code for linear regression of any
data, no matter how large the data set is. For this type of fit there is a
simple routine available in numpy, called [Link]. The code below
uses this routine for the data that we have specified above and running
the code will result is exactly the same coefficients and plots as we have
seen above.
1 import numpy as np
2 import [Link] as plt
3
In most cases the data that you need to fit is not described by a polynomial.
The routine optimize.curve_fit from scipy can be used to fit any function.
This function can be as complicated as we like and it does not have to
be a simple analytical function. As shown in the example below, we can
define our own function (lines 5-6) and then just use that function, on
line 12, to fit it to a certain data set. In this example the data point are
generated using the same function and on line 10 a random gaussian
noise is added. The latter is just for making the data a bit noisy for the
purpose of the example. The resulting (noisy) data with the exponential
fit are shown in Figure 6.3.
9 x = [Link](0, 4, 50)
10 y = func(x, 2.5, 1.3, 0.5)
11 yn = y + 0.2*[Link](size=len(x))
12
15 print(popt)
16
The examples shown this section have consisted of quite nice data that
accurately follows a certain trend. This is not always the case and the
correct selection of data requires some attention. One of the very general
issues that are true for all curve fits is the fact that each point weighs
equally in the fit. This means that if there is one point out of 10 that is
very far off the trend, for instance a linear trend, the result can be a very
bad linear fit. An example is shown in Figure 6.4 where one of the points
in the linear regression example from above is moved way off the linear
trend. This results in a line that tries to include this point, resulting in a
line that clearly does not describe the data very well at all! A solution to
this is data selection before the fitting procedure is performed.
A related issue may occur in the fitting of nonlinear functions, for instance
an exponential function. Since in almost all curve fitting routines the
sum of the square of the difference between the data and the fitted
curve is minimised, it is the absolute difference that plays a role. For
an exponential function of which the values can span several orders of
magnitude this may result in a poor correspondence in the tail of the
exponential, where the absolute differences may still be small but the
relative difference can be very big. As a result, after fitting, the exponential
fit will follow the high points with reasonable accuracy but plotting it Figure 6.4: Illustration of the effect of out-
liers. One point that is far off the linear
on a log scale typically shows quite severe deviations in the tail of the trend can make the resulting linear trend
exponential. very bad.
−𝛽
𝑘 𝐶𝑇 (𝑑) = 𝐴0 𝑒 𝑑 , (6.11)
which can easily be recast in the form of a linear equation, plotting the
natural logarithm of the rate against the inverse of the distance:
1
ln 𝑘 𝐶𝑇 = ln(𝐴0 ) − 𝛽 . (6.12)
𝑑
Fitting of this linear relation can easily be done, for instance by linear
regression as described in detail above. The result will in most case be a
more evenly weighted fit where the low points are better described than
when just fitting an exponential function directly. The latter is particularly
6 Curve Fitting 43
In this exercise we will use the definition of the Gaussian function that
we did in Chapter 6 and make a fit to some experimental data. Some one
has done some measurements, obtaining a series of 10 x-values and a
corresponding series of 10 y-values:
1 import numpy as np
2
1. Use the lines above in your code to generate two numpy arrays.
2. Make a scatter-plot that visualizes the data.
3. Copy the function that you defined for the Gaussian in Chapter 6
and add it to your code.
4. Use the [Link] function 0 𝑐𝑢𝑟𝑣𝑒 _ 𝑓 𝑖𝑡 0 and your Gaussian
function to make the best fit to the data-points and print the values
of the parameters a, b and c to the screen.
5. Add a smooth Gaussian curve (100 points) to the scatter-plot that
you made above and confirm that the fitted curve nicely follows the
experimental data. Make sure you make a pretty plot with legend,
axis labels, a title and different colors for the scatter plot and the
fitted curve.
Reading and analyzing data 7
The rigorous statistical analysis of large amounts of data is an important 7.1 Reading data from a file . . . 44
aspect science and engineering problems. Experimental results often 7.2 Statistics and histograms . . 47
comes in the form of a large data set that needs further processing to 7.3 Exercise: fitting many curves 49
7.4 Extra exercises . . . . . . . . . 50
extract meaningful information. This further processing can be done in
Process Data . . . . . . . . . . 50
many ways, for instance simply in programs such as Excel or Origin,
or specialized software packages for statistical analysis such as SPSS.
Powerful tools have also been developed to perform statistical analysis
directly within Python. This is for example the case in the module
[Link]() that contains a range of methods to process complex data.
In this chapter we will discuss some very basic aspects of analyzing
large data sets. The subjects treated will be very useful and often directly
applicable to student research projects.
Before data can be processed in a Python code, we first need to read the
data from a file, and later maybe write the resulting processed data back
into another file. Several options are available to read/write data from
file depending on its format. The most common file formats are space-
separated values, where each line contains several numbers separated by
a blank space, and character-separated values (csv) where the numbers
are separated by a character, typically a comma or semicolon. This format
is very popular to share data between different programs. Excel can, for
example, export spreadsheets in .csv format.
Space-Separated Values Data files that contain column of data where
the different value on each line are separated by a space can be directly
read into your Python code by using the function 𝑛𝑢𝑚𝑝 𝑦.𝑙𝑜𝑎𝑑𝑡𝑥𝑡().
Similarly, you can write data to a file with the function 𝑛𝑢𝑚𝑝 𝑦.𝑠 𝑎𝑣𝑒𝑡𝑥𝑡().
As an example, we would like to import the data shown in Fig. 7.1, into a
Python code:
This data file contains, as the first column, a list of numbers running
from 0 , 1 , 2...., 𝑁 . The other columns in the file are filled with a series of
random numbers. If the name of this file is ’[Link]’, this data can be
loaded into a numpy array using the function [Link]():
7 Reading and analyzing data 45
1 matrix_data = [Link](’[Link]’)
This function will read the data from the file and return a two-dimensional
numpy-array that is assigned to matrix_data. This matrix contains all the
numbers present in the file. We can now access the data in the file, for
instance to plot it in a graph by taking different ’slices’ of it. For example
the first column can be accessed using the indexing that indicates all
lines in the first index (:) while only selecting the column index 0:
1 first_col = matrix_data[:,0]
Note that in this case the variable name first_col points to the same
values in memory as the first column in the array column_data, it is not
a copy!
Similarly, the values contained in a two-dimensional numpy-array can be
written to a file using the function [Link]():
1 [Link](’my_file.dat’,mat_data)
This function will generates a file called my_file.dat that contains all the
elements of the matrix mat_data. The save and read commands have a lot
of other argument as for example the format used to print each number
(%f,%e,...), but the default options work very well in most cases.
Character-separated values. The space separated values may seem like
the most logical format, but character-separated values formatting is
used by many programs and experimental equipment. These files cannot
be imported using the [Link]() function as they contains commas. A
typical .csv file is shown in Fig. 7.2
General file format. It is always possible to manually load all the values
contained in a file using the standard input/output functionality of
Python. This can be be done by reading the file line by line and storing
the values contained on each line in a List or a [Link]. This method
can be used to read any type of files with any formatting (or no formatting
at all). We illustrate this approach in the snippet of code below that reads
a .csv file with a header.
1 import numpy as np
2
4 def read_csv(filename):
5
This examples requires much more explanation since the loading of the
numerical data in the file takes place in many steps. The function open()
opens the file called filename, where the option ’r’ indicates that the file
is for reading only. This is also the default option. Other options such
as ’w’ allows writing to a file. The function [Link]() reads all the
lines contained in the file one by one and stores them in the variable
data_string. This variable is a 1D array where each element correspond
to one line of characters. Hence, data_string[0] is the first line stored as a
string of characters and so on. The next important function that is used
in the code above is line = data_string[iL].split(’,’). This function takes
one individual line ( data_string[iL]) and splits it at each comma ,. For
example the string ’1.00, 2.00, 3.00’ will become a 1D vector containing
3 separate elements: [’1.00’, ’2.00’, ’3.00’ ]. All these elements are still
strings of characters and not numbers. To convert these text-strings into
numbers we use the function float(x) that converts a string into a floating
point number if possible.
The approach of reading a file line-by-line and processing the data in the
form we want is a bit more complicated than the much simpler numpy
functions used above but it does allow full control over the reading of
the data. In addition, it is sometimes impossible to use the [Link]() or
[Link]() and the method shown here is the only solution capable
7 Reading and analyzing data 47
Once we have loaded the data into a proper numpy array we can start
processing its elements. In the file presented above, the first column
represent the time of the experiment and the each following column
contains the value of certain physical quantity measured at these times.
We can therefore plot the time evolution of this data using matplotlib
as shown in Fig. 7.3.
1 import numpy as np
2 import [Link] as plt
3
4 DATA = [Link](’[Link]’)
5 X = DATA[:, 0]
6 Y = DATA[:, 1:]
7
11 [Link](’time’, fontsize=15)
12 [Link](’Measurements’, fontsize=15)
13 [Link]()
The data looks very noisy and not much information can be extracted
from inspection of the plots. We can just observe that the mean value of
the blue curve seems lower than that of the red ones. To obtain a more
quantitative analysis of the data we can analyze the data in some more
detail. Different modules such as Numpy and Scipy offers a vast library
of statistical function that can be accessed easily. For example the mean
values of the elements of a vector 𝑋 defined as
1 𝑁−
X1
𝜇= 𝑋[𝑛] (7.1)
𝑁 𝑛=0
s
𝑁−
1 X1 1 𝑁−
X1
var = (𝑋[𝑛] − 𝜇)2 𝜎= (𝑋[𝑛] − 𝜇)2 (7.2)
𝑛 𝑛=1 𝑛 𝑛=1
where 𝜇 is the mean value of the vector. These two quantities can be
calculated for the values in a numpy-array can be calculated using the
functions [Link]() and [Link]() or the equivalent syntax [Link](X) [Link](X).
These quantities define how the values are dispersed around the mean
values and this plays a key role in error analysis. Finally, the median
of a time series, i.e. the value that is superior to half of the number
contained in the series and inferior to the second half is returned by the
function [Link](X).The code below illustrates the basic usage of these
functions
1 nLine, nCol = [Link]
2 for iCol in range(nCol):
3 print ’\n === Data Column %02d’ %(iCol+1)
4 print("\tMean : %1.6f" %(Y[:,iCol].mean()))
5 print("\tMedian : %1.6f" %([Link](Y[:,iCol])))
6 print("\tMinimum : %1.6f" %(Y[:,iCol].min()))
7 print("\tMaximum : %1.6f" %(Y[:,iCol].max()))
8 print("\tVariance : %1.6f" %(Y[:,iCol].var()))
9 print("\tStd deviation : %1.6f" %(Y[:,iCol].std()))
This code will produce the values reported in the table below for the first
two columns of data.
The calculation of these quantities confirms our observations. The mean
values of the first two columns are -0.89 and 0.06, respectively. In addition,
the variance of the first time series is smaller than for the second indicating
a smaller spread of values around the mean. Finally, we can see that the
mean and median values are equal for the first time series but not for the
second, which indicates that the shapes of the distributions of random
values are quite different.
The same insights can also be obtained if we plot these data, but not as a
function of the time but as a histogram. Matplotlib offers an easy way to
construct such a histogram. Using the function [Link](X,bins=nbin,...)
will directly generate a histogram of the data contained in X using a
certain number of intervals (or bins) to compute the count. The use of
this function is illustrated in the code below.
1 import numpy as np
2 import [Link] as plt
3
9 # plot histograms
10 [Link](Y[:, 0], bins=50, facecolor=’#00AFFF’, alpha=0.5)
11 [Link](Y[:, 1], bins=50, facecolor=’#FF003F’, alpha=0.5)
12 [Link](’Count’, fontsize=15)
7 Reading and analyzing data 49
13 [Link](’Measurements’, fontsize=15)
14 [Link]()
In this code the function [Link]() accepts as the most important argument
a numpy array containing the data to be plotted. In addition, different
optional arguments can be supplied to specify the number of ’bins’ to be
used or the color of the histogram and its level of transparency (alpha).
This generates the plot in Fig. 7.4.
Plotting the histogram gives a quick insight in the main differences
between the two data sets. The blue distribution seems quite symmetric
around its mean value, which explains why its mean is identical to its
median. In contrast, the red distributions is not symmetric, leading to the
difference between its mean and median value. These two distributions
consequently look like a normal (Gaussian) and a Poisson distribution,
respectively.
𝑡
𝐼(𝑡) = 𝐴 · exp − (7.3)
𝜏𝑓 𝑙
where t is time, is the fluorescence life time and A is a pre-
exponential factor.
4. Use the [Link] function 0 𝑐𝑢𝑟𝑣𝑒 _ 𝑓 𝑖𝑡 0 (see Chapter 6) to fit
this one data set to the exponential function that you defined above.
Make a plot where you combine the scatter plot you made above
with the optimal fit.
7 Reading and analyzing data 50
Process Data
To facilitate your analysis the company has built 250 devices on three
different days. They’ve measured the light coming out of these devices
and gave you 3 space separated files ’exp_data_1.dat’, ’exp_data_2.dat’
and ’exp_data_3.dat’ that you can download on the Blackboard. Each
file contains 251 columns where the first column is the time (i.e. the
x-axis) and the 250 others are the light intensity coming out of the 250
devices built and tested that day. The company wants you to estimate
the distribution of decay parameters obtained for all these devices. To
analyze the data create a script called ’[Link]’. In this script you
will do for each file:
I Plot the histogram of the k-value. You should obtain a figure similar
to Fig. 7.6 (we show here the distributions obtained for the three
files.)
I Estimate the mean value and standard deviation of each of the
distribution.
𝑉 = 𝑣0 𝑣1 ... 𝑣𝑁
(8.1)
where 𝑣 0 is the first element of the vector, 𝑣 1 the second .... As we have
already seen such a numpy-array can be constructed using the numpy-
function-array() that takes a Python-list of numbers as an argument:
𝑎 0 ,0 𝑎 0 ,1 ... 𝑎0,𝑁
𝑎 1 ,0 𝑎 1 ,1 ... 𝑎1,𝑁 ®
© ª
𝐴 = . .. .. ®® (8.2)
.. . ... . ®
« 𝑎 𝑁 ,0 𝑎 𝑁 ,1 ... 𝑎 𝑁 ,𝑁 ¬
Where 𝐴 is the name of the matrix and 𝑎 𝑖,𝑗 its element on the 𝑖 -th row
and 𝑗 -th column. Matrices can be constructed using the numpy-function
array() that takes a list filled with lists as its argument in this case. For
example the code:
1 M = [Link]([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
1 2 3
M = 4 5 6® (8.3)
© ª
«7 8 9¬
8 Vectors and Matrices 54
Note that the argument given to the array() function is a list of which the
elements are lists themselves. Therefore, there is an extra set of square
brackets.
𝑣 contains the values of 𝑢 ranging from the first to the second index
(remember the index starts at 0). In general the statement:
1 v = u[a:b]
extracts the elements of 𝑢 ranging from 𝑎 upto (but not including) 𝑏 and
stores them in 𝑣 . If you omit the 𝑎 or 𝑏 , the default values 𝑎 = 0 and
𝑏 = 𝑁 will be used, where 𝑁 is the total length of the vector. We have
therefore:
1 u = [Link](10) #u = [0 1 2 3 4 5 6 7 8 9]
2 v = u[:4] #v = [0 1 2 3]
3 v = u[6:] #v = [6 7 8 9]
4 v = u[:] #v = [0 1 2 3 4 5 6 7 8 9]
5 v = u[1:-2] #v = [1 2 3 4 5 6 7]
On line 5 a negative index is used which means that the second element
from the end of the array is indicated. The index -1 corresponds to the
last element of the array, -2 the second last and so on.
refers to the element on the 2nd row and 2nd column of the matrix. To
extract a complete row of a complete column of a matrix we can use:
1 A = [Link](25)
2 mat = [Link](A,(5, 5))
3 col = mat[:, 2]
4 line = mat[3, :]
The first two lines of this code serve to generate a two-dimensional matrix.
First, we use the function arange() to create a linear array. This array is
transformed into a two-dimensional array using the function reshape()
that takes as its second argument a ’shape’ that consists of two numbers.
These two numbers indicating the shape have to be contained between
8 Vectors and Matrices 55
This statement extracts a 3-by3 sub-matrix containing only the rows and
column ranging between from the first (with index 0) to the third (with
index 3).
If the indexes of the elements to be extracted are not contiguous, the
numpy function 𝑖𝑥 can be used to extract them. This function constructs
a new matrix taking elements from an existing matrix where the ’mesh’
of elements has to be supplied as an argument. This mesh consist of two
Python lists that contain the row and column indexes to be extracted and
can be used as a ’mask’ to extract certain values to form a sub-matrix. An
example is shown in the code below.
1 mat = [Link](25).reshape(5, 5)
2 ind1 = [0, 2, 4]
3 ind2 = [0, 2, 4]
4 sub = mat[np.ix_(ind1, ind2)]
The first line of the code shows an alternative use of the reshape() function
in which is is directly appended to the arange function that generates a
linear array. The effect is exactly the same as in the example above where
the reshape function was called separately. On the second and third lines
of the code the row and column indexes are indicated in a Python list.
On line 4, these lists are passed as arguments to the [Link]() function
to generate the mesh of indexes to be assigned to the sub-matrix. The
codes eventually leads to the following matrices:
0 1 2 3 4
5 6 7 8 9® 0 2 4
© ª
mat = 10 11 12 13 14® −→ sub = 10 12 14® (8.4)
® © ª
15 16 17 18 19® «20 22 24¬
®
«20 21 22 23 24¬
The sub-matrix is generated by taking only the elements on the 0th, 2nd
and 4th row and 0th, 2nd and 4th [Link] this case the same indexes
have been chosen for the columns and rows, but they can be different
also, as long as they do not contain an index larger than the size of the
matrix.
It is finally possible to change the values of the elements contains in a
row, a column or a part of a matrix, for instance by a random number as
in the example below.
1 mat = [Link](25).reshape(5, 5).astype(float)
2
3 mat[:, 2] = [Link](5)
4 mat[1, :] = [Link](5)
5
6 ind2 = [0, 2, 4]
7 mat[np.ix_(ind1, ind2)] = [Link](3, 3)
8 Vectors and Matrices 56
On the first line we have used the function ’astype(float)’, which forces
Python to interpret the number in the matrix as floating point num-
bers rather integers, which is the default for [Link](). The function
[Link]() generates a series of random numbers between 0
and 1. The argument, in this case 5, indicates how many random numbers
are to be returned. The conversion of the elements in the matrix ’mat’
to floating point numbers is important here because otherwise, all the
random numbers would be truncated to make them an integer and hence
they would all be equal to 0.
The two components of the ’shape’ are pass with an extra set of brackets,
(𝑁 , 𝑀), to indicate that they together form a single argument! On the
last line of the code above we use the function [Link]() to generate
an identity matrix, which is a square matrix by definition and hence no
information other than the size is needed.
This code contains two loops: one loop inside another loop, or nested
loops. The first loop, called outer loop, runs over all the rows, while
the second one (the inner loop) runs over the columns. This code works
perfectly fine and gives full control over what happens, but in many cases
there are faster ways of doing the same thing by exploiting vectorization
in numpy routines. As an example, the matrix created above can also be
obtained by calculating the outer product of two vectors. If we have two
vectors u and v that are 𝑀 × 1 and 𝑁 × 1 vectors, then the outer product
of these two vectors is given by:
𝑢1
𝑇
©𝑢1 𝑣1 𝑢1 𝑣 2 𝑢1 𝑣 3
u ⊗ v = uv = 𝑢2 𝑣 1 𝑣2 𝑣 3 = 𝑢2 𝑣1 𝑢2 𝑣 2 𝑢2 𝑣 3 ®
(8.5)
ª
𝑢3
«𝑢3 𝑣1 𝑢3 𝑣2 𝑢3 𝑣 3 ¬
For two one dimensional vectors the function [Link]() has the same
effect as the function [Link]() that we have seen earlier.
This complex operation can be performed with the function 𝑛𝑢𝑚𝑝 𝑦.𝑚𝑎𝑡𝑚𝑢𝑙()
or 𝑛𝑢𝑚𝑝 𝑦.𝑑𝑜𝑡():
8 Vectors and Matrices 58
1 m = 10
2 n = 15
3 p = 5
4 A = [Link](n,m)
5 B = [Link](m,p)
6 C = [Link](A,B)
𝑢0 𝑎 0 ,0 𝑎 0 ,1 ... 𝑎0,𝑁 𝑣0
𝑢1 ® 𝑎1,0 𝑎 1 ,1 ... 𝑎1,𝑁 ® 𝑣1 ®
© ª © ª © ª
. ®= . .. .. ®® · .. ®® (8.8)
. ® . .
. ® . ... . ® . ®
«𝑢 𝑁 ¬ « 𝑎 𝑁 , 0 𝑎 𝑁 ,1 ... 𝑎 𝑁 ,𝑁 ¬ «𝑣 𝑁 ¬
where
𝑁
X
𝑢𝑖 = 𝑎 𝑖,𝑗 · 𝑣 𝑗 (8.9)
𝑗=0
where the elements of the vector u are calculated one by one. For each
element, the inner product of the 𝑖 -th row of A with v. This for-loop can
be included in a Python program as shown below.
8 Vectors and Matrices 59
1 import numpy as np
2
5 nLine = [Link][0]
6 nCol = [Link][1]
7 sizeV = [Link][0]
8
10 if nCol != sizeV:
11 print("Error : Size inconsistent")
12 return
13
14 u = [Link](nLine)
15 for iL in range(nLine):
16 u[iL] = [Link](A[iL, :], v)
17 return u
18
19 A = [Link](4, 5)
20 v = [Link](5)
21
22 u = mat_vect(A, v)
23 ucheck = [Link](A, v)
24
25 print([Link](u-ucheck))
𝑎 𝑏
Det(M) = | M | = = 𝑎·𝑑−𝑐·𝑏 (8.10)
𝑐 𝑑
𝑎 𝑏 𝑐
𝑒 𝑓 𝑑 𝑓 𝑑 𝑒
|M| = 𝑑 𝑒 𝑓 =𝑎· −𝑏· +𝑐· (8.11)
ℎ 𝑖 𝑔 𝑖 𝑔 ℎ
𝑔 ℎ 𝑖
the first element ( 𝑎 ) the minor matrix is the right-bottom 2x2 matrix that
remains, etc. Not that the sign in front of the co-factor alternates. Now
the calculation of the determinant of this 3x3 matrix has been reduced to
the calculation of determinants of three 2x2 matrices, which was defined
in Eq.8.10. The same approach can be used for matrices of higher order,
decomposing them step by step until they are given only in the form of
2x2 matrices. The general approach for the decomposition of matrices of
any order is given by the Laplace-expansion, in this case along the 𝑖 -th
row:
𝑛
(−1)𝑖+𝑗 𝑚 𝑖𝑗 𝐴 𝑖,𝑗
X
|M| = (8.12)
𝑖=1
The 𝑚 𝑖 𝑗 are element of the matrix M that in this case run along one of the
lines of the matrix. The matrix A is the (𝑛 − 1) x (𝑛 − 1) ’minor’ matrix
that does not contain the line and column of the current element 𝑚 𝑖 𝑗 .
Special attention should be paid to the factor (−1)𝑖+𝑗 that defines the sign
of the different contributions. As seen in the evaluation for the 3x3 matrix
the sign alternates along the line that is used for the decomposition. Note
that we use the common convention in mathematics here, which implies
that the indexing of matrices starts at 1. This is different from the numpy
indexing of arrays that starts at 0!
The evaluation of determinants of higher order square matrices quickly
becomes tedious and prone to small mistakes. However, it is done by a
systematic step-wise approach that is very suitable to do in a Python
code. In fact, it is an excellent example that can be solved by a so-called
recursive methods. In a recursive method we use a function that actually
calls itself, generally multiple times. This is exactly what we do in the
code below: it decomposes the matrix for which we want to calculate
the determinant into its minors and then we call the same function to
calculate the determinant of those minors! This sequence continues until
we arrive at a 2x2 matrix for which the determinant is easily calculated.
1 import numpy as np
2
3 def my_determinant(M):
4 rows = [Link][0]
5 columns = [Link][1]
6 # Start by checking whether the matrix is square
7 if rows != columns:
8 print(’this is not a square matrix’)
9 return(0)
10
15 # In all other cases: create submatrix, and add the terms up with the
16 # appropriate prefactor
17 # the function my_determinant is called recursively until we end up at the
18 # 2x2 matrix
19
20 else:
21 sum = 0
22 for i in range(columns):
23 submatrix = [Link]((rows-1, columns-1))
24 submatrix[0:, :i] = M[1:, :i]
25 submatrix[0:, i:] = M[1:, i+1:]
8 Vectors and Matrices 61
29
30
34 determinant = my_determinant(matrix)
35 print(’My determinant:’, determinant)
36
37 det = [Link](matrix)
38 print(’linalg:’, det)
6 C[iR,iC] = ......
7
8 return C
To test the function, generate two random matrices with the statements:
1 N, M, P = 100, 50, 150
2 A = [Link](N, M)
3 B = [Link](M, P)
This leads to the creation of two random matrices. You can check that
your function works correctly by comparing the result it gives to the
result provided by the standard numpy function to multiply matrices:
C = [Link](A,B). In addition you can time the performance of the function
with two nested loops, by using the function clock() from the module
time.
1 import time
2 start = [Link]()
3 # put the instruction you want here
4 end = [Link]()
5 print(’Calculation done in %f sec’ % (end-start))
𝑎1 𝑎2 ... 𝑎𝑛 𝑥1 𝑓1
𝑏 1 𝑏2 ... 𝑏 𝑛 ® 𝑥 2 ® 𝑓2 ®
© ª © ª © ª
. .. .. ®® · .. ®® = .. ®® (9.5)
. .
. .® . ® .®
«𝑧1 𝑧2 ... 𝑧 𝑛 ¬ « 𝑥 𝑛 ¬ « 𝑓𝑛 ¬
The first matrix contains all the constant multiplication factors. This
matrix is multiplied with the middle vector that holds the unknowns
and the right sides is a vector that holds the solution of each equation.
Such a problem is often written as A · x = f. Given a set of constants 𝑎 𝑖 ,
𝑏 𝑖 ... and 𝑓𝑖 , the goal is to find the values of all the different unknowns.
Similar linear equations can be written down for ethane and propane
and the overall system of linear equations can be written in a matrix form
A · x = f. This results in the matrix equation:
The system of equations 9.7 can be solved directly using a numpy function.
The numpy has a library (or sub-module) that contains many common
linear algebra functions. The most versatile sub-module is [Link]
that can be imported in a Python program in the usual way. In the code
below we rename the [Link] sub-module to ’nl’ while loading to
9 Systems of linear equations 65
reduce the amount of typing we have to do. We can then use function
[Link]() to solve systems of linear equations:
1 import numpy as np
2 import [Link] as nl
3
The most well-known 𝑑𝑖𝑟𝑒 𝑐𝑡 solution for solving sets of linear equations
is the Gauss-Jordan algorithm, also referred to as Gaussian elimination.
This approach relies on the use of an augmented or extended matrix,
made up of the matrix 𝐴 but to which the column vector 𝑓 is added
on the right side. The Gauss-Jordan algorithm transforms this matrix
in two steps to a matrix where all the elements of the original matrix
𝐴 are zero, except for the diagonals that become 1 after the procedure.
After this is achieved, the last column will contain the solution vector x.
This sequential transformation of the augmented matrix is schematically
indicated by:
9 Systems of linear equations 66
As seen in the above equation the algorithm takes place in two steps. First
we transform 𝐴 into an upper-triangular matrix, i.e. a matrix containing
non-zero elements only on and above its diagonal. In the second step,
we back-substitute the solution of the equation to obtain the solution
of the linear system. To perform these matrix transformations, only 3
operations, that do not change the solutions of the linear system are
allowed. These operations are:
This approach can be readily applied to the linear system in Eq. 9.7 for
the distillation column. The augmented matrix of this system reads:
We will refer to the first second and third row (line) of this matrix by 𝐿1 ,
𝐿2 and 𝐿3 , respectively. Our first goal it to eliminate the term 𝑏 1 = 0.1 in
𝐿2 . To do so, we perform the following operation: 𝐿2 = 𝐿2 − 00..91 𝐿1 . The
result of this operation is shown in the second matrix of Eq. 9.11, where
𝑏1 = 0. We can then eliminate the term 𝑐 2 = 0.2 from this new matrix by
performing 𝐿3 = 𝐿3 − 00..462
𝐿2 , leading to the last matrix in Eq. 9.11. Note
that 𝑐 1 happened to be already equal to 0, otherwise we should have also
made that zero by adding a certain factor time 𝐿1 , in addition to adding
a factor time 𝐿2 .
In Eq. 9.11 all floating point numbers have here been truncated to two
decimals to save space. Using the procedure above we have transformed
the matrix 𝐴 into an upper triangular matrix. We can now use the
so-called back-substitution approach to arrive at an ’identity matrix’
9 Systems of linear equations 67
containing only ones on the diagonal and zeros everywhere else. This
can be done systematically by starting from the bottom. We therefore
start with 𝐿3 = 𝐿3/0.61, leading directly to 𝑚3 = 1.15 as seen in the
first matrix in Eq. 9.12. We can use that new matrix to solve for 𝑚2 by
performing 𝐿2 = 1/0.46(𝐿2 − 0.18 ∗ 𝐿3 ). Finally we solve for 𝑚1 by doing:
𝐿1 = 1/0.90(𝐿1 − 0.30𝐿2 − 0.10𝐿3 ).
As expected, the last column now contains the solution of this linear
system. Although the method that is used by [Link]() is a bit more
complex than the Gauss-Jordan algorithm shown, here it relies on the
same general principle and of course leads to exactly the same solution.
These two stages of the solution are somewhat different and therefore the
implementation requires that we write two functions that each perform
one of these two operations.
A simple implementation of the first stage of the Gauss-Jordan algorithm
(i.e. obtaining an upper diagonal matrix) is shown in the code below.
1 import numpy as np
2
3 ########################################################
4 # Function to decompose the matrix
5 # in an upper triangular form
6 ########################################################
7
8 def LU(A,f):
9
23
The function LU() receives a matrix A and the right hand side vector ’f’ as
arguments. In the first few lines of code the dimensions of the matrix are
checked. Subsequently, a loop over all the columns is performed (with
exception of the last one) to eliminate all the terms below the diagonal
using the method discussed above. Note that we also make sure that
the diagonal element M[iC, iC] is not null as it would lead to an error
when we try to divide by M[iC, iC]. A more advanced version of the
Gauss-Jordan algorithm is required to treat such cases (interchanging
the different lines would solve this problem). Once all the terms are
eliminated we return the modified augmented or extended matrix M as
well as an integer value called ’info’, that equals 1 if the elimination was
successful and 0 otherwise.
The implementation of the second stage, the backward substitution, is
listed in the code below.
1
2 ########################################################
3 # Function to backsubstitute the results
4 # and get the final solution
5 ########################################################
6 def BS(M):
7
21
25 # loop over all the lines that are above this one
26 for iLL in range(iL-1, -1, -1):
27 M[iLL, :] -= M[iLL, iL]*M[iL, :]
28
29 info = 1
30 return M, info
31
32
33 ########################################################
34 ########################################################
17 # backsubstitute
18 M, info = BS(M)
19
20 print(M[:, -1])
Computational cost
3 import gauss_jordan as gj
4 import time
5 import [Link] as plt
6
10 # create list
11 cpu_time_numpy = []
12 cpu_time_mycode = []
13
21 # nummpy
22 t0 = [Link]()
23 [Link](A, f)
24 cpu_time_numpy.append([Link]()-t0)
25
26 # mycode
27 t0 = [Link]()
28 M,info = [Link](A, f)
29 M,info = [Link](M)
30 cpu_time_mycode.append([Link]()-t0)
31
The results of this code are shown in Fig. 9.2. As you can see finding the
solution of relatively big linear system is rather inexpensive. However,
the graph also shows clearly that [Link]()is considerably more
efficient than our simple Gauss-Jordan implementation. Numpy uses a
different algorithm to solve linear systems that is more complicated, but
the Gauss-Jordan algorithm is a very good illustration.
From these equations we can write for each of the values of 𝑚 in terms of
all the other values in the vector m. Of course at the start of the algorithm
we do not know any of the values in this vector yet, but we will guess
them at first. Therefore, we can now write these equations for each 𝑚 𝑖 in
terms of the values of the values of all other values in m as they were in
the previous cycle (or at the starting point in terms of the initial guess:
!
(𝑘+1) 1 X (𝑘)
𝑚𝑖 = 𝑓𝑖 − 𝑎 𝑖𝑗 𝑚 𝑗 (9.13)
𝑎 𝑖𝑖 𝑗≠𝑖
1
(1)
𝑚1 = (30.0 − 0.3 ∗ 20 − 0.1 ∗ 20) = 24.44 (9.14)
0.9
(1) 1
𝑚2 = (25.0 − 0.1 ∗ 20 − 0.2 ∗ 20) = 38.00 (9.15)
0.5
(1) 1
𝑚3 = (10.0 − 0.0 ∗ 20 − 0.2 ∗ 20) = 8.57 (9.16)
0.7
If we compare the the new value of the masses we can see that they are
collectively closer to the exact solution than our initial guess. We then
take these updated values and insert them as the new initial guess and
perform the same operation:
(2) 1
𝑚1 = (30.0 − 0.3 ∗ 38.00 − 0.1 ∗ 8.57) = 19.71 (9.17)
0.9
(2) 1
𝑚2 = (25.0 − 0.1 ∗ 24.44 − 0.2 ∗ 8.57) = 41.68 (9.18)
0.5
(2) 1
𝑚3 = (10.0 − 0.0 ∗ 24.44 − 0.2 ∗ 38.00) = 3.43 (9.19)
0.7
which gives updated values that are even closer to the exact solution. If
we keep iterating that way the solution will converge quickly to values
that are really close to the exact solution that we obtained from the
Gauss-Jordan method. In the exercise below we will write our own code
to implement the Jacobi algorithm and compare its performance to exact
solutions for different linear systems.
To ensure that the while-loop stops even if the solution is not converging,
which can sometimes happen, we typically impose a maximum number
of iterations 𝐼𝑇𝐸𝑅𝑀𝐴𝑋 and combine several conditions to stop the
𝑤 ℎ𝑖𝑙𝑒 loop:
1 while (res>tol) and (niter<ITER_MAX):
2 # block of instruction
3 # ...
4 res = .....
5 niter += 1
The loop will continue as long as (𝑟𝑒 𝑠 > 𝑡𝑜𝑙) and 𝑛𝑖𝑡𝑒𝑟 < 𝐼𝑇𝐸𝑅𝑀𝐴𝑋).
Do not forget to update 𝑟𝑒 𝑠 and increment 𝑛𝑖𝑡𝑒𝑟 within the loop otherwise
it will never stop! Test your Jacobi function to solve the linear system in
Eq. 9.7.
9 Systems of linear equations 74
Lake contamination
𝑘𝑔
180 = 𝑄 𝑆𝐻 𝐶𝑆 (9.20)
𝑦𝑟
Once you have written these equations rearrange them so that only term
without any unknown appears on the right-hand side. You can then solve
this system of equations with one of the methods seen in this chapter
and calculate the concentration of PCB contained in each lake.
An environmental organization is considering to build a bypass that
would go directly from lake Michigan to Lake Ontario with a flow rate of
20 𝑘𝑚 3 /𝑦𝑟 in order to reduce the concentration of PCB in lake Michigan.
Examine the effect of this bypass by rewriting the set of linear.
!
(𝑘+1) 1 X (𝑘+1)
X (𝑘)
𝑚𝑖 = 𝑓𝑖 − 𝑎 𝑖𝑗 𝑚 𝑗 − 𝑎 𝑖𝑗 𝑚 𝑗 (9.21)
𝑎 𝑖𝑖 𝑗<𝑖 𝑗>𝑖
The difference with the Jacobi method is that we always use the most
recent information we have of the solutions, even if they have not been
9 Systems of linear equations 75
(𝑘+1)
calculated for all values yet. For example, the calculation of 𝑚2 in
(𝑘)
the Jacobi method only depends on 𝑚1,3 . However, in the Gauss-Seidel
(𝑘+1) (𝑘+1) (𝑘)
method 𝑚2 is calculated from 𝑚1 (that we compute from 𝑚2,3 ) and
(𝑘)
𝑚3 . Implement your own version of the Gauss-Seidel algorithm and
compare its performance to the Jacobi-method. Note: the implementation
largely follows the Jacobi approach, but taking into account that the
solution-values from the current cycle are used where available.
Integration, Derivatives and
Non-linear systems
Numerical integration 10
The evaluation of the definite integral (or quadrature) of a function, as 10.1 Equal intervals . . . . . . . . . 77
shown in EQ. 10.1 is an operation that occurs in many problems, and Rectangle approximation . . 77
is a central issue in the integration of differential equations. Analytical Trapezoid approximation . . 79
Simpson’s rule . . . . . . . . . 80
solutions for the evaluation of such integral are sometimes available, but
Newton-Cotes formulas . . . . 81
in many cases they are difficult to derive, and often impossible. A specific
10.2 Gaussian quadrature . . . . . . 81
case where analytic integration is impossible is in cases where the data
10.3 Multiple Integrals . . . . . . 82
is obtained from experiments, i.e. only a series of tabulated values. The 10.4 Exercises . . . . . . . . . . . . . 83
numerical integration of functions is straightforward and in this chapter Composite methods . . . . . 83
we will discuss some common, efficient approaches for integration. Simpson’s rule . . . . . . . . . 85
∫ 𝑏
𝐼= 𝑓 (𝑥)𝑑𝑥 (10.1)
𝑎
In general, integrals can be calculated by summing up the areas of a
series of geometric shapes, such as rectangles or trapezes. The summation
is done in intervals and in principle an arbitrary accuracy can be reached
by making the size of the intervals smaller and smaller. In this chapter
we will consider two of such ’composite’ methods. We will first consider
approaches that use evenly spaced intervals, based on the so-called
Newton-Cotes formulas. Subsequently, we will consider more efficient
methods where the sizes of the intervals vary and are optimized to
give the best accuracy. Such methods, of which Gaussian quadrature is
an example, are more difficult to implement but several functions are
available in numpy to perform such integrations.
Rectangle approximation
∫ 𝑏 𝑁−
X1
𝑓 (𝑥)𝑑𝑥 ≈ ℎ 𝑓 (𝑥 𝑗 ) (10.2)
𝑎 𝑗=1
Three of the arguments of this function are optional: the ones that defines
the integration interval and the number of points in which the interval
in divided. For all the points in the interval, except for the last one, the
function adds up the area of all the rectangle, ℎ × 𝑓 (𝑥) and stores them in
the variable the variable 𝑟𝑒 𝑠 . An example of how this function is called
is given by:
1 A = int_rect([Link], a=0, b=[Link], N=1001)
∫ 𝑏 𝑁−
X1
𝑥 𝑗 + 𝑥 𝑗+1
𝑓 (𝑥)𝑑𝑥 ≈ ℎ 𝑓 (10.3)
𝑎 𝑗=1
2
10 Numerical integration 79
Trapezoid approximation
" #
𝑏
ℎ 𝑁−
X1
∫
𝑓 (𝑥)𝑑𝑥 ≈ 𝑓 (𝑥 𝑗 ) + 𝑓 (𝑥 𝑗+1 ) (10.4)
𝑎 2 𝑗=1
Simpson’s rule
The integration can now be performed for the two intervals at once and
it is easily derived that this integral is defined as:
𝑖+2
ℎ
∫
𝑓 (𝑥)𝑑𝑥 ≈ 𝑓 (𝑥 𝑖 ) + 4 𝑓 (𝑥 𝑖+1 ) + 𝑓 (𝑥 𝑖+2 )
(10.5)
𝑖 3
∫ 3
1
𝑑𝑥 = ln(3) (10.6)
1 𝑥
In Table 10.1 error for the numerical evaluation of this integral is given for
both the trapezoid rule and for Simpson’s rule. It is clear that Simpson’s
rule is by far superior and quickly approaches the numerical accuracy of
the computer itself. Importantly, the computational cost of Simpson’s rule
is similar as that for the trapezoidal approach, but for 𝑁 = 501 the error
is 5-6 orders of magnitude smaller! Therefore Simpson’s rule is generally
a high-accuracy method at a modest computational cost, making it a
method of great practical use.
10 Numerical integration 81
N Error trapezoid Error Simpson Table 10.1: Error in the evaluation of the
integral in Eq.10.6 for the trapezoidal ap-
5 1.8054e-02 1.3877e-03 proach and Simpson’s rule.
9 4.5984e-03 1.1306e-04
25 5.1401e-04 1.5619e-06
53 1.0956e-04 7.1788e-08
101 2.9628e-05 5.2624e-09
249 4.8175e-06 1.3923e-10
501 1.1852e-06 8.4284e-12
Name Formula
ℎ
𝑓 (𝑥 𝑖 ) + 𝑓 (𝑥 𝑖+1 )
Trapezoidal rule 2
ℎ
𝑓 4 𝑓 (𝑥 𝑖+1 ) + 𝑓 (𝑥 𝑖+2 )
Simpson’s rule 3 (𝑥 𝑖 ) +
3ℎ
𝑓 (𝑥 𝑖 ) + 3 𝑓 (𝑥 𝑖+1 ) + 3 𝑓 (𝑥 𝑖+2 ) + 𝑓 (𝑥 𝑖+3 )
Simpson’s 3/8 rule 8
2ℎ
7 𝑓 (𝑥 𝑖 ) + 32 𝑓 (𝑥 𝑖+1 ) + 12 𝑓 (𝑥 𝑖+2 ) + 32 𝑓 (𝑥 𝑖+3 ) + 7 𝑓 (𝑥 𝑖+4 )
Boole’s rule 45
Newton-Cotes formulas
The trapezoidal approach and Simpson’s rule are closely related method.
In the first the function over an interval is approximated as a straight
line, while in the second a quadratic function is used. We can extend
this sequence of methods to higher order polynomials. In every step an
additional interval is used to define the method and a higher polynomial is
obtained that runs through these points. Following this general approach
we obtain the (closed) Newton-Cotes formulas summarized in Table
10.2. The higher order methods are not often used, especially because
in modern computers it is easy to use a large number of intervals in the
evaluation. Generally, Simpson’s rule gives a good trade-off for accuracy
vs. computational cost.
∫ 1 𝑛
X
𝑓 (𝑥)𝑑𝑥 = 𝑤 𝑖 𝑓 (𝑥 𝑖 ) (10.7)
−1 𝑖=1
𝑏 1
𝑏−𝑎 𝑏−1 𝑎+𝑏
∫ ∫
𝑓 (𝑥)𝑑𝑥 = 𝑓( 𝑥+ )𝑑𝑥 (10.8)
𝑎 2 −1 2 2
𝑛
𝑏−𝑎 X 𝑏−1 1+𝑏
= 𝑤𝑖 𝑓 ( 𝑥𝑖 + ) (10.9)
2 𝑖=1
2 2
It can be shown that the the weights 𝑤 𝑖 are then given by:
2
𝑤𝑖 = (10.10)
(1 − 𝑥 𝑖 )2 [𝑃𝑛0 (𝑥 𝑖 )]2
5 def func(x):
6 return x*[Link](x)
7
12 a = 0.0
13 b = 2.0
14
15 exact = exactIntegral(a, b)
16 estimate = quad(func, a, b)
17
∫ 𝑏 ∫ 𝑑
𝐼= 𝑓 (𝑥, 𝑦)𝑑𝑦 𝑑𝑥 (10.11)
𝑎 𝑐
∫ 𝑑 𝑁𝑦
X
𝐼1 = 𝑓 (𝑥, 𝑦)𝑑𝑦 = ℎ 𝑦 𝑓 (𝑥, 𝑦 𝑗 ) (10.12)
𝑐 𝑗=0
𝑁𝑦
!
∫ 𝑏 X
𝐼 = ℎ𝑦 𝑓 (𝑥, 𝑦 𝑗 ) 𝑑𝑥
𝑎 𝑗=0
𝑁𝑦 ∫ 𝑏
X
= ℎ𝑦 𝑓 (𝑥, 𝑦 𝑗 )𝑑𝑥
𝑗=0 𝑎
𝑁𝑦
𝑁𝑥 X
X
= ℎ𝑥 ℎ𝑦 𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) (10.13)
𝑖=0 𝑗=0
The resulting expression for integral is now the sum of the volume of the
rectangular columns that spanning the 𝑥 𝑦 plane, where the height of the
column is taken to be the function value at the average coordinate in both
the 𝑥 - and the 𝑦 -direction. Similar expressions can be obtained for any of
the composite methods that we have used above in one direction.
10.4 Exercises
Composite methods
∫ 𝑏
𝐼= 𝑓 (𝑥)𝑑𝑥 (10.14)
𝑎
by the summing the area of a number of rectangles that spans the interval
to be integrated. The expression for the midpoint and trapeze methods
are given by
10 Numerical integration 84
𝑁−
X1
𝑥 𝑗 + 𝑥 𝑗+1
𝐼𝑚𝑖𝑑 = ℎ 𝑓 (10.15)
𝑗=1
2
" #
ℎ 𝑁−
X1
𝐼𝑡𝑟 𝑎𝑝 = 𝑓 (𝑥 𝑗 ) + 𝑓 (𝑥 𝑗+1 ) (10.16)
2 𝑗=1
where ℎ is the width of the rectangle. In the midpoint method the function
is evaluated in the middle of the 𝑗 -th rectangle (see Fig. 10.2). In contrast,
in the trapezoid method the average function values is taken of the two
points on either side of the interval (Fig. 10.3).
The two functions to be implemented should take as argument, the
beginning and end of the integration interval, 𝑎 and 𝑏 , the number of
integration points 𝑁 and the function to be integrated 𝑓 𝑢𝑛𝑐 . You will
compare the results provided by code with results given by the routine
[Link] presented in section 10.2.
𝑥
𝑓 (𝑥) = 2𝜋𝑥 2 sin(𝜋𝑥) exp (− ) (10.17)
2𝜋
between 0 and 0.1. Plot the difference in the Δ𝑚𝑖𝑑 = 𝐼𝑚𝑖𝑑 − 𝐼 𝑠𝑝 and
Δ𝑡𝑟 𝑎𝑝 = 𝐼𝑡𝑟 𝑎𝑝 − 𝐼 𝑠𝑝 (where 𝐼 𝑠𝑝 is the value of the integral obtained with
the quadrature of scipy) as a function of 𝑁 . As Δ𝑡𝑟𝑎𝑝 and Δ𝑚𝑖𝑑 can reach
very small values, use a loglog plot to visualize them:
1 [Link](....)
2 [Link]
∫ 𝑇2
Δ𝐻 = 𝐶 𝑝 (𝑇)𝑑𝑇 (10.18)
𝑇1
Simpson’s rule
𝑖+2
ℎ
∫
𝑓 (𝑥)𝑑𝑥 ≈ 𝑓 (𝑥 𝑖 ) + 4 𝑓 (𝑥 𝑖+1 ) + 𝑓 (𝑥 𝑖+2 )
(10.21)
𝑖 3
∫ 𝜋
1
𝑒 𝑥 sin(𝑥) = (1 + 𝑒 𝜋 ) (10.22)
0 2
𝑑 𝑓 (𝑥 𝑖 ) ℎ 2 𝑑2 𝑓 (𝑥 𝑖 )
𝑓 (𝑥 𝑖 + ℎ) = 𝑓 (𝑥 𝑖 ) + ℎ + +... (11.1)
𝑑𝑥 2 𝑑𝑥 2
𝑑 𝑓 (𝑥 𝑖 )
= 𝑓 (𝑥 𝑖 ) + ℎ + 𝑂(ℎ 2 ) (11.2)
𝑑𝑥
𝑑 𝑓 (𝑥 𝑖 ) 𝑓 (𝑥 𝑖 + ℎ) − 𝑓 (𝑥 𝑖 )
= + 𝑂(ℎ 2 ) (11.3)
𝑑𝑥 ℎ
𝑓 (𝑥 𝑖 + ℎ) − 𝑓 (𝑥 𝑖 )
≈ (11.4)
ℎ
𝑑 𝑓𝑖 𝑓𝑖+1 − 𝑓𝑖
= (11.5)
𝑑𝑥 ℎ
𝑑 𝑓 (𝑥 𝑖 )
𝑓 (𝑥 𝑖 − ℎ) = 𝑓 (𝑥 𝑖 ) − ℎ + 𝑂(ℎ 2 ) (11.6)
𝑑𝑥
𝑑 𝑓𝑖 𝑓𝑖 − 𝑓𝑖−1
= (11.7)
𝑑𝑥 ℎ
This backward approximation has the same 𝑂(ℎ) accuracy as the forward
approximation and its use is illustrated in Fig. 11.1b. Now we have two
equations for the derivative of the function in point 𝑥 𝑖 . We can combine
these two equations by adding them up:
𝑑 𝑓 (𝑥 𝑖 ) 𝑓𝑖+1 − 𝑓𝑖−1
= (11.9)
𝑑𝑥 2ℎ
This approximation is called the centered difference approximation and
it can be shown that its accuracy is on the order 𝑂(ℎ 2 ) instead of 𝑂(ℎ)
for the forward and backward approximation. This means that for a
given step size ℎ , the centered difference approximation will give much
11 Numerical differentiation 88
more accurate results than either the backward or the forward forward
approximation. The improved accuracy can be easily imagined when
looking at the illustration of the centered approximation in Fig. 11.2. The
red line through the two points has a slope that is much closer to the
slope of the function itself, as compared to the forward and backward
approximations. This is even true for relatively large values of ℎ .
Example
As a first step, we write a piece of code to calculate and plot the first
derivative of this function. The code below calculates the derivative using
the centered difference approximation.
1 import [Link] as plt
2 import numpy as np
3
16 # function to derive
17 def func(x):
18 return x*[Link](x)
19
25
The central part of this code is the use of the [Link]() function:
1 [Link](y,1)
This functions returns an array identical to 𝑦 , but with all the elements
have been ’rolled’ along the vector. This means that all values are moved
to the right by one position. The last value of the original array will ’roll
over’ to the first position. For example:
𝑦 = [1 , 2 , 3, 4, 5] −→ 𝑦roll = [5 , 1 , 2 , 3 , 4] (11.11)
The trick with the [Link]() function works well for points that
have neighboring points on both sides. The first and last element of the
derivative do not have neighbors on both sides and therefore cannot be
evaluated using the central difference. The only solution is to use the
forward and backward difference approximation for these two points,
which is done in the code above on two separate lines. The comparison
between the numerical and exact values of the derivative for the function
in Eq.11.10 is shown in Fig. 11.3.
3 def fun(x):
4 return x*[Link](x)
5
6 def exactDeriv(x):
7 return [Link](x) + x*[Link](x)
8
9 x0 = 1.0
10 exact = exactDeriv(x0)
11
14 h = 10**(-i)
15
16 fa = (fun(x0+h)-fun(x0))/h
17 ba = (fun(x0)-fun(x0-h))/h
18 ca = (fun(x0+h)-fun(x0-h))/2/h
19
20 print(" h = %1.4e" % h)
21
As can be seen in Fig. 11.4, the error decreases much faster for the centered
approximation than for the backward and forward approximation that
both have a similar accuracy.
In all three approaches used above for the first derivative, we started from
the Taylor expansion. This is a very common approach in the derivation
of numerical methods, especially for derivatives and it can also be used
to derive an expression for the second derivative. Start once again from
the Taylor approximation of a function around 𝑥 𝑖 ± ℎ , but now truncating
it after the term with the second derivatives:
𝑑 𝑓 (𝑥 𝑖 ) ℎ 2 𝑑2 𝑓 (𝑥 𝑖 )
𝑓 (𝑥 𝑖 + ℎ) = 𝑓 (𝑥 𝑖 ) + ℎ + + 𝑂(ℎ 3 ) (11.12)
𝑑𝑥 2 𝑑𝑥 2
𝑑 𝑓 (𝑥 𝑖 ) ℎ 2 𝑑2 𝑓 (𝑥 𝑖 )
𝑓 (𝑥 𝑖 − ℎ) = 𝑓 (𝑥 𝑖 ) − ℎ + − 𝑂(ℎ 3 ) (11.13)
𝑑𝑥 2 𝑑𝑥 2
These two expressions both contain the second derivative, and upon
adding them together the terms with the first derivative disappear:
𝑑2 𝑓 (𝑥 𝑖 )
𝑓 (𝑥 𝑖 + ℎ) + 𝑓 (𝑥 𝑖 − ℎ) = 2 𝑓 (𝑥 𝑖 ) + ℎ 2 + 𝑂(ℎ 4 ) (11.14)
𝑑𝑥 2
obtain:
𝑑 2 𝑓𝑖 𝑓𝑖+1 − 2 𝑓𝑖 + 𝑓𝑖−1
≈ (11.15)
𝑑𝑥 2 ℎ2
Again, the elements 𝑑 2 𝑓 [0] and 𝑑 2 𝑓 [−1] are not evaluated correctly.
To evaluate these elements we have to use the forward and backward
expression that in fact makes the second derivative of the first and last
points equal to their neighbours:
𝑑 2 𝑓𝑖 1
= 2 ( 𝑓𝑖+2 − 2 𝑓𝑖+1 + 𝑓𝑖 ) + 𝑂(ℎ) (11.16)
𝑑𝑥 2 ℎ
𝑑 2 𝑓𝑖 1
= 2 ( 𝑓𝑖 − 2 𝑓𝑖−1 + 𝑓𝑖−2 ) + 𝑂(ℎ) (11.17)
𝑑𝑥 2 ℎ
We have seen how to extract numerical expressions of the first and second
derivatives for one-dimensional functions. Similar techniques can be
used to derive the partial derivatives of a function of two variables 𝑓 (𝑥, 𝑦).
We will adopt the notation:
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = 𝑓𝑖𝑗 (11.18)
For such a function the grid we use has to extend along both the variable
axes. We will take an evenly spaced grid with an increment ℎ in the 𝑥
direction and 𝑘 in the 𝑦 direction. For the first derivatives the central
difference formula directly leads to:
𝑑 1
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = ( 𝑓𝑖+1,𝑗 − 𝑓𝑖−1,𝑗 ) (11.19)
𝑑𝑥 2ℎ
𝑑 1
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = ( 𝑓𝑖,𝑗+1 − 𝑓𝑖,𝑗−1 ) (11.20)
𝑑𝑦 2𝑘
The second partial derivatives are also given by the previously derived
expressions:
11 Numerical differentiation 92
𝑑2 1
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = ( 𝑓𝑖+1,𝑗 − 2 𝑓𝑖,𝑗 + 𝑓𝑖−1,𝑗 ) (11.21)
𝑑𝑥 2 ℎ2
𝑑2 1
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = ( 𝑓𝑖,𝑗+1 − 2 𝑓𝑖,𝑗 + 𝑓𝑖,𝑗−1 ) (11.22)
𝑑𝑦 2 𝑘2
𝑑2 𝑑 𝑑
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = 𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) (11.23)
𝑑𝑥𝑑𝑦 𝑑𝑥 𝑑𝑦
𝑑 1
= ( 𝑓𝑖,𝑗+1 − 𝑓𝑖,𝑗−1 ) (11.24)
𝑑𝑥 2 𝑘
1
𝑓𝑖+1,𝑗+1 − 𝑓𝑖−1,𝑗+1 − 𝑓𝑖+1,𝑗−1 + 𝑓𝑖−1,𝑗−1(11.25)
=
4ℎ 𝑘
11.4 Summary
There are many ways of deriving finite difference expressions for the
derivatives of a given function depending on the method used (forward,
backward, central) and the desired accuracy. In the table below, some
expressions are summarized for future reference.
11 Numerical differentiation 93
𝑑𝑦 𝑖 1 𝑑𝑦 𝑖 1
= (𝑦 𝑖+1 − 𝑦 𝑖 ) = (−𝑦 𝑖+2 + 4 𝑦 𝑖+1 − 3 𝑦 𝑖 )
𝑑𝑥 ℎ 𝑑𝑥 2ℎ
𝑑2 𝑦𝑖 1 𝑑 𝑦𝑖
2
1
= 2 (𝑦 𝑖+2 − 2 𝑦 𝑖+1 + 𝑦 𝑖 ) = 2 (−𝑦 𝑖+3 + 4 𝑦 𝑖+2 − 5 𝑦 𝑖+1 + 2 𝑦 𝑖 )
𝑑𝑥 2 ℎ 𝑑𝑥 2 ℎ
𝑑𝑦 𝑖 1 𝑑𝑦 𝑖 1
= (𝑦 𝑖 − 𝑦 𝑖−1 ) = (3 𝑦 𝑖 − 4 𝑦 𝑖−1 + 𝑦 𝑖−2 )
𝑑𝑥 ℎ 𝑑𝑥 2ℎ
𝑑2 𝑦𝑖 1 𝑑2 𝑦𝑖 1
= 2 (𝑦 𝑖 − 2 𝑦 𝑖−1 + 𝑦 𝑖−2 ) = 2 (2 𝑦 𝑖 − 5 𝑦 𝑖−1 + 4 𝑦 𝑖−2 − 𝑦 𝑖−3 )
𝑑𝑥 2 ℎ 𝑑𝑥 2 ℎ
𝑑 2 𝑓𝑖 𝑓𝑖+1 − 2 𝑓𝑖 + 𝑓𝑖−1
=
𝑑𝑥 2 ℎ2
𝑑2 1
𝑓 (𝑥 𝑖 , 𝑦 𝑗 ) = 𝑓𝑖+1,𝑗+1 − 𝑓𝑖−1,𝑗+1 − 𝑓𝑖+1,𝑗−1 + 𝑓𝑖−1,𝑗−1
𝑑𝑥𝑑𝑦 4ℎ 𝑘
11.5 Exercises
Second derivative
I Write a function that takes as its argument an array of values and
numerically evaluates the second derivative in each point. Make
sure you treat the points at either end of the interval correctly. All
the equations needed are given in section 11.2.
I To test your script use it to calculate the second derivative of the
function 𝑓 (𝑥) = 𝑥𝑠𝑖𝑛(𝑥) in the interval 0 to 4𝜋.
I Compare your result to the analytical solution of the problem and
plot in a single image the function, your solution of the second
derivative and the analytical solution, similar to Fig.11.3.
Non-linear equations 12
In chapter 9 we have discussed how systems of linear equation can be 12.1 The bisection method . . . . 94
solved using techniques from linear algebra, i.e. using matrix operations Recursion . . . . . . . . . . . . 95
such as Gauss-Jordan elimination or the Jacobi algorithm to find the Recursive bisection method 96
12.2 The Newton method . . . . . 97
N solutions to a set of N linear equations. However many problems
Two coupled equations . . . 99
in science and engineering involve non-linear equations, rather than
Example . . . . . . . . . . . . . 101
linear ones. Such nonlinear equations can contain powers of the variables
12.3 Reduced Newton Method . . 103
( 𝑥 3𝑖 ), products of different variables ( 𝑥 𝑖 𝑥 𝑗 ) or even more complicated Multiple equations . . . . . . 103
terms (cos(𝑥 𝑖 ), 𝑒 𝑥 𝑘 ...). Such non-linear equations or systems of multiple 12.4 Exercises . . . . . . . . . . . . . 104
non-linear equations can not be solved in the same straightforward way The bisection method . . . . 104
and different approaches are needed. In this chapter we will discuss two The Newton method . . . . . 104
general methods to find the solutions of non-linear equations numerically; Newton for two equations . 106
the bisection method and Newton’s method. 12.5 Extra exercises . . . . . . . . . 106
Gradient descent in 2D . . . 106
𝑓 (𝑥) = 4 𝑥 3 − 2 𝑥 2 + 4 𝑥 − 4 = 0 (12.1)
Solving this equation requires finding the values of 𝑥 such that 𝑓 (𝑥) = 0.
An intuitive algorithm by which we can find the root of an arbitrary
function is the bisection method. The algorithm locates the root of 𝑓 (𝑥)
in a predefined interval [a, b] where 𝑓 (𝑎) and 𝑓 (𝑏) have opposite signs
(i.e. the value of 𝑥 for which the function value is zero is between 𝑎 and
𝑏 . The method requires that there is only a single root between a and b.
This is illustrated in Fig.12.1.
The first step in the algorithm is to determine the middle of the interval,
between 𝑎 and 𝑏 . By calculation the function value, 𝑓 (𝑐) and comparing
it to 𝑓 (𝑎) and 𝑓 (𝑏) we can determine whether it is on the same side of
the ’root’ as either 𝑎 or 𝑏 . In the case of Fig. 12.1, 𝑓 (𝑐) has the same sign
as 𝑓 (𝑏) and is therefore on the same side of the root (remember that 𝑓 (𝑎)
and 𝑓 (𝑏) are required to have opposite signs). Therefore, we replace 𝑏 by
𝑐 and start from the new values of 𝑎 and 𝑏 .
The algorithm for the bisection method can be summarized as follow:
The method successively divides the search interval in two parts, select the
interval where the root is and repeat the process on the new interval until
the points of the interval are closer together than a predefined convergence
criterion. This algorithm can be implemented in a straightforward way
in a ’while’-loop. However, the structure of the algorithm is very suitable
for the use of a recursive approach that we have seen before when
determining the determinant of a matrix.
Recursion
10 factorial(5)
The code is first called with 𝑛 = 5 that leads to a call of the function
with 𝑛 = 4, which leads to a call for 𝑛 = 3, etc., until finally we obtain
a call with 𝑛 = 1. This calls return a value of 1 that is multiplied by 2
and returned as the return value of the call factorial(2). This value of 2 is
then multiplied by 3 and returned as return value of the factorial(3) and
so on. Finally, we obtain the correct value of factorial(5). The recursive
approach fully omits the need for a for-loop that we have used in Chapter
3 to calculate the factorial and implements the calculation of the factorial
in a single if-else statement.
11 print(’’)
12 print(’-’*20)
13 print(’Iteration %d/%d’ %(n, NMAX))
14
24
𝑓 (𝑥) = 𝑥 2 − 2𝑥 = 0 (12.2)
The starting point for deriving the Newton method is a Taylor expansion
as we have seen before already. Using the Taylor expansion we can
expand the function around a certain starting value 𝑥 0 in terms of the
derivatives of the function in 𝑥 0 :
(𝑥 − 𝑥0 )2 00
𝑓 (𝑥) = 𝑓 (𝑥0 ) + (𝑥 − 𝑥0 ) 𝑓 0(𝑥0 ) + 𝑓 (𝑥 0 ) + 𝑂(𝑥 − 𝑥0 )3 (12.3)
2!
where the last term tells us that the error is on the order of (𝑥 − 𝑥 0 )3 ).
We are interested in finding the value for which 𝑓 (𝑥) = 0 using the first
derivative of the function at 𝑥 0 . Therefore, we can truncate the Taylor
expansion to obtain a first order approximation:
0 = 𝑓 (𝑥 0 ) + (𝑥 − 𝑥 0 ) 𝑓 0(𝑥 0 ) (12.4)
12 Non-linear equations 98
This truncation means that for the time being, we assume that the
function is linear with a slope equal to the derivative of 𝑓 (𝑥). In this
linear approximation we can easily find the value 𝑥 for which it becomes
zero by rearranging it into:
𝑓 (𝑥0 )
𝑥 = 𝑥0 − (12.5)
𝑓 0(𝑥0 )
in Eq.12.5, the value of 𝑥 𝑛+1 is calculated using the the function value
and the derivative at 𝑥 𝑛 : 𝑓 (𝑥 𝑛 ) and 𝑓 0(𝑥 𝑛 ). For a simple function such as
the quadratic one in Eq.12.2 the derivative is easily obtained analytically,
but in general we can insert in Eq.12.5 the numerical value for the first
derivative:
𝑓 (𝑥0 ) 2ℎ
𝑥 = 𝑥0 − = 𝑥0 − 𝑓 (𝑥 0 ) (12.6)
𝑓 (𝑥 0 )
0 𝑓 (𝑥0 + ℎ) − 𝑓 (𝑥0 − ℎ)
In this case we have used the centered approximation for the derivative,
allowing us to obtain an accurate value without the need for analytical
derivatives. The Newton method for a single dimension is implemented
in the code below:
1 import numpy as np
2
11 h = [Link](’h’)
12 return (func(x+h)-func(x-h))/2./h
13
14 # initial guess
15 x0 = 5.0
16
𝑓1 (𝑥1 , 𝑥2 ) = 0 (12.7)
𝑓2 (𝑥1 , 𝑥2 ) = 0 (12.8)
We can again start from a Taylor series to expand each of these functions
(0) (0)
around an initial estimate of the solution indicated by 𝑥 1 and 𝑥 2
12 Non-linear equations 100
Since we are dealing with functions that depend on two variables, the
Taylor expansion includes partial derivative with respect to both of these.
The left-hand side of these equation can be set equal to zero since we are
looking for the root. If we truncate the Taylor series to the first derivative
terms and rearrange the equations we obtain:
(0) (0)
𝛿1 = 𝑥1 − 𝑥1 (12.11)
(0) (0)
𝛿2 = 𝑥2 − 𝑥2 (12.12)
𝜕 𝑓1 𝜕 𝑓1
! ! !
(0) (0)
|
𝜕𝑥 1 𝑥 (0)
|
𝜕𝑥2 𝑥 (0) 𝛿 𝑓
𝜕 𝑓2 𝜕 𝑓2 · 1(0) = − 1(0) (12.15)
|
𝜕𝑥 1 𝑥 (0)
|
𝜕𝑥2 𝑥 (0)
𝛿2 𝑓2
or even
𝕁 · 𝜹 = −f (12.16)
The matrix 𝕁 is the Jacobian matrix or first derivative matrix of the system.
The elements of the Jacobian can either be calculated analytically or
numerically for any well-behaved function. The numerical calculation of
partial first derivatives was discussed in Chapter 11. Once we can set up
the Jacobian, the matrix equation in Eq. 12.15 can be solved using any of
the techniques discussed in Chapter 9 for linear systems. This solution
(0) (0)
provides values for 𝛿 1 and 𝛿 2 that can be added to our initial guess of
the solution to obtain an improved guess of the solution:
12 Non-linear equations 101
Example
𝜕 𝑓1 𝜕 𝑓1
!
3 𝑥 12 2𝑥2
𝜕𝑥1 𝜕𝑥 2
𝕁= 𝜕 𝑓2 𝜕 𝑓2 = (12.21)
2𝑥1 −3𝑥 22
𝜕𝑥1 𝜕𝑥 2
The Newton algorithm is initiated by taking initial guess values for the
solutions as:
(0)
𝑥1 = 1 (12.22)
(0)
𝑥2 =2 (12.23)
!
(0)
𝛿
3 4 5
· 1(0) = − (12.24)
2 −12 𝛿2 − 7
Solving this matrix equation with, for example, the Jacobi iterative
method leads to:
(0)
𝛿1 = −0.7272727 (12.25)
(0)
𝛿2 = −0.70454545 (12.26)
We can now use these new values of the solution to recompute the
(1) (1)
Jacobian 𝕁 and the right-hand side vector f to obtain 𝛿 1 and 𝛿 2 and
then new a solution. The Newton method is implemented in the code
listed below.
1 import numpy as np
2
14 # initial guess
15 x0 = [Link]([1.0, 2.0])
16
17 print(’initial guess’)
18 print("x1 = %f \t x2 = %f" % (x0[0], x0[1]))
19 print(’’)
20 # initial convergence and tolerance
21 eps, tol = 1, 1E-6
22
45
50 else:
51 x0 = x1
52 niter += 1
As we can see the code converges towards the correct solution, i.e. 𝑥 1 = −1
𝑥2 = 1 in 613 iterations. As can be seen by printing the convergence criteria,
the error increases quite dramatically in the first iterations of the method.
This explains why so many iterations are required to converge toward
the solution. This also illustrates that the choice of the initial guess is
very important: if we choose values close to the actual solution we have a
good chance of fast convergence, but sometimes some tweaking of the
12 Non-linear equations 103
To limit the number of step taken during the optimization we can search
along the direction of the Newton step for a point where the error is less
than the error at the starting point. This can be done simply by adding
the following lines of code at line 32 of the code shown above
1 # uptade the solution
2 isearch = 0
3 max_search = 10
4 eps1 = 2*eps
5
8 # new solution
9 x1 = x0 + 0.5**(isearch)*delta
10
15 # increment
16 isearch += 1
17
Multiple equations
𝑓1 (𝑥1 , . . . 𝑥 𝑘 ) = 0 (12.30)
...... (12.31)
𝑓 𝑘 (𝑥1 , . . . 𝑥 𝑘 ) = 0 (12.32)
Using the same approach as above we can collect these 𝑘 equations and
12 Non-linear equations 104
𝜕 𝑓1 𝜕 𝑓1 𝜕 𝑓1
© 𝜕𝑥1 𝜕𝑥 2
... 𝜕𝑥 𝑘 ª 𝛿1 𝑓1
𝜕 𝑓2 𝜕 𝑓2 𝜕 𝑓2 ®
... 𝛿2 ® 𝑓2 ®
© ª © ª
𝜕𝑥1 𝜕𝑥 2 𝜕𝑥 𝑘 ®
.. .. .. .. ® · .. ® = − .. ® (12.33)
® ® ®
. . . . ® .® .®
𝛿 « 𝑓𝑘 ¬
𝜕𝑓 𝜕 𝑓𝑘 𝜕 𝑓𝑘
®
𝑘
... « 𝑘¬
« 𝜕𝑥1 𝜕𝑥 2 𝜕𝑥 𝑘 ¬
The elements of the Jacobian can be either calculated via their analytic
expressions or via the numerical expression of the first derivatives given in
Eq. 11.9. The Newton approach to solving systems of non-linear equations
nicely illustrates how the different methods that we have considered in
earlier chapters come together. For solving systems of multiple equation,
we require both the numerical calculation of derivative and the solution
of a linear system of equations collected in a matrix equation.
12.4 Exercises
In this chapter we have described the bisection method that can be used
to find the root of any function. In this exercise we will implement the bi-
section method. We briefly summarize the method again for convenience.
To find the root of a function 𝑓 (𝑥) in the interval [a, b] the bisection
method works in 3 steps:
𝑓 (𝑥) = 4𝑥 3 − 2 𝑥 2 + 4 𝑥 − 4 (12.34)
In this exercise we will implement the Newton method for the solution of
a non-linear equation and use the method to determine the cause of a gas
tank explosion. The Newton method has been presented in section 12.2
and a pseudo code is given on page 98. You can use this structure to
create your own implementation of the Newton’s method.
Part I : Implementation As a reminder, the Newton method allows to
find the solution of a nonlinear equation, such as 𝑓 (𝑥) = 𝑥 2 − 2 𝑥 = 0, in
an iterative approach. We therefore start with an initial guess 𝑥 0 . Most
12 Non-linear equations 105
likely, this guess does not satisfy 𝑓 (𝑥 0 ) = 0, i.e, it is not the solution we
are looking for. We therefore have to improve the quality of our solution.
To do that the Newton method relies of the iterative procedure
𝑓 (𝑥 𝑛 )
𝑥 𝑛+1 = 𝑥 𝑛 − (12.35)
𝑓 0(𝑥 𝑛 )
𝑥 𝑥
cos(𝑥) exp(− )− =0 (12.36)
2𝜋 𝜋
!
𝑎 𝑉
𝑃+ − 𝑏 − 𝑅𝑇 = 0 (12.37)
𝑉 2
𝑛
𝑛
where
27 𝑅 𝑇𝐶
2 2
1 𝑅𝑇𝐶
𝑎= 𝑏= (12.38)
64 𝑃𝑐 8 𝑃𝑐
12 Non-linear equations 106
Your firm has provided you with the last reading of the tank 𝑇 = 384 K
and P = 4891.3 kPa. The volume of the tank is 𝑉 = 0.15 𝑚 3 . In addition
the critical temperature and pressure of propane are 𝑇𝑐 = 369.9 K and
𝑃𝑐 = 4254.6 kPa. Finally the gas constant is 𝑅 = 0.008314 𝑚 3 𝑘𝑃𝑎/(𝑚𝑜𝑙.𝐾).
Write a Python code called ’[Link]’ and solve the van der Waals equation
to obtain the number of mole 𝑛 in the tank at the time of the explosion.
To solve this equation, you can use the Newton method implemented
above. Check that there is only one solution of this equation by plotting
the left side of the van der Waals equation as a function of 𝑛 .
The correct answer is 𝑛 = 546.51.
Gradient descent in 2D
𝑓1 (𝑥, 𝑦) = 𝑥2 + 4 ∗ 𝑦2 𝑥, 𝑦 ∈ [−(12.42)
3; 3]
𝑥2 𝑦 2
𝑓2 (𝑥, 𝑦) = − sin( − + 3) cos(2 𝑥 + 1 − 𝑒 𝑦 ) 𝑥, 𝑦 ∈ [−3; 1]
(12.43)
2 4
In each case you will start the search for different initial guess of the
solution ( 𝑥 0 ) and report on the minimum found by the gradient de-
scent approach. You will obtain figure similar to the one represented in
Fig. 12.4.
k1 k3
A+B↽ −⇀
−−−−− C + D −−−→ E
k2
𝑑𝐶 𝐴
= −𝑘1 𝐶 𝐴 𝐶 𝐵 + 𝑘2 𝐶 𝐶 𝐶 𝐷 (13.1)
𝑑𝑡
𝑑𝐶 𝐵
= −𝑘1 𝐶 𝐴 𝐶 𝐵 + 𝑘2 𝐶 𝐶 𝐶 𝐷 (13.2)
𝑑𝑡
𝑑𝐶 𝐶
= 𝑘1 𝐶𝐴 𝐶𝐵 − 𝑘2 𝐶𝐶 𝐶𝐷 − 𝑘3 𝐶𝐶 𝐶𝐷 (13.3)
𝑑𝑡
𝑑𝐶 𝐷
= 𝑘1 𝐶𝐴 𝐶𝐵 + 𝑘2 𝐶𝐶 𝐶𝐷 − 𝑘3 𝐶𝐶 𝐶𝐷 (13.4)
𝑑𝑡
𝑑𝐶𝐸
= +𝑘3 𝐶 𝐶 𝐶 𝐷 (13.5)
𝑑𝑡
𝑑𝑦
1st order = 𝑓 (𝑦) (13.6)
𝑑𝑡
𝑑2 𝑦 𝑑𝑦
2nd order +𝑦 = 𝑒 −𝑡 (13.7)
𝑑𝑡 2 𝑑𝑡
𝑑3 𝑦 𝑑2 𝑦 𝑑𝑦
3rd order + 𝑎 +𝑏 +𝑦=0 (13.8)
𝑑𝑡 3 𝑑𝑡 2 𝑑𝑡
The first of these equations is the general form of a first order differential
equation, containing only first derivative terms. Eq. 13.7 contains a
second derivative and is therefore of the 2nd order. In addition, one of the
terms 𝑒 −𝑡 depends explicitly on the independent variable. Any equation
where one term depends on the explicit variable is called non-autonomous.
Similarly, eq. 13.8 contains not only a second derivative but also a third
derivative and is therefore of the third order. Numerical integration of
ODEs is easier and computationally more efficient if the system is in a
canonical form: first-order only and autonomous. Usually, higher-order
and non-autonomous equation can be recast in a canonical form by
substitution, i.e. replacing the variables used by another.
Consider the third-order eq. 13.8. To transform this high-order ODE in a
series of coupled first-order ODEs we can introduce new variables as:
𝑑𝑦
𝑦1 ≡ (13.9)
𝑑𝑡
𝑑𝑦1 𝑑2 𝑦
𝑦2 ≡ = 2 (13.10)
𝑑𝑡 𝑑𝑡
𝑑𝑦2 𝑑3 𝑦
and therefore 𝑑𝑡 ≡ 𝑑𝑡 3
. Using these new variables Eq. 13.8 becomes:
𝑑𝑦
= 𝑦1 (13.11)
𝑑𝑡
𝑑𝑦1
= 𝑦2 (13.12)
𝑑𝑡
𝑑𝑦2
= −𝑎 𝑦2 − 𝑏 𝑦1 − 𝑦 (13.13)
𝑑𝑡
𝑑𝑦
𝑦1 ≡ (13.14)
𝑑𝑡
𝑦2 ≡ 𝑒 −𝑡 (13.15)
𝑑𝑦
= 𝑦1 (13.16)
𝑑𝑡
𝑑𝑦1 𝑑2 𝑦
≡ = −𝑦 𝑦1 + 𝑦2 (13.17)
𝑑𝑡 𝑑𝑡 2
𝑑𝑦2
≡ −𝑒 −𝑡 = −𝑦2 (13.18)
𝑑𝑡
𝑑2 𝑥
= −𝑔 (13.19)
𝑑𝑡 2
(13.20)
𝑑𝑥
= 𝑣 (13.21)
𝑑𝑡
𝑑𝑣
= −𝑔 (13.22)
𝑑𝑡
(13.23)
These two equations can be solved using the same two initial values
mentioned above.
13 Classification of ODEs and linear differential equations 112
k1 k3
A↽ −⇀
−−−−− B ↽ −⇀
−−−−− C
k2 k4
𝑑𝐶 𝐴 (𝑡)
= −𝑘 1 𝐶 𝐴 + 𝑘2 𝐶 𝐵 (13.24)
𝑑𝑡
𝑑𝐶 𝐵 (𝑡)
= 𝑘 1 𝐶 𝐴 + 𝑘4 𝐶 𝐶 − (𝑘2 + 𝑘 3 )𝐶 𝐵 (13.25)
𝑑𝑡
𝑑𝐶 𝐶 (𝑡)
= −𝑘 4 𝐶 𝐶 + 𝑘 3 𝐶 𝐵 (13.26)
𝑑𝑡
(13.27)
𝑑𝐶 𝐴
𝑑𝑡 ª
© 𝑑𝐶 −𝑘 1 𝑘2 0 𝐶𝐴
𝐵 ® = © 𝑘1 −(𝑘2 + 𝑘3 ) 𝑘4 ª® · © 𝐶 𝐵 ª® (13.28)
𝑑𝑡 ®
𝑑𝐶 𝐶 0 𝑘3 −𝑘4 ¬ «𝐶 𝐶 ¬
« 𝑑𝑡 ¬ «
or using a more compact notation
𝑑
C=K·C (13.29)
𝑑𝑡
To obtain an idea on what the solution of this matrix equation may look
like we will first recall the solution of a single first order linear differential
equation: equation of the type
𝑑
𝑦 = 𝑘𝑦 (13.30)
𝑑𝑡
with the initial condition 𝑦(𝑡0 = 0) = 𝑦0 . The solution for this equation
can be derived by separating the variables and integrating both sides of
the equation:
𝑦 𝑡
𝑑𝑦 𝑦
∫ ∫
= 𝑘𝑑𝑡 −→ 𝑙𝑛 = 𝑘𝑡 −→ 𝑦 = 𝑒 𝑘𝑡 𝑦0 (13.31)
𝑦0 𝑦 𝑡0 𝑦0
The solution of the matrix equation 13.29 takes a similar form when we
introduce the exponential of a matrix. The exponential of a matrix is
13 Classification of ODEs and linear differential equations 113
defined in such a way that the solution of the set of differential equation
assembled in the matrix equation 13.29 takes the form:
C = 𝑒 K 𝑡 C0 (13.32)
𝑒 K𝑡 = U𝑒 𝚲𝑡 U−1 (13.33)
𝑒 𝜆1 𝑡 0 0 ... 0
0 𝑒 𝜆2 𝑡 0 ... 0
© ª
®
0 0 𝑒 𝜆2 𝑡 ... 0
®
𝑒 𝚲𝑡 = ® (13.34)
.. .. .. .. .. ®
. . . . .
®
®
« 0 0 0 ... 𝑒 𝜆𝑁 𝑡 ¬
It is important to note here that the element by element exponentiation
only gives the correct matrix exponential for a diagonal matrix. For
non-diagonal matrices this is not the case! The calculation of the matrix
exponential in Eq. 13.33 can easily be done in a Python program by
first determining the eigenvectors and eigenvalues of the matrix, and
subsequently calculating the relevant matrix products. However, as can
be expected the module scipy contains a function [Link]() that
returns the exponential of a matrix that it is given as an argument. These
two approaches to calculate the matrix exponential are compared in the
code below:
1 import numpy as np
2 import [Link] as scln
3 import time
4
5 # create a matrix
6 N = 500
7 A = 0.1*[Link](N,N)
8 t = 0.01
9
26 else:
27 print(’error %1.6e’ %(error))
Not surprisingly, the use of the scipy function results in a much faster
evaluation of the matrix exponential; being almost ten times faster for a
4x4 matrix.
To see how this works in practice we will return to our system of first
order chemical reactions presented above. The solution Eq. 13.32 can be
used to calculate the evolution of the concentrations of 𝐴, 𝐵 and 𝐶 as a
function of time. As a starting point, we need to know the initial value of
each of the concentrations, collected in the vector C0 . We assume that
these concentrations are given by:
For this example the rate constants of all four equations are known:
𝐶 𝐴 (𝑡 𝑝 ) −1 0 0 1
𝐶 𝐵 (𝑡 𝑝 ) ® = exp 1 −2 3 ® 𝑡 𝑝 · 0®
© ª © ª
(13.37)
© ª
«𝐶 𝐶 (𝑡 𝑝 )¬ 2 −3¬ «0¬
0
«
𝑛
C(𝑡 𝑛 = 𝑛Δ𝑡) = 𝑒 K𝑛Δ𝑡 C0 = 𝑒 KΔ𝑡 C0 (13.38)
Therefore, we only need to calculate 𝑒 KΔ𝑡 and all multiples of the time
step only requires and additional multiplication. For a whole series of
time steps between 𝑡0 = 0Δ𝑡 and 𝑡 𝑓 = 𝑁𝑡𝑜𝑡 Δ𝑡 , we can therefore use the
recurrence relationship:
1 import numpy as np
2 import [Link] as plt
3 import [Link] as scln
4
5 # rate constants
6 k1 = 1
7 k2 = 0
8 k3 = 2
9 k4 = 3
10
11 # initial condition
12 C0 = [Link]([1,0,0])
13
14 # evolution time
15 tmax = 5
16 nT = 250
17 T = [Link](0,tmax,nT)
18 dT = T[1]-T[0]
19
31 for i in range(1,nT):
32 C0 = [Link](eKdt,C0)
33 [Link](C0)
34
35 #print C
36 C = [Link](C)
37 [Link](T,C[:,0],linewidth=2,label=’C_A’)
38 [Link](T,C[:,1],linewidth=2,label=’C_B’)
39 [Link](T,C[:,2],linewidth=2,label=’C_C’)
40
41 [Link](’time (min)’,fontsize=12)
42 [Link](’Concentration’,fontsize=12)
43 [Link](loc=1)
44 [Link]()
13.3 Exercises
k1 k2 k4
↽
A −−−→ B −
−−−⇀
−− C −−−→ D
k3
Initially, we only add the starting material 𝐴, implying that the concen-
tration 𝐶 𝐴 = 1.0, while the concentration of the other four species are
zero. The reaction rate constants are given by:
𝑘1 = 2.5min−1
𝑘2 = 0.5min−1
𝑘3 = 0.3min−1
𝑘4 = 0.1min−1
The goal of this exercise is to make a plot of the concentrations of the four
species as a function of time. To achieve this take the following steps:
1. Set up the system of differential equations that describe the reac-
tions.
2. Collect the differential equations in a matrix form.
3. Calculate the exponential of the matrix using the [Link]()
function.
4. Calculate the concentrations of 𝐴, 𝐵, 𝐶 and 𝐷 during the first 10
minutes in 250 steps using a for-loop.
13 Classification of ODEs and linear differential equations 117
𝑑y
= f(𝑡, y) (14.1)
𝑑𝑡
y(𝑡0 ) = y0 (14.2)
As a starting point we first consider the case where only one differ-
ential equation defines the system to illustrate the numerical solution
of non-linear differential equations in general. These approaches can
easily be extended to multiple simultaneous nonlinear ODEs. In the
numerical solution of nonlinear ODEs the goal is to obtain for instance
the concentrations of different reaction species as a function of time. In
such cases, ’time’, 𝑡 is called the independent variable. A first step is the
discretization of the independent variable (here 𝑡 ) in a series of intervals:
𝑡 = [𝑡0 , 𝑡1 , ...𝑡 𝑖 , ...𝑡 𝑁 ]. The solution is then obtained as a series of values
𝑦 = [𝑦0 , 𝑦1 , ...𝑦 𝑖 , ...𝑦 𝑁 ] and if the interval becomes small enough we get a
continuous solution for the system of differential equations. The solution
of the general differential equation above can be obtained by rearranging
the terms in equation 14.1 and integrating on both sides:
∫ 𝑦 𝑖+1 ∫ 𝑡 𝑖+1
𝑑𝑦 = 𝑓 (𝑡, 𝑦)𝑑𝑡 (14.3)
𝑦𝑖 𝑡𝑖
∫ 𝑡 𝑖+1
𝑦 𝑖+1 − 𝑦 𝑖 = 𝑓 (𝑡, 𝑦)𝑑𝑡 (14.4)
𝑡𝑖
∫ 𝑡 𝑖+1
𝑦 𝑖+1 = 𝑦 𝑖 + 𝑓 (𝑡, 𝑦)𝑑𝑡 (14.5)
𝑡𝑖
14 Initial value problems 119
∫ 𝑡 𝑖+1
𝑓 (𝑡, 𝑦)𝑑𝑡 = ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.6)
𝑡𝑖
𝑦 𝑖+1 = 𝑦 𝑖 + ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.7)
𝑡 𝑖+1
ℎ
∫
𝑓 (𝑡, 𝑦)𝑑𝑡 = 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) + 𝑓 (𝑡 𝑖+1 , 𝑦 𝑖+1 )
(14.8)
𝑡𝑖 2
ℎ
𝑦 𝑖+1 = 𝑦 𝑖 + 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) + 𝑓 (𝑡 𝑖+1 , 𝑦 𝑖+1 )
(14.9)
2
We can then use this prediction to obtain a corrected value for 𝑦 𝑖+1 :
ℎ
(𝑦 𝑖+1 )𝐶𝑜𝑟 = 𝑦 𝑖 + 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) + 𝑓 𝑡 𝑖+1 , (𝑦 𝑖+1 )𝑃𝑟
(14.11)
2
ℎ ℎ
𝑦 𝑖+1 = 𝑦𝑖 + 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) + 𝑓 (𝑡 𝑖+1 , 𝑦 𝑖+1 ) (14.12)
2 2
= 𝑦𝑖 + 𝑤1 𝑘1 + 𝑤2 𝑘2 (14.13)
1
𝑤1 = 𝑘1 = ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.14)
2
1
𝑤2 = 𝑘2 = ℎ 𝑓 (𝑡 𝑖 + ℎ, 𝑦 𝑖 + 𝑘 1 ) (14.15)
2
𝑚
X
𝑦 𝑖+1 = 𝑦 𝑖 + 𝑤𝑖 𝑘𝑖 (14.16)
𝑖=1
𝑘1 = ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.17)
𝑘2 = ℎ 𝑓 (𝑡 𝑖 + 𝑐 2 ℎ, 𝑦 𝑖 + 𝑎21 𝑘1 ) (14.18)
𝑘3 = ℎ 𝑓 (𝑡 𝑖 + 𝑐 3 ℎ, 𝑦 𝑖 + 𝑎31 𝑘1 + 𝑎32 𝑘2 ) (14.19)
... (14.20)
𝑘𝑚 = ℎ 𝑓 (𝑡 𝑖 + 𝑐 𝑚 ℎ, 𝑦 𝑖 + 𝑎 𝑚 1 𝑘1 + 𝑎 𝑚 2 𝑘2 + ... + 𝑎 𝑚,𝑚−1 𝑘 𝑚−1()14.21)
𝑝−1
X
𝑘 𝑝 = ℎ 𝑓 (𝑡 𝑖 + 𝑐 𝑝 ℎ, 𝑦 𝑖 + 𝑎 𝑝𝑙 𝑘 𝑙 ) (14.22)
𝑙=1
1
𝑦 𝑖+1 = 𝑦 𝑖 + (𝑘1 + 2 𝑘2 + 2 𝑘3 + 𝑘4 ) (14.23)
6
with:
14 Initial value problems 122
𝑘1 = ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.24)
ℎ 𝑘1
𝑘2 = ℎ 𝑓 𝑡𝑖 + , 𝑦𝑖 + (14.25)
2 2
ℎ 𝑘2
𝑘3 = ℎ 𝑓 𝑡𝑖 + , 𝑦𝑖 + (14.26)
2 2
𝑘4 ℎ 𝑓 𝑡 𝑖 + ℎ, 𝑦 𝑖 + 𝑘3
= (14.27)
This method is probably the most used technique in practice to solve ODEs
and it can be programmed in a relatively simple way. The implementation
of the fourth-order RK method is the subject of one of the exercises in
this chapter.
k
A −−−→ B
𝑑 𝑑
𝐶 𝐴 = −𝑘𝐶 𝐴 𝐶 𝐵 = 𝑘𝐶 𝐴 (14.28)
𝑑𝑡 𝑑𝑡
Using the initial conditions 𝐶 𝐴 (0) = 1 and 𝐶 𝐵 (0) = 0, the exact solutions
of these ODEs can be readily derived
12 # 4th order RK
13 def RK4(func, yi, ti, dt):
14 k1 = dt*func(yi,ti)
15 k2 = dt*func(yi+0.5*k1,ti+0.5*dt)
16 k3 = dt*func(yi+0.5*k2,ti+0.5*dt)
17 k4 = dt*func(yi+k3,ti+dt)
14 Initial value problems 123
20 # variation of CA
21 def df(ca,t):
22 return -k*ca
5 # variation of CA
6 def df(t,y,k=0.2):
7 return -k*y # return the value
8
9 # time argument
10 t_start = 0.0 # min
11 t_end = 30.0 # min
12
13 tmax,nT = 30.0,500
14 T = [Link](0,tmax,nT)
15
16 # initial condition
17 C0, T0 = 1.0, 0.0
18
19 # storage
20
25 [Link](result.t, result.y[0],linewidth=2)
26 [Link](’time (min)’,fontsize=12)
27 [Link](’Concentration’,fontsize=12)
28 [Link]()
14 Initial value problems 124
Simultaneous ODEs
We have so far only discussed the different integration methods for ODEs
by applying them to one equation at a time. The same methods can be
extended in a straightforward way to deal with systems coupled ODEs.
To illustrate this, we consider the following set of 𝑛 ODEs
𝑑𝑦1
= 𝑓1 (𝑡, 𝑦1 , 𝑦2 , ...𝑦𝑛 ) (14.30)
𝑑𝑡
𝑑𝑦2
= 𝑓2 (𝑡, 𝑦1 , 𝑦2 , ...𝑦𝑛 ) (14.31)
𝑑𝑡
... (14.32)
𝑑𝑦𝑛
= 𝑓𝑛 (𝑡, 𝑦1 , 𝑦2 , ...𝑦𝑛 ) (14.33)
𝑑𝑡
Any other methods presented above for a single ODE can be used to
solve a set of simultaneous ODEs and we will explore this in one of the
exercises at the end of this chapter.
𝑑𝐴
= 𝛼𝐴𝐵 − 𝛽𝐴 (14.35)
𝑑𝑡
where 𝛼 is the growth rate due to prey consumption and 𝛽 is the death
rate by overpopulation. The population change of prey-animals is given
by
𝑑𝐵
= 𝛾𝐵 − 𝜇𝐴𝐵 (14.36)
𝑑𝑡
where 𝛾 is the birth rate and 𝜇 the death rate due to the predator. We can
use any of the methods presented above to study the dynamics of these
populations. The code below shows how to do this using the fourth-order
14 Initial value problems 125
5 # 4th order RK
6 def RK4(func, yi, ti, dt):
7 k1 = dt*func(ti,yi)
8 k2 = dt*func(ti+0.5*dt,yi+0.5*k1)
9 k3 = dt*func(ti+0.5*dt,yi+0.5*k2)
10 k4 = dt*func(ti+dt,yi+k3)
11 return yi + 1./6*( k1+2*k2+2*k3+k4 )
12
13 # population change
14 def df(t,Y):
15 Yp = [Link](2)
16 Yp[0] = alpha*Y[0]*Y[1] - beta*Y[0]
17 Yp[1] = gamma*Y[1] - mu*Y[0]*Y[1]
18 return Yp
19
20 # time
21 tmax, nT = 10.0,500
22 T = [Link](0,tmax,nT)
23 dT = T[1]-T[0]
24 # parameters
25 alpha, beta, gamma, mu = 1.5,2.0,2.0,1.0
26 # initial population
27 A0, B0 = 2.0, 0.2
28 T0 = 0.0
29
30 # storage
31 POP = [Link]((nT,2))
32 POP[0,:] = [A0,B0]
33
38
39
40 # scipy approach
41 t_span = [Link]([0., tmax])
42
14.6 Exercises
Forward Euler
In this first exercise you will implement the forward Euler method for
the solution of ODE’s. As a reminder, the forward Euler method allows
finding the solution of the differential equation
𝑑𝑦
= 𝑓 (𝑡, 𝑦) (14.37)
𝑑𝑡
𝑦 𝑖+1 = 𝑦 𝑖 + ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 ) (14.38)
To implement this approach you will write a function, for example called
forwardEuler() that returns the value of 𝑦 𝑖+1 on taking the value of 𝑦 𝑖 as an
argument. Subsequently, this function can be called in a for loop running
over all the time steps desired. To test this program, calculate the solution
of the simple reaction
𝑘
𝐴→
− 𝐵 (14.39)
with the initial condition 𝐶 𝐴 (0) = 1 and 𝐶 𝐵 (0) = 0 (where 𝐶 𝑋 (𝑡) is the
concentration of 𝑋 at time 𝑡 ) and a value of 𝑘 = 0.2 min−1 . Simulate the
dynamics of this reaction over a period of 30 min.
The differential equations for the evolution of 𝐶 𝐴 and 𝐶 𝐵 can be easily
derived and solved analytically. You can therefore test your program by
comparing the solution it provides with the exact solution that you will
have derived analytically. Vary the number of time step you use in the
simulations of the reaction. Use between 10 and 500 and comment on the
accuracy of the method.
14 Initial value problems 127
1
𝑦 𝑖+1 = 𝑦 𝑖 + (𝑘1 + 2 𝑘2 + 2 𝑘3 + 𝑘4 ) (14.40)
6
with:
𝑘1 =
ℎ 𝑓 (𝑡 𝑖 , 𝑦 𝑖 )
𝑘2 =
ℎ 𝑘1
ℎ 𝑓 𝑡 𝑖 + , 𝑦𝑖 +
2 2
(14.41)
𝑘3 =
ℎ 𝑘2
ℎ 𝑓 𝑡 𝑖 + , 𝑦𝑖 +
2 2
𝑘4 =
ℎ 𝑓 𝑡 𝑖 + ℎ, 𝑦 𝑖 + 𝑘3
k1 k2 k4
↽
A −−−→ B −
−−−⇀
−− C −−−→ D
k3
𝑘1 = 2.5min−1
𝑘2 = 0.5min−1
𝑘3 = 0.3min−1
𝑘4 = 0.1min−1
14 Initial value problems 128
Just like in section 13.3, the goal of this exercise is to make a plot of the
concentrations of the four species as a function of time. To achieve this
take the following steps:
1. Write your Runge-Kutta routine or copy it from the exercise above.
2. Write a function that returns a vector containing the function values
calculated for the right-hand side of the differential equations (one
value for each differential equation)
3. Calculate the concentrations of 𝐴, 𝐵, 𝐶 and 𝐷 during the first 10
minutes in 250 steps using a for-loop.
4. Plot the concentrations as a function of time.
5. Compare the resulting graph with the one obtained from the
exponentiation solution in Chapter 3. It should be exactly the same!
You can again play around with the end-time to see if eventually all of
𝐴 is transformed into 𝐷 . Note: in this particular case we are dealing
with linear differential equations, which is particularly easy to solve
and the Runge-Kutta method will give an exact result, regardless of the
number of steps taken in the integration. This is in general not the case
for non-linear differential equations.
Epidemic propagation
You have been hired by the World Health Organization to estimate the
impact of a disease outbreak on a city of 50.000 people. The WHO has
developed a model, called the 𝑆𝐼𝑅 model to simulate the propagation of
the infection. Let 𝑆 be the number of healthy people that can contract the
virus, (𝑆 stands for susceptible here) 𝐼 is the number of infected people
and 𝑅 the number of people who have recovered from the infection and
are immune to it. The evolution of the populations 𝑆 , 𝐼 and 𝑅 is then
given by
𝑑𝑆
= −𝛾 · 𝑆 · 𝐼 (14.42)
𝑑𝑡
𝑑𝐼
= 𝛾·𝑆·𝐼−𝛼·𝐼 (14.43)
𝑑𝑡
𝑑𝑅
= 𝛼·𝐼 (14.44)
𝑑𝑡
where 𝛾 is the transmission rate between sick and healthy people and 𝛼 the
recovery rate. For the present case the initial conditions are 𝑆(0) = 50000,
𝐼(0) = 2 and 𝑅(0) = 0. The transmission rate (𝛾) has been estimated at
10−4 day−1 while the recovery rate is 0.5 day−1 . Simulate the dynamics
of the three different populations over a period of 30 days. Report the
maximum number of infected people that can be expected and when
this will happen.
To limit the effect of the outbreak, the WHO has proposed to quarantine
the population. They estimate that it will decreases the transmission rate
to 0.5 10−4 day−1 but will slow down recovery to 𝛼 = 0.1 day−1 . Before
implementing their plan they asked you to forecast the impact of the
14 Initial value problems 129
𝑑
𝑦1 = 𝑓1 (𝑥, 𝑦1 , 𝑦2 ) (15.1)
𝑑𝑥
𝑑
𝑦2 = 𝑓2 (𝑥, 𝑦1 , 𝑦2 )
𝑑𝑥
𝑦1 (𝑥 = 𝑥0 ) = 𝑦1 , 0 (15.2)
𝑦2 (𝑥 = 𝑥 𝑓 ) = 𝑦2 , 𝑓
This problem is graphically shown in Fig. 15.1 where for one of the
functions the initial values is known, while for the other the y-value is
known at the end of the interval. Boundary values problem occur in
many branches of chemical engineering. Several methods of different
complexity have been developed to solve such problems. In this chap-
ter, we will consider two of these: the shooting method and the finite
difference methods.
𝑦2 (𝑥0 ) = 𝛾 (15.3)
We can then solve the equations 15.1 as an initial value problem with the
initial condition 𝑦1 (𝑥 0 ) = 𝑦1,0 and 𝑦2 (𝑥 0 ) = 𝛾 . This can be done using any
of the methods presented in previous chapters. Since the value of 𝑦2 (𝑥 0 )
was a guess, it will not generally be correct and therefore it misses the
target at 𝑥 = 𝑥 𝑓 and the solution does not satisfy the boundary condition
𝑦2 (𝑥 = 𝑥 𝑓 ) = 𝑦2, 𝑓 . Therefore, we have to optimize our guess-value until
we obtain the correct boundary condition for 𝑦2 (𝑥 𝑓 ). To optimize the
value of 𝛾 in a systematic way we define a so-called cost-function:
The cost function quantifies how far the solution obtained with the
guessed value of 𝑦2 (𝑥 0 ) is from the boundary condition. The cost function
can be then be expanded in a Taylor series:
𝜕Φ
Φ(𝛾 + Δ𝛾) = Φ(𝛾) + Δ𝛾 + 𝑂(Δ𝛾 2 ) (15.5)
𝜕𝛾
𝜕Φ −Φ(𝛾)
0 = Φ(𝛾) + Δ𝛾 −→ Δ𝛾 = (15.6)
𝜕𝛾 𝜕Φ
𝜕𝛾
𝜕Φ 𝜕𝑦2 (𝑥 𝑓 ,𝛾)
Since 𝜕𝛾
= 𝜕𝛾
we finally obtain
𝜕𝑦2 (𝑥 𝑓 , 𝛾) −1
Δ𝛾 = − 𝑦2 (𝑥 𝑓 , 𝛾) − 𝑦2, 𝑓 (15.7)
𝜕𝛾
15 Boundary-value problems 132
𝜕𝑦2 (𝑥 𝑓 , 𝛾) 𝑦2 (𝑥 𝑓 , 𝛾𝑛+1 ) − 𝑦2 (𝑥 𝑓 , 𝛾𝑛 )
= (15.9)
𝜕𝛾 𝛾𝑛+1 − 𝛾𝑛
Application
𝑑2 𝑦
= 4(𝑦 − 𝑥) (15.10)
𝑑𝑥 2
𝑑𝑦1
= 𝑦2 (15.11)
𝑑𝑥
𝑑𝑦2
= 4(𝑦1 − 𝑥) (15.12)
𝑑𝑥
We now have two first order differential equation but only for the first
equation the initial value is known ( 𝑦1 (0) = 0). For the second equation
we do not have an initial value, but we do know the final value for the
first equation, 𝑦1 (1) = 2. Therefore, we need to find a suitable initial
value for 𝑦2 that satisfies the extra boundary condition, 𝑦1 (1) = 2. As
a starting point we make a guess and say: 𝑦2 (0) = 𝛾0 = 1.0. Using this
initial value, we can now solve the system of equations above as an initial
value problem (for example with 𝑠𝑐𝑖𝑝 𝑦.𝑖𝑛𝑡𝑒 𝑔𝑟 𝑎𝑡𝑒.𝑠𝑜𝑙𝑣 𝑖 𝑣𝑝 ) and check
whether the boundary condition 𝑦1 (1) = 2 is respected or not. As seen
15 Boundary-value problems 133
in Fig. 15.2, our initial guess misses the target and hits 𝑦1 (1 , 𝛾0 ) = 1
instead.
We therefore have to improve our guess. However, at this stage we only
know 1 value 𝛾0 and we cannot improve the value of 𝛾 in a systematic
way using Eq. 15.9. In order to evaluate the partial derivative in Eq. 15.7
we need an additional value and we therefore make a second guess,
𝑦2 (0) = 𝛾1 = 0.0. As shown in Fig. 15.2 the solution obtained with
this initial condition also misses the target (even more) and leads to
𝑦1 (1, 𝛾1 ) = −0.813. However we now have two values for and we can
use both these values to estimate how to update our value of 𝛾 using Eq.
15.7 and write
𝑦1 (𝑥 𝑓 , 𝛾1 ) − 𝑦1 (𝑥 𝑓 , 𝛾0 ) −1
Δ𝛾 = − 𝑦1 (𝑥 𝑓 , 𝛾1 ) − 𝑦1, 𝑓 (15.13)
𝛾1 − 𝛾0
𝑁 equations
𝑑𝑦 𝑖
= 𝑓 (𝑥, 𝑦1 , 𝑦2 , ...𝑦 𝑁 ) 𝑥 ∈ [𝑥 , 𝑥 𝑓 ] (15.14)
𝑑𝑥
15 Boundary-value problems 134
The boundary conditions are split between initial conditions for 𝑟 equa-
tions and at least 𝑁 − 𝑟 final conditions:
Following the principle we have outlined above for the shooting method
we will make an initial guess for the missing initial conditions:
(0)
𝑦 𝑗 (𝑥0 ) = 𝛾 𝑗 𝑗 = 𝑟 + 1 , ....𝑁 (15.17)
𝑦𝑟+1 (𝑥 𝑓 , 𝛾𝑟+1 )
−1 © 𝑦𝑟+1 (𝑥 𝑓 , 𝛾𝑟+2 )ª®
Δ𝜸 = J(𝑥 𝑓 , 𝜸)
· (15.19)
...
®
®
« 𝑦 𝑁 (𝑥 𝑓 , 𝛾 𝑁 ) ¬
𝑑2 𝑦 𝑑𝑦
+ 𝛼(𝑥) + 𝑦𝛽(𝑥) = 𝑓 (𝑥) (15.21)
𝑑𝑥 2 𝑑𝑥
𝑑𝑦 𝑦 𝑖+1 − 𝑦 𝑖−1
= (15.22)
𝑑𝑥 2ℎ
𝑑2 𝑦 𝑦 𝑖+1 − 2 𝑦 𝑖 + 𝑦 𝑖−1
= (15.23)
𝑑𝑥 2 ℎ2
ℎ
(𝑦 𝑖+1 − 2 𝑦 𝑖 + 𝑦 𝑖−1 ) + 𝛼(𝑥 𝑖 )(𝑦 𝑖+1 − 𝑦 𝑖−1 ) + ℎ 2 𝛽(𝑥 𝑖 )𝑦 𝑖 = ℎ 2 𝑓 (𝑥 𝑖 ) (15.24)
2
By writing this equation for each node and accounting for the boundary
conditions we obtain a system of linear equations given by:
𝑦0 = 𝑦 𝑎 (15.25)
ℎ ℎ
𝑦 𝑖+1 (1 + 𝛼(𝑥 𝑖 )) + 𝑦 𝑖 (ℎ 2 𝛽(𝑥 𝑖 ) − 2) + 𝑦 𝑖−1 (1 − 𝛼(𝑥 𝑖 )) = ℎ 2 (15.26)
𝑓 (𝑥 𝑖 )
2 2
𝑦 𝑁 = 𝑦𝑏 (15.27)
𝑎 00 𝑎01 ... 𝑎 0𝑁 𝑦0 𝑏0
𝑎10 𝑎11 ... 𝑎 1 𝑁 ® 𝑦1 ® 𝑏 1 ®
© ª © ª © ª
. .. .. .. ®® · .. ®® = .. ®® (15.28)
. .
. . . ® . ® . ®
«𝑎𝑁 0 𝑎𝑁1 ... 𝑎 𝑁 𝑁 ¬ « 𝑦 𝑁 ¬ «𝑏 𝑁 ¬
or more compactly
A·Y=B (15.29)
As an example, the code below uses the finite difference method to solve
the boundary value problem:
𝑑2 𝑦
+𝑦=0 (15.30)
𝑑𝑥 2
5 # Number of nodes
6 N = 10
7 X = [Link](0,[Link]/2.,N+1)
8 h = X[1]-X[0]
9
28 Y = solve(A,b)
29
30 [Link](X,Y,marker=’o’,color=’black’)
31 [Link](X,[Link](X)+[Link](X),linewidth=2,color=’blue’)
32 [Link](’X’)
33 [Link](’Y’)
34 [Link]()
15 Boundary-value problems 137
15.3 Exercises
In this first exercise you will implement the shooting method for the
solution of boundary value problems as described in section 15.1. As a
reminder, the shooting method treats a boundary value problem as an
initial value problem. We therefore make a guess of the initial values that
are not provided, solve the equations and compare the solution obtained
with the boundary conditions provided. This allows us to make a better
guess for the missing initial condition. We can therefore converge toward
the correct solution by iterating the process. A detailed description of
the method can be found in section 15.1 and 15.1.
To test your implementation you will compute the solution of the bound-
ary value problem given by Eq. 15.10, i.e.
𝑑2 𝑦
= 4(𝑦 − 𝑥) (15.31)
𝑑𝑥 2
𝑑𝑦1
= 𝑦2 (15.32)
𝑑𝑥
𝑑𝑦2
= 4(𝑦1 − 𝑥) (15.33)
𝑑𝑥
We now know the initial value for the first equation ( 𝑦1 (0) = 0) but
not for the second. We therefore have to guess it and improve our
guess following the procedure introduced in section 15.1. As shown in
15 Boundary-value problems 138
section 15.1, the shooting method should converge rapidly toward the
solution as represented in Fig. 15.2.
HINT: To help you in the implementation we provide below a possible
pseudo code for the shooting method
1 from [Link] import solve_ivp
2
10 ....
11
12 # boundary condition
13 y1_0,y1_f = 0.0, 2.0
14
1D transport
𝑑2 𝑣 Δ𝑃
𝜇 = ≡𝐺 (15.34)
𝑑𝑦 2 Δ𝑥
with the boundary conditions 𝑣(𝑦 = 0) = 0 and 𝑣(𝑦 = 𝐵) = 𝑉𝑢𝑝 . The goal
of this exercise is to compute the profile 𝑣(𝑦) using the expression for the
numerical evaluation of the second derivative and the techniques used
to solve linear systems of equations. Using the finite difference approach
we can rewrite the eq. 15.34 as
𝑣 𝑖+1 − 2𝑣 𝑖 + 𝑣 𝑖−1
𝜇 =𝐺 (15.35)
ℎ2
Note that G You can therefore rewrite the problem as a system of linear
equations
1 −2 1 ... ... 𝑎 0𝑁 𝑣0 𝑠0
0 1 −2 1 ... 𝑎 1𝑁 ® 𝑣 1 ® 𝑠 1 ®
© ª © ª © ª
. .. .. ®® .. ®® .. ®®
. .
. . ® . ® . ®
. .. .. ®® · .. ®® = .. ®® (15.36)
. .
. . ® . ® . ®
. .. .. ®® .. ®® .. ®®
. .
. . ® . ® . ®
«𝑎𝑁 0 𝑎𝑁2 ... 𝑎 𝑁 𝑁 ¬ «𝑣 𝑁 ¬ « 𝑠 𝑁 ¬
In this equation 𝑣 𝑖 is the velocity profile at the 𝑖 -th point. You must
determine the expression for the 𝑎 𝑖𝑖 and 𝑠 𝑖 . Once you have established
this system of linear equations, use one of the techniques presented in
chapter 9 (e.g. [Link]) to solve it and obtain the values of the
𝑣 𝑖 . Compare your numerical results to the analytic expression given by
𝑦 1 Δ𝑃 2
𝑣(𝑦) = 𝑉𝑢𝑝 + (𝑦 − 𝑦𝐵) (15.37)
𝐵 2𝜇 Δ𝑥
During the calculation you will take : 𝐵 = 5 · 10−3 m, 𝑉𝑢𝑝 = (1/6) · 10−4
m/s, 𝜇 = 10−3 . You will vary the dynamic pressure gradient from 0.0
to -0.06 Pa/m (take 5 points linearly spaced between these two values ).
Take about 50 to 100 points along the 𝑦 axis. You will then compare your
numerical results to the analytic solution. You should find variations of
the velocity profile similar to the ones represented in Fig. 15.4b.
The functions returns one random number in the open interval [0 , 1[. The
generation of random numbers is done by deterministic algorithm and
therefore they are not really random. However, good random number
generators lack any patterns in the series of numbers generated and can,
for all practical purposes be considered as truly random. In fact, it is
sometimes useful when coding to obtain the exact same sequence of
random numbers when running a code multiple times, especially for
debugging purposes. This reproducibility of the random numbers is the
result of the algorithms used, which generally produce numbers based
on a certain ’seed-number’. By default, this seed value is taken for the
system-time on the computer that runs the code but it is also possible to
imposing a given seed to the generator with the function:
1 r = [Link](13)
16 Monte Carlo schemes 142
The seed (here equals to 13) is an integer that entirely determines the full
sequence of random number generated. For example in the following
code:
1 [Link](13) # seed the generator
2 r1 = [Link](10) # creates 10 random numbers
3 r2 = [Link](10) # creates another 10 random numbers
4
1 𝑥2
𝑝(𝑥) = √ 𝑒 − 2 (16.1)
2𝜋
1 X = [Link](100)
4 N = 1E4
5 XYu = [Link](N,2)
6 XYn = [Link](N,2)
7
8 [Link](XYu[:,0],XYu[:,1],’o’,color=’#007FFF’)
9 [Link](XYu[:,0],20,color=’#007FFF’)
10
11 [Link](XYn[:,0],XYn[:,1],’o’,color=’#007FFF’)
12 [Link](XYn[:,0],20,color=’#007FFF’)
This code generates two arrays (each with two columns representing the
𝑥 and 𝑦 values) containing the position of 𝑁 = 1𝐸4 points that are either
distributed either uniformly or according to a normal distribution in the
𝑥, 𝑦 plane. These points are visualized both by plotting them directly
and by computing the histogram of their 𝑥 coordinate. The results of this
simple code is show in Fig. 16.1
By default [Link]() generates random numbers between
0 and 1, however, other distributions can sometimes be useful. As an
example the command:
1 X = [Link](min,max,N)
1 (𝑥−𝜇)2
−
𝑝(𝑥) = √ 𝑒 2𝜎2 (16.2)
2𝜋𝜎 2
The random number generators above all generate floating point numbers.
For some purposes it is of interest to generate a random sequence of
integers, for which the function [Link].random_integers() is
available:
1 X = [Link].random_integers(low,high,N)
also creates 𝑁 random integers but in the interval [𝑙𝑜𝑤, ℎ𝑖 𝑔 ℎ[ (i.e. the
high value is never generated). In addition to the random generators
contained in 𝑛𝑢𝑚𝑝 𝑦 the module 𝑠𝑐𝑖𝑝 𝑦.𝑠𝑡 𝑎𝑡𝑠 gives access to a large
number of distributions that may be useful for specific physical systems.
For example the code:
1 from [Link] import maxwell
2 r = [Link](size=1000)
One of the Monte Carlo schemes that we will consider is the use of random
numbers in the calculation of integrals. We have seen that integrals over
a certain integration interval can be calculated by dividing the interval in
16 Monte Carlo schemes 144
Standard MC integration
∫ 𝑏
𝐼= 𝑓 (𝑥)𝑑𝑥 (16.3)
𝑎
𝑁
𝑏−𝑎X
𝐼= 𝑓 (𝑥 𝑖 ) (16.4)
𝑁 𝑖=1
where the points 𝑥 𝑖 are chosen randomly instead of using a regular grid.
This is graphically shown in Fig.16.2. This approach to integration may
seem somewhat surprising, but in fact it turns out to be a very efficient
way of evaluation high-dimensional integrals.
To illustrate the Monte-Carlo integration method we will calculate the
∫2
integral 1 (2 + 3 𝑥)𝑑𝑥 , of which the exact value is 6.5. The code below
evaluates the integral using the Monte Carlo method:
1 import numpy as np
2
3 def f(x):
4 return 2.+3.*x
5
6 def MCInt(a,b,f,N):
7 X = [Link](a,b,N)
8 I = (float(b-a)/N) * [Link](f(X))
9 return I
10
11 a,b,N = 1.,2.,1000
12 I = MCInt(a,b,f,N)
13 print(I)
By executing this code we find a value of ≈ 6.5 ± 0.05 with only random
100 points. Of course the exact value changes each time we execute the
program as the sampling points are chosen randomly. We can add more
points to obtain a better estimate but for simple 1-dimensional integrals
the Monte-Carlo integration converges slowly and 1E6 points are required
to obtain an error below 10−3 . Therefore, for simple integrals the methods
16 Monte Carlo schemes 145
Throwing darts
𝑏
𝑁0
∫
𝑦= 𝑓 (𝑥)𝑑𝑥 ' × [𝑦0 (𝑏 − 𝑎)] (16.5)
𝑎 𝑁
In the Python code below this method is used to calculate the integral:
4
𝑥
∫
𝐼= 4 + sin(2 𝑥) exp( ) (16.6)
0 2
1 import numpy as np
2
3 def f(x):
4 return 4 + [Link](2*x)*[Link](0.5*x)
5
6 def MCInt(a,b,y0,f,N):
7 X = [Link](a,b,N)
8 Y = [Link](0,y0,N)
9 N0 = Y[Y<=f(X)].size
10 return float(N0)/N*y0*(b-a)
11
12 a,b,y0,N = 0.,4.,15.,10000
13 I = MCInt(a,b,y0,f,N)
14 print(I)
For 100 points this approach leads to 𝐼 ' 17.9 where the quadrature
method implemented in 𝑠𝑐𝑖𝑝 𝑦 leads to 𝐼 = 17.8365. The same method is
used in one of the exercises at the end of this chapter to determine the
value of 𝜋.
16 Monte Carlo schemes 146
Random Walk in 1D
4 # number of walker
5 M = 1000
6
10 # positions
11 X = [Link]((NT,M))
12
23 [Link]([Link](NT),[Link]([Link](X**2,1)),color=’blue’,linewidth=2)
24 [Link]([Link](NT),[Link](X**1,1),color=’black’,linewidth=2)
25 [Link]([Link](NT),[Link]([Link](0.,NT,NT))*dX,color=’#007FFF’,
linewidth=2)
26 [Link]([Link](NT),[Link](M),color=’#5F5F6F’,linewidth=2)
27 [Link]()
As we can see on the first panel of Fig. XX, every walker take a very
different trajectory than than the other ones. Some go left (negative
position), some right (positive values). However the mean value of the
16 Monte Carlo schemes 147
position (< 𝑋 >) is null over time since the probability to go left equals
the one to go right. However it can be shown that the average of the
square distance, i.e. < 𝑋 2 > increase linearly with the number of steps
𝑁 . Taking the square root of < 𝑋 2 > we obtain the root-mean-squared
distance which represent the average positive distance away from the 0.
We easily see that this distance increases with the square-root of 𝑁 :
√ √
< 𝑋 2 >(𝑁) = 𝑁Δ𝑋 (16.7)
As we can see our numerical simulations here with 1000 walkers gives
agrees very well with this theoretical results. Running the simulation with
even more walkers would lead to an even better agreement. Finally we
can note from the last panel that the distribution of final positions follows
a Gaussian distribution centered around 0. That also agree perfectly with
analytical analysis that are however beyond the scope of this chapter.
𝜕 𝜕2
𝑃(𝑥, 𝑡) = 𝐷 2 𝑃(𝑥, 𝑡) (16.11)
𝜕𝑡 𝜕𝑥
√
𝑥 𝑛 = 𝑥 𝑛−1 + Δ𝑡𝜇𝑥 𝑛−1 + 𝜎𝑥 𝑛−1 𝑟𝑛−1 Δ𝑡 (16.12)
where 𝑥 𝑛 is the stock price at time 𝑡 𝑛 , Δ𝑡 the time interval, 𝜇 is the growth
rate of the stock price, 𝜎 its volatility and 𝑟0 , 𝑟1 , . . . 𝑟𝑛−1 are normally
distributed random numbers with means 0 and unit standard deviation.
A common technique to evaluate the expected price of the stock is to
simulate 𝑁 realization of the equation above providing 𝑥 0 , 𝜇, 𝜎 and Δ𝑡 .
A few realizations, simulated with 𝑥 0 = 100, Δ𝑇 = 0.1, 𝜇 = 𝜎 = 0.01, are
shown in Fig. XX
From these realizations alone it is rather hard to make any prediction.
Running a large number of these simulations allows to get a better idea
of the price evolution. Three cases are reported in Fig. XX. As you can
see there depending on the growth rate and volatility it is a more or less
good idea to invest in the stock.
Random walk in 2D
Random walks can be generalized to the 2D. The process is very similar to
the 1D cases except that each walker can not only go left or right but also
16 Monte Carlo schemes 149
Consider the the 2D random walk presented above. During the simulation
the motion of each particle is uncorrelated with those of the others
particles. The probability for one particle to move in a given direction is
independent of the position of the other ones. As a consequence all the
possible configurations of 𝑁 particles are equiprobable. This is of course
not the case when the particles interact with each other. In presence of
interactions the particle will tend to adopt preferential configurations
that minimize the potential energy of the entire system. The simulation
of such system can be done with the Metropolis Monte-Carlo algorithm
16 Monte Carlo schemes 150
presented here.
Lattice Metropolis
For a given configuration of the atoms the total potential energy of the
system can be obtained by summing up all the interactions between the
pairs of atoms. We assume that in their initial configuration the atoms 𝐴
and 𝐵 occupied each half of the domain as represented in Fig. XX. We
want to obtain the equilibrium configuration for a given temperature 𝑇 .
The metropolis algorithm then proceed as follows:
(1) Consider one particle with a randomly chosen position on the lattice
𝑖, 𝑗
(2) Randomly chose one move for the atom at 𝑖, 𝑗 and perform that move
(3) Compute the energy change Δ𝐸 that is induced by this trial movement
(4) if Δ𝐸 < 0 the move is accepted and we return to (2)
(5) if Δ𝐸 > 0, a random number 𝑟 is generated between 0 and 1.
(6) if 𝑟 < exp(−Δ𝐸/𝑘 𝐵 𝑇), the move is accepted. Otherwise the move is
rejected. In both case we return to (2).
This very simple algorithm ensure that the final configuration (providing
enough Monte-Carlo moves have been performed) respect a Boltzman
distribution, i.e. is in equilibrium with the temperature imposed to the
system. We have reported in Fig. XX two final configurations obtained
after 10000 moves for two different temperature. During the simulations
we have set 𝐸rep = −𝐸att = 1.0. The thermal energy was then 𝑘 𝐵 𝑇 = 0.5
(low temperature) and 𝑘 𝐵 𝑇 = 10.0 (high temperature). As you can see on
the corresponding configuration the a very ordered phase is obtained at
low temperature as a very small number of unfavorable moves (i.e. move
that increases the total potential energy) are performed. On the contrary
at high temperature many unfavorable move are performed resulting in
a much more disordered phase.
Off-lattice Metropolis
two dimensional space. We assume that the particles interact with each
other via a Lennard-Jones potential
𝜎 12 𝜎 6
𝑈(𝑟) = 4𝜖 − (16.14)
𝑟 𝑟
As you can see the algorithm is almost identical to the previous case. The
principal difficulties arises in the calculation of the total energy that can
be extremely time consuming if all the pair wise interaction are calculated
16 Monte Carlo schemes 152
Approximating 𝜋 by MC
Stock Price
√
𝑥 𝑛 = 𝑥 𝑛−1 + Δ𝑡𝜇𝑥 𝑛−1 + 𝜎𝑥 𝑛−1 𝑟𝑛−1 Δ𝑡 (16.15)
Take an initial price of $50, and explore different values for the growth
rate and the volatility.
Mixing of particles
We here consider a box divided in two equal size part by a wall. One
half of the box contains 𝑁 particles that are uniformly distributed in a
random fashion. We now remove suddenly the wall and we want to study
how the molecule fills in the entire box. We can simulate this process
using a 2D random walk in a box. We set the dimension of the box to be
𝐿𝑋 = 1 ad 𝐿𝑌 = 1 for simplicity.
To initiate the simulations you will randomly place 𝑁 particles (N=100 is
a good way to start) in the region defined by the boundary [0; 12] × [0; 1].
Implement a 2D random walk dynamics to simulate the trajectory of
the particles. The particles cannot escape from the box. Hence a particle
at one edge of the domains cannot go in all 4 directions. Visualize the
results as an animation.
16 Monte Carlo schemes 153
K2 𝑡 2 K3 𝑡 3
𝑒 K𝑡 = I + K 𝑡 + + + ... (A.1)
2! 3!
We can use this expression to demonstrate that eq. 13.32 is a valid solution
of eq. 13.29
𝑑 𝑑 K𝑡
C = 𝑒 C0 (A.2)
𝑑𝑡 𝑑𝑡
𝑑 K2 𝑡 2 K3 𝑡 3
= I + K𝑡 + + + ... C0 (A.3)
𝑑𝑡 2! 3!
K3 𝑡 2
= K+K 𝑡+2
+ ... C0 (A.4)
2!
K2 𝑡 2
= K I + K𝑡 + + ... C0 (A.5)
2!
= K 𝑒 K𝑡 C0 = KC (A.6)
Exponential of a matrix
𝐴 = 𝑈Λ𝑈 −1 (A.7)
where the 𝑖 -th column of 𝑈 contains the 𝑖 -th eigenvector and where Λ
is a diagonal matrix containing the eigenvalues of 𝐴. The characteristic
equations 𝐴𝑋 = 𝑋Λ gives by right-multiplying by 𝑈 −1
1
𝑒𝐴 = 𝐼 + 𝑈Λ𝑈 −1 + 𝑈Λ2 𝑈 −1 + ... (A.10)
2!
1 2
= 𝑈 𝐼+Λ+ Λ + ... 𝑈 −1 (A.11)
2!
= 𝑈 𝑒 Λ 𝑈 −1 (A.12)
Index