F18XD Matlab 1

Maths for Engineers and Scientists 4: Matlab

Contents
1 Introduction to Matlab 1
1.1 What is it? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 On-line help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Using the Editor to create a new file . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 Matrix eigenvalues and eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5 Functions, for loops, if–then–else . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.6 Graphs and strings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.7 M-files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.8 Laplace and inverse Laplace transforms . . . . . . . . . . . . . . . . . . . . . . . . 9
1.9 The Matlab ODE solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Problem sets 11
2.1 Problem set 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Problem set 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Problem set 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Problem set 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5 Problem set 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1 Introduction to Matlab
These notes mostly refer to Matlab under Windows on the University’s desktop PC system.
The Univeristy has a comprehensive license for Matlab and you can download your own copy.
If you e-mail the IT Helpdesk ithelp@hw.ac.uk, they will send you installation and activation
instructions.
The following notes provide a quick summary of Matlab to allow a quick start. The booklet
An Introduction to Matlab by David F. Griffiths (downloadable through Vision) provides a
broader introduction to Matlab.

1.1 What is it?

Matlab is essentially a programming language and is similar to C, Pascal, Fortran and other
widely-used programming languages. (Technically it is a bit of a mixture between a compiler
and an interpreter.) It also has many built-in functions for finding approximate solutions of
mathematical problems, doing linear algebra etc. and producing good graphs of the results.
Matlab runs on most standard computer systems (PCs, Macs, Unix workstations etc.) and
looks and behaves “almost” the same on each one.
Matlab is widely used in higher education to teach linear algebra, numerical analysis, signal
processing and other topics because it is much easier to learn and to use than most other
programming languages. However, once you have mastered Matlab, you will find it easy to
learn any programming languages required for work in the real world1 or further study. Matlab
is increasingly being used in industry, both as a quick way to develop and test complex software,
and for the facilities provided by the many add-on “toolboxes” tailored for different applications.
1
real world = outside university
F18XD Matlab 2

These toolboxes include financial modelling, image enhancement, systems control, simulation
and optimisation.

1.2 On-line help

• If you know the name of the command you want to find out about, type doc *** in the
Matlab command window to get information on the Matlab command *** or M-file
***.m (here *** stands for any command name e.g. mean, std, sum etc. listed in Section
26 of the Introduction to Matlab).

1.3 Using the Editor to create a new file

• Launch Matlab.
• Click on New (at the top left of the Matlab command window) and select Script.
• The Matlab Editor/Debugger window should open, and type the line below into it
avec = [1, 2, 3]
• Now save the file by clicking on the Save option at the top left of the editor window.
Select Save As and you should see a dialogue box opening. You need to specify a name
for your file (call this one myprog.m) and put it in your home directory or in a folder in
your home directory.
• To test this out, click on the Run arrow/button in the editor window. You should get the
response

avec =

1 2 3

You might need to change the Matlab working directory to match the directory you
stored the file in.
• To test it out a different way, go back to the Matlab Command Window prompt » and
type myprog to see if you get the same output. You should get the same output.
• If it doesn’t work then go back and carefully follow the itemised steps above again. It is
very important that you learn how to save files before you create anything complicated.

Matrices
Matlab deals mainly with one type of object — rectangular matrices, whose entries can either
be numbers (real or complex), or text string characters. Matrices with only one row or column
can be treated as vectors, and 1 × 1 matrices as scalars. (See Sections 5, 8 and 15 of An
Introduction to Matlab.)
Type the lines
x = [-2, 4, -8]
y = [3, 5, 7]
z = [1; 2; 3]
a = 5
A = [-1 0; 3 1]
F18XD Matlab 3

to create two 1 × 3 (row) vectors (x and y), a column vector (z), a scalar (a) and a 2 × 2 matrix
(A). Note that unlike some programming languages the case of variable names does matter in
Matlab , so that a and A represent different quantities. Matrices are usually entered in for
loops (An Introduction to Matlab Section 18) or using colon notation. Note that ending a
line with a semicolon ; stops the answer being displayed. Try it out.
Individual matrix entries can be referenced by indices in brackets: e.g. typing

