MATLAB Review

26/09/2018

This is a MATLAB tutorial guide for the course Mathematical Modeling (MATH3290) in Chinese University of Hong Kong (CUHK). All the codes can be executed on MATLAB 2016a (or version above). It requires MATLAB 2016a (or version above) to open this live script.

There are many applications for MATLAB. In this guide, we focus on the following topics:

Basic operations in MATLAB;
Common built-in functions in MATLAB;
Flows of control, including if-else, for/while-loop;
Basic plotting: how to visualize your results;
Debugging and getting help.

Finally, we provide some (programming) practice problems related to our course.

Basic operations in MATLAB

Addition (or subtraction)

% Between scalars
3 + 5   

3.14 + 1.57

2/7 + 0.16
format short
2/7 + 0.16

% Between scalar and vector (or matrix)
1 + [1;1;1;1]

I = 1 + [1,0;0,1]
disp(I)

% Between vectors / matrice
[1;2;3;4] - [5;6;7;8]

new_vector = [1;2;3;4]

Remark

Each row of matrix end with the sign ";".
For subtraction, just replace "+" by "-" with exactly the same syntax.
It is useful to grasp part of the matrix using the sign ":" and the keyword "end".

A = magic(5)

A1 = A(2:3, 4:5)
disp(A1)

A2 = A(1:3,:)
disp(A2)

A3 = A(:,1:3)
disp(A3)

A4 = A(:,2:end)
disp(A4)

Multiplication

% Between scalars
2 * 4

3.14 * 1.57
format short
3.14 * 1.57
format long
ans

2.7 * 3/2
format short
ans

% Between scalar and vector (or matrix)
3.5 * [1;1;1;1]

pi * [1,2;3,4]
format long
pi
% use frpintf to print out more digits

% Between matrice (MATRIX-MULTIPLICATION)
[1 2 3] * [4;5;6]

[1 2 3; 4 5 6] * [7 8; 9 10; 11 12]

% Between matrice (ENTRY-WISE)
A = [1 2 3; 4 5 6]

B = [7 8 9; 10 11 12]

A .* B

a = A(1,:)'

b = A(2,:)'

sum(a.*b)

a'*b

Division (Inversion)

r1 = 2 / 7

r2 = 2 \ 7

A = [1 2; 4 5]
B = [7; 10]
x = A\B
% A x = B;

format rat
invA1 = A^(-1)
invA2 = inv(A)
disp(invA1)
disp(invA2)

Remark: We prefer to use x = A\B rather than x = inv(A)*B.

Power

p1 = 2^5
p2 = power(2,5)

A = [1 2; 4 5]
B = [3 4; 5 6]
result = A .^ B
disp(A)
disp(B)
disp(result)

a = 27^(1/3)
b = nthroot(27,3)
c = 16^(1/2)
d = sqrt(16)

Summary

Basic operations including: addition, subtraction, multiplication, division and power.
Display format in MATLAB: short, long, rat, bank, etc.
For some operations (e.g., multiplication or power) there are entry-wise versions .* or .^.
If you are doing the entry-wising operation, the two operants have to be same size.
One may solve Ax = b by using the operation "\" in MATLAB: x = A\b.

Common built-in functions in MATLAB

There are many useful built-in functions in MATLAB and we just name very few of them here, which are frequently used in this course.

zeros, ones, repmat, kron, linspace, diag, eye;
mean, norm, sum, abs, max, min, det, rank;
rand, randn, mod, floor, round, ceil, eig;
who, whos, clear, clc;
plot, scatter;
disp, save, load;

zero_mat = zeros(3,4);
disp(zero_mat)
ones_mat = ones(3,4)
disp(ones_mat)
A = ones(2,1);
A_rep = repmat(A,1,2)
disp(A)
disp(A_rep)
random_field = rand(3,3);
disp(random_field)

Users may want to define a function by themselves. There are two common ways to do so.

% Function handle
% function_name = @ (list of arguments) mathematical expressions
clear a b c;
quad_root = @(a,b,c) [(-b+sqrt(b^2-4*a*c))/(2*a) (-b-sqrt(b^2-4*a*c))/(2*a)];
r = quad_root(1,2,1)
disp(r)
p = quad_root(1,-3,2)
disp(p)

% Define a function in a separate m-file.
[x1,x2] = quad_solver(1,2,1)
% function [x1,x2] = quad_solver(a,b,c)
%     x1 = (-b+sqrt(b^2-4*a*c))/(2*a);
%     x2 = (-b-sqrt(b^2-4*a*c))/(2*a);
% end

Flows of control

If-else statement.

The logical expression can be some variables, inequality or equality. We list some symbols used in logical expressions.

<, less than;
>, greater than;
<=, less than or equal to;
>=, greater than or equal to;
==, equal to;
~=, not equal to.

% if-else
% Syntax
% if (logical expression)
%     program statements ...
% end
% if (logical expression)
%     program statements ...
% else
%     another program statements ...
% end
a = 0; b = 0; c = 1;
if (b^2 - 4*a*c) >= 0
    disp('This quad-eqn has two real roots.');
