# MATLAB Basics

*
Ville-Pekka Eronen,
Jarkko PeltomÃ¤ki
Department of Mathematics and Statistics, University of Turku
Last updated: 31.8.2017*

This document contains an introduction to the MATLAB software (later simply Matlab). The aim of the material is to give sufficient knowledge on the software to attend to Tero Aittokallio's course Scientific Computing; it is not meant to be an introduction to programming or problem solving with Matlab. The material is mainly for self-study, experimenting is arguably one of the best ways to learn specific software and programming.

## Contents

- Short History of Matlab
- Links to Other Materials
- The User Interface and the Interpreter
- m Files and On Writing Own Functions
- On Debugging
- Importing Data to Matlab
- Exporting Data from Matlab
- 2D Graphics
- 3D Graphics
- Other Types of Graphics
- Linear Algebra
- Numerical Integration and Differentiation
- Solving Differential Equations Numerically
- On Minimizing a Function
- Finding Roots
- Data Fitting
- On Interpolation

## Short History of Matlab

Matlab (stands for matrix laboratory) is a proprietary software that is meant for matrix-based numerical computing. It was developed in the 1970's by Cleve Moler who wanted his students too access Fortran-based linear algebra libraries without knowing Fortran itself. Later Moler commercialized the software and cofounded the MathWorks software company in 1984. The company actively develops Matlab even today: the 50th release version, dubbed R2017a, was released in March 2017. What makes Matlab special is that it is a pure matrix language; even real-valued variables are strictly 1-by-1 matrices. A significant portion of numerical tasks is naturally interpreted to the language of matrices, so Matlab, being based on the optimized LINPACK linear algebra library, is very efficient and popular numerical computation software. Its matrix operations are notably faster than, for instance, Mathematica's. Unlike Mathematica, Matlab cannot natively manipulate symbolic computations, even though the feature is available as an additional toolbox. Matlab is a proprietary software, and its licence is quite expensive. Its is, however, possible for students to get a cheaper licence. It is possible to use free software, such as Scilab or GNU Octave, to learn the principles of Matlab. These programs are not complete Matlab clones, but the basic programming covered in this document is very similar. Matlab is available for the usual platforms: Unix, Mac OS X, and Windows.

## Links to Other Materials

- Numerical Computing with MATLAB
- G. Recktenwald's course material
- M. Vuorinen's course material (in Finnish; the page contains additional links)
- Aalto University's short introduction to Matlab (in Finnish; the page contains additional links)

## The User Interface and the Interpreter

By starting the Matlab software, the following window opens:

In the middle of the above picture, there is the interpreter (command window) where you can interactively run Matlab commands. Displayed on the left, there is the current folder. On the right, there is the workspace that records the current variables, i.e., the state. Below it, we see the command history.

The interpreter's prompt >> is distinct from the printed lines not displaying it. A command is executed by pressing the ENTER button. For example:

x = 5 + 1

x = 6

phi = (1+sqrt(5))/2

phi = 1.6180

We observe that Matlab prints output after every command. This interactivity is often unnecessary, so printing can be suppressed by ending a command with a semicolon ;. Using a semicolon or a comma, several commands can be written on a single line, if desired. The percent sign % starts a comment. The syntax of Matlab is straightforward to anyone who has used some C-like programming language. We will not dwell on different aspects of the syntax here, it should become evident from the examples and the experiments by the reader. However, we will say this: Matlab is a dynamically typed language. Notice that by creating a new variable, it is possible to overwrite a built-in Matlab command. Defining a variable named `sin` prevents the usage of the function `sin`; this should be avoided.

More information on a Matlab command can be obtained by writing `help command` or by writing the command into the search box (top right in the above picture). By simply writing `help`, a list of available commands organized by topic is printed. The command `lookfor word` finds all commands that contain the given word on the first line of their help message (the first line contains a description for the whole function). In addition, the command `lookfor word -all` lists all commands whose whole help message contains the given word. The web page of Mathworks also provides a wealth of information.

**Exercise.** Make yourself comfortable with Matlab's user interface and try basic arithmetic operations in the interpreter. Use the `help` command to study the command `format`, and print the golden ratio with greater accuracy. Does the amount of displayed digits affect the precision used in the computations? Find out what `format compact` does.

By writing the above example in the interpreter, we see that the defined variables `x` and `phi` will stay in memory and are displayed in the workspace. The command `who` lists all defined variables, `clear` deletes all of them, and `clear x` removes the variable `x`. Try also the command `whos`. Like in a calculator, Matlab offers the special variable `ans` that contains the result of the previous command. The `clc` clears the text displayed in the interpreter.

The command `disp` is used for printing. For example:

disp(phi)

1.6180

Printing strings works as follows (in Matlab a string is indicated by single quotes):

```
disp('Hello world!')
```

Hello world!

In Matlab, strings are arrays of characters, so more complex output can be generated as follows (cf. the creation of vectors and matrices below).

disp(['The value of the variable phi is ' num2str(phi) '.'])

The value of the variable phi is 1.618.

The elementary functions, such as `abs`, `sin`, `cos`, `tan`, `exp`, `log`, `sinh`, are built in Matlab. More information can be obtained by typing `help elfun`. Special functions, like Bessel functions, are found by entering `help specfun`.

**Exercise.** Visit the Matlab documentation on page Complex Numbers and find out how complex numbers work in Matlab. Try calling elementary functions with complex parameters.

### Matrices

Like mentioned already, Matlab is a matrix language. The easiest way to define a matrix is to give its elements by rows:

A = [1 2 3; 4 5 6; 7 8 9]

A = 1 2 3 4 5 6 7 8 9

The definition can be split onto several rows. Then the semicolon is not needed.

A = [1 2 3 4 5 6]

A = 1 2 3 4 5 6

By writing three dots, a long row can be split onto two separate lines:

```
A = [1 2 3 4 5 6 ...
7 8 9
10 11 12 13 14 15 16 17 18]
```

A = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

The elements of a matrix can be separated by comma instead of a space. Naturally the elements can be defined as results of arithmetic operations.

A = [1, 2*sin(pi/4), exp(2)]

A = 1.0000 1.4142 7.3891

The command `length` returns the length of a vector:

length(A)

ans = 3