x(2)
x(1,2)
A(2,2)

results in 4, 4 (the first two commands refer to the same element of x), and 1. Note that a
matrix or vector will only accept positive integers as indices, so that B(12,197) is OK, whilst
B(0,1), B(2,-3), or B(1.3,1) are not.

Matrix operations
Most of matrix operations are listed in Section 15 of An Introduction to Matlab. Note that
many of these can be made to work entry-wise by typing a dot . before the operation — e.g.
A*A is the matrix A2 , whilst A.*A is the matrix with (i, j) entry equal to A(i, j)2 . Confirm this
in Matlab with the 2 × 2 matrix A below:

>> A = [-1 0 ; 3 1]
A =
-1 0
3 1

>> A*A
ans =
1 0
0 1

>> A.*A
ans =
1 0
9 1

See Section 12 and 15.9 of An Introduction to Matlab.

1.4 Matrix eigenvalues and eigenvectors

The number λ is an eignenvalue of the square matrix A if det(A − λI) = 0.
The vector v is an eignevector of the same matrix A associated with eigenvalue λ if

Av = λv.

There will be much more about these in the course. Meantime we just have a look to see how
to compute them in Matlab.
F18XD Matlab 4

To find eigenvalues of a matrix A we can use the eig command. e.g.

>> A = [-1 0 ; 3 1]

A =

-1 0
3 1

>> d=eig(A)

d =

1
-1

Vector d will contain the eigenvalues of matrix A.

To compute eigenevectors of the matrix A we type

>> A = [-1 0 ; 3 1]

A =

-1 0
3 1

>> [V,D]=eig(A)

V =

0 0.5547
1.0000 -0.8321

D =

1 0
0 -1

The columns of matrix V are eigenvectors of matrix A and the diagonal matrix D contains the
corresponding eigenvalues such that A*V = V*D.

1.5 Functions, for loops, if–then–else

Functions
These are introduced in Section 7 and listed in Section 26 of An Introduction to Matlab.
Note that many scalar functions apply entry wise to matrices or vectors. Try out the following
examples (remember that you can find out what the command *** does by typing doc ***).
F18XD Matlab 5

exp(x)
sy = sqrt(y)
abs(A)
floor(sy)
ceil(sy)
max(x)

for loops and colon notation

Type in the following lines

sum20 = 0
for i = 1:20
sum20 = sum20 + i
end

(the indentation isn’t necessary but makes the code clearer). This code:

• initialises the variable sum20 to zero;

• for each integer i between 1 and 20 performs the command sum20 = sum20 + i.

The command i = 1:20 is an example of colon notation. It means “all the integers between
1 and 20, starting with 1 and ending with 20”. If a, b and c are numbers, then a:b:c means “all
the numbers a, a+b, a+2b, . . . ending at or just before c”. The syntax is first:increment:last.
If you don’t include an increment (e.g. i=1:20 above), then it is assumed to equal 1.
Type in the following lines to get a better idea of what is going on.

1:2:10
5:-7:-30
14.32:-0.2:8.5
1:20
1.24:9.1

Why does the last example give a vector whose last element is 8.24 and not 9.1?
An Introduction to Matlab Sections 18 and 20 deal with for loops and related things in detail.

if . . . else statement
The if statement allows the execution of different commands depending on the value of some
logical test.

a = 3;
b = 2;
if (a<b)
j = -1;
else
j = 0;
end

What is the value of j in the above example?

F18XD Matlab 6

The code below produces a vector v = 1 2 2 2 2 3 3 3 3 3.

v = zeros(1,10);
for i = 1:10
if (i < 2)
v(i) = 1;
elseif (i > 1 && i <= 5)
v(i) = 2;
else
v(i) = 3;
end
end

