CS 101 - Introduction to Computing

Machine Problem 1
Spring 2012
Released: Monday, February 6, 2012
              Due Date: Friday, March 2, 2012, 9pm

No Group Work on MPs:

Copying code or allowing someone else to copy your code is cheating. All participants will receive a score of zero. So don't cheat! You are welcome to talk to the TAs and the professor during their office hours, or email us and ask questions, but the TAs should not be asked for help on the MP during lab hours.

Start on this MP as soon as possible after the exam, if not earlier. The TA will hold office hours on a regular basis, but the time you will have to wait for the TA's attention generally increases exponentially as the deadline approaches.

To handin your work follow these instructions.

1. Problem Statement

You are the engineer of a project that is in charge of building a suspension bridge over a river. In a nutshell, the bridge consists of a tower on the left bank, a tower on the right bank, and a bridge with a number of support (vertical) cables spanning the river. You are given the heights of the left and right towers. You are told how many support cables there are going to be across the bridge and at what distances they are from each other and from the towers. You therefore also know the total horizontal distance across the river. You are given the weights that are to be attached to each cable. You are constrained by the fact that you have only a certain amount of main (horizontal) cable to work with (L). Given these constraints, your goal is to calculate the optimum heights of the support cables so that you can construct the bridge with exactly the amount of cable that you have with you, and no more.
Look at the diagram in the next section for a better understanding of the system.

2. Math

To restate the elements of the system that you are given:

• The heights of the left tower: y_L (or y₀) and the right tower: y_R (or y_n+1)
  
• The number of support cables: n (n is 5 in the system above)
  
• The weights suspended from each of the support cables: W = [W₁; W₂; W₃; .. W_n] = [W₁; W₂; W₃; W₄; W₅]  - (W is a column matrix)
  
• The distances separating the towers and the cables: d = [d₁ d₂ d₃ .. d_n d_n+1] = [d₁ d₂ d₃ d₄ d₅ d₆]
• The amount of the main cable available to you: L
  

You're trying to calculate the heights of the support cables: [y₁ y₂ y₃ .. y_n] = [y₁ y₂ y₃ y₄ y₅]

If you resolve the forces acting on the system at the point where each support cable (in red) meets the main cable (in blue), you have a figure that looks like the following:

The weights suspended from the cables exert a force in the downward direction W_i. Since the system is in a stable state, we resolve the forces acting along the main cables into their vertical and horizontal components:

Since all the forces acting in the horizontal direction are equal, we can call all of them HF.

Also, from geometry (similar triangles):

Equating the forces acting in the vertical direction, we get:

Solving this system of equations, we arrive at the general equation for the system:

*

=

Proceeding from here as in Lab 4 , we arrive at the equation for the system in the matrix form: Dy = wp ; which is what we solve in the function find_xy.

3. Preparation

mp1_data

1. main
2. diff_L
3. find_xy (will be done in Lab 4)
4. xy_length
5. final_graphs

Note that these functions call each other a lot. So if you get an error while running the program, and MATLAB tells you that the error is at Line p, Column q, in file somefile.m, that needn't necessarily be the case. You might have to trace backward a little bit to find out where the error actually originated. Quite often, the error will be in the someotherfile.m that just called the function somefile, where the error was manifested. Sometimes you might have to go back a couple of steps. If you get totally lost, the best thing to do is start at the beginning, and go through all of your code, line by line and check everything.
Also, it goes without saying that the first line in all of your functions (aside from the comments that you write) must be the function definition. This isn't mentioned in the explanation of what goes into the functions, but duh!

It would also be a good idea to go through the entire MP before you start coding. Write down some ideas in pseudo-code on paper, or on this sheet if you've printed it out, before you start.

At the end of some of the functions, there is a small footnote titled Unit Test. These have instructions for testing that function in isolation. That means that if you test that function according to the instructions in that section, and something doesn't work right, then you should fix that before going on ahead.

When you begin coding, please note that you MUST put your name and netid, and a short (one or two lines will do) description of that function, at the beginning of the function, in the form of a comment. You must do this for each of the five files that you will be handing in.

4. Implementation

This is where you come in. You might want to write the functions in the order that they are explained below, rather than in the order that they're listed under the Preparation section, because of the manner in which the functions call each other.