On the other hand, the dimension of a matrix is given by the command `size`, try it! Matlab has several built-in commands for matrix generation. For example, the command `eye` returns an identity matrix:

eye(3)

ans = 1 0 0 0 1 0 0 0 1

**Exercise.** Use the `help` command, and get to know the commands `rand`, `ones`, and `zeros`.

The usual arithmetic operations +, -, *, /, \, ^, ', addition, subtraction, multiplication, division from right and left, exponentiation, transposition, are available for matrices. The programmer needs to make sure that the operated matrices are of the appropriate size. It is easier to get to know the operations by experimenting; we shall return to the division later when we discuss solving a system of equations. It is worth mentioning here that adding or subtracting a scalar has special syntax:

A = eye(3) A - 1

A = 1 0 0 0 1 0 0 0 1 ans = 0 -1 -1 -1 0 -1 -1 -1 0

We see that Matlab can automatically interpret the subtracted scalar to be multiplied by the matrix `ones(3)`. It is also worth mentioning, that ' is Hermitian transpose, that is, in addition to the usual transposition, also complex conjugation is applied. The usual transposition is available with the command `transpose` or with the notation .'.

Multiplication, exponentiation, and division also have elementwise versions. So for example, .* computes the elementwise product of two matrices instead of the matrix product. For example:

A = [1 2; 3 4], B = [1 1; 2 2] A .* B

A = 1 2 3 4 B = 1 1 2 2 ans = 1 2 6 8

Many Matlab functions, such as the elementary functions mentioned above, operate elementwise. For instance, the function `sin` computes the sine of every element of a matrix:

sin(eye(2))

ans = 0.8415 0 0 0.8415

### Colon Notation and Indexing

It is often necessary to produce vectors where the distance between two consecutive elements is constant. This is easily accomplished with the colon notation:

start:step:end

This produces a vector containing the numbers `start`, `start` + `step`, `start` + 2* `step`, ... that are smaller than or equal to `end`. For example:

1:0.1:2 1:0.3:2

ans = Columns 1 through 7 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 Columns 8 through 11 1.7000 1.8000 1.9000 2.0000 ans = 1.0000 1.3000 1.6000 1.9000

If `step` is omitted, then its value is assumed to be `1`. The following code produces the same result using the `for`-loop. This method should not be used, this is only for introducting the syntax of the `for`-loop and for showing how to add elements to vectors.

X = []; s = 1; e = 2; t = 0.1; for n = 0:(e-s)/t X = [X s + n*t]; end X

X = Columns 1 through 7 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 Columns 8 through 11 1.7000 1.8000 1.9000 2.0000

**Exercise.** Figure out how the command `linspace` works and how it is different from the colon notation. The command `logspace` might also be useful later.

The next code shows how to select specific elements, rows, and columns from a matrix.

A = [1 2 3 4; 5 6 7 8; 9 10 11 12] A(2,3) % The element at (2,3). A(8) % The element at (2,3) A(:,1) % 1st column A(1,:) % 1st row A(:) % The whole matrix as a column vector A(2:3,2:3) % Submatrix, rows 2-3, columns 2-3 A(:,[1,4]) % All rows, columns 1 and 4. B = logical([0 1 0 0; 1 0 1 0; 0 0 1 0]) % Creates a "logical" matrix which tells which elements of A are to be chosen. A(B) % Selects the elements of A whose corresponding positions in B contain 1, ordering is top to bottow, left to right.

A = 1 2 3 4 5 6 7 8 9 10 11 12 ans = 7 ans = 7 ans = 1 5 9 ans = 1 2 3 4 ans = 1 5 9 2 6 10 3 7 11 4 8 12 ans = 6 7 10 11 ans = 1 4 5 8 9 12 B = 0 1 0 0 1 0 1 0 0 0 1 0 ans = 5 2 7 11

A matrix can be analogously modified:

A(2,2:3) = [0 0] % Change the columns 2 and 3 on row 2. A(:,4) = [12; 8; 4] % Change 4th column.

A = 1 2 3 4 5 0 0 8 9 10 11 12 A = 1 2 3 12 5 0 0 8 9 10 11 4

**Exercise.** Define some matrix and add its first row multiplied by 3 to its second row. The colon notation thus allows to implement the elementary row operations for matrices. It is thus quite straightforward to program a program that puts a matrix into row echelon form.

### On Strings

Strings can be searched and edited with regular expressions. The next example finds the nonoverlapping occurrences of the string `aa` in a longer string.

regexp('abbabaabaaa', 'aa')

ans = 6 9

Naturally it is possible to write proper regular expressions. The next example finds "greedily" all occurrences of words that begin with k and end with o from a longer string.

s = 'koko kokoa koko kokko kokoon'; regexp(s, 'k[\w]*o') % returns the starting positions of the matching occurrences regexp(s, 'k[\w]*o', 'match') % returns the matching words themselves

ans = 1 6 12 17 23 ans = 'koko' 'koko' 'koko' 'kokko' 'kokoo'

Substrings are replaced as follows:

regexprep(s, 'k', 'o') % change letters k to letters o

ans = oooo ooooa oooo ooooo ooooon

Read more on regular expressions on the page Match regular expression (case sensitive).

## m Files and On Writing Own Functions

### Scripts

Often there is a need to save the typed commands to a file after an interactive session. This is accomplished by saving the code into a text file with extension .m. The code can be written with any editor, but it is very convenient to do it using Matlab's built-in editor. The editor can be opened by selecting `New -> Script` from the user interface. The same can be accomplished by typing `edit filename` into the interpreter. Let us copy the following code to the editor and save it as the file `testi.m` into the current folder (so that Matlab finds it, this is important).

x = randn(1, 100); m = mean(x) s = std(x)

After this, write `testi` in the prompt (you can also click `Run` in the editor).

testi

m = 0.0312 s = 0.9653

The result is thus the same as if the contents of the file `testi.m` were written in the interpreter. Notice that exactly the same output was produced and after the execution the defined variables `x`, `m`, and `s` stayed in the workspace. By using m files, it is thus possible to initialize variables with proper values in the beginning of an interactive session, for instance. It is possible in an m file to call the code of another m file in the exact same way. Remember that the contents of a m file can be commented using the comment sign %; this sould always be done. If you select some lines of code in the editor, their contents can be executed by pressing the F9 button.