Section 20.2 of An Introduction to Matlab explains the if statement in more detail.

1.6 Graphs and strings

Plotting graphs
The simplest plot command has the form plot(x,y), where x and y are vectors of the same
length, n say. It plots a graph of the points (x1 , y1 ), (x2 , y2 ), . . . , (xn , yn ). For example, type
in the following lines to get a plot of sin t against t.

t = [0:0.01:2*pi];
v = sin(t);
plot(t,v)

Section 13 of An Introduction to Matlab describes how to make multi-plots, 3D plots and how
to change the axes scaling and line styles. Axis labels and a title are essential to tell people
what is being plotted, and can be added using the xlabel, ylabel and title commands; e.g.
try out the command ylabel(’sin t’)

Text strings
The argument of the above ylabel command is an example of a text string — i.e. a collection
of characters surrounded by single quotes. Another example is
st = ’my name is Frankenstein’
Note that both quotes are the same. Commands like disp, input and error use text strings.
Another useful command is num2str, which converts a number into a string. As an example of
how to use it and a slightly more complicated string expression, type
title([’plot of sin t up to ’, num2str(2*pi)])
to label your graph.
Note that in the above example using title, the string has more than one element. If more than
one string element is used (typical when mixing words and numbers), then the string elements
to be joined together must be enclosed in square brackets and separated by commas.

1.7 M-files
Matlab can be used in two main ways. You can either type instructions directly at the prompt
or make up a file or files containing the instructions you want Matlab to execute and feed this
file (called an M-file) to Matlab. An M-file is a computer program written in Matlab’s own
F18XD Matlab 7

language. Typing commands at the prompt is fine if you want to execute single commands, but
errors usually creep into the typing of longer expressions and sequences of operations and it is a
pain having to retype things used over and over again. In general, you should use M-files, and
the mechanism for creating them is described in Section 1.3 above.
There are two different kinds of M-files, script files and function files, which act in slightly
different ways. Script M-files act as a shorthand for a sequence of Matlab commands and are
a bit easier to use than function M-files. We discuss both types of M-files below, and Sections
10 and 22 of An Introduction to Matlab also describe them in detail.

A sample script M-file

Suppose we want to write an M-file to find the sum of the first 20 integers, instead of typing the
commands at the Matlab prompt. The following file which we shall call S20.m (note that all
M-files end in .m) does what is required. Follow the instructions in Section 1.3 and create this
example file.

% S20 Find the sum of the first 20 integers

%
% Computes and outputs the sum 1 + 2 + ... + 20
% in the variable sum20.

% Note that this line is ignored.

sum20 = 0; %initialise sum20 to zero

for i = 1:20
sum20 = sum20 + i;
end
disp([’The value of sum20 is ’,num2str(sum20)])

Typing the name S20 at the Matlab prompt will execute the M-file and should result in the
line

The value of sum20 is 210

Note that if you miss the semicolon from the line sum20 = sum20 + i then you will get a whole
lot of unwanted output.
The semicolon ; suppresses the output.
Create the file and try it out (make sure that you have learnt how to Save files first). Typing
doc S20 will print out the lines at the top of the file that start with a % (the comment lines),
up to the first line that doesn’t start with a %. This is a very useful way of letting yourself (and
anyone else who uses your Matlab programs) know what a particular file does. It is good style
to always start an M-file with a comment line. Matlab ignores any commands that follow a %,
and so comments can be inserted in the body of the program after a %.
F18XD Matlab 8

Sample function M-files

The following file (called mkvec.m) is a simple function M-file.

function v = mkvec(x)
%MKVEC Construct the vector [1, x, x^2, x^3] from x
% mkvec(x) returns the vector [1, x, x^2, x^3]

v = [1, x, x^2, x^3];

It contains an input argument (x), the 4–component vector output argument v and a function
statement (the first line). Write and save the file and then type mkvec(2) at the command line.
You should get the result

ans =

1 2 4 8

