First, focus on only one weight W_i in the system:

This is the equation that we solved in Lab 4.

find_xy is the function that we wrote in Lab 4. For a detailed description of the function find_xy, refer to the text in Lab 4.


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    7.5000   15.0000   22.5000   30.0000   37.5000   45.0000   52.5000   60.0000   67.5000   75.0000

y =

   70.0000   56.6500   45.8000   37.4500   31.6000   28.2500   27.4000   29.0500   33.2000   39.8500   50.0000

Step 1 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 has two input parameters x and y, which hold vectors,  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. See lecture 3 slides 9-18.

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

Step 2 diff_L

The function diff_L that we'll be using here is a lot like the function diff_L that you saw in lecture 6 slides 3-13.
The purpose of this function is to compute the difference, delta_L, between the target length L and the total length of the springs running from the left tower to the right tower. The input argument to this function is the guesstimate of the horizontal force, HF, and the return value is delta_L, the error in your approximation of the length of springs needed for you bridge. Your function should perform the followings:

1. You need to create a variable called Data. Use the global command.
   
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.
   
3. Call the function find_xy using the horizontal force HF. A call of find_xy returns two vectors of coordinates of the topmost points of the vertical cables (including the towers). Store them in two vectors; x and y.
   
4. Use these vectors, x and y as input arguments and call xy_length, to find the sum of the lengths of the springs between the left tower and the right tower. Store the return value of xy_length in a variable; calc_length.
   
5. Use the vectors x and y to make a plot of this iteration. Use blue solid lines with circle 'o' markers. Use Matlab help to find how to use the MarkerEdgeColor, MarkerFaceColor and Markersize properties in the plot function. Use 'k' for MarkerEdgeColor, 'g' for MarkerFaceColor and 10 for MarkerSize. 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. Your plot command will look like this:
   
   handle = plot(x,y,'b-o',  ____________________________);
   
   
   
6. Use the pause command for a tenth of a second.Use Matlab help to learn about the pause command.
   
7. Update the global variable Data that you're using to store the horizontal force , the calculated length and the plot handle for each of the iterations. Basically, we want to append three values (HF , calc_length and handle ) as the next row of Data during each iteration (i.e. during each call to this function). Each of the three values should be added as a new row to Data. Take the Data that existed from the previous iteration, and adds the HF , calc_length  and handle from this iteration to it. This will create a timeline of horizontal forces and lengths. We'll use this later to plot how we approached the optimal answer.
   The variable Data should therefore be a matrix with i rows and 3 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 springs lengths that you calculate at each iteration and the third element in each row( the entire third column) will contain the plot handle at each iteration.
   
   After i iterations, Data should have this structure:
   

       Data = [ HF1  calc_length1  handle1; 
                HF2  calc_length2  handle2;
                HF3  calc_length3  handle3;
                   ... ...
                HFi  calc_lengthi  handlei]
   

   
   
8. We're interested in getting the difference between the calc_length that we need for all the springs, and the target length given to us (L), down to zero. So subtract calc_length from L and make this the return value delta_L of the diff_L function. This delta_L represents how far away you are from the target value L.
   
   Unit Test: Call diff_L with the initial value of 20 like so;
   
   >>b = diff_L(20)
   
   You should get the answer as:
   
     b =
   
       6.3740
   
   Run the mp1_demo program to check your plot.
   
   

Step 3 main

This is the function that you will actually call when you run the program as a whole. The function accepts one input parameter, HF, and it returns two vectors x and y, and two values, true_HF and max_F. The function definition line should look like this:

function [x,y,true_HF,max_F] = 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 , true_HF, max_F]=main (10)

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 springs. We're doing this so that at the end of our program we can plot a graph that tells us how we converged from our initial guess HF 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.
   
2. 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. )
   
3. Add commands that closes all the open figure windows, and clear the command window. 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
   clc
   
   
4. Use fzero to find a root (or a zero) of diff_L using an initial value of HF, and store the value returned in a variable true_HF.
   