An m file prepared in this fashion is called a script, and it cannot take any parameters when called. The next subsection tells that m files may contain functions, which naturally can take parameters.

Matlab can be used to produce documents that contain code, text, and other markup and that can be compiled into a HTML or PDF file. This document itself is a single m file that is compiled into a HTML file using the `publish` command. More information on this on the page Publishing Markup.

### Function m Files

Function definitions are written to m files with the following syntax:

```
function y = avg(x)
y = sum(x)/length(x);
```

In order for Matlab to recognize the m file as a function definition, the keyword function needs to be found on the first line of the m file. The m file containing the function must have the same name as the function (here `avg.m`) so that Matlab finds the desired function. Notice that the commands within a function should end with a semicolon; otherwise lots of unnecessary output might be produced. Remember that, as is usual with common programming languages, variable definitions inside a function definition are local and are not available after the execution of the function.

In the above function, `y` is the value returned by the function, and it is set on the file line of the function definition. The name of the function is `avg` and its parameter's name is `x`. A function may take several parameters: `avg(x,y)` etc.

A function may also return several values. Then the syntax is the following:

```
function [a,b] = foobar(x)
a = sum(x);
b = sum(x)/length(x);
```

This function is called as follows:

[a,b] = foobar(x)

If the user does not care for all the values, then only the first one can be saved:

a = foobar(x)

The preceding bracket notation can be analogously be used to pick one or more return values of a function. If a function returns nothing, the return value can be omitted:

function func(x) ...

Matlab supports optional parameters, that is, not all parameters need to be specified. Then it is up to the function if this raises an error or if it modifies the behavior of the function in some way (cf. More on Programming Functions below).

A single m file may contain the definition of two or more functions, but these extra functions are local to the m file. Only the file's "main" function may call these auxiliary functions. In this case, all function definitions must be terminated with the keyword `end`; otherwise it would be unclear where each function definition ends.

### Anonymous Functions

Several Matlab functions take functions as parameters, like the function `fzero` that finds a root of a function. Let us suppose that the file `fun.m` has the following contents:

```
function y = fun(x)
y = 1./(exp(x) + sin(x)) - x.^2;
```

The function defined in the file thus computes the value of the function . Let us find its root in the interval :

fzero(@fun, [-0.5, 1])

ans = 0.6350

The significance of the @ sign is that it creates a reference (function handle) to the function `fun`. References can be saved into variables, so the following code works as expected:

f = @fun; fzero(f, [-0.5, 1])

ans = 0.6350

It is tedious to write a separate m file for each little function, so Matlab provides syntax for defining a so-called anonymous function. The following code does the same as the preceding code without the need to write a separate m file:

f = @(x) 1./(exp(x) + sin(x)) - x.^2; fzero(f, [-0.5, 1])

ans = 0.6350

An anonymous function can also depend on a parameter:

a = pi; f = @(x) 1./(exp(x) + sin(x)) - a*x.^2; fzero(f, [-0.5, 1])

ans = 0.4090

Read more on anonymous functions here.

### The Location of m Files

By default Matlab searches for scripts and functions from the directories listed by the `path` command. The first directory in the list is the current working directory. The function `addpath` can be used to add more directories to the list. The directories matter especially when several m files with the same name are specified. In this case, Matlab picks the first m file it encounters when searching the directories given by `addpath` in order. The command `which` tells which m file Matlab attempts to use. Note that an m file should not have the same file as a built-in Matlab command as the builtins are always called first.

### Matlab Toolboxes

Matlab provides several toolboxes. These are typically collections of functions, that is, collections of m files, but they can also contain applications and data.

Toolboxes available during a session are found with the command `ver`. By executing the command, we learn, for instance, that the "optimization toolbox" is available. Thus all functions provided by the toolbox are directly available. This page lists additional purchasable toolboxes. Matlab Central provides third-party toolboxes; consult their readme files to see how to install them. If a toolbox is simply a directory containing m files, it is sufficient to add it using the `addpath` command. If the directory `G:\ELBM` contains m files, they are made available as follows:

help ELBM % Matlab does not know this function. addpath('G:\ELBM') help ELBM % Now it does. path

### More on Programming Functions

Matlab naturally has the usual control structures `if`, `for`, `while`, etc. We will not look at these in this document; for more, see the page Control Flow. In some of control structures it is useful to use logical operations. These are found with command `help ops`. The implementation of several Matlab commands are viewable using the `type` command. Try `type('quad')`.

Next we introduce a function that computes an approximation of pi based on the series (very slow convergence). The behavior of this function depends on the number of parameters given; it also prints an error if a parameter is invalid.

function p = pi_approx(iterations) % If the user did not specify the number of iterations, do 1000 iterations. if nargin() ~= 1 iterations = 1000; end if iterations <= 0 error('The number of iterations must be a positive integer.') end s = 0; for n = 0:(iterations-1) if mod(n, 2) == 0 s = s + 1/(2*n + 1); else s = s - 1/(2*n + 1); end end p = 4*s

The command `nargout` returns the number of return values requested by the caller. This is sometimes good to know if some return value is expensive to compute; if it is not requested, it is useless to compute it.

Acquiring input from the user is performed with the command `input`. This command takes as a parameter a string to be displayed when the user is prompted; it returns the value inputted by the user. The function below is the above function modified; now it asks the user to input the number of iterations.

function p = pi_approx_inter() iterations=input('How many iterations? '); if ~isnumeric(iterations) error('The number of iterations must be numeric') end if iterations <= 0 error('The number of iterations must be a positive integer.') end s = 0; for n = 0:(iterations-1) if mod(n, 2) == 0 s = s + 1/(2*n + 1); else s = s - 1/(2*n + 1); end end p = 4*s

As Matlab is dynamically typed, the functions shold be programmed to accept different types of parameters. For instance, a function computing the value of a real function at a given point might as well accept a vector or matrix of points as its input. Accomplishing this in Matlab is often straightforward due to the elementwise operations on matrices. The next function computes the value of the sin function, using the Taylor series, for all elements of the input matrix (or vector, or number).

