Project 2 Eigenvalues in Plankton Blooms

Abstract

Finding eigenvalues of large matrices is a problem that comes up frequently in applications. The strategy of nding such that the equation (A I)v = 0 is nonsingular becomes too costly as the size of n grows. We must resort to other methods for computation of eigenvalues in large systems. In this project we analyze two methods: QR Algorithm and the Power Method. QR iterations takes the original matrix A and nds a sequence of matrices similar to A that converges to an upper triangular matrix. Once we have a the upper triangular matrix T that is similar (similar means the two matrices have the same eigenvalues) to A, then we are able to simply read o the eigenvalues from the diagonal of T. Our power method is di erent in the sense that it is used to compute the extreme eigenvalues i.e., the largest and smallest eigenvalues of a matrix. We use a model describing the size of plankton population as a case study of how the eigenvalue problem can arise in applications. In this project we derive a condition based on the smallest eigenvalue of the system for the population to grow.

    • Introduction

The motivation for this project comes from the growth-di usion equation of plankton populations. Phy-toplankton provide 50% of the world’s oxygen and allow some of the atmosphere’s carbon dioxide to be absorbed into the ocean. It is thus important for us to know when a plankton bloom will be growing or decaying. This problem was studied in an important paper of Kierstead and Slobodkin [1]. The simplest growth model of population of plankton is the following [1]

dNdt = (K    )N


where K is the growth rate, N is the rate of population loss by di usion, and N is the population size. For the simple growth model, growth of a plankton population will occur if K > . However, since the

1
population is modeled with di usion, we instead use the di usion term in the simple growth model arriving at a partial di erential equation that describes growth in space and time

@c
@2c



= D

+ Kc
(1)


2


@t
@x



for the concentration c(x; t) of phytoplankton. Here, we sum the di usion term with the term Kc which represents the linear growth term of the plankton bloom. The concentration is assumed to be zero at the ends of an interval of length l, so we have the boundary conditions u(0; t) = u(l; t) = 0.

The di usion equation is very well understood and easy to work with, so we make the transformation

c(x; t) = u(x; t) eKt
(2)
and get the di usion equation

@2u



@u

= D


(3)

@t

@x2


Introduction Questions

1. Show that the substitution of (2) into (1) reduces the model to (3).

    • The Exact Solution

To solve equation (3) when the di usion coe cient D is constant, we will use separation of variables. We assume a solution of the form u(x; t) = X(x)T (t). Substituting into equation (3) leads to the two equations

T0+ DT =0
(4)
X00+ X=0
(5)
where  are our eigenvalues and the solution to the di usion equation. The boundary conditions are
X(0) = 0 = X(l)

This gives us the eigenvalues and eigenfunctions

n;  Xn(x) = sin  p

x  ;  Tn(t) = e
D nt

n


Each product Xn(x)Tn(t) is a solution of (3), and satis es the boundary conditions. The general solution of the di usion equation (3) is obtained by taking the weighted sum of all these solutions XnTn:

1
    • p

u(x; t) =    An sin(    nx)e D nt
n=1

The concentration is, from (2),

1




X
p

x)e(K  D n)t

