1. health and fitness
2. disease
3. diabetes

MSc course Fluvial systems (GEO4-4436):
MATLAB MANUAL
Marion Tissier
Maarten Kleinhans
August 30, 2014
Contents
1 Scripts, matrices and syntax
1.1 Getting started . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.1.2 Creating you own work-directory . . . . . . . . . . . . . . . . . . . .
1.1.3 Introductory exercises . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Script files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Arrays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.2 Manipulating arrays . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.3 Operations on arrays . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3.4 Matlab built-in functions for arrays . . . . . . . . . . . . . . . . . . .
1.4 Visualising data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Final assignment : Analysis of the Rhine flow discharge measured at Lobith
2 Tests, loops and functions
2.1 Relational and logical operators
.
.
.
.
.
.
.
.
.
.
.
.
5
5
5
6
6
10
13
13
14
17
19
23
29
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
33
1
.
.
.
.
.
.
.
.
.
.
.
–
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
part I
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2.2
.
.
.
.
.
.
37
37
39
44
47
50
3 Numerical modelling
3.1 A very simple morphodynamic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Flow computation: backwater equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Morphodynamic model with flow computation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
55
58
60
4 Creative Matlab exercise
65
2.3
2.4
Flow control . . . . . . . . . . . . . . .
2.2.1 if-statements . . . . . . . . . .
2.2.2 for-loops . . . . . . . . . . . . .
2.2.3 while-loops . . . . . . . . . . .
User-defined functions . . . . . . . . . .
Final assignment : Analysis of the Rhine
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . .
flow discharge measured at Lobith
2
.
.
.
.
.
–
. . . .
. . . .
. . . .
. . . .
. . . .
part II
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Begin at the beginning
• Each chapter of this manual corresponds to one of the 4-hour MATLAB practicals. Below you can find a short
description of their content.
First practical How to use the main features of Matlab (manipulation of arrays of data and use of script files).
First data analysis and visualisation.
Second practical How to write more complex programs using tests, loops. How to create your own functions
and use them in a script file to analyse data.
Third practical Introduction to the modelling of physical processes. During this practical you will build a
simple morphodynamic model.
Fourth practical Creative assignment, where you should be able to answer a research question using Matlab.
• During the practicals, you will have to work in couples. Change partner every practical.
• The Matlab assignments should be handed in before Friday 4:00 pm in the week that the computer practical
took place. The different files (Matlab code and answers to the questions) should be combined in a zip-file and
sent to a.w.baar@uu.nl. The answers to the questions should be written as comments in the script files or (if
specified) in a separate word document.
3
The evaluation of the codes will be based on three main criteria :
- does the program run?
- does it do the right thing?
- can we read/understand it easily? (= is it commented enough?)
• You will often get error messages while running your Matlab programs. Do not panic and carefully read the
displayed messages. Solving them without help is part of the learning process! The devil is in the details: most
errors are in syntax. To understand what is going wrong with the larger code and data matrices, it may help to
play with very small and simple matrices first.
4
Chapter 1
Scripts, matrices and syntax
1.1
Getting started
To start MATLAB, navigate into the program list and click on ”Matlab.exe”.
1.1.1
Overview
Your Matlab window should roughly look as in Figure 1.1. If it does not, click on ”View” in the menu bar and select
”Desktop layout”, ”Default”. Four main divisions appear in the window :
- ”Current folder” : displays the content of the current work-directory,
- ”Workspace” : lists the variables assigned by the user,
- ”Command history” : shows the last commands used,
5
- ”Command window” : where the commands are typed and the outputs displayed.
In the following, we will mainly work in the ”Command window”, highlighted in Figure 1.1.
1.1.2
Creating you own work-directory
The current directory is displayed on the top part of the Matlab window (see Fig. 1.1).
We advise you to start every new practical with creating a new work-directory. To do that, click the ”browse for
folder” button (see Fig. 1.1) right to the ”Current directory” drop-down list on the Matlab button-bar. You can then
browse through the directories and create a new folder at the desired location using the ”New directory” button. For
instance, create a first folder called Matlab_practicals. In this folder create a folder called Practical1 in which
your will be working on today.
You can also use this button to move to an already existing directory.
1.1.3
Introductory exercises
We will perform in this section some introductory exercises where Matlab is basically used as an ”advanced” calculator.
• Type in the command window a = 32.8 and then press ”enter”. The following lines appear:
>> a = 32.8
a =
32.8000
With this command, we defined a new variable a and assigned the value 32.8 to it.
• Now let’s type:
>> b = 3*a^2
and then press ”enter”. We defined a new variable, b, equal to 3 times the value of a squared. Its value is displayed
in the command window. The basic operations and their symbols are presented in Table 1.1.
6
Figure 1.1: Organisation of the Matlab window
7
Operation
Addition
Substraction
Multiplication
Division
Exponentiation
Symbol
+
−
∗
/
ˆ(e.g. 2ˆ3 = 8)
Table 1.1: Basic operations and associated symbols
• Type now:
>> whos
and then press ”enter”. The following lines are displayed in the Matlab command window:
Name
Size
Bytes Class
a
1x1
8 double
b
1x1
8 double
This command displays the different variables which have been previously defined (a and b), and gives their
main characteristics (see presentation of arrays, §1.3, for more explanations). These variables are now listed in
the ”Workspace” sub-window, and can be visualised by double-clicking on them.
• Use the command clear, and then the command whos again. What is displayed? Any interpretation?
NB: To get more information about the function clear, type help clear and press ”enter” in the command
window.
• Type now:
>> c=a*30/200
and then press ”enter”. We are defining a new variable c, as function of a, which has been deleted previously
(clear). The command cannot be executed since a is not defined anymore. An error message, in red, appears
in the command window.
8
Function
sin(), cos(), tan()
asin(), acos(), atan()
abs()
sqrt()
round()
ceil() (or floor())
log(), log10()
exp()
Description
sine, cosine and tangent (arguments in radians)
inverse sine, cosine and tangent
absolute value
root square (equivalent to ˆ(0.5))
round towards nearest integer
round towards plus (or minus) infinity
natural logarithm, logarithm base 10
exponential
Table 1.2: List of some mathematical functions
• Type:
>> a =-4, c=a*30/200
and then press ”enter”. Two commands are written in the same line, separated by a comma. They are executed
successively by Matlab. a is first defined and then used in the calculation of c.
• Type now the same command, but replace the comma by a semicolon. What is the difference? And if you add
a semicolon at the end of the line?
• Type now:
>> c_bis =abs(c)
In this command, we applied the Matlab function abs() to the variable c. The result of this operation is a second
variable, c_bis, which is the absolute value of c. c is called input argument of the function abs. c_bis is the
output argument.
Numerous elementary mathematical functions are already implemented in Matlab. The table 1.2 presents a
selection of the most commonly-used elementary functions. Type help elfun in the command window for an
extensive list of these functions.
9
Remarks
- To obtain help in Matlab:
help function name displays information about the function directly in the command window;
doc function name opens a separate window where a more detailed help is displayed. The explanations usually
include a few examples, which are worth studying carefully;
lookfor keyword looks for all the Matlab functions related to the keyword;
- The commands you’ve recently entered can be retrieved by pressing the up-arrow key; They are also displayed
in the ”Command history” window.
- Variable names must start by a letter, and can be followed by character chain, including letters, numbers and
any symbols except those used as mathematical operators by Matlab (+,-,*,...) or punctuation signs. Variable
names can contain up to 31 characters; When no variable name is specified, a temporary variable name, ans, is
generated by default.
- Try to use the name of a function to name one of your variables.
1.2
Script files
A script, also called program, is a succession of Matlab commands written in a separate text file. Using script files
allows you to organize and save series of commands. These scripts can be edited/extended/corrected anytime you
want, and executed several times. Writing the commands into a script is necessary when you perform advanced data
analysis, since it allows you to check/improve your methodology, but also to be sure to perform the exact same data
processing for several datasets for instance.
The scripts are user-created text files which can be written using Matlab editor. Type:
>> edit
10
The Matlab editor window opens. In this window, you can write series of commands and save them. The resulting
script is saved as *.m file. The ”new script” button (see Fig. 1.1) can also be used to open the editor.
It is important to add comment lines to your scripts to make them easier to understand. This is particularly
useful when scripts are re-used long after they have been created. In a script file, anything written to the right of a
percentage symbol is considered as a comment and will be ignored during the execution. In the Matlab editor, the
comments appear in a different color (green by default).
It is recommended to always give a similar structure to the script files. An example of structure for a ”classic”
script file is given below:
% Header
Describe the script in a few words (objective, date of creation, author)
% Initialization
Clear your workspace using the command clear
Close all the figure windows (close all)
Define the main variables, load data files used in the script
% Calculations
Execute the calculations
% Output and visualisation
Commands to display the results, plot graphics, save variables
A simple example of script can be find below:
11
% Example of a very simple Matlab script
% M. Tissier, 8-8-2012
% Initialization
clear
displacement = 12;
time = 10;
% displacement (m)
% duration (s)
% Calculation
v = displacement/time;
% migration speed (m/s)
% Output
disp([’The migration speed is ’ num2str(v) ’m/s’])
Once the script has been saved, you can run it by typing its name (without the ”.m” extension) in the ”Command
window” and then press ”enter”. When a script file is executed, Matlab execute the commands in the order they are
written.
EXERCISE 1 The aim is to calculate flow velocity and water depth in a river from known variables and parameters.
(Find out through a web search what the difference is between ”variable” and ”parameter”.)
Use the law of Chezy (see Kleinhans (2005) paper for equations). The flow resistance parameter for this river is
C = 44 m0.5 /s, the bankful discharge Qbf = 2500 m3 /s, the channel width W = 500 m and the channel gradient
S = 1.6 × 10−4 m/m.
The commands will be written/organised in a small script file called Exercise1.m, saved in your current directory.
Do not forget to comment the file, and give meaningful names to the variables.
HINT: 1.6 × 10−4 can be written in Matlab as 1.6e-4.
12
1.3
1.3.1
Arrays
Definition
The array is the fundamental form that Matlab uses to store and manipulate data. This way of storing data allows
for efficient calculations and analysis over entire datasets.
An array is basically a list of elements (numbers or characters, see below) organised in rows and/or columns.
Different types of arrays can be defined:
- The simplest array, or scalar, contains only one value (1 row, 1 column).
- The vector is a one-dimensional array which is basically a row or column of values.
- Two-dimensional arrays are also called matrices, and are characterized by their amount of rows and columns.
- Multi-dimensional arrays (with 3 or more dimensions) can also be defined.
These arrays can contain different type of data, for instance:
- Integers : whole numbers (-10 -4 1 0 3 12)
- Floating numbers : fractional numbers (3.271 -2345.17...)
- Characters : text ASCII character (a b z / ] ?)
In the two first cases, they are called double arrays, and in the third, character array.
In this section, we will learn how to create arrays in Matlab, extract the relevant information from one of them,
and perform calculations with these arrays.
13
1.3.2
Manipulating arrays
Creating vectors and matrices
• Type:
>> M = [1 2 5 10; 3 4 -2 7]
and then press ”enter”. A 2-by-4 matrix, called M, has been created and is now displayed in the command window. When defining an array other than a scalar, brackets [ ] should be used. The elements belonging to a same
row are separated by a ”space” (or a comma), while semicolons are used to skip to the next line of the array.
• Use the whos command, and check the information concerning the new variable M. The size of the array is
displayed. It gives 2 × 4, which means that M has 2 rows and 4 columns. This information can also be retrived
using the function size():
>> size(M)
which returns a 1 × 2 array containing the number of rows and columns of M.
• Type now
>> u=[2; 4; 6; 3]
What kind of array do you obtain? Type now [n_r n_c]=size(u); and check the values of n_r and n_c. What
do they contain?
• Vectors can also be created in the following way. Type
>> k = [1:1:6], l = [1:2:6], m=[12:-3:0]
and analyse the outputs. Here vectors are defined in the following form: [begin:interval:end], where
interval can be any positive or negative number.
Remark: When interval is equal to 1, it does not need to be specified. 1:6 returns the same vector as 1:1:6.
• Similar commands can be applied to create matrices. Type:
>> Z = [1:9; 3:11]
and then press ”enter”.
14
• Type now
>> Z = [1:8; 3:11]
Why does it produce an error?
Concatenation It can be useful to combine several related arrays in one unique array, for instance when two datasets
are related and should be treated together. This can be done thanks to the operation of concatenation (=joining the
arrays one after the other).
• Type for instance:
>> a = [1:3]
>> b = [4:6]
>> c = [a b]
The space is used in this case to concatenate the arrays horizontally (alternatively, you can use the function
horzcat(), e.g. c = horzcat(a,b)).
• Type now
>> d=[a ; b]
In this case the arrays have been concatenated vertically (see also vertcat()).
• Concatenation can be used to add an extra column (or line) to a matrix. Type
>> e = [10;20]
>> f = [d e]
and analyse the outputs.
• Vertical (horizontal) concatenations can only be performed if the arrays have the same number of columns (rows).
Type:
>> [f; a]
Do you understand the error message?
15
Manipulating arrays by their subscripts Elements of arrays are usually identified by their position (row and
column numbers), also called subscripts. The following table shows how the different elements of a 3 × 2 matrix can
be identified using their subscripts:
(1,1)
(2,1)
(3,1)
(1,2)
(2,2)
(3,2)
• Type for instance
>> M(2,4)
It returns the element located at the 2nd line, 4th column. Type now
>> M(3,4)
Do you understand the error message?
• This notation can be used to alter a value of an array. Type for instance:
>> M(1,2)=0
• Several elements can be selected using a vector instead of a scalar subscript. Type for instance
>> M(2,[1 3])
The first and third elements of the second row of the matrix M are displayed.
• To create a vector containing the entire row, type for instance:
>> vectM = M(2,:)
Similarly M(:,2) would refer to the second column of M.
• Replace the elements of the first column of M by 10.
Manipulating arrays by their indices
16
• For a vector, it is easier to define the position of an element using only one number, called its index, which is
its position in the row or the column. For instance, vectM(2) is equivalent to vectM(1,2). Several elements of
vectM can also be selected using a vector containing the relevant indices. Execute for instance:
>> vectM([1:3])
• Similar conventions can be used for matrices. The elements are then numbered column by column, starting from
the top left element and progressing downward. For a matrix of dimension 3 × 2 (note that the index of a given
element will depend on the size of the matrix), we have:
(1)
(2)
(3)
(4)
(5)
(6)
Type the following commands and compare the output:
>> M(5)
>> M(1,3)
• How would you call the element M(1,4) using indices?
1.3.3
Operations on arrays
Operations between a scalar and an array
• Multiplication by a scalar:
All the elements of the array are multiplied by the scalar. Type for instance :
>> 2*M
and check the output. Remark: 2*M is equivalent to M*2.
17
• Addition (subtraction) with a scalar:
Type
>> M-2
The scalar is added (subtracted) from each element of M.
Operations between two matrices
• Addition (subtraction) of two matrices. These operations can only be applied to matrices of identical size. The
resulting matrix is obtained by adding (subtracting) their corresponding elements. Example:
>> M-f
• Multiplication of two matrices. This operation is executed according to the rules of linear algebra. Let’s define
two matrices A and B such as:


