DAY21
Thursday, 28 June 2012


Today is a LAB day.

Your project information is due today. You must sent a title and a one sentence description to Detelina by email.

We have a LAB #6 to be done during class time, for credit. When you have completed the lab, please let Detelina know so that she can approve and record your work.

Detelina will not be available for office hours this coming Monday. She will make it up with office hours on Tuesday, at the usual time.


LAB#6:

Rabbits and foxes live on an island. The foxes hunt the rabbits, but the rabbits simply eat grass. Over a year, the rabbit population would tend to increase at the rate of 2 new rabbits for every existing rabbit, while the fox population would tend to decrease at a yearly rate of 10 deaths for every existing fox.

However, if any rabbit meets any fox over the year, this decreases that rabbit's chance of living and improves the fox's chance. There are R * F possible encounters of this kind, and we will estimate the relative cost of such an encounter to the rabbit population by - 0.004, and the benefit to the fox population as + 0.003.

Under these assumptions, the fox and rabbit populations may be estimated by the differential equations:

        dR/dt =    2 * R - 0.004 * R * F
        dF/dt = - 10 * F + 0.003 * R * F
      

We can try to approximate these equations by using a small time step DT, so that dR/dt is approximated by ( R(new) - R(old) ) / DT, with a similar approximation for dF/dt. If we think about storing the successive values of T, R and F as arrays, then here is what we are doing:

        T(i) = i * DT
        R(i) = R(i-1) + DT * (    2 * R(i-1) - 0.004  * R(i-1) * F(i-1) )
        F(i) = F(i-1) + DT * ( - 10 * F(i-1) + 0.003 * R(i-1) * F(i-1) )
      

We start the problem at time T(0) = 0, with rabbit population R(0) = 4,000, and fox population F(0) = 100. Now we want to take N equal steps of size DT, to get us from time 0 to time 5 years, using the approximate equations to fill in the values of our arrays.

Write a program that takes N steps to go from time = 0 to time = 5, filling in the values of the vectors T, R and F. Try using a value of N that is 10,000. Think about it: this means your arrays should be declared to have size N+1!

After you have computed T, R, and F, print all three arrays. Print T(i), R(i), and F(i) on one line, so that we can feed the data to gnuplot. Create a data file with the command

        ./a.out > data.txt
      

Using gnuplot, make the following two plots from your data:

        1) T versus R and T versus F
        2) R versus F
      

When you are done, show Detelina your results.


HOMEWORK #6 (must be turned in by next Thursday):

Homework 6.1: Every year, 1/10 of the people outside of California move in; in the same year, 2/10 of the people inside California move out!. Suppose this statement is true. Let's ask ourselves what would happen if this statement remains true for the forseeable future. Assume the population of California is 37 million people, and that the US population is 300 million. Let A be the matrix

        8/10  1/10
        2/10  9/10
      
and let p be the vector
        37
       263
      
Then p[0] is the California population and p[1] is the non-California population. Moreover, the populations next year will be
        pnew[0] = 8/10 p[0] + 1/10 p[1]   0.8 Californians stay put, 0.1 NonC people move in
        pnew[1] = 2/10 p[0] + 9/10 p[1]   0.2 Californians move out, 0.9 NonC people stay put.
      

If the same pattern holds year after year, we can model the population over 20 years by:

        Set matrix A, and initialize vector p
        for i = 1 to 20
          pnew = A * p
          p = pnew
        end
      
where we must use the temporary vector pnew to collect our new results.

Write such a program, which models the exchange in population between California and the rest of the US. Print your population estimates p[0] and p[1] for each of the 20 years. Your results, if correct, will show a definite trend or pattern.

Homework 6.2:

If we have an m1 by n1 matrix A, and an m2 by n2 matrix B, the we are allowed to compute the product C=A*B, but only if n1 = m2; when the multiplication is allowe, the result C will have dimension m1 by n2. (In the most common case, all the dimensions are equal.) We are allowed to multiply a 3x2 matrix by a 2x7 matrix, in which case we will end up with a 3x7 result.

To compute the product C = A*B, we need a triple for loop:

        for ( i = 0; i < m1; i++ )
        {
          for ( k = 0; k < n2; k++ )
          {
            c[i][k] = 0.0;
            for ( j = 0; j < n1; j++ )
            {
              c[i][k] = c[i][k] + a[i][j] * b[j][k];
            }
          }
        }
      

Suppose the entries of the matrices A and B are defined by:

        A[I][J] = min ( I, J ) + 1;
        B[I][J] = I - J;
      
I am using C notation here, so I and J indexes start at 0, and thus the first entry A[0][0] = 1 and B[0][0] = 0.

Write a program which sets up a matrix A of dimension 4x3, and a matrix B of dimension 3x5, computes the matrix product C=A*B, and prints out C.

Homework 6.3: Suppose we are given the following matrix A:

        -2  3  0  1
         3 -6  3  0
         0  3 -6  3
         1  0  3 -2
      
and right hand side vector b:
        18
         6
        24
        12
      
One way to solve a (small) linear system A*x=b is to compute the inverse matrix (call it "C") and solve for x by multiplying x=C*b. For our matrix A, the inverse C is defined by the following formula
        C[I][J] = abs ( I - J ) / 6
      
Write a program which defines the right hand side vector b, evaluates the inverse matrix C and computes the solution x by multiplying x=C*b. Print your result x. (If you want to check your work, your result x is correct if A*x=b).

Homework 6.4: (Graduate students only!) Let A be the following 4x4 matrix:

        2 -1  0  0
       -1  2 -1  0
        0 -1  2 -1
        0  0 -1  2
      
This is an example of the NxN second difference matrix which is used to approximate second derivatives. For any N, the inverse of this second difference matrix can be computed by:
        C[I][J] = (I+1) * (N-J) / (N+1)   if I <= J
        C[I][J] = (J+1) * (N-I) / (N+1)   if I >  J
      
(Note that I have already adjusted this formula for C indexing. Therefore, for N=4, C[0][0] should be 4/5.) Write a program which evaluates C for the case N = 4; compute the product A*C and print it out. (Your result is correct if it is the 4x4 Identity matrix.)

For each of your homework programs, turn in a copy of the program, and a copy of the output.


Programs we might discuss:


Last revised on 28 June 2012.