function y = sin_taylor(A, iterations) if nargin() ~= 2 iterations = 10; end if iterations <= 0 error('The number of iterations must be a positive integer.') end y = 0; for n = 0:(iterations-1) if mod(n, 2) == 0 y = y + (1/factorial(2*n+1))*A.^(2*n+1); else y = y - (1/factorial(2*n+1))*A.^(2*n+1); end end

With the default number of iterations, the result agrees with the output of Matlab's `sin` function.

sin_taylor([0 0.1; 0.5 1]), sin([0 0.1; 0.5 1])

ans = 0 0.0998 0.4794 0.8415 ans = 0 0.0998 0.4794 0.8415

**Exercise.** Program a function that takes as its input a square matrix `A` (check that the user does not give an inappropriate matrix) and computes the matrix with the series using as many terms as the user desires.

## On Debugging

If the code contains an error, Matlab prints a reason for it. Sometimes Matlab prints the line number where the error occured. Click the link provided by Matlab to get to edit the errorneous line. An error message can also be studied with the command `lasterror`.

A=zeros(3); A(1,4); sin('a'); fzero('a'); err = lasterror err.message err.stack err.identifier

The code can be debugged and single-stepped by using Matlab's debugger. The debugger is activated by adding a breakpoint into the code. This can be done either by clicking the dash next to a line number (it turns into a red ball) or by the `dbstop` command. Let us try this in the case of the following code from Matti Vuorinen's course (lecture notes p. 64):

dbstop in e310.m at 15 e310

This sets a breakpoint at line 15. When executed, Matlab goes into debugging mode, which is recognizable from the letter K next to the prompt. Matlab commands can be used normally in debug mode, and the code can be single-stepped. The next line is executed with the `dbstep` command. The command `dbcont` resumes the execution until the next breakpoint. The debug mode can be exited with the command `dbquit`; this terminates the execution of the code.

who n dbstep who dbstep who dbcont n dbquit n

By using the debugger, it is thus possible to check that the variables have the values they should have.

Sometimes valid code executes slowly. With the commands `tic` and `toc` you can measure how long the code spends time in a section of code.

tic for i=1:10000 v(i)=i; end; toc tic v=1:10000; toc

Elapsed time is 0.001604 seconds. Elapsed time is 0.000280 seconds.

More time-related functions are found by typing `help timefun`.

## Importing Data to Matlab

Matlab has many commands for importing date such as `load`, `fscanf`, `textscan`, `importdata`, and `xlsread`. More functions are found from the page Supported File Formats for Import and Export.

Let us first look at the command `load`. It can load mat files and ASCII files containing numerical values. Mat is Matlab's own format, which is used to save workspace variables with the `save` command.

m = 1; M = [1,2; 3,4]; save('test.mat', 'M'); % Save matrix M. who; % Look which variables occupy workspace. clear; who; % Now the workspace is empty. load('test.mat'); who; % M loaded.

If you do not specify the variables to be saved, the `save` command saves the state of the whole workspace. If no filename is specified, the default name `matlab.mat` will be used.

The file plotdata.txt contains two columns of numerical values. These can be loaded as a matrix by typing

```
M = load('plotdata.txt');
```

Occasionally the data does not consist of only numerical values. Then other functions can be used to import the data. Let us next look at importing xlsx data with the `xlsread` command.

The data consists of three sheets. In the first sheet there is a simple table. Let us first check how we can read this by typing `help xlsread` in Matlab.

Read the table with the following command

```
M = xlsread('xlsxtest.xlsx')
```

We observe that no strings are read and the numerical data is stored in the variable `M`. The strings can be saved to the variable R as follows:

```
[M,R,U] = xlsread('xlsxtest.xlsx')
```

The variable `U` contains "raw data". Its class is cell array and each cell includes one read element. In general a cell can contain any data type. The data contained in a cell can be extracted by using braces (look the exqample below) or by using the function `cell2mat`.

class(U) class(U(1,2)) class(U{1,2}) N=cell2mat(U) N=cell2mat(U(2:4,2:5)) sum(sum(abs(N-M)))

Let us move on to the next sheet. In this data there are two missing values and one NA value. In addition, there is a new column of strings. Let see how `xlsread` will deal with this.

```
[M,R] = xlsread('xlsxtest.xlsx',2)
```

We observe that the missing values and the NA value has been replaced by NaN (double), which stands for "not a number". Furthermore, the string column was correctly read in the variable R.

In the third sheet there are decimal commas. These are nontrivial in the sense that Matlab uses points to indicate decimals.

```
[M,R] = xlsread('xlsxtest.xlsx',3)
```

We observe that `xlsread` was able to change decimal commas to dots correctly.

Exercise: Try to read the following datas by using `importdata` function. (spoiler:problems ahead)

The import wizard (File->Import Data...) of Matlab uses the function `importdata`. As can be seen from the exercise it is not always straightforward to use and sometimes it is necessary to use other functions.

One easy way to import the datas in the exercise is to read them in excel, save an xlsx-file and use the function `xlsread`. Another way is to use the function `textscan` or `fscanf`. For example, the data in the file "dotdata" can be read as follows:

fo=fopen('G:\matlab\dotdata.txt'); % open the file. header=fscanf(fo,'%s %s %s %s',4); % variable for the header row. m=textscan(fo,'%s %f %f %f %f'); % The rest of the rows are of the same form. fclose(fo); % close the file. M=[m{2},m{3},m{4},m{5}]; % Creating the numeric matrix.

We can also change the decimal commas into decimal dots in the file "commadata.txt" with `fscanf`.

fo=fopen('G:\matlab\commadata.txt'); commadata=fscanf(fo,'%c') % '%c' preserves empty spaces fclose(fo); dotdata=regexprep(commadata,',','.') % Replace commas by dots in the string pilkkudata. fc=fopen('G:\matlab\commadatamod.txt','w'); % open and, if necessary, create the file "pilkkudatamuok". fprintf(fc,'%c',dotdata); fclose(fc);

More help for importing data can be found on page https://mathworks.com/help/matlab/import_export/recommended-methods-for-importing-data.html.

## Exporting Data from Matlab

For most of the commands to import data there are corresponding commands to export data. An example of this correspondence are the `save` and `write` commands. Detailed information on commands to export data is found from the two preceding links. Let us as an example save the matrix