1. find_xy
   find_xy is the function that we wrote in Lab 4. Flashback: we created the function find_xy that would simulate one run of the algorithm to find the optimum weights. By calling this function over and over again from within the function diff_L, and with the help of the functions you see here, we should be able to find the optimum heights of the cables such that they fit our predefined total length.
   For a detailed description of the function find_xy, refer to the text in Lab 4.
   Note: You can use the same code that we wrote in Lab 4. However, if you inserted a plot statement into the function find_xy in the lab, then take it out. The diff_L function will do the plotting of the cables for you.
   
   Unit test: If you did Lab 4, then this function should be fine. But to test it again, if you call it from the MATLAB prompt like so:
   
   >>[x,y] = find_xy(15)
   
   Then you should get the following values for x and y:
   

   x = 
     0  5.3931  10.3278  14.6056  20.0877  25.0048  30.1779
   y =
     30.0000  23.4883  19.8001  18.1607  18.3307  20.4408  25.0000
   

2. xy_length
   The purpose of the function xy_length is to calculate the sum of the distances between the points contained in the two vectors that are given to it as input. The function takes as its input two vectors x and y, and its return value is a single variable (let's call it total), that represents the sum of the Euclidean distances between each consecutive pair of points. Let's look at how we would approach this.
   Given two points with coordinates (x₁,y₁) and (x₂,y₂), the geometric (Euclidean) distance between them is given by the formula
   

   

   Therefore there are two things you need to note when you're writing your function xy_length. One is that you need to use a function that, starting with the vectors x and y of x-coordinates and y-coordinates, gives you the differences between consecutive values. The second is that you need to take the sum of this result over all such differences, so that you get the total length of the cable segments that you would need. (This is what you will return, and therefore what you will compare with L in the function diff_L.)

   Set the value of the return variable total to be the sum of the distances between each pair of successive cables (left tower to leftmost cable, between all the cables and rightmost cable to right tower).

   Unit Test: Calculate the distance between the points (0,0), (3,4) and (6,8):
   
   >> a = xy_length([0,3,6] , [0,4,8])
   
   You should get the answer as:
   

   a =
     10
   

3. diff_L
   
   The purpose of this function is to compute the difference delta between the target length L and the total length of the cables running from the left tower to each of the support cables, to the right tower. The input argument to this function is the guesstimate of the horizontal force, HF, and the return value is delta, the error in your approximation of the length of cables needed for you bridge.
   
   
   1. Make sure you allow global access to the Data vector. This function will create this vector by appending values to it during each iteration. For a detailed description, see 3.g to understand what the vector Data would look like, and see 4.a to understand why we want it to be global.
   2. Call the mp1_data.m file as a script. Though it would work the same, please don't copy all the code from mp1_data.m and paste it here. As long as you've downloaded the file to your home directory, you can run it like a script. It reads in the input data for you.
   3. Call the function find_xy on the horizontal force that you passed in. Each iteration of find_xy gives you two vectors of coordinates of the topmost points of the vertical cables (and the towers). Store them in two vectors; say x and y.
   4. Use these vectors to call xy_length, to find the sum of the lengths of the cables between the left tower, the vertical cables and the right tower. Store the return value of xy_length in a variable; say c1.
   5. Use the vectors x and y to make a plot of this iteration. Use blue solid lines with '*' markers. Since diff_L is called over and over, you'll see a visualization of the progression from the initial system based on the estimated horizontal force to the correct answer based on the amount of cable you have.
   6. Since there's no major processing involved in this MP, the visualization that was talked about in part 3.e will flash before you too fast for you to make any sense of it. To avoid this, slow it down by inserting a pause of a tenth of a second.
   7. Update the global variable Data that you're using to store the horizontal force and the calculated length for each of the iterations. Basically, append (attach at the end) a vector with two elements HF and c1, to the already existing global vector Data. When you do the update, make sure that each new iteration adds data to a new row, and that each row in the Data vector has the horizontal force and the length for one iteration. What this does is it takes the Data vector that existed from the previous iteration, and adds the data from this iteration to it. This kind of creates a timeline of horizontal forces and lengths. We'll use this later to plot how we approached the optimal answer.
      The Data vector should therefore be a vector with i rows and 2 columns; where i is the number of iterations we perform. The first element in each row (the entire first column) will contain the horizontal forces that you calculate at each iteration, and the second element in each row (the entire second column) will contain the lengths that you calculate at each iteration.
   8. We're interested in getting the difference between the current length that we need for all the cables (c1), and the amount of cable given to us (L), down to zero. So subtract c1 from L and make this the return value delta of the diff_L function. This delta represents how far away you are from the final answer.
      

   Unit Test: Call diff_L with the initial value of 15 like so;
   
   >>b = diff_L(15)
   
   You should get the answer as:
   
   

   b =
     3.0722
   

   While the plot should look like this:

   
    

4. main
   This is the function that you will actually call when you run the program as a whole, to see if everything is working. The function accepts one input parameter, HF and it returns two vectors x and y. The function definition line should look like this:
   

   function [x,y] = main(HF)

   That is, when you type in the following at the command prompt, the program should run and give you the output you expect, plus the graphs.
   
   >> [x,y]=main(15)
   
   
   1. You need to define a vector called Data to be global, to archive the horizontal force and the error in approximation of the total length of cables required. We're doing this so that at the end of everything we can plot a graph that tells us how we converged from our initial guess to the right answer. The reason we're making this variable global is that we want to be able to access it from inside a number of functions, and we don't want to go through the pain of passing it about as an argument and returning it and what-not. Just accessing one global structure is more elegant in this case, and much simpler! To enable a function to access and perhaps update this vector (in spite of not having been passed it as an argument), we must explicitly tell MATLAB that the data structure is global. Every function that wants to access the Data vector has to do this. Think of it as adding someone to your Buddy List on an Instant Messenger like MSN or AIM. Both of you have to add the other to your list so that you can communicate seamlessly. Same way, both (or all) functions have to say that they are associated with the vector Data. The functions are not communicating with each other in this way, but each of them is communicating with Data.
      Initialize the vector Data to be an empty vector. (This is just to start off with; we'll add history data to it as we go along.)
      
      
   2. Add a command that deletes all the open figure windows. This won't make any sense the first time you run it of course, but as you go from iteration to iteration, you want to minimize the clutter on your desktop.
      

      close all

   3. Find a function that, when passed the name of another function somefunc as a string, and an initial guess someguess, finds a root (or a zero) of somefunc near someguess, and returns that root. Give this function diff_L in place of somefunc, HF in place of someguess, and store the value returned in a variable true_HF.
      
      
   4. Call the function final_graphs to plot your progression from the initial guess to the optimal solution. Use the command,
      handles = final_graphs();
      If you want to test everything you've written till this point, before going on to write final_graphs, you can do so, you just won't be able to see all of the figures that you'll need in the end. Just leave this line out for the time being, and insert it when you have written final_graphs.m.
      
      
   5. Just so that you can present the final answer to the user, call the function find_xy with true_HF as it's argument, and store the values returned in vectors x and y.

5. final_graphs
   Since we're applying an iterative process to get to the optimal solution, we'd be interested in knowing how we approached that solution. That's what this function is going to do. It does not accept any input parameters and it has no return values. The function definition line should look like this:

   function handles = final_graphs

   1. final_graphs is definitely going to want to access the information in our global vector Data, so make Data visible by declaring it as global here too.
      
      
   2. Open a new figure window.
      

      figure

   3. Read in the data from mp1_data.m, the same way as above.
      
      
   4. Split up the vector Data into two vectors Forces and Lengths as follows:
      Create a vector Forces that has all the forces that were in Data. To do this, you need to access the first element in every row of Data.
      Create a vector Lengths that has all the lengths that were in Data. To do this, you need to access the second element in every row of Data.
      Yes, there is a simple way to do this!
      
      
   5. Plot a graph of Lengths vs. Forces (FYI: When you say 'plot y vs. x', you mean plot(x,y) ). Use red for the color specification, and a * for the symbol of the points.Use the command,
      handles(1) = plot(...);
      
      
   6. Set hold on.
      
      
   7. Create a vector called True_Length as follows:
      True_Length has one row, and the number of columns is equal to the size of the Forces vector.
      Each element of True_Length has its value equal to L; the allowed total length of the cables (which is in mp1_data.m and which you can access because you read it in.)
      
      
   8. Now plot a graph of True_Length vs. Forces, using blue for the color specification, and leave the line solid (the default setting).
      Adjust the axes (use the axis function) so that you can view the picture in better perspective; make the range of the X axis from the minimum Force to the maximum Force, and the range of the Y axis from the minimum Length to the maximum Length. Use the command.
      handles(2) = plot(...);
      Set the title of the graph to be 'Lengths vs. Forces'. Label the X axis as 'Forces'. Label the Y axis as 'Lengths'. Display a legend for the graph, showing which colors signify what quantity, and position it inside the axes boundary, obscuring as few points as possible ( use the legend function, with 0 as the value of the last argument).
      
      
   9. Open a second figure window, and plot a graph of the Lengths, in red with a * as the symbol. Use the command,
      handles(3) = plot(...);
      Set hold on, and then plot a graph of the Forces, in green with a * as the symbol. Use the command,
      handles(4) = plot(...);
      Plot the graph True_Length in blue, with the default solid line. Use the command,
      handles(5) = plot(...);
      Let the title of this graph be 'Length vs. iteration'. Label the X axis as 'iteration'. Label the Y axis as 'Lengths'. Display the legend of the graph as describing what each of the three colors mean, and again position it inside the axes boundary, obscuring as few points as possible( use the legend function, with 0 as the value of the last argument).
      

   Your should get the following plots (in addition to the plot by diff_L) when you run main with the initial guess of 20 for the horizontal force (HF):

   

   
    

Extra credit (5 points): Improving find_xy using spdiags. It's All So Empty...

1. Introduction to spdiags and Sparse Matrices

Now we're going to perform a little optimization. Notice that the coefficient matrix we create to solve the system of equations in find_xy can have a lot of zeroes if there is a large number of weights. This can bog down calculations, because MATLAB will spend the time to use each of those zeroes, even though they don't do anything, and they also take up space because MATLAB stores them as part of the matrix. The solution to this is to use something called a sparse matrix. Sparse matrices store only the non-zero elements; all other elements are assumed to be zero. This saves on storage space for matrices with lots of zeroes, and it also speeds up computation, since MATLAB will simply ignore the non-existent zeroes in the sparse matrix when doing matrix multiplication.

To create sparse matrices, we'll use the command "spdiags", which will replace diag. spdiags is a tad more complicated (but more powerful) than diag. The format for using spdiags is spdiags(D, c, m, n), where:

• D: matrix containing the diagonal columns.
• c: "control vector" of diagonal positions. c has as many elements as D has columns.
• m: number of rows in resulting matrix.
• n: number of columns in resulting matrix.

That is, spdiags will take each column of the D matrix and place it on the diagonal indicated by the control vector in an m x n sized matrix. Here's an example of how it works. First, we construct a matrix with the diagonal columns:

>> D = [1:5; 6:10; 11:15]'

D =

     1     6    11
     2     7    12
     3     8    13
     4     9    14
     5    10    15

Note the three columns of D. We're going to use the columns as diagonals in a 5 x 5 matrix:

>> A = spdiags(D, [0 1 -2], 5, 5)

A =

   (1,1)        1
   (3,1)       11
   (1,2)        7
   (2,2)        2
   (4,2)       12
   (2,3)        8
   (3,3)        3
   (5,3)       13
   (3,4)        9
   (4,4)        4
   (4,5)       10
   (5,5)        5

Figure 4.

Given D, c, m, and n, the resulting matrix A from spdiags is of m x n size. Each column of D is used to create a new diagonal. The first column is on the main diagonal of A, the second column is on the first upper diagonal, and the third column is on the second lower diagonal. How did they get there? Note the control vector c: 0 for the first value, corresponding to the main diagonal, 1 for the second value (first upper diagonal), and -2 for the last value (second lower diagonal). Thus, the c vector controls which diagonal the columns of D appear in A.

But wait, you may (should) be saying. How do the diagonals get created if the column of D isn't the same size as the diagonal it's being used for? How did the second column of D end up on the first upper diagonal of A if the column has five elements to the diagonal's four? The solution is that MATLAB will drop values from the D columns to make the diagonal fit. In the case of upper diagonals, the values are dropped starting from the top of the column matrix. Thus, the 6 from the second column of D was discarded, and 7 through 10 were used to make the diagonal. Similarly, in the case of lower diagonals, values are dropped starting from the bottom of the column matrix. For the second lower diagonal, 14 and 15 were dropped, while 11 through 13 were used.

2. Implementation

Now let's put all that newfound knowledge of spdiags to good use. Create a new function called "find_xy2" and store it in a file called "find_xy2.m". This function will take the same parameters as find_xy. You can copy most of the code over from find_xy. Instead of using multiple "diag" calls to create your coefficients matrix, now you will use a single call to "spdiags" to create a sparse version of your matrix. That is, the only changes that need to be made in your code for find_xy is to change how you create the matrix A (see section 2 Math) . You will need to make your own columns matrix D and control vector c, then feed them to spdiags along with the appropriate values for the matrix dimensions to get the same coefficients matrix as before, but in sparse form.

Enjoy!