else
    disp('This quad-eqn has two imaginary roots.');
end
a = 1; b = 0; c = 1;
if (b^2 - 4*a*c) >= 0
    disp('This quad-eqn has two real roots.');
else
    disp('This quad-eqn has two imaginary roots.');
end

For-loop statement

We use for-loop in MATLAB in order to execute some commands specific number of times. Usually, some_vector may be 1:1:n, where n is a large integer.

% For-loop

% for idx = some_vector

% program statements ...

% end

% Example: Forming Hilbert matrix without built-in function.

% Plotting the condition number of Hilbert matrix against its size.

% cond(H) = || H ||_2 * ||H^(-1)||_2

% || H ||_2 = max_{x neq 0 in R^n} ||Hx||_2 / ||x||_2.

cond_H = [];

for sz = 1:1:10

H = zeros(sz,sz);

for m = 1:sz

for n = m:sz

H(m,n) = 1/(m+n-1);

H(n,m) = H(m,n);

end

cond_H = [cond_H cond(H)];

end

scatter(1:10,log10(cond_H))

title('Condition number of Hilbert matrix','FontSize',16)

xlabel('Size of matrix','FontSize',14)

ylabel('Log10 of Condition number','FontSize',14);

While-loop statement

We use while-loop to execute (repeatedly) when the condition is true.

If the experssion is true, the program statements will be excuted;
After finishing all statements, the expression on top will be determined again;
The program will keep being run if expression is true, otherwise terminate.

% While-loop
% while expression
%    program statements ...
%    ... ...
% end
% For example, calculating sin(x) by power series up to certain number of power.
x = 2; s = 1; max_power = 10;
% sin(x) \approx x/1! - x^3/3! + x^5/5! - x^7/7! + x^9/9! 
format short;
res = 0;
while s <= max_power
    res = res + (-1)^((s-1)/2)*(x.^s)/(factorial(s));
    s = s+2;
end
s
res
error = abs(res - sin(x))

Basic plotting

Visualizing your results is important. A graph is better than a thousand of words. One may sketch some figures by MATLAB in assignments if necessary.

First, we look at a typical example of mathematical modeling. Let

be the population of mouse after n years and

the population fo owl after n years. The ecologist has suggested the following model

We want to know if two species can coexist in the habitat and if the outcome is sensitive to the starting population. Given

and

. How does the population of two species change after

years?

M = 200; O = 150;
n = 40;
record = zeros(2,n+1);
record(:,1) = [M;O];
for i=2:n+1
    M_new = 1.2*M - 0.001*O*M;
    O_new = 0.7*O + 0.002*O*M;
    record(:,i) = [M_new;O_new];
    M = M_new;
    O = O_new;
end

How can one show the results? Use the built-in function plot to sketch the figures. Further, one can set different kind of properties of the figure such as title, label and legend. One may also save the figure as a fig.file using the built-in function saveas.

fh = figure;

h = plot(1:(n+1),record(1,:),1:(n+1),record(2,:));

set(h(1),'LineWidth',2); set(h(2),'LineWidth',2);

set(h(1),'Marker','*'); set(h(2),'Marker','s');

set(h(1),'Color','b'); set(h(2),'Color','r');

h_l = legend('Mice', 'Owl');

set(h_l, 'Location', 'NorthWest');

title('Owls vs Mice','FontSize',16);

xlabel('Number of years','FontSize',14);

ylabel('Population','FontSize',14);

% saveas(fh,'assignment_plot','fig');

Application: Cubic splline

In this section, we investigate in MATLAB the calculation of a cubic spline, which is a commonly used mathematical tool in our course.

Consider the following data set

and find the natural cubic spline

that satisfies the data above. Assume that

Next, list all the conditions

satisfies

Solve the above linear system to obtain the coefficients

We formulate the above linear system in MATLAB.

x = [3 4 5]; y = [12 17 38];
sz = size(x,2); % sz = 3;
b_c = zeros(4*(sz-1),1);
b_c([1;2;4]) = y; b_c(3) = y(2); % b_c = [12; 17; 17; 38; 0;0;0;0]
A_c = zeros(4*(sz-1),4*(sz-1));
A_c(1,1:4) = x(1).^(3:-1:0); % 3:-1:0 = [3 2 1 0];
A_c(2,1:4) = x(2).^(3:-1:0);
A_c(3,5:end) = x(2).^(3:-1:0);
A_c(4,5:end) = x(3).^(3:-1:0);
A_c(5,1:3) = (3:-1:1).*(x(2).^(2:-1:0));
A_c(5,5:end-1) = -A_c(5,1:3);
A_c(6,1:2) = ([6:-4:2]).*(x(2).^(1:-1:0));
A_c(6,5:6) = -A_c(6,1:2);
A_c(7,1:2) = (6:-4:2).*(x(1).^(1:-1:0));
A_c(8,5:6) = (6:-4:2).*(x(3).^(1:-1:0));
x_c = A_c \ b_c; % A_c * x_c = b_c;