M = [1,2,3; 4,5,6; 7,8,9]

as an excel file. The right syntax is found with the command `help xlswrite`.

xlswrite('xlswritetest.xlsx', M);

**Exercise.** Using the command `fprintf`, save the above matrix `M` into a text file in such a way that the columns have headers "one", "two", and "three".

## 2D Graphics

### Plot

The command `plot` draws the given points on the plane and connects the points by a line (by default). For example:

x = [0 1 2]; y = [1 -1 0.5]; plot(x, y);

Matlab sets the x-axis and y-axis limits based on the plotted data. This can be adjusted as follows:

```
axis([-1,1,-2,2]) % x-axis -1 - 1, y-axis -2 - 2
```

The scale of the axes is controlled with the commands `axis equal` (same scale on all axis) and `axis square` (axes are equally long on the screen). Resetting is done by typing `axis normal`. Notice that these commands need to be given after the `plot` command (Matlab cannot know the data before it has drawn it). The command `grid` displays a grid on the plot. Try it!

The commands `axis` and `grid` implicitly modify the current figure. Matlab may have several figures open at the same time. A new figure is opened with the command `figure`; they are referred to by a running sequence of integers. The command `figure(n)` changes the current figure, which means that the commands `axis`, `grid`, etc. operate on another figure. In the next example, we draw the sine and cosine functions in different figures with different scales.

clear; x = linspace(0, 2*pi); y = sin(x); z = cos(x); figure; % Create new figure, we are now modifying implicitly figure no. 1. plot(x, y); axis equal; figure; % Modifying implicitly figure no. 2. plot(x, z, 'x'); % Draw only data points, not connecting lines. axis([0, 2*pi, -2, 2]); figure(1); % Modifying again figure no. 1. grid; figure(2);

The figure can be saved in common formats with the command `save`. For instance, `print('image', '-dpdf')` saves the current figure as a PDF file (saved in the working directory as the file `image.pdf`). The parameter `-deps` saves the figure as EPS. Saving is also accomplished by selecting `File -> Save As` from the menu of the figure window.

It is possible to draw several plots in the same figure. Calling the command `plot` erases the previous picture. To prevent this, use the command `hold on`. Now all calls to `plot` draw on top of what was previously displayed until the command `hold off` is used. The `plot` command can handle several data sets at once:

plot(x1, y1, 'string1', x2, y2, 'string2', x3, y3, 'string3', ...)

The given strings (LineSpec) control the line color and shape. Try the values '--or' and ':^g'. More information on this on the page Plot.

**Exercise.** Give a vector of complex numbers for the command `plot`. Using the above link, find out how to draw the data points without the connecting line segments.

### Fplot

The preceding `plot` command drew the data points given by two vectors and connected them by line segments. Indeed, this command merely draws the given data, which does not need to represent a function. If you want to plot a built-in Matlab function or a function you have written yourself, the command `fplot` can be used instead.

fplot(@(x) x.^2 .* sin(x), [0, 4*pi]);

Now Matlab automatically handles the discretization of the interval.

Plotting parametrized functions also works:

fplot(@(t) cos(3*t), @(t) sin(2*t));

Naturally, the usage of polar coordinates is possible:

t = linspace(0,2*pi); polar(t, sin(2*t).^2);

### Multiple Plots in the Same Figure

The command `subplot` divides the figure into multiple subplots. Calling `subplot(m, n, p)` divides the figure into m*n subplots that are arranged as a m*n grid. The parameter selects the pth subplot as the current figure (numbering is by rows). In the next example, we create two horizontally adjacent plots, draw to them, and also set titles for the subplots (study the command `title` in more details).

close all; % Close all figure windows. subplot(1, 2, 1); fplot(@(x) sin(x), [0, 2*pi]); title('Sine Curve'); subplot(1, 2, 2); fplot(@erf, [-1,1]); title('Something Else');

### Miscellanous

Here is one more example on adding text to plots, adjusting line width etc. Try it out!

figure() x = linspace(-2*pi, 2*pi); y1 = sin(x); y2 = cos(x); plot(x, y1, x, y2) title('Title'); legend('y = sin(x)', 'y = cos(x)'); xlabel('text below the x-axis alle'); ylabel('text next to the y-axis'); axis([-5, 5, -1.5, 1.5]); set(gca, 'FontSize', 15); % gca refers to the axes of the current figure set(gca, 'FontWeight', 'bold');

## 3D Graphics

The command `plot3` is analogous to the command `plot`:

x = [-1:0.5:1]; y = [2 2 2 3 3]; z = [1:5]; plot3(x, y, z); xlabel('x'); ylabel('y'); zlabel('z'); grid; axis([-2,2,2,3,0,5]);

Plotting a space curve also works in the same way:

fplot3(@(t) sin(t), @(t) cos(t), @(t) t, [0, 8*pi]);

Let us next consider drawing a surface (that is, a plot of a function from the plane to the line). This is accomplished with the `surf` command. Let us first look at the command `meshgrid`. This command takes as its parameters two vectors `x` and `y` and forms two matrices `X` and `Y`. These matrices correspond to the lattice points defined by the vectors `x` and `y`.

x = [0:4]; y = x; [X, Y] = meshgrid(x, y)

X = 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 Y = 0 0 0 0 0 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4

Now the elements `X(i,j)` and `Y(i,j)` at `(i,j)` form the lattice point determined by `x(j)` and `y(i)`. It is easy to evaluate a function in these lattice points. When we provide the parameters `X`, `Y`, and the evaluated values to the command `surf`, Matlab draws the plane pieces connected by the lattice points. Let us as an example draw the plot of the function .

```
x = linspace(-10,10); y = x; % the limits of x-axis and y-axis are the same
[X, Y] = meshgrid(x, y);
d = sqrt(X.^2 + Y.^2);
Z = sin(d)./d;
surf(X, Y, Z);
```

**Exercise.** Draw the preceding function (or your personal favorite function). Open the Tools menu of the figure window, and try out especially the options Rotate 3D and Data Cursor.

**Exercise.** Use the preceding values, and call `contour(X, Y, Z)`. What is the meaning of the displayed plot? Consider also the command `surfc`.

## Other Types of Graphics