c(x; t) =
An sin(
n

(6)
n=1

The constants An are determined by the initial conditions c(x; 0). The population will grow if the growth rate K > 1, where 1 is the smallest eigenvalue known to be
1 =


2



(7)

l


We will show how discretization of the di usion equation (3) leads to this result.

2
Section 2: Questions

    1. Solve the ODEs obtained by the separation of variables (4) and (5) subject to the boundary conditions to obtain the eigenvalues and eigenfunctions stated above.

    2. Explain why taking a weighted sum of solutions is also a solution to the di erential equation. (Hint: can you write (3) as a linear operator)

    3. Use (2) to derive the solution for the concentration from the solution of u(x; t).

    4. Explain why the plankton population will grow if the growth rate K > D 1.

    • Discretization

In this section, we will show how discretization of the eigenvalue problem leads to our result. We can discretize the operator X00, to obtain the matrix eigenvalue problem, which approximates the eigenvalue problem (4). Discretize the interval (0; l) into n + 1 parts of length h = l=(n + 1). Then the second centered di erence approximation of the second derivative is [2]:

X00(x)
X(x  h)  2X(x) + X(x + h)



h2


Let xk = X(hk). Then the eigenvalue problem (4) is approximated by the matrix-vector eigenvalue problem


2
1   2

1








3 2
x2
3

2 x2
3



6
2
1











7

x1



x1



















6

7

6

7

1


1
.
2

.
.
1
.
.





x3



x3



6



.
.


.


.


7

.



.



















6
.
7
=
6
.
7
(8)

2
6













7
6

7

6

7

h

6





..
.
..
.
..
.
7 6
.
7

6
.
7



6













7
6
.
7

6
.
7



6








.




7

.



.




6









..
2
1
7
6
xn  1
7

6
xn  1
7



6












7 6

7

6

7



6













7
6

7

6

7



6











1   2
7 6
xn
7

6
xn
7



6













7
6

7

6

7



4













5
4

5

4

5

Note that the equation is for the interior points, as X(0) = 0 = X(l), corresponding to x0 = 0 = xn+1. Let
Bn be the n  n \
1; 2;  1" matrix above. Then the eigenvalue problem can be written as










+ 1

2






















n

Bn x =  x
















l








We de ne An
n+1

2
Bn: For simplicity assume that l = 1 for the remainder of the project. The eigenvalues


l




of this matrix
approximate the eigenvalues for the plankton model. As n
! 1
the approximations become


























increasingly accurate.





















Section 3: Questions



















1. The derivative of a single valued function f at a point x can be de ned as either








f0(x) = lim
f(x + h)
f(x)

f(x + h)

f(x)
(9)














h













h
!
0



h












or as




































f(x)
f(x
h)

f(x)  f(x
h)










f0(x) = lim






:
(10)









h





h











h
0























!


















Where the approximation is increasingly accurate for small h. The rst approximation is known as a forward di erence approximation and the second is known as a backward di erence. Derive the approximation for the second derivative by applying (9) to X0 and (10) to your approximation for X0.

3
    2. Derive the system of equations that leads to the matrix equation. To show this explicitly write out the equations for

X00= X
at the points x1; x2; x3 and recognize the pattern.

3. Write a program that creates the An matrix.

    • QR Algorithm

We now use the QR Algorithm to nd the eigenvalues of several An. We extract the smallest eigenvalue and compare it to the actual value 1 = l 2. Before applying the algorithm we study how the algorithm works. To understand this algorithm we study Gram-Schmidt algorithm and the QR factorization.

Gram-Schmidt algorithm

Given a set of linearly independent vectors fvigni=1 it is possible to create an orthonormal basis fuigni=1. To illustrate this let us analyze low-dimensional cases and build up an algorithm. If we have one vector in our set fug, then to construct an orthonormal basis we can set
v
    • = jjvjj

For the case when we have two vectors in our set fv1; v2g we start o    by setting
v1
u1 =

jjv1jj

We now have to nd a vector that is orthogonal to v1. We can write a vector in terms of its it’s components parallel and perpendicular to the line spanned by u1.

v2 = v2jj + v2?

Then the vector
v2? = v2   projspanfu1g(v2) = v2   (u1  v2)u1

is orthogonal to u1. We divide v2? by its length to get the second vector u2 of the orthonormal basis.
v?
u2 =    2

jjv2?jj

For general number of vectors we have the following algorithm



n  1 (uk; vn)

u1
= v1
;  un = vn





uk




uk
jj
2




X jj





k=1

where we are subtracting o the projection of vn on to the subspace spanned by fui gni=11 to nd a vec-tor orthogonal to the previous vectors. The nal step is to divide each vector by its norm to obtain an orthonormal basis.

    1. Write a program that takes in a matrix of size 3 3 and uses the Gram-Schmidt algorithm to nd an orthonormal basis for the column space of the matrix. Apply this program to A3. Verify your algorithm works by showing the dot product (ui; uj) = 0 or numerical zero for i 6= j. Also, check that each vector has length 1.

4
QR Factorization

The Gram-Schmidt process takes in our \old" basis and produces a \new" orthonormal basis. This transformation can be described by a change of basis matrix R. We can write this as


v1   : : :
vn  = un : : : un  R
Here Q is an orthogonal matrix by
construction and R is an upper triangular matrix.




1. Use the    le cgs.m to obtain the QR factorization of An matrices for size 3    3 and 4    4.

QR Algorithm in Practice

The following algorithm is used to compute the eigenvalues of a matrix A.

Step 0: Set A0 = A
Step k: Set Ak = QkRk then set Ak+1 = RkQk

Here the matrix Ak+1 is similar to the matrix Ak and under certain conditions Ak converge to an upper triangular matrix (or diagonal if A is symmetric).

    1. Show that the matrix Ak+1 is similar to the matrix Ak. Recall a matrix A is similar to a matrix B if they have the same eigenvalues.

    2. Apply the QR Algorithm 1000 times to An of size 10 10, 20 20 and 100 100 give the smallest eigenvalue for each (smallest element on the diagonal) and compute the error between this value and 1. (Hint: you can extract the diagonal of a matrix A using the diag command)

The Power Method

We now analyze the extreme eigenvalues via the power method. We begin by giving a description. If you assume that there is an eigenbasis (there is for our problem). Then we can write any vector as

x = c1v1 + c2v2 +    + cnvn

for ci 2 R.
If we apply A to this equation we obtain the equation

Ax = c1Av1 + c2Av2 +    + cnAvn

and applying the eigenvalue relationship (Avi =    ivi) we obtain the equation

Ax = c1  1v1 + c2  2v2 +    + cn  nvn:

We assume that the largest eigenvalue is distinct from the others i.e., j 1j > j 2j : : : j nj. Then we can factor out 1.



Ax =  1  c1v1
+ c2  1 v2 +   + cn  1 vn







2



n


Repeated applications of A give us the equation
















Anx =  1n  c1v1 + c2

1
n
v2
+   + cn    1
n











vn






2








n


As n ! 1 we have
j

n

















will go to zero for j 6= 1.











1












    1. You are given the code that gives you the largest eigenvalue using this method power method.m. Use this code to nd the smallest eigenvalues of An of size 10 10, 20 20 and 100 100 compute the error between this value and 1. (Hint: what is the relationship between the eigenvalues of a matrix A and A 1)


5
Comparison

Compare the two methods. Which one did you like more?

    • References

References

    [1] H. Kierstead and L. Slobodkin. The size of water masses containing plankton blooms. J. Marine Res., 12:141{147, 1953.

    [2] David Kincaid and Ward Cheney. Numerical analysis. American Mathematical Society, Providence, RI, third edition, 2002. Mathematics of scienti c computing.

    [3] Nasa.gov. Chlorophyll a concentration, 2003. [Online; accessed April 03, 2017].

    • Project Guidelines

Your group will submit your project write-up on Canvas to the appropriate \Project Assignment" (you can nd these under the \assignments" tab in Canvas). Adhere to the following guidelines:

Do not put o nding a group (you must works in groups of 2-3). You should have a group set-up within one week of the project assignment due date.

Submit your project in a pdf format and submit ALL code used for your project (.nb les for Mathematica, .m les for Matlab, .py or ipynb for python). Code in Matlab, Mathematica, Maple, R, C, C++, Python or Julia is acceptable (Matlab is recommended). Code in Microsoft Word or Excel (or any other spreadsheet program) is not acceptable. All other languages need instructor permission (please ask as soon as possible).

Code may be included in the appendix if you wish. DO NOT submit anything on Canvas as a .zip le. Contents of .zip les will not necessarily be graded.

Have only ONE group member submit the project. Having multiple people in your group submit the project to Canvas will result in multiple grades, and we will take the LOWEST one.

Include the names and recitation section numbers of all group members working on the project on the cover page of the report.

When you submit the report to Canvas, please include each group member’s information (name, student number, and section number) in the comments. This allows us to quickly search for a student’s report.

Your report needs to accurately and consistently describe the steps you took to answer the posed questions.

This report should have the look and feel of a technical paper. Presentation and clarity are very important.

Here are some important items to remember:

Remember: you are to submit a complete report for this project. Documents submitted with numbered responses will be severely penalized.

Labs must be typed, including all equations (part of your learning experience is to learn how to use an equation editor). An exception can be made for lengthy calculations in the appendix, which may be hand written (as long as they are neat and clear), and minor labels on plots, arrows in the text, and a few subscripts.

6
Write your report in an organized and logical fashion. Section headers such as Introduction, Back-ground, Problem Statement, Calculations, Results, Conclusion, : : : are not mandatory but are highly recommended. They not only help you write your report, but help the reader navigate your paper.

Start with an introduction that describes what you will discuss in the body of the report. A brief summary of important concepts used in your discussion could be helpful here as well. Always introduce relevant equations that will be used or discussed in the report.

Always include units in your answers

DO include and label any plot that supports your conclusions or gives you insight in your investigations (these should be found in the body of your report). However, DO NOT use screen-shots of your  gures.

DO NOT include printouts of computer software screens or code snippets. You simply need to state which software you used in each step and what it did for you.

DO include the results of any calculations in the main body of your report. DO NOT include lengthy calculations in the main body of your paper (calculations should be included in a labeled appendix and should be referred to in the main body).

Your report does not have to be long. You need quality, not quantity, of work. Do not omit any important piece of information, but do not feel obligated to add any extras.

Summarize what you have accomplished in the conclusion. No new information or new results should appear in your conclusion. You should only review the highlights of what your wrote about in the body of the report.




































7