University of Waterloo
Faculty of Engineering
Department of Electrical and Computer Engineering
ECE 204A
Pre-Laboratory 6
Prepared by
Surname/Last Name, Legal Given/First Name(s)
UW Student ID Number: 2NNNNNNN
UW User ID: uwuserid @uwaterloo.ca
2A Electrical/Computer Engineering
18 September 2015
6.1 Use format long
6.1a Write a function that takes a vector x = (x1, x2)T as a parameter and returns the vector (cos(x1 – x2), sin(x2 – x1))T. Enter your function here:
function [y] = f61( x )
y =
return;
end
Starting with the initial vector x0 = (0, 0)T, perform 100 iterations of x = f61(x). What is the value of x100, the 100th iteration? It should be close to (0.3070, –0.9517)T.
Copy and paste your Matlab commands and the output here.
6.1b Assign the value of the 100th iteration to the variable x100 and calculate:
x101 = f61(x100)
x102 = f61(x101)
Does it appear that the values x100, x101 and x102 are the same?
Yes or No
6.1c Now, look at the values using format hex:
format hex
x100
x101
x102
format long
Are x100 and x101 or x101 and x102 the same vectors?
Yes or No
Are x100 and x102 the same vectors?
Yes or No
6.1d In Laboratory 4 for Math 211, you implemented the function fp.m . Modify this function to work with a vector-valued function by replacing the call to abs to one with a call to norm. Save this as the function fpv and enter your implementation here:
function [x, k] = fpv( f, x0, eps_step, N_max )
% Enter your implementation here
end
6.1e Again with x0 = (0, 0)T and the function f61, how many iterations are required for the fixed-point interaction to converge using estep = 10–5 and Nmax = 100?
Copy and paste your Matlab commands and the output here.
6.1f What is the norm ||x100 – xfpv||2 where xfpv is the point returned by the function fpv with the arguments from 6.1e and x100 is the previously calculated value after 100 iterations. Recall that || . ||2 is the 2-norm (see help norm).
Copy and paste your Matlab commands and the output here.
6.1g Given the matrix
M = [1 2 3; 4 5 6; 7 8 9];
what is the Matlab cod to extract the 2nd row of this matrix?
Copy and paste your Matlab commands and the output here.
6.1hGiven the matrix and vector
M = [1 2 3; 4 5 6; 7 8 9];
x = [5 7 3]';
y = [0 0 0]';
what is the Matlab command that will assign to the 1st entry of the vector y, the value of the first row of M multiplied by the vector x? At the end of your operation, y should be the vector
[? 0 0]' where ? is the dot-product of the first row of M and the column vector x.
Copy and paste your Matlab commands and the output here.
6.1iThe operation
y = M*x;
calculates the result of multiplying the vector x by the matrix M. Write a function matrix_vector_multiply which does the same thing, but does so by:
1. Creating a new vector of zeros called y = (yk) which has the same dimension as the row dimension of M,
2. For each entry yk where k goes from 1 to the row dimension of M, assign to yk the value found by multiplying the kth row of M with the column vector x.
function [y] = matrix_vector_multiply( M, x )
% Enter your implementation here
return;
end
6.1j What command accesses the row dimension of a matrix or vector? What command accesses the column dimension? See help size.
Copy and paste your Matlab commands and the output here.
6.1kRead the comments on and save the function jacobi in a file jacobi.m and execute the following:
> M = [5 3 -1; -2 6 1; 2 -1 7];
> x = [1 1 1]';
> M \ x
> jacobi( M, x, 1e-5, 100 )
> jacobi( M, x, 1e-10, 100 )
> jacobi( M, x, 1e-15, 100 )
Copy and paste your Matlab commands and the output here.
6.2 Use format long
6.2a Approximate the slope of the sine function at x = 1 by attempting the following:
cos(1) % the actual value of the derivative
h = 1e-5;
(sin(1 + h) - sin(1 - h))/(2*h)
h = 1e-10;
(sin(1 + h) - sin(1 - h))/(2*h)
Copy and paste your Matlab commands and the output here.
6.2b What happens if you attempt to calculate
(sin(1 + h) - sin(1 - h))/2*h
Your explanation in English.
6.2c Given two vectors, how do you find the element-wise multiplication of the two vectors?
x = [1 2 3]';
y = [4 5 7]';
Copy and paste your Matlab commands that produces [4 10 21]' using element-wise multiplication here.
6.2d Implement the following two functions in Matlab: First, implement the function
,
saved in f62.m and a second function stored in a file df62.m which evaluates the function
You will note that the second function, df62, is the derivative of the function f62. Make sure that your function works for a vector input.
function [y] = f62( x )
% Enter your implementation here
end
function [dy] = df62( x )
% Enter your implementation here
end
Execute the following two Matlab commands:
format short
f62( 0:3 )
df62( 0:3 )
format long
Copy and paste your Matlab commands and the output here.
6.2e Rather than always calculating the derivative by having to type
(f62(1 + h) - f62(1 - h))/(2*h)
every time we want to approximate a derivative, it would be easier to create a function which does this for us!
The formula described above is the centred divided difference approximation of a derivative because the point x0 was in the centre of the two values x0 – h and x0 + h. For this, we have would have to create a function Dc (derivative-centred) that takes three arguments: f, x and h and would assign to the return value
.
function dy = Dc( f, x, h )
% Enter your implementation here
return;
end
You can test your function:
h = 1e-5;
Dc( @sin, 1, h ) % these two should be the same
(sin(1 + h) - sin(1 - h))/(2*h)
Dc( @f62, 1, h ) % as should these two
(f62(1 + h) - f62(1 - h))/(2*h)
Copy and paste your Matlab commands and the output here.
You will note that we pass @sin and @f62. The @ operator is close to the unary operator in C++: create a reference or pointer to the function; however, in Matlab, the terminology is function handle. See help function_handle if you are interested in more information.
6.2f Section 6.2.1 of the laboratory describes how we can approximate not the derivative, but the error of the derivative based on the 3rd derivative of the function. Of course, in reality, we don’t have access to the third derivative; however, what this does is it gives us confidence in our approximations.
We will explore this here:
The third derivative of sin(x) is –cos(x) and thus,
h = 1e-4;
approx = Dc( @sin, 1, h )
exact = cos( 1 )
Should approximate the derivative of sin(x) at x = 1. The correct answer is, of course, cos(1). The laboratory describes how we can approximate the error by the formula
.
actual_error = exact - approx
approx_error = -(-cos(1))/factorial(3)*h^2
What is the relative error of the approximation of the error? Recall that the relative error is where x is the actual value and a is the approximation.
Copy and paste your Matlab commands and the output here.
This is a good opportunity to learn command completion in Matlab. Suppose you don’t want to type abs( (actual_error - approx_error)/actual_error ) every time you want to find the relative approximation. Every time you start to type a word, say act, you can hit the Tab key and you will be able to select from a menu all possible variables starting with act. If you hit the Tab key and there is only one unique variable, it completes it. For example, type actu and then hit the Tab key to complete it.
To demonstrate this is not a fluke, d3f62 calculates the 3rd derivative of the function f63.
function [d3y] = d3f62( x )
d3y = 3*exp(-x) - x*exp(-x);
end
and demonstrate that your code works correctly by approximating the derivative at
x = 0.42 and finding the actual error and our approximation of the error.
h = 1e-4;
approx = Dc( @f62, 0.42, h )
exact = df62( 0.42 )
actual_error = exact - approx
approx_error = -d3f62(0.42)/factorial(3)*h^2
What is the relative error of the approximation approx and what is the relative error of approx_error.
Copy and paste your Matlab commands and the output here.
6.2g One nice property of the formula is that if the error is O(h2), then halving h should reduce the error by a factor of 4. You may see this because and you will recall in ECE 250, if you double the size of a problem that has a O(n2) run time, it takes approximately four times as long to run.
h = 1e-4;
a1 = Dc( @f62, 0.42, h )
a2 = Dc( @f62, 0.42, h/2 )
exact = df62( 0.42 )
err1 = exact - a1
err2 = exact - a2
By referring to the two errors, explain why this example supports the previous argument.
Your answer here.
6.2h Section 6.2.2 of the web site describes how the centred divided difference method fails if h is too small. Demonstrate this by executing
h = 2^-29; % approximately 1e-9
approx = Dc( @f62, 0.42, h )
exact = df62( 0.42 )
actual_error = exact - approx
approx_error = -d3f62(0.42)/factorial(3)*h^2
What is the power of 10 of the actual error and what is the power of 10 of the expected approximation of the error.
Your answer here.
6.2i Now, use format hex and examine both the variables approx and exact. Describe in words the difference between the approximation of the derivative and the exact value of the derivative using df62.
format hex
approx
exact
format long
Describe the difference here.
6.2j In Section 6.2.3, we describe Richardson extrapolation. To begin, since you have already implemented Dc.m, save the following as
function dy = D1st( f, x, h )
dy = (f(x + h) - f(x))/h;
return;
end
This is the formula you learned for the derivative
which you learned in 1st year (hence the name). This approximation is not a good as the centred-divided difference formula, as we can take the Taylor series
.
It is quite easy to solve for to get
.
What this says is that using this formula, if we halve h, we only reduce the error by two (not four). Execute the following and explain how this demonstrates the concept explained above.
h = 1e-4;
a1 = D1st( @f62, 0.42, h )
a2 = D1st( @f62, 0.42, h/2 )
exact = df62( 0.42 )
err1 = exact - a1
err2 = exact - a2
Your answer here.
6.2k It may be useful to create a file called richardson11.m. You will need this in the laboratory.
10