Here we provide a cursory glance on how to visualize data as histograms, pie diagrams etc. See the page Pie Charts, Bar Plots, and Histograms for completeness.

Let us generate a thousand normally distributed numbers, and draw a histogram.

x = randn(1000, 1); histogram(x);

Let us draw a pie diagram out of the numbers given by the vector `x` (the areas rotate counter clockwise starting at noon).

x = [0.3 5 2.3 0.88 1]; pie(x)

## Linear Algebra

Let us consider solving a linear system of equations. Here we will be using Matlab's operations / and \, that is, division from the right and the left. The system is solved by multiplying from the left with the inverse of A. In Matlab, this is accomplished by "dividing" by A from the left:

A = [1 2 3; 0 -4 1; 0 3 -1]; y = [1 2 3]'; % column vector x = A\y A*x % check for validity

x = 65 -5 -18 ans = 1 2 3

Naturally, right division solves the system . This division operation in Matlab is not the same as multiplying with the inverse, but it has to do with solving the system numerically, cf. `mldivide` and `mrdivide`. The division operation gives an answer ever in the case that the system is overdetermined:

A = [1 2 3; 0 -4 1; 0 3 -1; 0 2 1]; y = [1 2 3 4]'; x = A\y

x = -6.7419 0.6774 2.1290

Here Matlab finds the solution in the sense of least squares.

The coefficient matrix should preferably be invertible. Let us have look what happens if it is not.

```
A = [1 -2 4; 0 0 6; 0 0 1];
det(A) % the determinant of A is 0
y = [1 2 3]';
A\y
```

ans = 0 Warning: Matrix is singular to working precision. ans = -Inf -Inf 3

As we can see from the output, Matlab gives a nonsensible result, which hints that the system does not have a solution. If the system had a solution, it would belong to the column space of the matrix A. By adding the vector `y` as a column of the matrix A, we see that the rank increases:

rank(A) rank([A y])

ans = 2 ans = 3

Thus `y` does not belong to the column space of A, so the system is contradictory. If a system has a solution when is not invertible, then a solution can be found using the Moore-Penrose inverse:

A = [1 3 7; -1 4 4; 1 10 18]; y = [5 2 12]'; x = pinv(A)*y A*x

x = 0.3850 -0.1103 0.7066 ans = 5.0000 2.0000 12.0000

If an exact solution does not exist, the above method yields a solution in the sense of least squares. For more on solving systems of linear equations, see here.

The inverse of a matrix is given by the command `inv`:

A = [1 1; 0 2]; B = [1 0; 5 2]; B*inv(A) B/A % this also works ans*A % check for validity

ans = 1.0000 -0.5000 5.0000 -1.5000 ans = 1.0000 -0.5000 5.0000 -1.5000 ans = 1 0 5 2

Above we introduced the commands `det` and `rank` that compute the determinant and the rank of a matrix. Eigenvalues and eigenvectors are found with the command `eig`:

A = [0 -6 -1; 6 2 -16; -5 20 -10]; eig(A) % find only eigenvalues (faster) [V, D] = eig(A) % finds eigenvectors (the columns of the result) ja eigenvalues as a diagonal matrix V*D*inv(V) % this gives back the matrix A (diagonal representation)

ans = -3.0710 + 0.0000i -2.4645 +17.6008i -2.4645 -17.6008i V = -0.8326 + 0.0000i 0.2003 - 0.1394i 0.2003 + 0.1394i -0.3553 + 0.0000i -0.2110 - 0.6447i -0.2110 + 0.6447i -0.4248 + 0.0000i -0.6930 + 0.0000i -0.6930 + 0.0000i D = -3.0710 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -2.4645 +17.6008i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -2.4645 -17.6008i ans = 0.0000 + 0.0000i -6.0000 + 0.0000i -1.0000 + 0.0000i 6.0000 - 0.0000i 2.0000 + 0.0000i -16.0000 + 0.0000i -5.0000 - 0.0000i 20.0000 + 0.0000i -10.0000 + 0.0000i

If the matrix is non diagonalizable, its Scur decomposition (unitary conjugation to an upper triangular matrix) is obtained with the command `schur`:

A = [13 8 8; -1 7 -2; -1 -2 7]; [U, T] = schur(A)

U = 0.9428 -0.3333 0 -0.2357 -0.6667 -0.7071 -0.2357 -0.6667 0.7071 T = 9.0000 -12.7279 -0.0000 0 9.0000 -0.0000 0 0 9.0000

Other typical decompositions can be found on the page Matrix Decomposition.

## Numerical Integration and Differentiation

The command `diff` returns the differences of the elements of a vector:

x = [4 3 0 -1 5]; diff(x)

ans = -1 -3 -1 6

Thus this command can be used to numerically differentiate (naively):

step = 0.1; % step size x = 0:step:5; f = x.^2 + 3*x; df = diff(f)/step; % difference quatient plot(x, f, x(1:length(df)), df); % the vector df is one element shorter

It is more straightforward to use the `gradient` command that can also compute the gradient of a function of several variables.

gradient(f, step);

Numerical integration is done with the `integrate` command. This command is only available after Matlab version R2012A. With an old version, the `trapz` command below needs to be used.

integral(@(x) x.^2, 0, 4)

ans = 21.3333

**Exercise.** Look up the help on the command `integral2` and integrate the function over the area. .

The command `integral` can only integrate a Matlab function, it does not accept a data vector as an input. The integral of an unknown function represented only as a data vector can be computed with the command `trapz` (trapezoidal rule). Let us apply the command to compute the integral of the function over the interval .

```
step = 0.1;
x = 0:step:4;
f = x.^2; % generate the data
trapz(x, f)
```

ans = 21.3400

## Solving Differential Equations Numerically

Next we consider how to use Matlab to solve differential equations or systems of differential equations. The command used in this document to solve differential equations is `ode45`, which is recommended as the default option. Matlab has several other solvers for different uses, cf. Choose and ODE Solver.

Matlab only solves ODE's of the form , higher-order equations need to by transformed into this form (cf. lecture notes on differential equations or the above link). The command `ode45` takes as its parametrs the function on the right side of the above equation, the integration interval, and an initial value. If we want to find a numerical solution to the ODE in the time interval with initial value , then we need to write