1 3
1 0 2
A=
and B =  4 1  .
5 10 7
2 2
The element (i, j) of the matrix C = A ∗ B is defined as :
C(i, j) =
n
X
A(i, k) ∗ B(k, j),
k=1
where n is the number of columns of A (=number if rows of B). C is therefore equal to:
1∗1+0∗4+2∗2 1∗3+0∗1+2∗2
5 7
C=
=
5 ∗ 1 + 10 ∗ 4 + 7 ∗ 2 5 ∗ 3 + 10 ∗ 1 + 7 ∗ 2
59 39
18
Operation
between
two matrices
Addition
Matlab expression
A+B
Multiplication
A*B
Element-by-element multiplication
Element-by-element division
Condition
size of A = size of B
number of columns of A = number of rows of B
A.*B
size of A = size of B
A./B
size of A = size of B
Table 1.3: Summary of the main operations between two matrices A and B. The requirements in term of matrix size
are also indicated.
Element-by-element operations
• To perform a multiplication element-by-element, it is necessary to use the operator ”.*” instead of ”*”. Type
for instance
>> M.*f
and examine the result. Type now the same command but without the dot. Do you understand the error message? (see conditions in Table 1.3). Similarly, element-by-element division is obtained using the operator ”./”,
and element-by-element exponentiation using ”.^”.
• These operations can be combined. Type for instance:
>> MM = vectM.^2 + vectM
1.3.4
Matlab built-in functions for arrays
Elementary mathematical functions The functions defined in Table 1.2 can also be applied to arrays. In this
case, the function is applied to each matrix element individually, and returns an array of the same size as the input
19
argument. Compare the outputs of the two following commands:
>> M/3
>> round(M/3)
Transposition Enter the following command
>> v = transpose(u)
and compare v to u. Now look at transpose(M) and compare it with M. Do you understand how the transposition
function works?
transpose() changes the rows of the matrix into columns, and vice versa. As a consequence, a matrix of size n × m
will turn into a matrix of size m × n containing the same elements organised in a different way. This operation can
also be performed using the symbol ’ just after the name of the array. u’ is therefore equivalent to transpose(u).
Creation of arrays The functions zeros() and ones() can be used to create a matrix. Type for instance:
>> C = ones(3,5)
Do you understand the meaning of the arguments of this function? As expected, the function zeros() is similar to
ones() but returns a matrix filled with zeros.
• With the help of the functions introduced previously , create a matrix of size 20 × 30 containing only 5s.
Data analysis The main Matlab functions available to perform data analysis are summarized in Table 1.4. Some
examples of the use of these functions are given below.
• Type the following commands:
>> mean(vectM), median(vectM), sum(vectM)
Each of these functions, applied to the vector vectM, returns a scalar, which is equal to (respectively) the averaged
value of the vector , its median value and to the sum of its elements.
• Type:
>> vmax =max(vectM)
20
vmax is the largest value of the vector vectM.
In some practical cases, it is important to know when the maximum value has been reached, that is to say what
is the position of the largest value in the original vector. In order to get this information, the function max (or
min) can be used with two output arguments. Type
>> [vmax ind]=max(vectM)
The second argument (ind) is the index of the maximum value in the vector. The function min works in a similar
way.
• Type now:
>> s = sort(vectM)
Compare s and vectM.
• These functions can also be applied to matrices. In this case, they handle each column separately. The output
is therefore a vector of length the amount of columns. Type for instance:
>> mean(M)
and analyse the result. When M is a matrix, mean(M) is a row vector containing the mean value of each column.
• It is generally possible to specify the direction along which the computations have to be performed as an additional input argument. Type for instance:
>> mean(M,1)
>> mean(M,2)
and compare the results. Do you understand the meaning of the second input argument?
The functions median, min, max and sum work in a similar way. You can check the Matlab help to find out how
to choose the direction of computation for these functions.
• An alternative approach is to apply the functions after transposition of the matrix. Type for instance:
>> mean(M’)
and compare with mean(M,2).
21
max
mean
median
min
std
var
prctile
sum
sort
Maximum value
Average or mean value
Median value
Smallest value
Standard deviation
Variance
Percentile values (the percent of interest should be specified)
Sum of the elements
Arrange the elements in ascending order
Table 1.4: Some important Matlab functions for data analysis
EXERCISE 2 Make a new file called Exercise2.m in your work-directory, where the commands needed to solve
the following assignment will be written.
1. Create the following matrix:

