Project 2 Eigenvalues in Plankton Blooms

Description

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]

^dN_dt = (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

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 X_n(x)T_n(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 X_nT_n:

1

• _p

u(x; t) = A_n 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 A_n are determined by the initial conditions c(x; 0). The population will grow if the growth rate K > ₁, where ₁ 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 ₁.

• Discretization

In this section, we will show how discretization of the eigenvalue problem leads to our result. We can discretize the operator X⁰⁰, 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 x_k = 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 x₀ = 0 = x_n+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 X⁰ and (10) to your approximation for X⁰.

3

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

X⁰⁰= X

at the points x₁; x₂; x₃ and recognize the pattern.

3. Write a program that creates the A_n matrix.

• QR Algorithm

We now use the QR Algorithm to nd the eigenvalues of several A_n. We extract the smallest eigenvalue and compare it to the actual value ₁ = _l ². 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 fv_igⁿ_i=1 it is possible to create an orthonormal basis fu_igⁿ_i=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 fv₁; v₂g we start o by setting

v₁

^{u1 =}

jjv₁jj

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

v₂ = v₂^jj + v₂^?

Then the vector

^v2^{? = v}2 ^projspanfu₁g^(v2^{) = v}2 ^(u1 ^v2^)u1

is orthogonal to u₁. We divide v₂^? by its length to get the second vector u₂ of the orthonormal basis.

v^?

u₂ = ²

jjv₂?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 v_n on to the subspace spanned by fu_i gⁿ_i=1¹ 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 A₃. Verify your algorithm works by showing the dot product (u_i; u_j) = 0 or numerical zero for i 6= j. Also, check that each vector has length 1.

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 A_n 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 A₀ = A

Step k: Set A_k = Q_kR_k then set A_k+1 = R_kQ_k

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

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

2. Apply the QR Algorithm 1000 times to A_n 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 ₁. (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 = c₁v₁ + c₂v₂ + + c_nv_n

for c_i 2 R.

If we apply A to this equation we obtain the equation

Ax = c₁Av₁ + c₂Av₂ + + c_nAv_n

and applying the eigenvalue relationship (Av_i = _iv_i) we obtain the equation

Ax = c_{1 1}v₁ + c_{2 2}v₂ + + c_{n n}v_n:

We assume that the largest eigenvalue is distinct from the others i.e., j ₁j > j ₂j : : : j _nj. Then we can factor out ₁.

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 A_n of size 10 10, 20 20 and 100 100 compute the error between this value and ₁. (Hint: what is the relationship between the eigenvalues of a matrix A and A ¹)

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.

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