You can use the M-file to create named vectors — for example if you type w = mkvec(6) then
you obtain

w =

1 6 36 216

instead.
Function M-files can have more than one input argument — e.g. the file mkvec2.m below:

function v = mkvec2(x,n)
%MKVEC2 Construct the vector [1, x, x^2, ..., x^n] from x and n
% mkvec2(x,n) returns the vector [1, x, x^2, ..., x^n]

for k = 0:n
v(k+1) = x^k;
end

and more than one output argument — e.g. the file mkvec3.m below:

function [v, len] = mkvec3(x,n)

%MKVEC3 Construct the vector [1, x, x^2, ..., x^n] from x and n
% [v, len] = mkvec3(x,n) returns the vector
% v = [1, x, x^2, ..., x^n] and its length len.

for k = 0:n
v(k+1) = x^k;
end
len = length(v);

You can also use dot operations to avoid loops in the previous example, e.g. as in the file
mkvec3a.m below:
F18XD Matlab 9

function [v, len] = mkvec3a(x,n)

%MKVEC3A Construct the vector [1, x, x^2, ..., x^n] from x and n
% [v, len] = mkvec3(x,n) returns the vector
% v = [1, x, x^2, ..., x^n] and its length len.
% Does not use loops

k = 0:n;
v = x.^k;
len = length(v);
See Section 22 of An Introduction to Matlab for more details.

1.8 Laplace and inverse Laplace transforms

To compute the Laplace and Inverse Laplace transforms we can use the functions laplace and
ilaplace from the Matlab Symbolic Toolbox
To comute a Laplace transform of function f (t) = t we type at the Matlab prompt
>> syms t
>> f=t;
>> laplace(f)

ans = 1/s^2
To compute an inverse Laplace transform of function F (s) = 1/s2 we type
>> syms s;
>> F = 1/s^2;
>> ilaplace(F)

ans = t
It is also possible to use the package Maple to compute Laplace and inverse Laplace transforms.
Maple is also avaialble on a comprehensive campus license. However, let us stick to Matlab
here.

1.9 The Matlab ODE solvers

The Matlab ODE solvers are written to solve systems of DE’s of the form
dx1
= f1 (t, x1 , x2 , x3 , . . . )
dt
dx2
= f2 (t, x1 , x2 , x3 , . . . )
dt
dx3
= f3 (t, x1 , x2 , x3 , . . . )
dt
...

i.e.
dx
= F (t, x).
dt
There are several built-in ODE solvers from which we can chose, e.g., ode45, ode23, ode113,
ode15s, ode23s, ode23tb (each solver has its advantages and disadvantages).
F18XD Matlab 10

Below is a sample code for the solution of the system of DEs

    
d x1 −1 −1 x1
=
dt x2 1 −2 x2

using the ode45 solver.

function main
% Use ODE45 to solve
% dx_dt(1) = -1*x(1)-1*x(2)
% dx_dt(2) = 1*x(1) -2*x(2)

options=odeset(’RelTol’,1e-6); %set error tolerance for the ODE solver

xinit = [1;1]; %initial condition

tspan = [0,5]; %time interval

[t,X] = ode45(@F,tspan,xinit,options); %call the solver

%plot the results

figure
hold on
plot(t,X(:,1));plot(t,X(:,2),’:’)
legend(’x1’,’x2’);ylabel(’x’);xlabel(’t’)
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rhs function x’=F(t,x)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx_dt]= F(t,x)
M = [-1,-1; 1,-2];
dx_dt = M*x;
return
F18XD Matlab 11

2 Problem sets
2.1 Problem set 1
1. (a) What is the semi-colon for?
(b) Generate the row vector x = (0.6, 0.7, . . . , 1.4, 1.5) using 2 colons and three numbers.
(c) Find the transpose of the vector x in (b).
(d) If A and B are n × n matrices, what is the difference between A ∗ B and A. ∗ B? Try
it out for some 2 × 2 matrices that you choose yourself.

