%!TEX output_directory=tex_output
% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
%
% This template is to be used with the class file "simiode.cls".
%
% The latest version of this class and template files
% can be found at
%
% in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with TEACHER VERSION, but also full author identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%*********** Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{xcolor}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\usepackage{verbatim} % for verbatim code environ
\usepackage{xspace} % for smart spacing in macros
\usepackage[justification=centering]{subcaption}
\renewcommand{\baselinestretch}{1.3}
% macros
\newcommand{\M}{\textit{Mathematica}\xspace}
\newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\bv}[1]{\mathbf{#1}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\bF}{\mathbf{F}}
\newcommand{\bG}{\mathbf{G}}
\newcommand{\br}{\mathbf{r}}
\newcommand{\bq}{\mathbf{q}}
% INSTRUCTOR COEFFICIENTS
\newcommand{\version}{I\xspace}
\newcommand{\period}{25\xspace}
\newcommand{\totalmoose}{1,600\xspace}
\begin{document}
\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote.
% See above for instructions on toggling from Teacher Version to Student Version.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\
Interdisciplinary Modeling on Isle Royale}
% header titles left / right
\markboth{Interdisciplinary Modeling on Isle Royale}{Interdisciplinary Modeling on Isle Royale}
\author{Steven Morse, Brian Allen, Stanley Florkowski\\
United States Military Academy\\
West Point, NY 10996 USA}
\makeStitlePDFLaTex
\begin{abstract}
The primary aim of this project is to draw a connection between differential equations and vector calculus, using population ecology modeling as a vehicle. This setting allows us to also employ multivariable optimization as a means of model fitting and multivariable integration in the context of density functions. The project is designed for advanced first-year undergraduates in a multivariable and vector calculus course, with a background in differential equations.
\end{abstract}
\section*{SCENARIO DESCRIPTION}
\subsection*{Introduction}
This is a group assignment. Your primary submission is a hard copy report that should not exceed four (4) pages (not including appendices and title page). This report will include graphs reflecting your work.
You should include any additional supporting graphs, equations, or computations in the appendices. An appendix is NOT a place to put printed computer code. If you think your code is important, you may include it in an appendix, but you should ensure it includes annotations describing what you are doing. You will also digitally submit \M and/or \textit{Excel} files you develop via email. If you think a computation or equation is important for the reader, please ensure it is in the main report or appendices and not just in your digital file.
Ensure all work is logical, neat, and organized. Submit a cover sheet and document any assistance you receive. Doing the computations correctly is important, but analyzing, communicating, and reflecting the relevance of your results is also critical.
In addition you will give a 10 minute presentation with time afterward for questions to your instructor reporting your results. Your presentation will be no more than five (5) slides (one slide per task, plus any background/introductory material).
\paragraph{Goals:}
\begin{enumerate}
\item Students clearly communicate their mathematical model and solution methods in written form.
\item Students explore aspects of the problem from a different disciplinary perspective to supplement their mathematical modeling.
\item Students develop and solve a mathematical model for a real world problem using multivariable calculus.
\item Students confront the ambiguity of problem solving and understand the necessity of making assumptions in formulating a mathematical model and the impact these assumptions have on their solutions.
\item Students apply technology learned in the course to analyze and solve a multivariable calculus model and develop an appreciation for how technology enhances their problem solving capabilities.
\item Students connect the concepts learned in multivariable and vector calculus to the knowledge previously obtained in modeling with differential equations.
\item Students practice collaboration on a technical project.
\end{enumerate}
%%%%%%%%%%%%%
% PROJECT
%%%%%%%%%%%%%
\section*{Background}
\textbf{The following scenario is based on recent, real events.} (See \cite{Milot}.)
Isle Royale National Park (IRNP) is an island in Lake Superior, near Northern Michigan, which is 45 mi long and 9 mi wide. There is a predator-prey relationship of moose and wolf on the island. Researchers have monitored this dynamic since the 1950s, observing (until recently) relatively stable population trajectories.
However, the wolf population dropped down to only two wolves in January 2018 (male and female), resulting in an explosion of the moose population. This is having negative second order effects on the Park: overpopulation of moose has led to overfeeding, the vegetation change has affected stream flows, which affected riverine ecological balance and even infrastructure damage on the island.
Researchers working with the Isle would like to see if they can model the past dynamics of the system. If so, they will assess whether they can return the population to equilibrium through careful human interventions (e.g., bringing more wolves to the island, controlled hunting).
You are an engineer officer assigned to the U.S. Army Corps of Engineers (USACE) at the District Headquarters in Detroit, MI. The Federal Parks Service has requested USACE aid in assessing the damage to the park, which includes road and facilities damage. The District Commander assigns you to temporary duty (TDY) as the liaison officer (LNO) between USACE and the IRNP.
While not part of your official duties, you volunteer to collaborate with the team of researchers who are trying to model the population dynamics of the park.
%%%%%%%%%%%%%
% TASK 1
\section*{TASK 1: Exploring the Predator-Prey Model}
We would like to model the population dynamic of moose and wolves on the Isle with a predator-prey model,\footnote{Reference Zill \cite{Zill} Chapter 10 (p. 109, 417-418).}
\begin{align}
\begin{split}
x' &= P(x,y) = -ax+bxy, \\
y' &= Q(x,y) = dy-cxy
\end{split}
\end{align}
where $x(t)$ represents the size of the predator population at time $t$, $y(t)$ represents the size of the prey population at time $t$, and $a,b,c,$ and $d$ are parameters of the model. Recall $d$ is the growth rate of the prey, $a$ is the death rate of the predator, and $b$ and $c$ measure the effect of interactions between the two species.
Define
\begin{equation}
\bG = \langle P, \ Q \rangle, \qquad \text{and} \qquad \br(t) = \langle x(t), \ y(t) \rangle,
\end{equation}
Here $\br(t)$ is a vector equation of a space curve representing a solution trajectory of the system $\bG$. Note also that $\br'(t)=\bG$ since
\begin{equation}
\br'(t) = \langle x'(t), y'(t) \rangle = \langle P, Q \rangle = \bG.
\end{equation}
In Task 1, we will reformulate solutions to the predator-prey model in terms of level sets of some underlying potential function. This will provide a vehicle for estimating parameters of the model using data, in Task 2.
First let's gain some intuition about the problem. Imagine a vector field $\bG$ of a typical predator-prey system, and form a new vector field $\bF$ which is orthogonal to $\bG$ at any given $(x,y)$. If there exists some potential function $f$ such that $\bF = \nabla f$, then it appears that the level curves of $f$ would correspond to solution curves of the system $\bG$. (See Figure \ref{fig:figure1}.)
Note: There is no need to use \M in this Task, however, you are welcome to employ it as you see fit to check any of the simple algebraic manipulations.
\begin{enumerate}
\item Let's explore some properties of $\bF$. Define
\begin{equation}
\bF = h(x,y) \ \langle Q, \ -P\rangle
\end{equation}
where $h(x,y)$ is some arbitrary nonzero function.
\begin{enumerate}
\item Show that $\bG=\langle P, \ Q\rangle$ is not conservative.
\item Show that, with $h(x,y)=1$, $\bF$ is also not conservative.
\item Show that $\bF$ is \emph{orthogonal} to $\bG$ for any $h(x,y)$.
\item Finally, show that in order for $\bF$ to be \emph{conservative}, $h(x,y)$ must solve the following \emph{partial differential equation}:
\begin{equation}
\pd{}{x}\big(hP\big) + \pd{}{y}\big(hQ\big) = 0
\label{eq:T1concondition}
\end{equation}
\end{enumerate}
\begin{figure}[hbt]
\centering
\begin{subfigure}[t]{0.4\linewidth}
\includegraphics[width=\linewidth]{fig1a.pdf}
\caption{Example Predator-prey system ($\bG$) and solution curve ($\br(t)$).}
\end{subfigure} \quad
\begin{subfigure}[t]{0.4\linewidth}
\includegraphics[width=\linewidth]{fig1b.pdf}
\caption{Example gradient field ($\bF=\nabla f$) and level curve $f(x,y)=k$.}
\end{subfigure}
\caption{Intuition: the predator-prey system vector field (left) is orthogonal to the gradient vector field of some potential function (right). Then a solution curve of the predator-prey system corresponds to a level curve of the potential function.}
\label{fig:figure1}
\end{figure}
% \newpage
\item We will now give mathematical justification of our observation in Figure \ref{fig:figure1}. We know $\br'(t)=\bG$. Let $\bq(t) = h(x,y) \ \langle y(t), -x(t)\rangle$, and we also have $\bq'(t) = \bF$.
\begin{enumerate}
\item Show that $\br'(t)$ is orthogonal to $\bq'(t)$, for any $t$ and any $h(x,y)$.
\item Assume that $f(x,y)$ exists such that $\bF = \nabla f$. Show that
\begin{equation}
\frac{d}{dt}f\big(x(t), \ y(t)\big) = 0
\label{eq:T1levelcurvedef}
\end{equation}
(Hint: Begin by expanding the left side of this equation using the multivariable chain rule. Then observe this results in a dot product of two quantities for which you have already shown a key property.)
\end{enumerate}
In other words, \eqref{eq:T1levelcurvedef} tells us the change in $f$ is zero along the solution trajectory given by $\br(t)$ --- the definition of a level curve!
\item We could attempt to solve the partial differential equation in \eqref{eq:T1concondition} directly, but instead, let's try a different attack. Following the derivation in \cite{Zill} (Zill, Chapter 10), we form a new differential equation in terms of $dy/dx$ and perform a separation of variables. We get the expression:
\begin{equation}
d \ln x + a \ln y - cx - by = k
\label{eq:T1levelcurve}
\end{equation}
where $a,b,c$, and $d$ are parameters of the original model. Here, $k$ is a constant that is determined by initial conditions $(x_0, y_0)$ --- i.e. $k = d \ln x_0 + a\ln y_0 - cx_0 - by_0$.
Define
\begin{equation}
f(x,y) = d \ln x + a \ln y - cx - by
\end{equation}
and note that $f(x,y)=k$ corresponds to a particular level curve of $f(x,y)$.
\begin{enumerate}
\item Show that $\bF = \nabla f$. (This will involve identifying $h(x,y)$.)
\item Verify that $h(x,y)$ satisfies \eqref{eq:T1concondition}.
\end{enumerate}
\item Summarize your findings in this Task. What relationship have we shown, and how did we show it? Looking toward the next Task, estimating the parameters $a,b,c$ and $d$, why might this be useful?
\end{enumerate}
%%%%%%%%%%%%%
% TASK 2
\section*{TASK 2: Estimating parameters based on data}
Researchers working with IRNP have collected data on the moose and wolf populations over the last $N=61$ years, in particular 1956-2017.\footnote{We will use a slight modification to this data which is given in the \M notebook issued with this project.} Each datapoint is a pair $(x_i, y_i)$, with the wolf ($x$) and moose ($y$) population for the $i$-th year. The observations are averaged over several park rangers' counts at the beginning of each year, and the moose are counted by 10's. For example, $(x_1, \ y_1) = (21, 53.4)$ would represent there were 21 wolves and 534 moose observed in 1956.
In Task 2 we will estimate parameters of a predator-prey model that models this data.
For convenience define $w=\ln x$ and $z=\ln y$, and we can now rewrite Equation \eqref{eq:T1levelcurve} as
\begin{equation}
dw-cx-by+az=k
\label{eq:T2hyperplane}
\end{equation}
Notice this represents an equation of a plane in 4-dimensional $(w,x,y,z)$-space (i.e. a ``hyperplane'').\footnote{You may also interpret it as a 4-dimensional level surface where the 5-dimensional function $f(w,x,y,z)=dw-cx-by+az$ equals $k$. Say that five times fast!}
Our data represents 61 points in this 4-dimensional space, with coordinates $(w_i, x_i, y_i, z_i)$ for the $i$-th datapoint. For example, a year-end reading of 21 wolves and 534 moose represents the datapoint
\begin{equation}
(w,\ x, \ y, \ z) = \big(\ln(21), \ 21, \ 53.4, \ \ln(53.4)\big) = (3.044, \ 21, \ 53.4, \ 3.98)
\end{equation}
(Remember that moose are in 10's.)
Unfortunately there is no set of parameters $a,b,c,d$, and $k$ such that Equation \eqref{eq:T2hyperplane} is fulfilled for every single datapoint! (This is similar to the idea of linear regression with two-dimensional data: the datapoints appear to form a line, but not a perfect line. Now, the datapoints appear to form a plane, but not a perfect plane.)
\footnote{Note: the typical regression method is \emph{ordinary least squares}. In this Task, we are performing \emph{total least squares}, or orthogonal regression, also known as error-in-variables regression. (See \cite{Fuller}.)}
\begin{figure}[hbt]
\centering
\begin{subfigure}[t]{0.25\linewidth}
\includegraphics[width=\linewidth]{2d_tls_2.pdf}
\caption{Fitting a line (2-D).}
\end{subfigure} \qquad
\begin{subfigure}[t]{0.25\linewidth}
\includegraphics[width=\linewidth]{3d_tls_2.pdf}
\caption{Fitting a plane (3-D).}
\end{subfigure}
\caption{Choosing parameters for the line (2-D) or plane (3-D) that minimizes the orthogonal distance between datapoints and the plane. In Task 2 we are dealing with a hyperplane (4D).}
\label{fig:figure2}
\end{figure}
So, we would like to choose parameters $(a,b,c,d)$ and $k$ for our hyperplane that make it fit the data as ``best as possible''. See Figure \ref{fig:figure2}. Recall that the distance of a point $(w_0, x_0, y_0, z_0)$ to the hyperplane in \eqref{eq:T2hyperplane} is\footnote{Reference \cite{Stewart} Stewart, Section 12.5, p. 829-830.}
\begin{equation}
\text{dist} = \frac{|dw_0-cx_0-by_0+az_0-k|}{\sqrt{a^2+b^2+c^2+d^2}} \quad \Rightarrow \quad
\text{dist}^2 = \frac{(dw_0-cx_0-by_0+az_0-k)^2}{a^2+b^2+c^2+d^2}
\end{equation}
Let's select parameters $a,b,c,d$, and $k$ that \textbf{minimize the sum of the squared distances from each point to the hyperplane}. That is, we'd like to determine $a,b,c,d$ and $k$ which solve the optimization problem
\begin{equation}
\text{minimize} \ \sum_{i=1}^N \frac{(dw_i-cx_i-by_i+az_i-k)^2}{a^2+b^2+c^2+d^2}.
\label{eq:T2illdefined}
\end{equation}
Notice this optimization problem is ill-defined: we could just send $a$, $b$, $c$, and $d$ toward infinity, the denominator of \eqref{eq:T2illdefined} would get arbitrarily big, and the objective function would go to zero.\footnote{This is clearer if we peek ahead at \eqref{eq:T2optunscaled} and \eqref{eq:T2opt}, and note that if $s^2$ is allowed to go to infinity, or equivalently if we drop the constraint in \eqref{eq:T2opt}, the objective value will be (trivially) zero.}
We know the parameters must be finite, so the sum of their squares must also be finite --- let's define this sum to be a constant $s^2$. We can use this fact to \emph{constrain} our optimization, giving:
\begin{align}
\begin{split}
\text{minimize} \quad &\sum_{i=1}^N \frac{(dw_i+az_i-cx_i-by_i-k)^2}{s^2} \\
\text{subject to} \quad & a^2+b^2+c^2+d^2 = s^2
\label{eq:T2optunscaled}
\end{split}
\end{align}
which we may solve using the technique of Lagrange multipliers.
But this requires knowing $s$, which we don't know! We can avoid this issue by dividing through by $s^2$ and creating the equivalent optimization problem
\begin{align}
\begin{split}
\text{minimize} \quad &\sum_{i=1}^N (\hat{d} w_i + \hat{a} z_i - \hat{c} x_i - \hat{b} y_i - \hat{k})^2 \\
\text{subject to} \quad & \hat{a}^2+\hat{b}^2+\hat{c}^2+\hat{d}^2 = 1
\label{eq:T2opt}
\end{split}
\end{align}
where
\begin{equation}
\hat{a}=\frac{a}{s}, \quad \hat{b}=\frac{b}{s}, \quad \hat{c}=\frac{c}{s}, \quad \hat{d}=\frac{d}{s}, \quad \text{and} \quad \hat{k}=\frac{k}{s}.
\end{equation}
\begin{enumerate}
\item Solve the constrained optimization problem in \eqref{eq:T2opt} for our data, using \M. What are the optimal scaled parameters $\hat{a}, \hat{b}, \hat{c}, \hat{d}$, and $\hat{k}$?
\item We really want the unscaled parameters $a,b,c,d$, and $k$. To do this, we use three facts:
\begin{itemize}
\item Recall from your study of systems of ODEs that we can closely estimate the period $T$ of the system using the eigenvalues of the Jacobian near the central critical point. (Refer to \cite{Zill} Zill Chapter 10.) Specifically, $T=\frac{2\pi}{\sqrt{ad}}$.
\item We know, through observation of the real populations, that the period $T$ is equal to about \period years.
\item Finally, recall that $\hat{a}\hat{d} = ad/s^2$.
\end{itemize}
Use these facts to find an equation for $s$ in terms of the period and the product $\hat{a}\hat{d}$. Then use this value of $s$ to recover $a,b,c,d$, and $k$.
% Use these facts to recover $(a,b,c,d)$.
\item Using \verb|NDSolve|, with your estimated model parameters, plot the trajectories of $x(t)$ and $y(t)$. Overlay the actual data with these plots. Make a qualitative assessment of your model against the real data.
\end{enumerate}
%%%%%%%%%%%%%
% TASK 3
\section*{TASK 3: Intervention Strategy}
Michigan Technological University's lead wolf researcher, Rolf Peterson, has recommended resorting to human intervention strategies to restore balance in the wolf-moose balance on Isle Royale. Specifically, he proposes relocating two wolfpacks from the nearby forests in Canada to the Isle.
You propose to assist this decision process by modeling the density of moose on the Isle. Throughout this task, we will place a coordinate axes on the Isle with origin in the center, the main body of the island running along the $x$-axis, and the scale in miles. (See Figure \ref{fig:isleaxes}.)
The Isle Royale National Park knows the moose tend to herd together in two distinct groups, North and South, based on their typical feeding patterns. Under our coordinate axes scheme, the observed geographic centers of the Moose Group North and Moose Group South are
\begin{equation}
\text{North:} \ \ (10,1), \qquad \text{South:} \ \ (-16,-1)
\label{eq:T3geocenters}
\end{equation}
There are a total of \textbf{\totalmoose moose} on the Isle at a \textbf{ratio of 3:2} between the North and South populations.
\begin{figure}[bt]
\centering
\includegraphics[width=0.6\linewidth]{fig_isle_axes3.pdf}
\caption{Isle Royale with coordinate axes. (\emph{Image: Google Maps.})}
\label{fig:isleaxes}
\end{figure}
\begin{enumerate}
\item First consider a density function of the form:\footnote{Although we are not using it to model probability, this density function is in the family of \emph{bivariate Gaussian probability density functions}, i.e. the normal distribution. The univariate version of this density function is sometimes called the ``bell curve.''}
\begin{equation}
P(x,y) = \frac{1}{25\pi} \text{exp}\Big\{-\frac{1}{25}\Big((x-m)^2 + (y-n)^2\Big) \Big\}
\label{eq:T3density}
\end{equation}
with parameters $m$ and $n$. Notation: $\text{exp}\{x\}=e^x$. Let's familiarize ourself with the behavior of $P(x,y)$ on the domain of all of $\R^2$.
\begin{enumerate}
\item Using \M, to what value does $P(x,y)$ integrate, for arbitrary parameters $m, n$? (Bonus: show how to solve this by hand, by switching to polar coordinates.)
\item\label{q:T3com} Using \M, what is the center of mass $(\bar{x}, \bar{y})$ of $P(x,y)$, for arbitrary parameters $m, n$? Why is this property of $P(x,y)$ useful?
\end{enumerate}
\end{enumerate}
Now we'll use a combination of density functions of the form $P(x,y)$ to model the moose population on the Isle. Let's model each moose subpopulation (North and South) with its own density function and associated parameters, that is
\begin{equation}
P_{\text{north}} (x,y) \quad \text{with} \quad m_{\text{north}}, \ n_{\text{north}}
\qquad \qquad
P_{\text{south}} (x,y) \quad \text{with} \quad m_{\text{south}}, \ n_{\text{south}}.
\label{eq:T3dfuns}
\end{equation}
Choose these 4 parameters based on the observed geographic subpopulation centers and your observation in Question \ref{q:T3com}. Let the units of $P_{\text{north}}$ and $P_{\text{south}}$ be moose per sq. miles (the same scale as the coordinate axes in Figure \ref{fig:isleaxes}).
Then create a new, combined density function $\rho(x,y)$ for the whole Isle:
\begin{equation}
\rho(x,y) = K \Big( 0.6 \ P_{\text{north}} (x,y) + 0.4 \ P_{\text{south}} (x,y)\Big).
\end{equation}
where $K$ is a scaling parameter. Note the units of $\rho(x,y)$ are also moose per sq. miles.
\begin{enumerate}
\setcounter{enumi}{1}
\item What is the purpose of the ``weights,'' 0.6 and 0.4, in $\rho(x,y)$, in relation to the distribution of the moose groups on the island?
\item Finally, we need to incorporate the shape of the island into our model.
\begin{enumerate}
\item Based on Figure \ref{fig:isleaxes}, let's use an ellipse to model the boundary of the island. Recall that an ellipse has the form $\frac{x^2}{a^2} + \frac{y^2}{b^2} = c^2$. Choose parameters $a,b,c$ for an ellipse that closely model the island and visualize your choice.
\item Integrate $\rho(x,y)$ over your chosen model of the island boundary, and choose an appropriate value for $K$ given the total moose population. You will likely need to use \verb|NIntegrate| to evaluate the double integral.
% For example,
% \begin{equation}
% \int_a^b \int_{g_1(x)}^{g_2(x)} f(x,y) \ dy dx \ \ \Rightarrow \ \
% \verb|NIntegrate[f[x,y], {x,a,b}, {y,g1[x],g2[x]}]|
% \end{equation}
% Notice the inside integral's limits should be listed second in this command.
\item Visualize $\rho(x,y)$ over your chosen model of the domain using \verb|ContourPlot|. The commands \verb|RegionFunction| and \verb|AspectRatio| will be useful. For example,
\begin{verbatim}
ContourPlot[x Cos[y], {x,-5,5}, {y,-2,2}, AspectRatio->2/5,
RegionFunction->Function[{x,y}, (x/4)^2+(y/2)^2 <= 1]]
\end{verbatim}
plots the function $x \cos(y)$ over the interior of an ellipse while maintaining the difference in scale in the $x$ and $y$ directions. (Bonus: try using \verb|Plot3D| with \verb|BoxRatios|.)
\end{enumerate}
\item Calculate the overall center of mass of the total moose population of the entire island, using $\rho(x,y)$. Include this location, along with the North and South Group centers of mass, on your visualization of the island-wide moose density above.
\item Make a recommendation to Rolf Peterson and the Isle Royale National Park as to where to place the two Canadian wolfpacks. Does your recommendation change if we are only able to acquire one wolfpack? Use your recommendation(s) as a new initial condition for the predator-prey system, and find the resulting solution curves. Note this is now a projection of future population trajectories. Use the results to quantify the effects of your recommendation.
\end{enumerate}
\begin{thebibliography}{98}
\bibitem{Fuller} Fuller, Wayne A. 1987. \textit{Measurement Error Models.} New York: Wiley.
\bibitem{Milot} Milot, Christine. 2017. Two wolves survive in world's longest running predator-prey study. \emph{Science}. \href{http://www.sciencemag.org/news/2017/04/two-wolves-survive-world-s-longest-running-predator-prey-study}{http://www.sciencemag.org/news/2017/04/two-wolves-survive-world-s-longest-running-predator-prey-study}. Accessed 19 May 2018.
\bibitem{IRNP} The Population Biology of Isle Royale Wolves and Moose. 2017. \href{http://www.isleroyalewolf.org/data/data/home.html}{http://www.isleroyalewolf.org/data/data/home.html}. Accessed 19 May 2018.
\bibitem{Stewart} Stewart, James. 2016. \textit{Calculus, 8th Edition}. Boston: Cengage Learning.
\bibitem{Zill} Zill, Dennis G. and Warren S. Wright. 2013. \textit{Differential Equations with Boundary Value Problems, 8th Edition}. Boston: Cengage Learning.
\end{thebibliography}
\teacher{
\section*{NOTES FOR TEACHERS}
\paragraph{Prerequisites.} Students must have some background in differential equations to be able to engage with the predator-prey application. The project uses some introductory vector calculus, but this is confined to Task 1 and could be omitted for a non-vector calculus course (see below on Task 1: Modification). The project makes heavy use of technology, in particular \M, although it is possible to modify this to a different platform (MATLAB, Python).
\paragraph{Timeline.} The project is designed to have three assessment points, in order: an In-Progress Review (IPR), a formal presentation, and a written report. The IPR is informal and designed to prevent any major derailment or procrastination; the presentation is formal, but not heavily weighted in the final grade; and the writeup is the clean, final product. As originally administered this took place over about 2.5 weeks of a 4.5 credit hour, semester-long (15 week) course. We dedicated 6 of 56 total lesson periods to the project. We operated on the expectation that each lesson corresponds to 3 hours of student time, giving a total of 18 hours for the project, per member of the 2-3 person groups.
% (A faculty member took the project ``for time'' beforehand, and it took him four hours to complete the analysis and another 1.5 for writeup. This satisfied our rule-of-thumb for assessments of a 1-to-3 ratio.)
\paragraph{Preparatory activities.} We incorporated some of the more challenging concepts from the project in smaller assignments earlier in the course. This helped set foundations, both conceptually and technically (i.e. coding). Specifically: we introduced the connections of Task 1 in an earlier lab assignment with the simple system $x' = -y$, $y'=x$, in an observational way without any rigorous mathematics. We also introduced parameter estimation in an earlier lab using ordinary least squares and the Cobb-Douglas production model (based on Stewart Chapter 14.1, Exercise 81).
\paragraph{Task 1.}
\begin{enumerate}
\item Ideally, students answer Question 4 (``summarize your findings'') by recognizing they are establishing a rigorous connection between the system of ODEs and a multivariable potential function, thus providing us a vehicle for estimating parameters in Task 2. Realistically, many will ``miss the forest for the trees'' and it is important to emphasize this connection.
\item \textbf{Modification.} In truth, the derivation of the implicit function $f(x,y)=k$ (i.e. \eqref{eq:T1levelcurve}) from the phase-plane method is enough to proceed to Task 2, and to shorten the project and/or eliminate the need for vector calculus, Task 1 could be completely replaced by a reference to this derivation. (However, this does not establish the stronger connection between the level sets of this function and solution curves of the ODE.)
\end{enumerate}
\paragraph{Task 2.}
\begin{enumerate}
\item Students may miss the fact that we are optimizing in the parameter space, and the datapoints $w,x,y,z$ are fixed. (This is the reverse of typical textbook multivariable optimization problems where $w,x,y,z$ are the variables and $a,b,c,d$ are fixed coefficients.)
\item Students may accidentally omit the constraint, leading to parameter estimates of zero (or numerically close to zero). This, of course, is the whole purpose of the constraint, and provides an opportunity to revisit the construction of the problem.
\item The parameter estimates $\hat{a}, \hat{b}, \hat{c}, \hat{d}, \hat{k}$ provide an acceptable approximation of the data \emph{in the phase plane}, that is, without regard to time. The re-scaling just stretches/shrinks the solution curves in \emph{time} and makes the approximation match in a $x$ or $y$ vs $t$ plot. This makes sense since our objective function was based in the phase-plane, and so we need to use external knowledge related to periodicity (i.e. the period $T$) to scale the model correctly. Get students to notice this effect by asking, what happens when we change the estimated period $T$?
\end{enumerate}
\paragraph{Task 3.}
\begin{enumerate}
\item The property of a bivariate Gaussian density function (and of course, many other density functions) that one can ``program in'' the center of mass may be an unfamiliar concept to students and is worth highlighting.
\item For brevity, Task 3 doesn't explore the actual modeling result of constructing the island density function, but this is a possible direction: for example, the estimated density of moose in a particular area of interest on the island.
\end{enumerate}
\paragraph{Other information.}
\begin{enumerate}
\item The data and solution provided with this modeling scenario is based closely on the \emph{real data} on Isle Royale, available publicly from \cite{IRNP}. However, for the project, we used a modification of this data with fewer outliers and cleaner population trajectories. Our rationale here was that use of the actual data would shift too much focus to statistical analysis, and increase the difficulty unnecessarily by introducing the extra step of downloading, cleaning, and loading the data into \M.
\item This project was a capstone assignment for the two-semester undergraduate Advanced Mathematics Program at the U.S. Military Academy (West Point), designed for first year students who display strong abilities in mathematics and validate the core mathematics program, and consisting of a semester in modeling with differential equations and a semester in multivariable and vector calculus.
\end{enumerate}
}
\end{document}