[t, y] = ode45(@(t, y) 3*t*y, [0 1], 1);

The command returns two vectors: the solution `y` as a function of the time `t`. The integral curve can thus be immediately drawn:

```
plot(t, y, '-o')
```

The general solution to the above ODE is . Thus in this case Matlab's solution looks correct. The command `ode45` uses the Runge-Kutta method to solve the ODE.

**Exercise.** Draw the integral curves of the solutions to the ODE with different values for the parameters and . Use an anonymous function that depends on a parameter.

Let us consider next solving a system of differential equations and use as an example solving the second order ODE . Using the substitution , , we see that we need to solve the system , . For this system, let us write the following fuction (an anonymous function would also do):

function r = dy_en(t, y) % Observe that the parameter y is a vector of two elements because the % system has two equations. The command od45 requires that the return % value is a column vector, so we cannot simply write % 'function [y1, y2] = dy(t, y)' as the function definition. r = zeros(2, 1); r(1) = y(2); r(2) = -sin(y(1)); end

The solution is obtained analogously:

[t, y] = ode45(@dy_en, [0 5], [0 1]); plot(t, y(:,1), t, y(:,2));

## On Minimizing a Function

Commands suitable for minimizing a function can be found from the optimization toolbox. Its documentation is available here or by typing "optimization toolbox" into the help.

Let us seek a minimum of the function . At least the commands `fminsearch` and `fminbnd` of the toolbox are suitable for the task. Documentation tells us that the syntax for `fminsearch` is the following:

[x, fval, exitflag, output] = fminsearch(fun, x0, options);

The first parameter is the target function, and the second one is the initial search point. Using the parameter `options`, it is possible to control how much information is printed.

opt = optimset('Display', 'iter'); % Set options to display informations on iterations. sincossum = @(x) sin(x)+cos(x); [x, fval, eflag, output] = fminsearch(sincossum, 0, opt); sincossum(x) output.funcCount

Iteration Func-count min f(x) Procedure 0 1 1 1 2 1 initial simplex 2 4 0.9995 expand 3 6 0.998499 expand 4 8 0.996494 expand 5 10 0.992472 expand 6 12 0.98438 expand 7 14 0.968009 expand 8 16 0.934527 expand 9 18 0.864728 expand 10 20 0.714808 expand 11 22 0.382525 expand 12 24 -0.333554 expand 13 26 -1.34737 expand 14 28 -1.38509 contract outside 15 30 -1.41225 contract inside 16 32 -1.41225 contract inside 17 34 -1.41412 contract inside 18 36 -1.41412 contract inside 19 38 -1.4142 contract inside 20 40 -1.41421 contract inside 21 42 -1.41421 contract inside 22 44 -1.41421 contract inside 23 46 -1.41421 contract inside 24 48 -1.41421 contract inside 25 50 -1.41421 contract inside 26 52 -1.41421 contract inside 27 54 -1.41421 contract inside Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 ans = -1.4142 ans = 54

In the preceding example, the class of the variables `opt` and `output` is `struct`. A struct is essentialy a named list of elements. Above `output.funcCount` prints those elements of the variable `output` that are under the name `funcCount`. In this particular example, there was only one element, but in general there can be several. In the general case, referring to a single value works as below.

structure = struct('names', {'x','y';'z','w'}, 'values', {1,2;3,4}); class(structure) structure.names structure(2).names

ans = struct ans = x ans = z ans = y ans = w ans = z

The command `fminbound` is used almost in the same way, this time instead of specifying an initial point, we need to specify an interval the minimum is searched for.

[x, fval, eflag, output] = fminbnd(sincossum, 0,2*pi, opt); sincossum(x) output.funcCount

Func-count x f(x) Procedure 1 2.39996 -0.0618786 initial 2 3.88322 -1.41286 golden 3 4.79993 -0.908745 golden 4 3.88982 -1.41324 parabolic 5 3.92941 -1.41421 parabolic 6 3.927 -1.41421 parabolic 7 3.92696 -1.41421 parabolic 8 3.92703 -1.41421 parabolic Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 ans = -1.4142 ans = 8

## Finding Roots

Finding a root of a one-variable function is accomplished with the `fzero` command. Its parameters are a reference to a function and an initial point. This command was earlier introduced in the section on anonymous functions.

```
fsini = @(x) sin(x);
[x, fval, ef, output] = fzero(fsini, 2)
fsini(x) % evaluate and check
```

x = 3.1416 fval = 1.2246e-16 ef = 1 output = intervaliterations: 10 iterations: 6 funcCount: 27 algorithm: 'bisection, interpolation' message: 'Zero found in the interval [0.72, 3.28]' ans = 1.2246e-16

Instead of specifying an initial point, an interval containing the root can be specified; see help.

The command `fzero` seeks only one root. All roots of a polynomial function can be found with the command `roots`.

a = [1, 2, 3, 4, 5, 6]; % specify the coefficients, x^5 + 2*x^4 + ... + 6 x = roots(a) polyval(a, x) % evaluate the polynomial at the found points

x = 0.5517 + 1.2533i 0.5517 - 1.2533i -1.4918 + 0.0000i -0.8058 + 1.2229i -0.8058 - 1.2229i ans = 1.0e-13 * 0.3464 + 0.1821i 0.3464 - 0.1821i -0.1332 + 0.0000i 0.1421 + 0.1998i 0.1421 - 0.1998i

**Exercise.** Find some root of the above polynomial using the command `fzero`.

The command `fsolve` is used to find a common root of several functions, that is, it can be used to solve a nonlinear system of equations.

nlineq = @(x) [sum(x.^2)-4, [-1,2]*transpose(x.^2)]; opt = optimset('Display', 'iter'); x = fsolve(nlineq, [1, 1], opt) nlineq(x) % check

Norm of First-order Trust-region Iteration Func-count f(x) step optimality radius 0 3 5 6 1 1 6 0.929784 0.849837 4.99 1 2 9 0.00256326 0.189768 0.235 1 3 12 2.83323e-08 0.01091 0.000777 1 4 15 3.53384e-18 3.64426e-05 8.68e-09 1 Equation solved. fsolve completed because the vector of function values is near zero as measured by the default value of the function tolerance, and the problem appears regular as measured by the gradient. x = 1.6330 1.1547 ans = 1.0e-08 * 0.1329 -0.1329

