% Oscillating Chemical Reactions Lab, Math 230, St Olaf College
% Matthew Wright
% Spring 2016
% based on the lab by Robert L. Devaney, http://math.bu.edu/people/bob/MA226/chem-osc-lab.html
\documentclass[10pt]{article}
\usepackage{amsmath,amssymb,amsthm}
\usepackage[usenames,dvipsnames]{color}
\usepackage{graphics}
\usepackage{framed}
\usepackage{fullpage}
\usepackage{enumitem}
\usepackage{tikz}
\usepackage{todonotes}
\usepackage{listings}
\lstset{basicstyle=\ttfamily\small,breaklines=true}
\parindent 0.0in
\newenvironment{solution}{\footnotesize \color{blue}}{\normalsize \color{black}}
\newcommand{\dps}{\displaystyle}
%%%%%%%%%% THE ASSIGNMENT %%%%%%%%%%
\title{Oscillating Chemical Reactions Lab}
\makeatletter
\let\thetitle\@title
\makeatother
\begin{document}
\pagestyle{empty}
\begin{center}
{\Large \textbf{\thetitle}} \\
Math 230 \\
due Tuesday, May 17 at 4pm
\end{center}
Certain chemical reactions oscillate rather than approach equilibrium. In this lab you will investigate a system of differential equations that models such a reaction. The system depends on a parameter, $B$, and you are to describe qualitatively what happens when $B$ passes through a certain threshold value. There are two parts to this lab. After making the appropriate computations in Part I, use technology to view the phase plane for the system and report on the qualitative behavior of the system in Part II.
\vspace{11pt}
Suppose that iodide and chlorine dioxide are each poured into a reaction vessel at a certain rate and are mixed thoroughly, causing a reaction to occur. The excess solution is removed so that the amount of reactant in the vessel remains constant. Let $x$ be the amount of iodide and $y$ the amount of chlorine dioxide in the vessel at time $t$. A mathematical model for a chlorine dioxide-iodide open system reaction is provided by the following system of (nonlinear) differential equations:
\begin{align*}
\frac{dx}{dt} &= -x + 10 - \frac{4xy}{1 + x^2} \\
\frac{dy}{dt} &= B\left( x - \frac{xy}{1 + x^2} \right)
\end{align*}
The quantity $B$ is a parameter that we will let vary between $1$ and $6$.
\vspace{11pt}
{\bf Instructions:} Complete each of the numbered items below.
Write your solutions clearly and neatly (typing is preferred).
Make sure that you clearly state the solution to each numbered item!
If you turn in a \emph{Mathematica} notebook, include section headers and explanations so that your notebook is easy to read. Pretend that you are writing for someone who does not understand \emph{Mathematica} code!
\subsection*{Part I: Calculations}
Do the following calculations first on scratch paper, and then write your solutions clearly and neatly.
\begin{enumerate}[label=\bf{\arabic*.}]
\item Find all of the equilibrium points for this system for all values of the parameter $B$.
\emph{Hint:} You should find exactly one equilibrium point for each value of $B$, and this equilibrium point should not depend on $B$. If this is not the case, your ship is sinking and you had better call a lifeguard!
\item Compute the Jacobian matrix of the system at a general point $(x, y)$. \emph{Be careful here:} If your matrix is wrong, your ship has sunk and you can't go on.
\item Now compute the Jacobian matrix at the equilibrium point you found in 1. This matrix should depend on $B$. Then find the eigenvalues of this matrix. These should also depend on $B$.
\item Classify these equilibria as spiral sinks, real sources, saddles, or whatever. Your answer should include a range of $B$ values for each type that occurs. Remember, we are assuming that $1 < B < 6$.
\item Next, answer the following question: At what value of $B$ does a bifurcation occur, and what happens at this bifurcation?
\item Sketch the curve in the trace-determinant plane given by this matrix as $B$ varies.
\item Sketch the nullclines for this system and indicate the direction of the vector field in the regions between these nullclines. How do these nullclines depend on the parameter $B$? It suffices to restrict attention here to the region $0 < x < 5$, $3 < y < 8$.
\end{enumerate}
\subsection*{Part II: Observations}
Use technology to help you answer the following questions. For example, you could use \emph{Mathematica} to numerically solve the system with various initial conditions and display solution curves in the phase plane. See sample \emph{Mathematica} code in the reference section below. Explain your answers clearly.
\begin{enumerate}[label=\bf{\arabic*.},resume]
\item Discuss the behavior of solutions of the system for the parameter value $B = 6$. Answer the question: What happens to the chemicals in the vessel as time passes?
\item Discuss the behavior of solutions of the system for the parameter value $B = 1$. Answer the question: What happens to the chemicals in the vessel as time passes?
\item As the parameter $B$ varies between $1$ and $6$, we know that there is a bifurcation (a qualitative change in the solutions) at some $B$ value, and from the linearization in Part I, we know what to expect there.
Write a few sentences explaining what this bifurcation entails, both mathematically and chemically. That is, as the parameter $B$ passes through the bifurcation value, explain what happens in the phase plane and what this means in terms of the chemicals involved in the reaction.
\end{enumerate}
\vspace{11pt}
\subsection*{Reference: Phase Plots in \em{Mathematica}}
The following code will help you plot the phase plane, with vector streams and a sample solution on the same diagram.
\begin{itemize}
\item Create a stream plot and save it as {\tt splot} (assume {\tt B} has a value between $1$ and $6$):
\begin{lstlisting}
splot = StreamPlot[{-x + 10 - 4x*y/(1 + x^2), B(x - x*y/(1 + x^2))},
{x, 0, 5}, {y, 3, 8}]
\end{lstlisting}
\item Numerically solve the system with initial condition $x(0) = \mathtt{r}$ and $y(0) = \mathtt{s}$; save the solution as {\tt sol}:
\begin{lstlisting}
eqns = {x'[t] == -x[t] + 10 - 4 x[t] y[t]/(1 + (x[t])^2),
y'[t] == B (x[t] - x[t] y[t]/(1 + (x[t])^2)),
x[0] == r, y[0] == s};
sol = NDSolve[eqns, {x, y}, {t, 0, 20}]
\end{lstlisting}
\item Create a parametric plot of the solution; save it as {\tt pplot}:
\begin{lstlisting}
pplot = ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 20},
PlotStyle -> Red]
\end{lstlisting}
\item Plot the stream plot and parametric plot together:
\begin{lstlisting}
Show[splot, pplot]
\end{lstlisting}
\end{itemize}
\vfill
{\scriptsize based on a lab by Robert L.\ Devaney: \texttt{http://math.bu.edu/people/bob/MA226/chem-osc-lab.html}}
\end{document}