2. Use the colon notation to generate the vector x = (0.6, 0.7, . . . , 4.4, 4.5).

(a) Find the sum of all its elements.

(b) Find its length.
(c) Find the sum of the squares of all its elements.

Use length and sum and see Section 21.2 of An Introduction to Matlab.

3. Generate the vector x = (0, 0.1, 0.2 . . . , 10).

(a) Find the vector s whose jth element is the sine of the jth element of x.
(b) Plot vector s against x using the line style “stars” and put appropriate labels on the
graph.

See Section 13 of An Introduction to Matlab.

4. (a) Enter the matrix  

1 1/2 1/3
A =  1/2 1/3 1/4  .
1/3 1/4 1/5

(b) Enter a vector b = (3, 1.9167, 1.4333)T (i.e. b is a column vector).

(c) Find the solution x to the linear system Ax = b (hint: use \).
Note that the matrix A is the 3 × 3 Hilbert matrix, and can be defined in Matlab more
simply by using the command A = hilb(3).

5. What effect do the commands format short e, format long e and format compact have
on Matlab output ? Try it out by typing pi after the various format commands. (See
help format).

6. Create

(a) a 9 × 9 identity matrix,

(b) a 7 × 9 zero matrix,
(c) a 4 × 6 matrix of ones.

Use (eye, zeros, ones).

F18XD Matlab 12

7. (a) Use the editor to create a function M-file called anysump.m which takes as input the
positive integer N , and as output the result p of the summation
N
X
p= j −2 .
j=1

(b) Compute from the Matlab command window, the value of the sum for N = 5, 10, 15, 20.

8. (a) Follow the examples in Section 13 of An Introduction to Matlab to plot f (x) =

x2 + x − 1 with a “dashed” line style. Use x = −5 : 0.05 : 5 for the x coordinates.
(b) Add labels to the x and y axes and add your name as the title of the graph.
(c) Annotate the graph to show which curve is being depicted. (legend, see the Lecture
notes)

9. (a) Follow the examples in Section 13 of An Introduction to Matlab to plot f (x) =

x + x2 − x3 /10 − 10 with a “dashed” line style and g(x) = cos(x) with a solid line in
the same graph. Use x = 0 : 0.05 : 10 for the x coordinates.
(b) Add labels to the x and y axes and add your name as the title of the graph.
(c) Annotate the graph to show which curve is which and place the text “f=g here” near
where the curves cross. (legend, text, see the Lecture notes)

10. Write a function M-file that takes as input the positive integer n and number r, and
outputs the sum
Xn
kr .
k=1

Use it to check the identity

n
X n (n + 1) (2n + 1)
k2 =
6
k=1

for various values of n.

11. Plot an (x, y) graph for x ∈ [−5, 5] when y = sin(x) + 1 − x. Use the information from the
graph to help the built-in function fzero find the root of the equation sin(x) + 1 − x = 0.
F18XD Matlab 13

2.2 Problem set 2

1. Produce 4 separate graphs against x for x ∈ [−π, π] (labelled appropriately), of the func-
tions (a)-(d) below all in the same Figure window (subplot).

(a) x3 /20 + x2 /10,

(b) cos x,
(c) sin2 x,
2
(d) e−x .

2. (a) Write a function M-file using for loops to define the N × N tridiagonal matrix A and
vector b given by
 
β γ 0 ... 0  
. sin(1)

 α β γ . 
.   sin(2) 
 
A= 0
 . . . . . .  , b=
 
···

 ,
. . . 0  
 .. sin(N −1) 
  

 . α β γ  sin(N )
0 ··· 0 α β

which will work for any integer N ≥ 3. Make sure that b is a column vector. The
matrix A is called tridiagonal because all its entries are zero apart from those on or
immediately next to the diagonal.
(b) For α = 1, β = 4, γ = −1, find the solution x of Ax = b in the case N = 4, and
find the 42nd element of the solution x in the case N = 199. (Use the backslash \
command.)
(c) Use the command eig to find the eigenvalues and eigenvectors of the matrix A with
N = 5, α = γ = 1, β = 2.

