UNB/ CS/ David Bremner/ teaching/ cs2613/ labs/ Lab 22

Working with grids

Time
15 minutes
Activity
Demo
• Octave has as special notation for ranges that we can use to initialize arrays.

``````  r = [-1:0.5:1]
``````
• The builting function meshgrid can be used to take Cartesian Products of two arrays.

``````  [X,Y] = meshgrid(r,r)
``````
• Despite what best practice is in languages like Java, Octave is quite fond of parallel arrays. In this example, the `(r(i),r(j))` element of the cartesian product is `[X(j,i),Y(j,i)]`. Since we just care about hitting all combinations, we won't worry too much about flipping indices

• For each point in our grid, we can compute a height as follows

``````  Z = X.^2 - Y.^2
``````
• We can then plot a surface using those heights with

``````  surf(X,Y,Z)
``````
• Things are very blocky here because we don't have many points. This is so we can understand the intermediate matrices, and how `surf` works. We'll use more points to give a smoother appearence in what follows.

Defining and testing a function

Time
20 minutes
Activity
Small Groups

In the rest of this lab we want to use Octave to explore the function

``````    δ(β, a, b) = βa·b - (β-1)b·b
``````

Here `β` is a scalar (a number) and `a` and `b` are vectors of length 2. We curious what the zeros of this function are for fixed `β` and `a`. One approach is to try plotting that surface as a function of `b`. We'll start by defining and testing a function for `beta`

• Translate the mathematical function `δ` into an Octave function.
```    ## compute βa·b - (β-1)b·b
## a, b are assumed to be column vectors, beta is a scalar
function ret = delta(beta, a, b)

endfunction
```
• Complete the tests according to given comments. Use `rand` to generate `a` and `b`. You will also have to look up the documentation for `assert` for how have a floating point error tolerance.
```    %% delta(0,a,b) = |b|²
%!test
%!
%!
%!

%% delta(1,a,b) = a.b
%!test
%!
%!
%!

%% delta(2,a,b) = 2a.b - |b|²
%!test
%!
%!
%!
```

Making 3D plots

Time
20 minutes
Activity
Small Groups
• Write one or more for loops to initialize `Z`
```    a = [4;4];
beta = 7.5;

%% Generating vectors
range = [-4:0.1:8];

% Compute cartesian product (grid)
[X,Y] = meshgrid(range,range);

surf(X,Y,Z);
```
• The resulting plot should look something like the following. • Try replacing `surf` with `contourf` to get a different view of the surface. What can you guess about the zero set (where the surface intersects the `z=0` plane).

Using `arrayfun`

Time
25 minutes
Activity
Small Groups
• Use an anonymous function to define `f` so that the following octave code produces the same plot as before.
```    a = [4;4];
beta = 7.5;
range = [-4:0.1:8];
[X Y] = meshgrid(range,range);

f=
Z=arrayfun(f,X,Y);

surf(X,Y,Z);
```
• Measure the speedup from the `for` loop version. You probably want to comment out the actual plotting step and just measure the matrix computation.

• How sensitive are both versions to the number of elements in `range`? If you double the size (half the step size), does the time double?

Fully vectorizing

Time
25 minutes
Activity
Small Groups

If we expand out a single call to `delta` in our for loop, we get

``````Z(i,j) = beta * (a(1)*X(i,j) + a(2)*Y(i,j)) - (beta-1) * (X(i,j)^2 + Y(i,j)^2)
``````
• Use this equation to complete the vectorized function `arraydelta` that satisfies the following test.
```    function ret = arraydelta(beta,a,X,Y)

ret =
endfunction

%!test
%! a = rand(2,1)
%! X = [2,2;2,2]
%! Y = [-1,-1;-1,-1]
%! assert (arraydelta(0, a, X,Y),[5,5; 5,5],eps)
```
• Use the anonymous function trick from above to complete the following test
```    %!test
%! a = rand(2,1);
%! beta = rand*6;
%! X = rand(10,10);
%! Y = rand(10,10);
%! assert (arraydelta(beta, a, X,Y), arrayfun(                            , X, Y), eps)
```
• Finally measure the speedup for new version. How sensitive is it to the size of `range`?
```    a = [4;4];
beta = 7.5;

range = [-4:0.1:8];
[X Y] = meshgrid(range,range);

Z=arraydelta(beta,a,X,Y);

surf(X,Y,Z);
```