## Data Fitting

Let us see next how to fit a curve into the given data. A good way to find suitable commands is to use a search engine or by typing `lookfor fit` to list all commands whose description contains the word 'fit'. Out of the listed commands, at least `polyfit` ja `lsqcurvefit` seem suitable for curve fitting. Let us generate some data and try these commands out.

x = 5*rand(200, 1); % 200 uniformly distributed points in the interval [0,5]. y = 2*x.^2 - 5*x + 1 + 0.5*randn(200, 1); % 2nd order polynomial 2x^2-5x+1 and noise. pcoef = polyfit(x,y,2) % Try fitting 2nd order polynomial obtaining coefficients. xsorted = sort(x); % Order the points for drawing; Matlab draws the connecting % line segments in the order given by the data points. ysorted = polyval(pcoef, xsorted); % Evaluate the polynomial given by the coefficients at the sorted points. plot(xsorted, ysorted, '-', x, y, '.'); title('Fitting data with polyfit', 'Color', 'k','Fontsize', 15); plottext=sprintf('Fitted curve: \ny=%3.2fx^{2}%3.2fx+%3.2f', pcoef); text(0.5,5,plottext);

pcoef = 1.9925 -4.9474 0.9577

Let us do the same fitting with the command `lsqcurvefit`.

fitf = @(a,x) a(1)*x.^2 + a(2)*x + a(3); pcoef = lsqcurvefit(fitf, [0 0 0], x, y)

Local minimum found. Optimization completed because the size of the gradient is less than the default value of the optimality tolerance. pcoef = 1.9925 -4.9474 0.9577

We observe that this command found exactly the same polynomial as `polyfit`. This is expected as both commands use the method of least squares. Let us next fit a curve `a+b*exp(c*x)` to new data. The command `polyfit` cannot be used, but using the command `lsqcurvefit` is as straightforward as above. Let us also have a look what else the command `lsqcurvefit` has to offer.

x = 5*rand(200, 1); y = 2*exp(0.4*x) + 1 + 0.5*randn(200, 1); fitf = @(a, x) a(1) + a(2)*exp(a(3)*x); opt = optimset('display', 'iter'); [fcoef,resnorm,res,exitflag,output] = lsqcurvefit(fitf, [0,0,0], x, y, [], [], opt); xsorted = sort(x); ysorted = fcoef(1) + fcoef(2)*exp(fcoef(3) * xsorted); plot(xsorted, ysorted, '-', x, y, '.'); title('Fitting data with lsqcurvefit','Fontsize', 15); plottext=sprintf('Fitted curve: \ny=%3.2f+%3.2fe^{%3.2fx}', fcoef); text(0.5,8,plottext); figure resnorm sum(res.^2) plot(x, res, '.'); title('Residuals','Fontsize', 15); exitflag output.funcCount % the number of evaluations of the function output.algorithm % the algorithm used

Norm of First-order Iteration Func-count f(x) step optimality 0 4 13443.6 1.47e+03 1 8 2616.96 5.20255 3.75e+03 2 12 635.566 4.22004 1.06e+03 3 16 245.512 1.22828 4.51e+03 4 20 83.6733 1.0153 2.53e+03 5 24 39.9543 0.085077 104 6 28 39.891 0.00629217 0.00452 7 32 39.891 3.16258e-05 5.22e-06 Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the default value of the function tolerance. resnorm = 39.8910 ans = 39.8910 exitflag = 3 ans = 32 ans = trust-region-reflective

The command `lsqcurvefit` can also be used to fit functions of several variables. The convergence could be slow, and you might need to increase the number of iterations.

w = 5*rand(200,1); z = 5*rand(200,1); y = 2*exp(0.4*w) + exp(0.2*z) + 1 + 0.1*randn(200,1); fitf = @(a,x) a(1) + a(2)*exp(a(3)*x(:,1)) + a(4)*exp(a(5)*x(:,2)); fcoef = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y); opt = optimset('MaxFunEvals', 30000, 'MaxIter', 5000); [fcoef,resnorm,res,ef,output] = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y, [], [], opt); output.funcCount output.iterations

Local minimum found. Optimization completed because the size of the gradient is less than the default value of the optimality tolerance. Local minimum found. Optimization completed because the size of the gradient is less than the default value of the optimality tolerance. ans = 108 ans = 17

When modelling statistical data it may be worthwile to look at functions in `Statistics toolbox` such as `glmfit` and `mvregress`. Let us solve the first example of the section with the function `glmfit`.

x = 5*rand(200,1); y = 2*x.^2 - 5*x + 1 + 0.5*randn(200,1); X = [x, x.^2]; % Model matrix. Not necessary to write the intercept term. pcoef = glmfit(X, y, 'normal') xsorted = sort(x); ysorted = polyval(flipud(pcoef), xsorted); plot(xsorted, ysorted, '-', x, y, '.'); title('Fitting data with glmfit', 'Color', 'k','Fontsize', 15); plottext=sprintf('Fitted curve: \ny=%3.2fx^{2}%3.2fx+%3.2f', flipud(pcoef)); text(0.5,5,plottext);

pcoef = 1.1679 -5.1365 2.0287

## On Interpolation

Interpolation is related to curve fitting. The point of interpolation is to estimate the values of a function between two data points; the interpolated curve will go through all data points unlike in curve fitting. Consider the following data.

x = 1:8; y = exp(0.2*x) + 0.2*transpose(rand(8, 1));

A naive way to find an interpolating function is to fit a polynomial of sufficiently high degree. However, a high-degree polynomial may oscillate significantly between two data points. It is better to use splines.

coef = polyfit(x, y, 7); xx = 0.6:0.2:8.4; yy = polyval(coef, xx); plot(x, y, 'k.', xx, yy, 'b-'); hold; yy = interp1(x, y, xx, 'spline'); % found using lookfor interpolation plot(xx, yy, 'r-'); legend('data', 'polynomial', 'spline', 'Location', 'NorthWest'); title('Interpolation', 'Color', 'k','Fontsize', 15); hold off;

Current plot held

**Exercise.** Do the above spline interpolation using the command `fit`.