3. So far we have seen two ways of solving Ax = b for x where A is an N × N matrix and b
an N vector:

Method Matlab commands

Gaussian elimination x = A\b;
Full inversion Ai = inv(A); x = Ai*b;

We would like to test which method takes the least cpu time.

(a) Write an function M-file where the input is the size N and the output is the time
(cputime) in seconds it takes by each of the two methods to compute the solution of
Ax = b for randomly chosen A and b (rand).
(b) Then write a second M-file for depicting the two resulting times for N = 100, 150, 200, . . . , 800.
Label the graph appropriately.

4. Newton’s method for finding a solution of the equation f (x) = 0 is the following iteration.
Step 1 Make a guess at the solution (call it x1 ),
Step 2 calculate the sequence of values x2 , x3 , . . . from
f (xk )
xk+1 = xk − for k = 1, 2, 3 . . . ,
f 0 (xk )
F18XD Matlab 14

Step 3 stop when |xk+1 − xk | < TOL, where TOL is a small number supplied by the user
or if we reach k = 100 without finding a solution.

(a) Write a function M-file with input: TOL, initial guess x1 and output: the Newton
method approximation of the solution of cos(x) = x, and the number of steps taken
to get there. It will make life easier for you if you write the M-file in such a way that
the equation to be solved can be changed easily.
(b) Type format long e to show numbers to full precision before the next part.
(c) Find the root of cos(x) = x using the built-in function fzero first of all, then apply
your Newton function to the equation cos(x) = x with TOL= 10−6 and try two
different starting guesses x1 = 1.4 (which works) and x1 = −1.4 (which doesn’t).
Does you solution agree with that from fzero? Now try a wider range of starting
guesses −2 ≤ x1 ≤ 2. Is there a pattern in which starting guesses work and which do
not? (for, if, break).
F18XD Matlab 15

2.3 Problem set 3

1. (a) Follow the example in Section 23 of An Introduction to Matlab to produce a mesh
plot of z = x sech(x2 + 2y 2 ) over the square [−2, 2] × [−2, 2] in the subplot(2,2,1).
(b) Now produce a contour plot of the function (contour) in subplot(2,2,2) and mark
it with a star at its maximum.

2. The derivative of a well-behaved function f (x) with respect to x can approximated by

various finite difference approximations:

f (x + h) − f (x) f (x + h) − f (x − h)
(F) : f 0 (x) ≈ , (C) : f 0 (x) ≈ .
h 2h
The labels F and C stand for “forward” and “central” respectively. The number h is usually
relatively small, and the approximations should, in theory, get better as h decreases.

(a) Write a function M-file that takes as input x and h, and outputs the forward and
central approximate derivatives of f (x) = sin−1 (x).
(b) Set x = 1/2 and plot the modulus of the error in the results against h for h =
10−10 , 10−9 , . . . , 10−2 . Which gives better results and why do they stop getting better
as h decreases? Use the loglog function to plot, and label the graph appropriately.
F18XD Matlab 16

2.4 Problem set 4

1. Use Matlab (or Maple) to compute Laplace transforms from the Table of Laplace trans-
forms from the lecture notes and from the Tutorial exercises.

2. Use Matlab (or Maple) to find inverse Laplace transforms of the functions from the
Tutorial exercises.

3. Solve the LCR circuit example problem from the Lecture notes using the Matlab built-in
ODE solvers and compare the computed solution to the exact solution from the lecture
notes.

4. Use the Matlab build-in ODE solver to solve the following system of ODEs for t ∈ (0, 5)

dx1
= x1 + 2x2 ,
dt
dx2
= 2x1 − 2x2 ,
dt
x1 (0) = 2 x2 (0) = 1, ,

Compare the computed solution with the exact solution from Section 3.5 of the Lecture
notes.