5. Call the function find_xy with true_HF as it's argument, and store the values returned in vectors x and y.
   
6. Call the function final_graphs to plot your progression from the initial guess to the optimal solution. 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.
   
7. Call the max_force function (you will write in step 5).
   
   max_F = max_force(y, true_HF);
   
   

Step 4 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, but it returns a list of two handles of plots for grading. The function definition line should look like this:

function handle = 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 the matrix 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 (i.e. the first column).
   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 (i.e. the second column).
   Yes, there is a simple way to do this!
   
5. Plot a graph of Lengths vs. Forces (FYI: When we say 'plot y vs. x', we mean plot(x,y) ). Use red for the color specification, and a * for the symbol of the points.
   
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 length of the Forces vector.
   Each element of True_Length has its value equal to L; the allowed total length of the springs (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 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.
   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. Here is how the version of the legend command that you need to use works: suppose you currently have a figure window opened where you have plotted two graphs, of vectors y and z (on the vertical axis), say versus the same vector x (on the horizontal axis), and suppose the plot of y versus x has a different format (color-linestyle-marker) than that of z versus x; then legend('this is y','this is z') , or simply legend('y','z') , has the effect of adding a legend box associating the first string argument (such as 'this is y' or 'y') to the format of the first plot and the second string argument ('this is z' or 'z') to the format of the second plot, so that one who looks at the figure knows that, for instance, the red stars represent values of y (versus x) and the blue circles values of z (versus x); and this works for any number of plots on a figure window, not just two. Like many Matlab commands, legend may take further optional input parameters (as pairs of strings indicating the type of argument and its value), which come after the parameters representing strings to be placed in the legend box. For obscuring as few points as possible on the figure, you need to use the pair 'Location' - 'Best'; for example, in the situation above, you would write
   legend('y','z','Location','Best').
   
9. Open a second figure window.
   
10. Plot a graph of the Lengths, in red with a * as the symbol. Save the handle of this plot as handle(1).
    
11. Set hold on.
    
12. Plot a graph of the Forces, in green with a * as the symbol. Save the handle of this plot as handle(2).
13. Plot the graph True_Length in blue, with the default solid line.
    Let the title of this graph be 'Lengths 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.
    
    

Step 5 max_force

The max_force function will determine the maximum magnitude of the F_i for any of the springs. This function should have two input parameters: h and HF . See the equation (8) in the Math section 3 above.

1. The function definition line for max_force is:

     function max_F = max_force(h,HF)

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.

3. Use the Matlab function namde max to find the maximum force and assign the result to the output variable named max_F .

Unit Test:


>> [x,y,true_HF,max_F] = main(10)
x =
         0    7.5000   15.0000   22.5000   30.0000   37.5000   45.0000   52.5000   60.0000   67.5000   75.0000
y =
   70.0000   57.6780   47.6296   39.8548   34.3535   31.1258   30.1716   31.4910   35.0840   40.9505   50.0000

true_HF =
   16.4939

max_F =
   31.7234

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

How did the columns of D turn into the diagonals of A? Let's break it down some more, using the same data as above:

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 3 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.

To test your new function by running the mp1_check code with option 7.

1. Edit the mp1_data.m file and change N = 9; to N=1000; Save the file.
   
2. Edit the main function and comment out the old line of code where you called find_xy and then write a new line of code calling find_xy2. Save the file.
   
3. Edit the diff_L function and comment out the old line of code where you called find_xy and then write a new line of code calling find_xy2. Save the file.
4. Test your program by typing the following:
   >> [x,y,true_HF,max_F]=main(1500);
   
1. At the Matlab prompt type the following,
   >> true_HF
   true_HF =
      1.6223e+03
   
   >> max_F =
      3.3514e+03
   
   

If you cannot get find_xy to work correctly then go back and change main and diff_L so that they call find_xy. Otherwise if both find_xy2 and find_xy work correctly, you don't have to change main or diff_L. The TAs may use a different mp1_data.m file so it doesn't matter what values you have in this file when you handin your work.

Part 5: Grading Scheme