16
 5
A=
 9
4
3
10
6
15

2 13
11 8 

7 12 
14 1
2. Define a vector A1 containing the first row of A. Define a scalar, sumA1, equal to the sum of its elements.
3. Create now a vector sumRows
P4 where each element contains the sum of the elements of one of the row of A. For
instance, sumRows(1) = k=1 A(k, 1). Can you do it in only one line of command?
4. Create a vector containing the sum of the column of A. What do you notice?
22
1.4
Visualising data
How to make a simple 2D-plot : example
• Define the following vectors:
x = -pi:.1:pi;
y = sin(x);
y2 = cos(x);
• Type
figure,plot(x,y)
The first command, figure, opens a new figure window, called Figure 1 by default since no other figure has
been opened previously. The plot command is then executed in the newly created figure.
• Type now
plot(x,y2)
The new plot has replaced the previous one in the figure window. To display the second graph on top of the
first one, the command hold on should to used. Type:
figure,plot(x,y)
hold on,plot(x,y2,’r’)
A new figure, called Figure 2 by default, is created, and contains the two plots on top of each other. For the
second plot, the line color is specified using the additional argument ’r’, with r standing for red.
• Title, labels on the x- and y-axis and legends can now be defined. Type the following commands one by one and
analyse their effects on the figure:
title(’Example plot’);
xlabel(’Radians’);
ylabel(’Function Value’);
legend(’sine’,’cosine’);
23
• To specify the limits of the x- and y-axis, type: for instance
axis([-pi pi -2 2])
Do you understand the meaning of these 4 numbers?
• The figure command can also be used with an input argument which is a positive integer, as for instance
figure(1). The number corresponds to the number of the corresponding figure window. If the Figure 1 does
not exist yet, this command will create it. If this figure already exists, it becomes the active figure window. This
means that new plot commands will be executed in this figure. Type for instance:
figure(1), hold on, plot(x,-y,’.’)
Because of the argument ’.’, the vectors are plotted as a cloud of points.
• plot can also be used without specifying the x argument. Type: figure,plot(y)
The elements are now plotted as a function of their position in the vector: y(1) at x = 1, y(2) at x = 2, etc.
• Multiple plots can also be displayed in one figure. This can be done using the command subplot(a,b,c) where
a, b and c are integers, and c is less than or equal to the product a*b. This command divides the figure in a*b
sub-figures, organized in a rows and b columns. To obtain 2 figures on top of each other, the following lines
should be used:
figure
subplot(2,1,1);plot(x,y)
subplot(2,1,2);plot(x,y2,’g’)
• It is possible to display on the figure the coordinates of a given data point using the ”data cursor” button
presented in Figure 1.2. Select this button and click on a data point.
• You can now close all the figures using the command:
close all
The function plot is the basic Matlab function for 2D plots. Numerous options are available to customize your
plots. Some examples are given below:
24
- plot(x,y,’r’): red solid line
- plot(x,y,’--r’): red dashed line
- plot(x,y,’--r’,’linewidth’,2): thick red dashed line
- plot(x,y,’*g’): the points are marked with green asterisks (no line)
- plot(x,y,’+-r’): solid red line with points marked with a plus sign.
The tables below show some of the colors and markers available. We advise you to look at the Matlab help (type
doc plot for instance) for a more extensive description of the plotting options. You can also navigate in the help
window to learn more about the numerous ways to visualise data in Matlab.
red
green
r
g
plus sign
asterisk
Color specifiers
blue b black k
cyan c yellow y
+
*
Marker specifiers
square
s
triangles v, < or >
magenta
white
circle
cross
m
w
o
x
EXERCISE 3 This exercise is based on the dataset contained in the file MPMtransportdata.xls. Meyer-Peter
and Mueller (1948) (abbreviated MPM) graphically report gravel transport data and derive their famous empirical
bedload transport predictor from this dataset. Wong and Parker (2006) recovered the original data and re-analysed
it to find that the original fit by MPM was wrong. Here you will plot the data in nondimensional form and calculate
a predicted transport. The Kleinhans 2005 paper contains all the equations necessary to solve this exercise. Remark:
We assume that the density of the water is 1000 kg/m3 .
25
Figure 1.2: Figure window
Remarks :
• It is recommended to test parts of your script regularly. An easy way to do it is to copy-paste these parts into
the command window and then press ”enter”. Alternatively, you can select the part of the text you want to run
and then press ”F9”.
• Always terminate the command by a semicolon if you do not want the result to be displayed in your command
window. This is important when working with large datasets.
1. Make a new script file called Exercise3.m in which the series of commands necessary to solve this assignment
will be written. Do not forget to clear your workspace and close the figures at the beginning of your script.
2. Download the file MPMtransportdata.xls and save it in your current work-directory. You can now load the
data using the following line:
26
dataMPM = xlsread(’MPMtransportdata.xls’);
Thanks to the function xlsread(), Excel sheets can be read and loaded in your workspace. The output of the
function is an array called dataMPM which contains the same information as the Excel sheet, but without the
description of the data (only numbers). Check its number of rows and columns and compare a few values with
the initial Excel sheet.
3. Define vectors containing the channel width W, the depth h, the slope S, the median grain size D50, the specific
gravity of the sediment s and the sediment transport rate qs.
4. From these vectors, calculate the total shear stress. The result will be a new vector, called tau for instance.
5. Calculate the Shields parameter using the shear stress previously defined.
6. Calculate the Einstein parameter.
7. Plot the sediment transport rate (y-axis) as a function of the shear stress (x-axis).
Plot the data in log-scales for both the x- and y- axis. To do that, the function plot() can be replaced by
loglog(). loglog is basically the same as the plot function, but with logarithmic axis, and can be used with
the same input arguments (see Matlab help for more details). Give an appropriate title to the plot, as well as
names to the x- and y-axis.
HINT: To write Greek letters in the text of a Matlab figure (axis labels, title, legend, etc.), use the following
characters:
α
β
\alpha
\beta
τ
Φ
θ
\tau
\Phi
\epsilon
\theta
φ
λ
\phi
\lambda
To write subscripts or superscripts, use the symbols _ and ^. For example, to display qs (m2 /s), you should
write q_s (m^2/s), and to display D50 (2 characters as subscript), use D_{50}.
27
8. Plot in a new figure the data in a nondimensional form, i.e. the Einstein parameter (y-axis) as a function of the
Shields parameter (x-axis). Give an appropriate title to the plot, as well as names to axis. Compare with the
previous figure : does the adimensionalization seem efficient?
9. Plot the MPM prediction as a additional red line in the second figure. To plot this line, first define two vectors
x and y such as x contains a range of values reachable by the Shields parameter (θ) and y = 8(x − θcr )1.5 (see
eq. (30) in Kleinhans 2005). We use θcr = 0.047.
10. In the same figure, plot a vertical line corresponding to x = θcr .
11. To further analyse the data, we will now remake the figures of the questions 7 and 8 using the function scatter().
scatter is a function which allows you to plot the data as circles of different sizes and colors depending on their
characteristics. Check the Matlab help to understand how to use the function scatter() and what are its input
arguments. The lines below are an example of use of the function:
figure; scatter(tau,qs,10,s,’filled’);
set(gca,’xscale’,’log’,’yscale’,’log’);
colorbar;
Execute these lines one-by-one and analyse the outputs. How many colors are displayed in the figure and why?
do you understand why the size of the circles does not vary?
12. Create a new figure, which will be divided in two subplots.
- In the top subplot of the figure, plot sediment transport (y-axis) as a function of shear stress (x-axis), with
colors depending on the grain size D∗ (Bonnefille number).
- On the bottom part, plot the same figure but in its nondimensional form. Add colorbars to both plots.
For a better distribution of the colors, it is recommended to use the logarithm of the Bonnefille number as the
”color vector” argument of the function scatter.
13. Look at the data plotted in non dimensional form. One set of points seems to follow a different trend than the
others. To which range of D∗ does it correspond?
28
1.5
Final assignment : Analysis of the Rhine flow discharge measured
at Lobith – part I
In this exercise, we will perform some analysis on a dataset containing the flow discharge of the Rhine measured at
Lobith (where the Rhine enters the Netherlands) for the period 1901-2000 AD. The dataset is given in the text-file
Qlobith19012000.txt as an array. Each element of this array is the discharge in m3 /s for 1 day, and each column
corresponds to one year (first column = 1901, last column = 2000).
Preliminary analysis
1. Create a new script file called Lobith1.m in which the assignments of the exercise below are executed. Start the
script by clearing your workspace and closing the windows.
2. Load the new dataset. As the data are stored as a text (ASCII) file, the following command should be used:
Q=load(’Qlobith19012000.txt’);
Check the size of the array Q. Do you understand how it is organised?
3. The data corresponding to the 29th of February are included for each year, leading to a number of lines of 366.
Display the 20 first elements of the row corresponding to the 29th of February. What do you notice?
Remark: NaN stands for ”not a number”. It is often used to represent missing values in datasets.
4. To obtain an overview of the dataset, use the following command:
figure,plot(Q)
In this command, the function plot is applied to a matrix instead of a vector. plot handles each column (and
therefore each year of data) separately. As a result, each line appearing in the figure corresponds to the evolution
of the discharge for a given year. Give a title to the figure and names to the x- and y-axis.
5. Zoom on the figure around the 29th of February (you can use the ”zoom” buttons presented in Figure 1.2 or
redefine the x-axis limits using xlim()). How are the NaN-values handled by the function plot?
29
6. Create a vector called year containing the years when the data were collected, as well as a vector days containing
the day numbers (1=first of January).
Analysis of one year of data
7. Create a new vector called Lobith1916 containing the flow measured at Lobith in 1916.
8. Determine the mean and median discharges for this year. Compute also the maximum and minimum values of
the discharge. When were these values reached?
9. Compute the mean discharge for the following year (1917). What do you obtain?
To understand better what is happening, check the results of the following operations:
>> NaN + 10
>> 20 * NaN
This explains why we cannot compute directly the mean value for a year without 29th of February.
10. To ignore the NaN values in the computations, alternative functions can be used, such as nanmean, nansum, etc.
Use the Matlab help for a description of these functions and compute the mean discharge in 1917.
Statistics over the entire data set
11. Build a vector of size 1 × 366, where each element is the discharge for a given day averaged over 100 years.
12. Plot this vector as a function of the time. When are the highest mean discharges occurring (which month?)?
Plot in the same figure the maximum et minimum values for each day (use different type of lines to differentiate
the curves). Add names to the axis and a legend to the graph.
13. Create a vector containing the annual flood discharges for the 100 years of data.
14. Compute the mean annual flood discharge over the entire period.
30
15. Compute the median value of the annual flood data as well as the 75th, 90th, 95th percentiles (see http://en.wikipedia.org/wik
for definition).
16. We will now use the annual flood discharge data to build a flood-frequency curve, which is a graph showing the
relationship between flood magnitude and their recurrence interval for a specified site.
The recurrence interval (T r) of a flood is a statistical measure of how often a flood of a given magnitude (QT r )
is likely to be equalled or exceeded. For instance, the ”fifty-year flood” (Q50 ) is one which will, on the average,
be equalled or exceeded once in any fifty-year period.
Follow the following procedure to compute the recurrence interval:
- Create a vector QT r containing the annual flood discharges sorted in ascending order.
- Define a vector R of same length as QT r containing the ”ranks” of the discharges. The largest discharge of
QT r should have a rank equal to 1, and the smallest a rank equal to N, where N is the number of years for
which we have flood data.
- The recurrence interval T r (in years) can then be computed using the formula:
T r = (N + 1)/R.
You can now plot QT r versus T r : this is the flood-frequency curve. Use logarithmic axis for a better display of
the data. Add a title and labels to the axis.
17. Read in the figure what is the magnitude of the fifty-year flood (use the data cursor button indicated Fig. 1.2
to help you). What about the 100-year flood? (the answers should be written as a comment in your script file).
Remark Numerous Matlab functions are available to study the distribution of data and the occurrence of a given
event. Histograms can be built using the function hist, the distribution of the data can be compared with different
probability distribution with the help of probplot, etc.
31
PRODUCTS TO HAND IN Exercises 1 to 3 and final assignment (Section 1.5). Please note that only the
relevant information should be displayed while running the scripts ⇒ use semicolons to avoid the display of endless
lists of numbers!
32
Chapter 2
Tests, loops and functions
2.1
Relational and logical operators
Relational operators A relational operator (see Table 2.1) compares two numbers (or arrays) to determine if a
comparison statement (e.g., a < b) is true or false.
• If the comparison statement is false, it returns a zero. Type for instance:
>> 3<2
If the statement is true, such as
>> 3<6
it returns 1.
• Such comparisons are generally performed on arrays and can be very useful to analyse datasets. Define an array
A such as:
>> A = [1:2:10;10:-2:1]
33
and then type
>> E = (A<5)
This command returns an array E which has the same size as A, but contains only ones and zeros. Each element
of A is compared to 5, and, depending on the outcome of each comparison, the elements of E are filled with 1 or
0.
• Use the function sum to determine how many elements of A are smaller than 5.
This method can be very useful to characterize large datasets where it is not possible to count one-by-one the
elements satisfying a certain condition.
• Relational operators can also be used to compare two arrays. Type:
>> B = [3 7 1 8 3; 1 45 0 2 9]
>> F = (A<B)
and analyse the result. The comparison between A and B is done element-by-element. The results are stored in
a third array, F, which has the same size as A and B. Each element of F is equal to 1 if the comparison is true at
this location, and equal to 0 otherwise.
• These operations can be performed using any operator described in Table 2.1. Type for instance:
>> G=(A==B)
Do you understand the output?
• Arrays such as E, F and G, resulting from relational and logical operations, are called ”boolean” or ”logical”
arrays. They can be used to select elements from another array. Type:
>> a = A(E)
A new vector, a, is defined, and contains the elements of A corresponding to the locations where the ”logical”
matrix E was equal to one. It therefore contains the elements of A smaller than 5. The matrix E can be used to
extract elements from any other matrix, as long as all the matrices have the same size. Type for instance:
>> b = B(E)
and interpret the result.
34
Description
Less than
More than
Less than or equal to
More than or equal to
Equal to
Not equal to
Symbol
<
>
<=
>=
==
∼=
Table 2.1: Description of the different relational operators and their symbols in Matlab
Of note, it would have been possible to write directly (without defining separately the logical array E):
>> a = A(A<5)
>> b = B(A<5)
Logical operators
2.2).
Comparison statements can be combined using the logical operators AND and OR (see Table
• Type for instance:
>> H=(B>2)&(B<10)
H contains ones only if the elements of B are larger than 2 and smaller than 10. The two comparison statements
should be true at the same time to return a 1.
• When the operator OR is used, the output is true if either one, or both, comparison statements are true. Replace
the & in the previous command by the OR symbol |. Did you expect the output?
• A third logical operator, NOT, can be used. Applied to a logical array, it returns its opposite. Type for instance:
∼ H, and compare with H.
35
Description
OR
AND
NOT
Symbol
|
&
∼
Table 2.2: Logical operators
Use of the function find The function find looks for the non-zeros elements of an array, and returns their indices.
If no elements are found, it returns an empty matrix.
• Type the following commands:
>> u=[32 6 4 9 -6 14]
>> v=(u>7)
>> ind = find(v)
find returns a vector, ind, containing the indices of the non-zeros values in v, that is to say the indices of the
elements of u larger than 7. Again, the two last lines could have been written in one time as:
>> ind = find(u>7)
• Type now:
>> ind2 = find(A<5)
We obtain the indices of the elements of A smaller than 5 (see definitions of indices for matrices page 16). Compute now A(ind2) and compare it with A(E). Do you understand the result?
EXERCISE 1 Open a new script file where the following assignments will be solved. The matrix A should be defined
at the beginning of the script.
1. Compute the average value of the elements of A larger than 2 and smaller than 6.
2. Replace the elements of A smaller than 3 by 0.
36
3. Replace the elements of A larger or equal to 9 by 2 times their value.
2.2
Flow control
In a simple program, the commands are executed one after the other, independently from the output of the previous
one. However, it is sometimes necessary to include some sort of decision making in the program. For instance, we can
decide to compute the sediment transport using different predictors depending on the sediment characteristics and
flow conditions. The program should therefore be able to choose the appropriate calculation and skip the others. It
can also be needed to repeat a sequence of commands several times, until a certain condition is reached for instance.
The order in which the commands or group of commands are executed, also called the flow of the program, can be
controlled in several ways in Matlab. In this section, the three basic flow control methods will be presented.
2.2.1
if-statements
Purpose if-statements are used when a decision must be taken depending on the outcome of a comparison statement (also called conditional expression). If the comparison is true, then the group of commands are executed, and if
it is not, it is skipped.
Format The most simple if-statement is the if-end construction, as in the following example:
x=2;
if x<5
disp(’x is less than 5’)
end
• Execute this small program and interpret the result.
• Replace x=2 by x=5 and re-run the program.
37
When the program has to choose between two alternatives, a if-else-end-constructions can be used. Read carefully
the piece of code presented below (please note that in this example Theta is a scalar):
if Theta<Theta_cr
%no sediment transport
Phi = 0;
else
%sediment transport computation
Phi = 8*(Theta - Theta_cr)^1.5;
end
Based on the outcome of the comparison statement Theta<Theta_cr, either the commands between the if and else
or the commands between else and end are executed. This small program is therefore simply expressing the fact
that if the Shields parameter Theta is below a critical value, there is no sediment transport, and that otherwise the
sediment transport is calculated using the MPM predictor.
The general structure of the if-statement is presented below:
38
if expression1
group of commands evaluated
if expression1 is true
elseif expression2
group of commands evaluated
if expression2 is true
elseif expression3
group of commands evaluated
if expression3 is true
. . .
else
group of commands evaluated if none
of the previous expressions is true
end
Only the commands associated with the first true expression are evaluated, the others are skipped. The conditional
expressions can be defined using the operators presented in §2.1.
2.2.2
for-loops
Purpose for-loops are used to repeat a command or a group of commands for a fixed (predetermined) number of
times.
Format
The general structure of a for-loop is presented below:
for index = start:increment:end
commands to be executed in the for-loop
end
39
Examples
• Let’s start with a simple example to understand how these loops are working. Run the following lines and
interpret the results:
for j=1:8
disp(’Hello’);
end
If necessary, check the Matlab help to learn more about the disp function.
• Loops can be for instance used to define arrays. Run now the following script and examine carefully its results:
clear
for j=1:9
s(j) = j^2;
end
s
The loop is executed for each element of j, which starts at one and increases by 1 until it reaches 9. At each
step, the square j is computed and stored in jth element of the s. In other words, this small script computes
the squares of the numbers from 1 to 9.
• Note that this for-loop could be replaced by the following lines:
clear
j=[1:9];
s=j.^2
These two ways of solving the problem are equivalent. However, to reduce the computational time it is recommended to avoid loops when it is possible.
• Loops can also be nested, as in the following example, where we define a new matrix B and fill it in element-byelement:
40
clear
B=zeros(2,3);
for i=1:2
for j=1:3
B(i,j) = i*j
end
end
In this case, every time i increases by one, the nested loop is executed 4 times (from j = 1 to j = 4). Run this
script to check in which order the matrix is filled in. Do you understand why?
Remark In the previous example, we first defined B as a 2 × 3 matrix of zeros, and then filled it in within
the loop. This initialising step is not compulsory in Matlab (comment the second line and re-run the script),
but it makes your program more efficient. Indeed, it is a way to pre-allocate memory: Matlab ”knows” from
the beginning the size of the matrix you are building up, and does not need to re-size it at each iteration. The
program runs faster (important when working with large arrays!).
• for-loops can be combined with if-statements, as in the following example
clear
s = [2 1 -4 23 0];
for i=1:5
if s(i)>=2
s(i) = 3*s(i);
end
end
s
Run this script and interpret the output.
In this case, the combined use of the for-loop and if-statement could have been replaced by the use of relational
41
operators. Do you understand how?
Some remarks about for-loops
- The first time through the loop, the variable is assigned the first value in the vector [start:increment:end].
For the second time through, MATLAB automatically assigns to the variable the second value in the vector, etc.
- The number of times a loop is executed equals the number of values in the vector.
- After the loop is finished executing, the loop variable still exists in memory and its value is the last value used
inside the loop.
EXERCISE 2 This exercise deals with the computation of the evolution of Carbon 14 in a sample due to radioactive
decay. Each thousand of years a fraction α of C14 is lost. This means that, calling C0 the initial mass of C14, the
mass remaining in the sample will be:
- C1 = C0 (1 − α) in 1 thousand years,
- C2 = C1 (1 − α) in 2 thousand years,
- C3 = C2 (1 − α) in 3 thousand years, etc.
We want to calculate the quantity of C14 still in the sample after 25000 years. Write a script which creates a vector
C such as its first element C(1) is the initial amount of U238, C(2) the amount after 1 thousand years, etc. The last
element of the vector should contain the amount of C14 after 25 thousand years. You can use the following script
lines as a basis and complete them. Data: C0 = 0.5 kg; α = 0.117;
42
% Initialization
???
???
% Calculation of the radioactive decay
for i=2:26
??? % Computation of the i^th element of the vector C
end
% Definition of a time vector in thousand years
t = (0:25);
% Visualisation
figure; plot(t,C/C(1));
ylabel(’C(t)/C(0)’); xlabel(’t (10^3 years)’);
EXERCISE 3 Let’s come back to exercise 3 of the last chapter (page 25) concerning the analysis of the dataset
from Meyer-Peter and Mueller. We observed that part of the data was following a different trend than the others after
non-dimensionalisation. A close look at it shows that these data points correspond to experiments performed using
a very fine sediment (D50 = 0.3402 mm). In this case, bedforms probably developed, and the use of a grain-related
shear stress instead of the total shear stress is more appropriate and might improve the prediction.
1. Copy the data file MPMtransportdata.xls in your current work directory and open a new script file.
2. Create a new vector tau2 containing the total shear stress for data such as D50 ≥ 1 mm and the grain-related
shear stress otherwise (D50 < 1 mm). You can use parts of the script written during the first tutorial.
Remarks:
- this question can be solved using if-statements and loops or relational/logical operators only. Choose the
43
most convenient way for you!
- in Matlab, the function calculating the base-10 logarithm is log10(). log() is the natural logarithm.
3. Calculate the associated Shields parameters.
4. Plot the data (using the function scatter) as the Einstein parameter as a function of the newly computed
Shields parameter. Compare with the figure created during the previous tutorial.
5. Additional question for those who found the exercise easy. Now also remove wall friction and explain how you did
it. See the Van Rijn fluid mechanics book or the paper by Wong and Parker 2006, J. of Hydraulic Engineering,
doi 10.1061/(ASCE)0733-9429(2006)132:11(1159).
2.2.3
while-loops
Purpose while-loops are used to repeat a command, or a group of commands, an indefinite number of times until
the condition specified by while is no longer satisfied.
Format
The general structure for a while-loop is the following:
while expression
group of commands repeated as long as
expression is true
end
When the program reaches the first line of the block, the conditional expression is checked. If it is true, Matlab
executes the group of commands. Then Matlab jumps back to the while command and check again the expression: if
it is still true, the commands are executed again. The looping process will stop when the expression becomes false.
44
Examples
• Execute the following script and interpret the output:
clear
x=1;
while x<14
x=x+4
end
How many times are the commands executed and why?
• Execute now:
x=1;
while x<14
x=x-4
end
Do you understand what is happening?
Remark: To terminate a computation, press Ctrl+C. This command can be very useful in case of programming
errors leading to never-ending loops.
Some remarks on the use of while-loops
- The variables used in expression must be defined BEFORE the ’while’ loop, so you can use it to enter the loop
- At least one variable used in expression must change INSIDE the while loop, or you will never exit the loop.
- After the loop is finished executing, the loop variable(s) still exist in memory and its value is the last value the
variable had in the while loop.
45
EXERCISE 4 In this exercise, we aim at characterizing the flow in typical experimental conditions. The channel
is 0.3 m wide and has a slope of 0.006 m/m. The Nikuradse roughness length is 0.0013 m. We try to answer to the
following questions : how deep will the water be when the water flow is 1 × 10−3 m3 /s? what will then be the flow
velocity?
This problem will be solved using an iteration process. We start with a first guess of the water depth, and correct
it at each loop to make it more accurate. The looping process will stop when the desired accuracy has been reached.
• Create a new script where this problem will be solved. Some help concerning the structure of the program can
be found below:
Initialisation
- Make a first inaccurate guess to define h (should be such as
0 < h < hmax , with hmax = 3 cm)
- Define also htemp to allow for the first loop to start, e.g.
htemp = 0.9 ∗ h.
As long as the absolute value of (h−htemp ) is larger than 1×10−6
- New guess for h : the ”real” depth should be contained between
h and htemp . We try the value (h+htemp )/2 (other options are possible!).
- Computation of u
- Computation of htemp
End of the loop
A good estimation of h has been found. Calculate the associate flow
velocity.
Some additional comments Each iteration starts with defining a new, presumably more accurate, value of the
flow depth h. We then compute the flow velocity u associated to h using the Chezy equation. This flow velocity
is then used to compute a second estimation of the water depth called htemp , based on the mass conservation
46
equation. h and htemp are then compared. They should be nearly equal if our last estimation of h was good
enough.
• To check your code, run the script with different guesses for the initial water depth. If the while loop is correctly
coded, it should always return the same final value for the depth h.
2.3
User-defined functions
You have already encountered functions that are built in to MATLAB (see Tables 1.2 and 1.4 for instance). As you
start to write your own programs in MATLAB you will need to create your own functions that take one or more input
arguments, operate on them, and return the result (output).
Function files are simply text files with a *.m extension, and can be written, as script files, using Matlab editor.
They are organized in the following way:
- The first line of a function m-file must always be of the following form:
function [ output variables ] = function name( input variables )
This line always starts with the word function. Following that, the output variables are enclosed in square
brackets [ ]. Then, the characters function_name define the name which will be used to call the function. The
input variables are finally specified into brackets.
NB: The function function_name should be saved in a file called function_name.m.
- The definition line of a function file must be followed by a block of comments, containing at least the following
fields: a short description of the function aim as well as definitions of the input and output variables (see example
below).
- The body of the function can then start. It is made of a succession of commands which ultimately defines the
output variables.
47
• Let’s start with an example of a simple function. Copy the following lines in a new m-file and save it in your
current directory as cylinder_charact.m:
function [S,V] = cylinder_charact(R,h)
% [S,V] = cylinder_charact(R,h);
% Computation of the volume and surface area of a cylinder
% input R radius (m)
%
h cylinder height (m)
% output
%
S surface area (m2)
%
V volume (m3)
a_lateral = 2*pi*R*h;
a_circles = 2*pi*R^2;
S = a_lateral + a_circles;
V = pi*R^2*h;
In this example, a function called cylinder_charact has been defined. The definition line is followed by a block
of comment lines, which are meant to provide information about the function and its arguments. These lines
are important since they allow the function file to be referenced by the Matlab help.
• Type now in the command window
>> help cylinder_charact
The block of comments following the definition line is displayed.
• Clear your workspace and type the following commands:
>> [S1,V1] = cylinder_charact(2,3);
The function cylinder_charact() is called with the input arguments 2 and 3. It returns the surface area and
volume of a cylinder of radius 2 m and height 3 m. These values are assigned to the variables S1 and V1.
• Type whos in the command window. The intermediate variables (=variables not defined as outputs), such as
48
a_lateral, a_circles only exist within a function and do not appear in the workspace. The use of functions
is therefore important since it avoids to overload the workspace with variables which are not going to be used
later.
• These functions can be used in a script file, as any built in function in Matlab. Modify the function file so as it
allows for non-scalar input arguments. You should then be able to run the following script:
% Computation of the volume and surface area
% of cylinders of increasing radius
% M. Tissier, 3-09-12
% Initialisation
clear
h=3;
% constant height of the cylinders (m)
r=(1:6); % increasing radius of the cylinders (m)
% Computation of the volumes S_tot and areas V_tot
% of the cylinders
[S_tot,V_tot] = cylinder_charact(r,h);
figure; plot(r,S_tot,’o-’);
xlabel(’radius (m)’);ylabel(’surface area (m^2)’)
Some (additional) reasons to use functions in your scripts
- Script files in which all the calculations are entered in one unique file can become too long and difficult to
understand. By using functions files, you can organize better your program, and make it easier to read and
understand.
- A function can be called several times within a script with different input arguments: you do not need to write
similar calculations over and over again...
49
EXERCISE 5 We now continue from exercise 3 of Chapter 1. Here, the objective is to compare the measured
sediment transport to the predictions given by the MPM and EH formulations, using several user-defined functions.
1. Create a function calculating the (total) Shields parameter. The input parameters should be the channel slope,
the densities of the water and the density of the sediment, the water depth and the median grain diameter.
2. Create two functions which calculate the non-dimensional sediment transport according to MPM and Engelund
and Hansen (EH) relations. The input arguments for MPM should be the Shields parameter and its critical
value. Choose the adequate arguments for EH.
HINT: For the MPM predictor, do not forget to account for the fact that sediment is only transported if the
Shields parameter is higher than a critical value (θcr = 0.047).
3. Create a function converting the non-dimensional transport into the dimensional one (in m2 /s).
4. Write a script which will call the functions defined above and ultimately plot the measured sediment transport
against the predicted one for both MPM and EH equations. Add a line corresponding to x = y on the figure.
Label the x- and y-axis and add a legend.
2.4
Final assignment : Analysis of the Rhine flow discharge measured
at Lobith – part II
Morphologically relevant statistics of the discharge dataset
1. Calculate the maximal and minimal daily discharge for the whole dataset.
2. Calculate the percentage of days of year 1916 where Q exceeds the mean annual flood discharge calculated in
Section 1.5.
50
3. Create a function which takes as an input argument a vector and returns its number of non-NaN elements (HINT:
check Matlab help for isnan and length).
4. Compute the same percentage for each year of data and store the results in a vector. Remark: For some years,
there is one day less than for others. Use the function defined in question 3 to take it into account while
calculating the percentages.
5. Plot the percentage of days where Q exceeds the threshold value as a function of the year (give a title and names
to the axis).
6. Compute finally the same percentage but for the whole data set (over the 100 years of data).
Computation of the discharge in the channel during floods
7. We consider here that the bankfull discharge, Qbankf ull , can be approximated by two times the mean value of
Q over the entire dataset. Calculate Qbankf ull .
8. Calculate the number of days where the bankfull discharge has been exceeded over the 100 years of data.
When the total discharge Q gets larger than the bankfull discharge Qbf , the river floods. In this case, the discharge
really flowing in the channel and contributing to the channel bed shear stress is less than the total discharge Q:
Qchannel = Qbankf ull + Qoverbchannel (see Figure 2.1 for definitions). The following question aims at calculating
Qchannel , which is the relevant discharge for sediment transport calculations.
9. Define a matrix Qchannel of the same size as Q. The elements of this matrix should be equal to those of Q if
Q < Qbf , and equal to Qbankf ull + Qoverbchannel otherwise.
To check your calculations, compare Qchannel and Q. Does it match your expectations?
HINT: Qoverbchannel can be expressed as a function Qoverbank , W and Wf p . To find this relation, relate the
cross-sectional area Aoverbchannel to Aoverbank . As the depth-averaged flow velocity is considered to be constant
51
Figure 2.1: Skematic view of the profile of the Rhine at Lobith during a flood event. Wf p is the total width of the
flood plain, W the channel width defined earlier. The total discharge Q in the river is equal to Qbankf ull + Qoverbank .
The part of the overbank discharge flowing over the channel is called Qoverbchannel .
52
over the cross-section, you can then obtain a relation between the discharges Qoverbchannel and Qoverbank .
Data: W = 550 m; Wf p = 3000 m;
Sediment transport calculations
Additional data: Channel slope S = 1.5 × 10−4 m/m; Median diameter of the sediment D50 = 2 × 10−3 m; Density
of the sediment ρs = 2650 kg/m3 ; Water density ρw = 1000 kg/m3 ;
10. Calculate the water level from the discharge using the QH (stage discharge) relation given below:
H = (384.73 ln(Q) − 1949.3)/100.
H is the elevation in meters above mean sea level and Q the total discharge.
11. Calculate the water depth h knowing that the average bed elevation is 3 meters above sea level.
12. Apply the Engelund Hansen (1967) relation to compute the sediment transport rate associated to the channel
discharge Qchannel . Use the functions defined earlier to calculate the sediment transport.
13. Calculate the sediment discharge Qs in m3 /s. Plot Qs as a function of Qchannel for the whole dataset. Choose
the best representation (logarithmic axis or not).
14. Calculate the average volume (m3 ) transported through Lobith per year. To do that, you should first compute
the total volume of sediment transported through Lobith for each year of the dataset, and then average these
values.
15. Determine what value of flow discharge transports this same average sediment volume. This is a representative
discharge for the morphodynamics and could be used in a numerical model for long-term morphological prediction instead of the full hydrograph, which is computationally more expensive.
16. Compare this representative discharge to the mean flow discharge at Lobith (calculated over the entire period
covered with the dataset).
53
PRODUCTS TO HAND IN
Exercises 1 to 5 and final assignment.
54
Chapter 3
Numerical modelling
The answers to the questions of this practical should be organised in a short word document, in which you should
include the relevant Matlab figures.
Remark: To save a figure, use the Matlab function print or click on ”File” and then ”Save as” in the figure window
(see Figure 1.2).
3.1
A very simple morphodynamic model
We consider in this section a one-dimensional river with an initial constant slope S0 = 1/10. We want to compute
the evolution of the river bed in space and time using very basic physical laws. Similar approaches are used in
cellular automats of rivers and in so-called reduced-complexity models (try a search in Google Scholar to see debates
and controversies). In this model, we do not consider the influence of the water motion on sediment transport: the
sediment transport rate qt is assumed to depend linearly on the local bed slope S. More precisely, we assume that
qt = SedP rop S, with SedP rop a proportionality coefficient.
55
To solve this problem numerically, it needs to be discretized in space and in time. The reach of the river under
consideration, of length L, is divided into Nx sub-reaches of length dx = L/Nx . These sub-reaches are separated by
Nx + 1 nodes, of coordinates xi = (i − 1)dx, with i integer varying between 1 and Nx + 1. To each of these nodes are
associated a bed elevation and a sediment transport rate called ηi and qti (i is the node number). The bed elevation
ηin at time tn is then used to compute the bed elevation ηin+1 at time tn+1 = tn + dt (concept of time discretization).
This should be done in two main steps (see also Fig. 3.1):
- the sediment transport rate qt is first computed as a function of the local slope at each node (proportionality
coefficient given in the Matlab script). HINT : the local slope, defined as minus the gradient of the bed elevation
∂η
(S = − ∂x
), can be computed using the built in function gradient (see Matlab help).
- to compute the bed evolution, the sediment conservation equation (also called Exner equation) is then used. It
can be written as:
∂qt
∂η
=−
,
∂t
∂x
that is to say, in a discretized form:
n
ηin+1 − ηin
∂qt
=
,
dt
∂x i
and therefore
ηin+1 = ηin −
∂qt
∂x
n
dt.
i
The gradient of sediment transport should be calculated and then used to update the bed elevation.
• Download the script file very_simplified_model.m. It contains the structure of the model, and some commands
to visualise the results.
• Add in the initialisation part of the script the definition of the bed elevation vector eta at t=0. The initial bed
elevation decreases linearly with a slope S0 . The bed elevation at x = L (downstream boundary) should be zero.
Data: L = 10 m and S0 = 1/10.
56
Figure 3.1: Organization of our very simple morphodynamic model (no computation of the flow)
• Complete the time-loop using the explanations given above.
• Run the model and analyse the results.
• Let’s now modify the upstream boundary condition. We set the sediment transport rate qt to zero at the
upstream boundary (node 1, x = 0). This means that there is no incoming sediment flux in our model. Run the
code and interpret the results.
• Create a vector which will contain the time-evolution of the bed elevation at the first node. This vector should
be filled element-by-element during the time loop. Plot this vector as a function of time. What kind of evolution
do you obtain?
57
3.2
Flow computation: backwater equations
To model the river in a more realistic way, it is necessary to compute the flow in the channel at each iteration, and
then use it to compute the sediment transport and the bed evolution (see Fig. 3.3). In this section, we present the
backwater equations, which will be used in Section 3.3 to compute the flow in a more complex morphodynamic model.
More precisely, the flow in the channel will be computed using the backwater formulation for steady, gradually
varied flow under the constraint of constant water discharge per unit width qw and an imposed water elevation ξd at
the downstream boundary. Let H denote depth, ξ the water surface elevation, η the bed elevation (ξ = H + η, see
also Fig. 3.2), kc the composite roughness height and S the bed slope.
The backwater equation writes:
dH
S − Cf Fr 2
(3.1)
=
dx
1 − Fr 2
where S is the local bed slope and Fr the Froude number defined as Fr 2 =
coefficient Cf is assumed to follow a Manning-Strickler resistance relation:
Cf −1/2 = αr
H
kc
U2
gH
with U =
qw
H .
The bed friction
1/6
(3.2)
with αr = 8.1.
How to solve this equation
depth at the ith node):
The backwater equation is first re-written in the discrete form (Hi being the water
Hi+1 − Hi
= Fback (H),
dx
with Fback a function defined as
Fback (H) =
S − Cf Fr 2
1 − Fr 2
58
(3.3)
that is to say, using the definitions of Cf and Fr ,
Fback (H) =
ηi+1 −ηi
dx
−
1
αr 2
1−
H
kc
−1/3
2
qw
gH 3
2
qw
gH 3
.
(3.4)
Rearranging equation (3.3), we get:
Hi = Hi+1 − Fback (H)dx.
(3.5)
In practice, the resolution is performed in two steps in the code using a predictor-corrector scheme:
- a first estimation of the water depth Hpred i is computed:
Hpred i = Hi+1 − Fback (Hi+1 ) dx
- Hpred i is then used to obtain a more precise prediction of the water depth (corrector step):
Hi = Hi+1 −
1
(Fback (Hpred i ) + Fback (Hi+1 )) .
2
The water depth Hi at the node i can be therefore computed using the water depth Hi+1 one node downstream.
The equations are solved starting from the downstream boundary (node Nx + 1), where the water depth H(Nx + 1)
is known:
H(Nx + 1)
=
ξ(Nx + 1) − η(Nx + 1)
= ξd − η(Nx + 1),
with ξd the water surface elevation imposed by the user.
• Download the Matlab program called backwater_equation.m as well as the function F_back.
59
(3.6)
(3.7)
• Identify the different steps in the computation of the water depth. Why is the loop starting at Nx + 1 and not
at 1 as usual?
• Run the script and analyse the output.
• The normal water depth Hn and critical depth Hc are computed at the beginning of the script (for explanations
about these depths, look at the part relative to the backwater equations in the powerpoint presentation Chapter
5 of Gary Parker’s e-book). Add on the figure two lines corresponding to H = Hc and H = Hn. Modify the
legend.
• Run the code for different values of ξd . Different behaviours can be observed depending on the value of Hd =
ξd +η(N x+1). Make a representative figure for each type of observed behaviour (to add to your word document)
and comment on your findings. Discuss in particular what is happening for values less than the critical depth.
• Calculate now a second vector H2, which contains the values of water depth obtained without using a corrector
step. Plot the results with and without using a corrector on the figure and compare them. This comparison
should be performed for a fine spatial grid (the one given initially for instance) and a very large one. Do you
understand why predictor-corrector schemes are used?
3.3
Morphodynamic model with flow computation
The model studied in this section is based on the same formulation as the spreadsheet RTe-bookAdDegBW.xls from
Gary Parker e-book. It allows to compute the variations of the bed level η in space and time for different upstream
and downstream conditions (variables in red in Fig. 3.2). This resolution is performed in 3 steps, summarized in Fig.
3.3:
- the first part computes the water depth and flow velocity, based on the backwater equations presented in Section
3.2,
- the sediment transport is then computed using the flow characteristics determined in the previous step,
60
- finally, Exner equations are used to calculate the bed elevation from the sediment transport fluxes.
• The folder full model contains the different m-files needed for this assignment.
• A detailed description of the model is given in chapter 20 of the e-book (a word file, not the ppt).Download this
document and use it to identify the different parts of the script model.m. Look through the different functions
called in this script and try to understand what is going on.
• Run the script. Why doesn’t the bed change?
• Implement the sediment transport block in the model. Change the definition of qt and remove the command
qt_u=0;, which sets the sediment transport at the upstream boundary to zero. Calculate the shear stress from
the flow velocity and friction; not from the depth slope product.
• Run the script again, analyse the output. Where is sediment deposited and why?
• What is happening when you raise the sea water level while keeping the sediment input constant? Why?
• What is happening if you decrease the sediment input? Why?
• In the Rhine, both happened at the same time. What happens then?
• Play with the model!
61
REMARK Your model might produce incorrect results for some values of grid size and time step. More precisely, if
the hump on the bed migrates more than one grid step dx during one time step dt, the model becomes unstable. This
gives the following stability condition (also known as the Courant-Friedrichs-Lewy (CFL) condition):
c
dt
≤ 1,
dx
where c is the characteristic celerity of a disturbance in bed level. In practice, it means that if you increase significantly
dt (to decrease your computation time for instance) you may have to increase dx to keep the model stable. However,
keep in mind that increasing dx decreases the spacial resolution of your model.
62
Figure 3.2: Initial condition and variable definitions
63
Figure 3.3: Organization of the morphodynamic model introduced in Section 3.3
64
Chapter 4
Creative Matlab exercise
The last (creative) exercise gives you some freedom of subject and method and allows you to distinguish yourself (earn
higher grades). Furthermore, you can perhaps use the resulting code in the delta project later in this course. The idea
is that you conceive your own researchable (quantitative) question, do the research using matlab, and present it in a
powerpoint to the group.
Steps of the empirical research cycle, that you will go through and that you will include in the five minute powerpoint
presentation, are:
1. a researchable question
2. hypothesis (based on literature, provide references)
3. method
4. results
65
5. discussion (in view of hypothesis)
6. conclusions
7. new research questions (at least one).
Before you start on the research propose your research question to the supervisors early Monday afternoon (22
September), or send your research question by email to m.g.kleinhans@uu.nl on Monday evening 22 September the
latest. Example research questions have been mentioned during class.
The creative exercise will be presented using powerpoint in max 5 minutes per person (see schedule in study guide
for date, time and location). Bring your powerpoint file on a USB stick.
You may do the exercise on your own or in student pairs. If you work in pairs there should be two (related) research
questions and investigations, and two powerpoints of 5 minutes each because each should present his/her own work.
Indicate in the code (as comments) the contribution of each partner.
The individual creative exercise is worth 20% and the other three exercises are worth 5% each of the final grade.
Grades will be individual and will be based on the quality of all steps in the empirical cycle plus the presentation
quality plus, very importantly, creativity! Grades for having all steps done well and the presentation itself can be
troll (2), not yet good enough (4), just good enough (6), excellent (8) and outstanding (10) and the final grade is the
numeric average of these, where steps are whimsically weighed by staff. The Matlab code is graded on the same scale
separately on three criteria: does it work, does it do what it is supposed to do, and can we, experienced yet mortal
and impatient lecturers, understand it. Have fun!
66