5. Use the Matlab build-in ODE solver to solve the following system of ODEs for t ∈ (0, 5)

dx1
= −4x2 + δ(t − 3) ,
dt
dx2
= x1 ,
dt
x1 (0) = 0 x2 (0) = 1, ,

where δ(t) is the Delta function (Hint: replace the delta function by an approximation
δ̃(t) = 1/ε for t ∈ (0, ε) for a sufficiently small , e.g.  = 0.01). Compute for different
values of  and compare the computed solution with the exact solution (cf., Lecture notes
Figure 3.2) 
cos (2t) t < 3,
x2 (t) =
cos (2t) + 0.5 sin (2t − 6) t ≥ 3.
F18XD Matlab 17

2.5 Problem set 5

1. The Midpoint rule gives a relatively reliable
Rb way to approximate integrals which cannot
easily be worked out exactly. To find a f (x) dx, the interval a ≤ x ≤ b is split into N
strips of equal width, and the function is assumed to be constant over the strip making
it easy to integrate. The contributions from all the strips are added up to give the final
approximation. The approximation formula that results is

Zb N
X
I≡ f (x) dx ≈ h f (xj ) ≡ IN
a j=1

where the width of the strips is h = (b − a)/N and the mid points of the strips are
xj = a + h(j − 21 ). The error in the approximation IN is roughly one third of the difference
between IN/2 and IN .
R2
(a) Write an M-file to approximate the integral 1 sin(x3 ) dx with N = 8, 16, 32, . . . , 1024
strips.
(b) Tabulate the result to full precision and make an estimate of the error in the results.

2. One of the best known examples in dynamical systems and chaos theory is the Logistic
Map. It is defined by xk+1 = λxk (1 − xk ) for k = 1, 2, . . . with starting value x1 ∈ (0, 1).
The constant 0 < λ ≤ 4 controls the type of behaviour observed.

(a) Write a function M-file with inputs λ (controlling parameter), N (number of steps)
and x1 (starting value), and output the vector (x1 , x2 , . . . , xN ).
(b) Use this to plot xk vs k for various values of λ and x1 in four horizontal strips on the
same figure. (subplot).
(c) Try the case λ = 3.7, N = 100 with x1 = 0.7, 0.7 + 10−9 , 0.7 + 10−6 , 0.7 + 10−3 and
comment on how many steps it takes the results to deviate from the x1 = 0.7 case.
(d) Try the case x1 = 0.7, N = 500 with λ = 3.3, 3.5, 3.739, 3.835 and comment on the
pattern the solution settles into. (It might help to list the last 10 numbers in the
sequence in each case.)

3. Euler’s method for solving the first order differential equation

dy
= f (y, t)
dt
for 0 ≤ t ≤ l subject to the initial condition y(0) = a.
Step 1 Fix a small value h > 0.
Step 2 Set t0 = 0 and tn = tn−1 + h.
Step 3 Find an approximation of y step-by-step from the recursion

yn+1 = yn + hf (yn , tn )

with initial value y0 = a.

If h is sufficiently small, this often produces a good approximation of the solution to the
above first order initial value problem.
F18XD Matlab 18

(a) Write an M-function with inputs the parameter h and length of the interval l, and
output the two vectors t and y obtained by applying Euler’s method to the initial
value problem
dy
= (1 − t)y, y(0) = 1.
dt
(for end).
2
(b) The exact solution of this ODE can actually be found analytically, y(t) = et−t /2 .
Plot a graph comparing the exact and the approximate solutions for 0 ≤ t ≤ 5. Use
a blue filled line for the exact solution and a green dashed line for the approximated
one. Label the graph appropriately. (plot).
(c) Use one of the built-in Matlab ODE solvers (e.g., ode15s) to solve the above ODE
and compare the resulting solution to the one obtained by the Euler’s method and
to the exact solution.

Dugald Duncan, 2018.

Based on previous versions by
Ľubomír Baňas, Lyonell Boulton and Dugald Duncan