Then, we obtain the coefficients x_c:

and

One can sketch the graph of the cubic spline

. Note that it is differentiable, especially at

xd1 = 3:0.01:4; xd2 = 4:0.01:5;

% xd1 = [3 3.01 3.02 3.03 ... 4]

% xd2 = [4 4.01 4.02 4.03 ... 5]

yd1 = x_c(1)*(xd1.^3) + x_c(2)*(xd1.^2) + x_c(3)*(xd1) + x_c(4)*(xd1.^0);

yd2 = x_c(5)*(xd2.^3) + x_c(6)*(xd2.^2) + x_c(7)*(xd2) + x_c(8)*(xd2.^0);

plot(xd1,yd1,'b-','LineWidth',1.0);

hold on;

plot(xd2,yd2,'m-','LineWidth',1.0);

scatter(x,y,'ro','filled');

legend('S_1(x)','S_2(x)','Location','northwest');

title('Figure of S(x)','FontSize',16);

xlabel('x','FontSize',14);

ylabel('y','FontSize',14);

Debugging and getting help

How to debug?

One may usually have different kinds of bugs or error in the program when coding. The bugs may come from the error of syntax, which is relatively easy to find out. Another type of bugs may be logically wrong and it has correct syntax, leading to unexpected results, and this is usually harder to figure out.

Let us look at a simple example: quadsolver() and we test two cases

[x1,x2] = quadsolver1(1,-2,1);

This is good to know. Next, we try another case.

[x1,x2] = quadsolver1(1,0,1)

Why would it lead to the unexpected result? It should be the case that

This quad-eqn has two imaginary roots.

Go down to the script of the function of quadsolver and see what happens. One may set breakpoint in the program (if under the usual script environment) so that one may see how the program runs step by step.

One may define another function called check_quadsolver and verify the correctness of our code.

d = check_quadsolver(1,0,1)

After that, one may fix the code suitably and now check the final result using the function quadsolver2.

[x1,x2] = quadsolver2(1,0,1)

Problem fixed.

Note that the logical sign ">=" or "<=" only makes use of the real part of a complex number and one should avoid to compare the amount between a complex number and another quantity, no matter this quantity is complex or not.

There is no the concept of "order" in complex field.

You can do it in real field. (5>1, -5 < -1.35 ...)

Get help from MATLAB

The most efficient way to find a solution for a specific MATLAB issue is to use the keyword "help" or "doc" in MATLAB. For example, to explore the properties of the function "plot", one can type the following command to get some help.

help plot

Another way is to search the answer from internet, especially from Google or Stackexchange. Try to use accurate keywords to represent your problem. Note that, before asking any question, one should make the question as clear as possible and let people find it easier to provide the answer.

Practice problems

Problem 1

Consider the ordinary differential equation with initial condition

Here,

is given. One may use Euler's method to numerical solve this ODE and this scheme is as follows:

where

is a given step size of time,

and

In this problem, assume that the function and initial conditions are

Take

and

and write a program to solve the initial value problem. Record all the values of

and plot

against

. Also, compare the numerical solution to the exact solution of the initial value problem

Hint: To compare the accuracy, one may define a functional handle

g(t) = t*log(t) + 2*t and the vector variable error = abs(g(tn) - yn), where yn is the numerical approximation.

Problem 2

Consider the unconstrained minimization problem

where f is a differentiable (convex) function. Use classical gradient method to solve the problem. Choose an initial guess $x^0 = 0$ and repeat

Note that

is a descent direction and

is a step size. To determine

, using the line search strategy as follows: initialize

at some fixed

, repeat

until the following condition holds

where

and

and user-defined. Your task is to implement the gradient method with the following objective function

Plot a graph with the function value

against the iteration index k.

Appendix: definition of the funtions in this script.

function [x1,x2] = quadsolver1(a,b,c)
    d = sqrt(b^2-4*a*c);
    if (d >= 0)
        disp('This quad-eqn has two real roots.');
    else
        disp('This quad-eqn has two imaginary roots.');
    end
    x1 = (-b +d )/ (2*a);
    x2 = (-b -d )/ (2*a);    
end
function d = check_quadsolver(a,b,c)
    d = sqrt(b^2-4*a*c);
    disp('The delta is ');
    disp(d);
    s2 = sprintf('The logical expression is %d.', (d>=0));
    disp(s2);
end
function [x1,x2] = quadsolver2(a,b,c)
    d = b^2 - 4*a*c;
    if (d >= 0)
        disp('This quad-eqn has two real roots.');
    else
        disp('This quad-eqn has two imaginary roots.');
    end
    x1 = (-b + sqrt(d) )/ (2*a);
    x2 = (-b - sqrt(d) )/ (2*a);    
end
function [x1,x2] = quad_solver(a,b,c)
    d = sqrt(b^2-4*a*c);
    x1 = (-b +d )/ (2*a);
    x2 = (-b -d )/ (2*a);    
end        