\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts,amssymb,graphicx,amsthm,html,draftcopy}

\newtheorem{theorem}{Theorem}
\theoremstyle{definition}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\newtheorem{problem}{Problem}
\newtheorem{example}{Example}
\newtheorem{exercise}{Exercise}

\def\Z{\mathbb{Z}}
\def\N{\mathbb{N}}

\renewcommand{\theenumi}{\roman{enumi}}
\renewcommand{\labelenumi}{(\theenumi)}

\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\i}{\mathrm{i}}

\begin{document}

\title{General Mathematics and Computational Science I}
\author{Difference Equations}
\date{Marcel Oliver}
\maketitle

\begin{abstract}
These notes provide a first introduction to the qualitative theory of
difference equations.  We introduce the concept of equilibrium points,
periodic points, and stability to derive a general description of the
dynamics of a solution to a difference equation even in situations
where a complete closed-form solution does not exist.  These concepts
are illustrated by several examples from modeling real-world systems.
To a large part, this material is distilled from the book by Elaydi
\cite{E96}, which contains a wealth of more advanced topics and
further examples.
\end{abstract}

%\tableofcontents

\section{Introduction}

In general, a \emph{difference equation} of \emph{order} $k$ is a
relation of the form
\begin{equation}
  \label{e.difference}
  x_{n+1} = f(x_n, \dots, x_{n-k+1})
\end{equation}
together with $k$ starting values $x_0, \dots, x_{k-1}$.  A sequence
of real numbers $x_n$ with $n = 0, 1, \dots$ is called the
\emph{solution} to the difference equation if it satisfies
\eqref{e.difference} as well as the given starting values.  The set
$\{x_n \colon n=0,1,\dots \}$ is called an \emph{orbit}.

A difference equation is called \emph{linear} if it is of the form
\begin{gather}
  x_{n+1} = a_0 (n) \, x_n + a_1 (n) \, x_{n-1} + \dots + 
            a_{k-1}(n) \, x_{n-k+1} + b (n) \,.
  \label{e.inhomogeneous}
\end{gather}
If the $a_i$ do not depend on $n$, we speak of a \emph{constant
coefficient} difference equation; if $b(n)=0$ we say that the equation
is \emph{homogeneous}. 

We have already encountered difference equations before.  For example,
the Fibonacci sequence
\begin{subequations}
\begin{gather}
  a_{n+1} = a_n + a_{n-1} \,, \\
  a_0 = a_1 = 1 \,,
\end{gather}
\end{subequations}
is a linear, second-order, constant coefficient, homogeneous
difference equation.  We already know that such equations can be
solved by using the method of generating functions.

The main point here is to develop techniques that remain applicable
when closed-form solutions are not available.  However, we will
frequently use solvable equations as examples.  A relatively simple
class of such equations is first order linear difference equations,
which are discussed next.


\section{First order linear difference equations}
\label{s.first-order}

\begin{problem} \label{p1}
To treat a chronic disease, a dose of $d\,\mathrm{mg}$ of a drug is
administered starting from day.  Each day, a fraction of $p$ of
the drug is eliminated from the body.  What is the amount of drug in
the body after a very long time?
\end{problem}

Let $a_n$ denote the amount of drug in the body at day $n$.  Suppose
the treatment starts at day $n=1$.  Then $a_0=0$ and
\begin{equation}
  a_{n+1} = (1-p) \, a_n + d \,.
\end{equation}
This is a first order linear difference equation which can be solved
as follows.  On day $n$, the dose which was administered on day $k$ has
been in the body for $n-k$ days.  Of this dose, the amount remaining
after each day is $1-p$ times the amount that was present the previous
day.  So the amount left of the dose at day $n$ is $(1-p)^{n-k}\,d$.
Thus,
\begin{align}
  a_n & = \sum_{k=1}^n (1-p)^{n-k} \, d
          \notag \\
      & = d \sum_{j=0}^{n-1} (1-p)^j 
          \notag \\
      & = d \, \frac{1-(1-p)^{n}}p \,,
\end{align}
where, in the last step, we used the formula for the geometric series
which we encountered before.  In particular, as $n$ becomes large,
$a_n$ converges to the \emph{equilibrium} amount of drug in the body
\begin{equation}
  \lim_{n\to\infty} a_n = \frac{d}p \,.
\end{equation}

\bigskip

The same reasoning can be applied to general first order linear
difference equations with nonconstant coefficients,
\begin{subequations}
  \label{e.gen-first-order}
\begin{gather}
  x_{n+1} = a_n \, x_n + b_n \,, \\
  x_0 = c \,.
\end{gather}
\end{subequations}
By tracing out the same steps as before, we obtain the following
general form of the solution.

\begin{theorem}
The solution to the general first order difference equation
\eqref{e.gen-first-order} is given by
\begin{equation}
  x_n = c \prod_{j=0}^{n-1} a_j
      + \sum_{k=0}^{n-1} b_k \prod_{j=k+1}^{n-1} a_j \,.
  \label{e.general}
\end{equation}
\end{theorem}

\begin{remark}
Note that the sum and product in \eqref{e.general} are nested, i.e.\
the expression reads
\begin{equation}
  \sum_{k=0}^{n-1} b_k \prod_{j=k+1}^{n-1} a_j
  = b_0 \, a_1 \cdots a_{n-1} + \dots 
  + b_{n-2} \, a_{n-1} + b_{n-1} \,. 
\end{equation}
\end{remark}

\begin{proof}
We prove this formula by induction.  For $n=0$, the formula clearly
reproduces the initial value $x_0=c$.  Let us now assume that the
formula holds for some $n\geq 0$.  Then
\begin{align}
  x_{n+1} & = a_n \, x_n + b_n 
              \notag \\
          & = a_n \, \biggl(
                       c \prod_{j=0}^{n-1} a_j + 
                       \sum_{k=0}^{n-1} b_k \prod_{j=k+1}^{n-1} a_j 
                     \biggr) + b_n
              \notag \\
          & = c \prod_{j=0}^{n} a_j
              + \sum_{k=0}^{n-1} b_k \prod_{j=k+1}^{n} a_j
              + b_n
              \notag \\
          & = c \prod_{j=0}^{n} a_j 
              + \sum_{k=0}^{n} b_k \prod_{j=k+1}^{n} a_j \,.              
\end{align}
This completes the inductive step.
\end{proof}

\begin{example} \label{ex1}
The linear, constant coefficient equation
\begin{equation}
  x_{n+1} = a \, x_n + b \,,
\end{equation}
an instance of which we encountered already in Problem~\ref{p1}, has
solution
\begin{align}
  x_n & = x_0 \prod_{j=0}^{n-1} a 
          + \sum_{k=0}^{n-1} b \prod_{j=k+1}^{n-1} a 
          \notag \\
      & = a^n \, x_0 + b \sum_{k=0}^{n-1} a^{n-1-k}
          \notag \\
      & = a^n \, x_0 + b \sum_{j=0}^{n-1} a^{j}
          \notag \\
      & = \begin{cases}
            x_0 + n \, b & \text{if } a=1 \\
            a^n \, x_0 + b \, \dfrac{a^n-1}{a-1} 
            & \text{if } a \neq 1 \,,
          \end{cases}
          \label{e.ex1}
\end{align}
where the first step is a direct use of the general expression
\eqref{e.general} and the last step is again due to the formula for
the geometric series.
\end{example}

\begin{example}
The difference equation
\begin{gather}
  x_{n+1} = (n+1) \, x_n + 2^n \, (n+1)! \,, \\
  x_0 = 1
\end{gather}
has solution
\begin{align}
  x_n & = \prod_{j=0}^{n-1} (j+1) 
          + \sum_{k=0}^{n-1} 2^k \, (k+1)! 
            \prod_{j=k+1}^{n-1} (j+1) 
          \notag \\
      & = n! + \sum_{k=0}^{n-1} 2^k \, (k+1)! \,
               \frac{n!}{(k+1)!}
          \notag \\
      & = n! + n! \, \frac{2^n-1}{2-1} 
          \notag \\
      & = n! \, 2^n \,.
\end{align}
\end{example}


\section{Stability of equilibrium points}
\label{s.stability}

It is rare for a nonlinear, even first order, difference equation to
have a closed-form solution.  We thus seek means to qualitatively
describe properties of the solution that can be deduced from the
equation itself independent of whether a closed-form solution is
available.  Often, the index $n$ plays the role of a time-like
quantity in the modeling regime and the difference equation describes
how a dependent quantity changes in time.  In this context, difference
equations are often referred to as \emph{discrete dynamical systems}.
The adjective ``discrete'' refers the time-like quantity being an
integer index, not a real number.

This section introduces the important concept of equilibrium points
and discusses their stability.  We motivate the theory using an
(explicitly solvable) example from economics.

\begin{problem} \label{p2}
A commodity is freely traded.  The suppliers of this commodity
increase production whenever the price is high and decrease production
when the price is low.  Buyers, on the other hand, react adversely to
changes in price.  Moreover, due to the time it takes a supplier to
complete a unit of the product, changes in production always lag
behind changes in price.  Under the hypothesis that the market price
is the price at which supply equals demand, does the market adjust to
a stable equilibrium?
\end{problem}

Note that this problem is not fully specified yet, and that we also
have not yet defined the precise meaning of ``stable'' and
``equilibrium''.  To proceed further, we make the simplest modeling
assumptions that capture the essence of this problem.  

First, we assume discrete market periods, labeled with integer $n$.
Let $p_n$ denote the price of the commodity, $d_n$ the demand, and
$s_n$ the supply in period $n$.  Second, we assume that the
price-demand and price-supply relationships are linear, i.e.\ 
\begin{gather}
  d_n = - m_{\mathrm{d}} \, p_n + b_{\mathrm{d}}
  \label{e.p-d}
\intertext{and}
  s_{n+1} = m_{\mathrm{s}} \, p_n + b_{\mathrm{s}} \,,
  \label{e.p-s}
\end{gather}
where $m_{\mathrm{d}}$ and $m_{\mathrm{s}}$ are positive.  The delay
in the reaction of the supply side to changes in price is manifest in
\eqref{e.p-s} through the dependence of the supply on the price of the
\emph{previous} market period.  Finally, the market price hypothesis
is expressed as
\begin{gather}
  d_n = s_n \,.
  \label{e.d-s}
\end{gather}
By combining \eqref{e.p-d}, \eqref{e.p-s}, and \eqref{e.d-s}, we obtain
a first order linear difference equation for $p_n$, namely
\begin{gather}
  p_{n+1} = - \frac{m_{\mathrm{s}}}{m_{\mathrm{d}}} \, p_n 
            + \frac{b_{\mathrm{d}}-b_{\mathrm{s}}}{m_{\mathrm{d}}}
              \,. 
  \label{e.pn+1}
\end{gather}
An \emph{equilibrium price} $p^*$ is a price that does not change from
one period to the next.  Thus, at the equilibrium, we must have
\begin{gather}
  p^* = - \frac{m_{\mathrm{s}}}{m_{\mathrm{d}}} \, p^*
        + \frac{b_{\mathrm{d}}-b_{\mathrm{s}}}{m_{\mathrm{d}}} 
\intertext{and therefore}
  p^* = \frac{b_{\mathrm{d}}-b_{\mathrm{s}}}%
             {m_{\mathrm{d}}-m_{\mathrm{s}}} \,.
\end{gather}

\begin{figure}[tb]
\centering
  \includegraphics[width=10cm]{market-astable.pdf}
\caption{Asymptotically stable equilibrium price.} 
\label{f.c1}
\end{figure}

\begin{figure}[tb]
\centering
  \includegraphics[width=10cm]{market-stable.pdf}
\caption{Stable, but not asymptotically stable equilibrium price.} 
\label{f.c2}
\end{figure}

\begin{figure}[tb]
\centering
  \includegraphics[width=10cm]{market-unstable.pdf}
\caption{Unstable equilibrium price.} 
\label{f.c3}
\end{figure}

The behavior of this model can be visualized in a so-called
\emph{cobweb diagram}.  The difference equation \eqref{e.pn+1} is of
the form $p_{n+1} = f(p_n)$.  The procedure is as follows.
\begin{itemize}
\item Draw the graph of $f$ into an $x$-$y$ coordinate plane.
\item Choose an initial value $p_0$ on the $x$-axis.
\item Find the first iterate by moving vertically until intersecting
the graph of $f$.
\item Find the next iterate by moving horizontally until hitting the
line $x=y$, then vertically until intersecting the graph of $f$.
Repeat. 
\end{itemize}

Figures~\ref{f.c1}--\ref{f.c3} show the resulting cobweb diagrams for
three different values of the ratio
${m_{\mathrm{s}}}/{m_{\mathrm{d}}}$.  We can distinguish three
regimes. 

\begin{enumerate}
\item When ${m_{\mathrm{s}}}/{m_{\mathrm{d}}}<1$, the slope of the
graph of $f$ has magnitude less than one.  The distance from the
equilibrium price $p^*$ decreases in each market period, see
Figure~\ref{f.c1}.  Such a market is called \emph{stable}.

(Mathematically, we speak of \emph{asymptotic stability} whenever $p_n
\to p^*$ as $n \to \infty$; see Definition~\ref{d.stability} below.)

\item When ${m_{\mathrm{s}}}/{m_{\mathrm{d}}}=1$, the slope of the
graph of $f$ has magnitude one.  The distance from the equilibrium
point neither decreases nor increases; $p_n$ oscillates between two
values, see Figure~\ref{f.c2}.  This market is still called stable.

(Mathematically, we say that the sequence $p_n$ is a $2$-cycle.  The
equilibrium price $p^*$ is called stable, but not asymptotically
stable.)

\item When ${m_{\mathrm{s}}}/{m_{\mathrm{d}}}=1$, the slope of the
graph of $f$ has magnitude greater than one.  The distance from the
equilibrium price $p^*$ increases in each market period, see
Figure~\ref{f.c3}.  Such a market is called \emph{unstable}.

\end{enumerate}
This result can be summarized as follows.

\begin{theorem}[Cobweb Theorem of Economics]
A market is stable if and only if the suppliers are no more sensitive
to price than the consumers.
\end{theorem}

\begin{remark}
We could have arrived at this conclusion more directly by noting that
\eqref{e.pn+1} is a linear, constant coefficient difference equation.
Therefore, its solution is covered by Example~\ref{ex1} with
$a=-{m_{\mathrm{s}}}/{m_{\mathrm{d}}}$ and $b= (b_{\mathrm{d}} -
b_{\mathrm{s}})/{m_{\mathrm{d}}}$.  The explicit solution
\eqref{e.ex1} shows that $p_n$ remains bounded as $n \to \infty$ if
and only if $\lvert a \rvert \leq 1$, i.e.\ if ${m_{\mathrm{s}}} \leq
{m_{\mathrm{d}}}$.  However, the pedagogical value of going the
graphical route is that it suggests how stability may be defined and
verified in cases where no closed form solution is available.  This is
the subject of the remainder of this section.
\end{remark}

We now proceed to discussing stability for general first order
difference equations of the form
\begin{equation}
  x_{n+1} = f(x_n) \,.
  \label{e.general2}
\end{equation}

\begin{definition} \label{d.stability}
A number $x^*$ is called an \emph{equilibrium point} or \emph{fixed
point} of the difference equation \eqref{e.general2} if
\begin{equation}
  x^* = f(x^*) \,.
  \label{e.simple-iteration}
\end{equation}
\end{definition}

\begin{definition}
An equilibrium point $x^*$ of the difference equation
\eqref{e.general2} is called
\begin{enumerate}
\item \emph{stable} if for every $\varepsilon > 0$ there exists some
$\delta > 0$ such that $\lvert x_n - x^* \rvert \leq \varepsilon$ for
every $n \geq 1$ whenever $\lvert x_0 - x^* \rvert \leq \delta$.

\item \emph{attracting} if there exists some $\eta > 0$ such that
\begin{equation}
  \lim_{n\to\infty} x_n = x^*
\end{equation}
whenever $\lvert x^0 - x^* \rvert \leq \eta$.

\item \emph{asymptotically stable} if it is stable and attracting.
\end{enumerate}
\end{definition}

In other words, an equilibrium point is stable if a solution remains
nearby whenever it starts sufficiently close to the equilibrium point.
Asymptotic stability is a stronger property.  A solution that starts
sufficiently close to an asymptotically stable equilibrium point will
not only remain nearby, but will actually be ``attracted'' by it.

\begin{remark}
An equilibrium point can be attracting without being stable.  An
example is the difference equation
\begin{equation}
  x_{n+1} =
  \begin{cases}
    2 \, x_n & \text{if } \lvert x_n \rvert < 1 \,, \\
    0        & \text{if } \lvert x_n \rvert \geq 1 \,.
  \end{cases}
\end{equation}
Here $x^*=0$ is clearly an attracting equilibrium point, but the
condition for stability fails.  (Working out the details of this
argument is a good exercise for those of you interested in Analysis;
for our present purposes this is more a mathematical curiosity.)
\end{remark}

In many situations, asymptotic stability is easy to verify.  A first,
simple result is the following.

\begin{theorem} \label{t.stability} 
Let $x^*$ be an equilibrium point of the difference equation
\eqref{e.general2} with continuously differentiable right hand side
$f$.
\begin{enumerate}
\item If $\lvert f'(x^*) \rvert < 1$, then $x^*$ is an asymptotically
stable equilibrium point.
\item If $\lvert f'(x^*) \rvert > 1$, then $x^*$ is an unstable
equilibrium point.
\end{enumerate}
\end{theorem}

\begin{proof}
First, consider the case when $\lvert f'(x^*) \rvert < 1$.  Since $f'$
is continuous, there exist $\eta>0$ and $L<1$ such that $\lvert
f'(x) \rvert < L$ for all $x \in [x^*-\eta, x^*+\eta]$.  Let us
assume that $x_0,\dots,x_n$ lie in this $\eta$-neighborhood of $x^*$.
Then
\begin{align}
  \lvert x^* - x_{n+1} \rvert
  & = \lvert f(x^*) - f(x_n) \rvert \notag \\
  & = \lvert f'(\xi_n) \rvert \, \lvert x^* - x_n \rvert \notag \\
  & \leq L \, \lvert x^* - x_n \rvert \,,
  \label{e.one-step}
\end{align}
where the second step is an application of the mean value theorem
which asserts the existence of some $\xi_n \in [x^*,x_n]$ such that
equality holds.  Since $L<1$, estimate \eqref{e.one-step} shows that
$x_{n+1}$ is also contained in the same $\eta$-neighborhood of $x^*$
and we can iterate the argument to conclude that
\begin{equation}
  \lvert x^* - x_n \rvert \leq L^n \, \lvert x^* - x_0 \rvert 
  \to 0
\end{equation}
as $n \to \infty$.  This proves part (i).

If, on the other hand, $\lvert f'(x^*) \rvert > 1$, by continuity of
$f'$ there exist $\eta>0$ and $L>1$ such that $\lvert f'(x^*) \rvert >
L$ for all $x \in [x^*-\eta, x^*+\eta]$.  Take any $x_0$ in this
$\eta$-neighborhood of $x^*$.  By re-tracing the argument above, we
find that 
\begin{equation}
  \lvert x^* - x_n \rvert \geq L^n \, \lvert x^* - x_0 \rvert 
\end{equation}
so long as $x_n$ remains in the $\eta$-neighborhood of $x^*$.
However, since $L>1$, the right hand side grows without bounds as
$n\to\infty$.  This implies that there exists an $n \in \N$ such that
$x_n \notin [x^*-\eta, x^*+\eta]$.  Thus, $x^*$ is not stable and we
have also proved part (ii).
\end{proof}

\begin{remark}
When $\lvert f'(x^*) \rvert = 1$, this simple test is inconclusive and
a more sophisticated analysis is required.  Some theorems based on
higher derivatives of $f$ are given in \cite{E96}.  However, there is
no easy criterion that is sufficient and necessary for stability; if
standard results fail, equations have to be investigated on a
case-by-case basis.
\end{remark}

\begin{example}
Re-considering the market model from Problem~\ref{p2} in the framework
of Theorem~\ref{t.stability}, we have, from \eqref{e.pn+1}, that
\begin{gather}
  f(x) = - \frac{m_{\mathrm{s}}}{m_{\mathrm{d}}} \, x
         + \frac{b_{\mathrm{d}}-b_{\mathrm{s}}}{m_{\mathrm{d}}}
\intertext{and therefore}
  f'(x) = - \frac{m_{\mathrm{s}}}{m_{\mathrm{d}}} \,.
\end{gather}
Thus, by Theorem~\ref{t.stability}, the equilibrium point $p^*$ is
asymptotically stable if $m_{\mathrm{s}} < m_{\mathrm{d}}$ and
unstable if $m_{\mathrm{s}} > m_{\mathrm{d}}$.  This coincides with
our previous analysis.
\end{example}

\begin{example}
Consider the difference equation
\begin{equation}
  x_{n+1} = \cos x_n \,.
\end{equation}
Since $f(x)=\cos x$ is decreasing from $1$ to $0$ on the interval
$[0,\pi/2]$ and $g(x)=x$ is increasing from $0$ to $\pi/2$ on the same
interval, the graphs of the two functions must intersect at some point
$x^* \in (0,\pi/2)$.  Similar arguments show that this is the only
point of intersection.  Thus, the difference equation has $x^*$ as its
only equilibrium point.  Since $f'(x)=\sin x$, we must have $\lvert
f'(x^*) \rvert < 1$ and we conclude that the equilibrium point is
asymptotically stable.  (You can try this out on your calculator by
repeatedly hitting the $\cos$-key on any input.)
\end{example}

\section{Periodic Points and Cycles}

In Problem~\ref{p2} we have already encountered the situation where
a solution to a difference equation alternates between two values.
This concept easily generalizes to the situation when the solution
returns to the starting point after cycling through a finite number of
points.  As in the previous section, we consider first order
difference equations of the form
\begin{equation}
  x_{n+1} = f(x_n) \,.
  \label{e.simple2}
\end{equation}

\begin{definition}
We say that $a$ is a $k$-\emph{periodic point} if it is an equilibrium
point of the difference equation 
\begin{equation}
  y_{n+1} = f^k(y_n) \,.
  \label{e.kth-iterate}
\end{equation}
\end{definition}

Note that we write $f^k$ to denote the $k$th composition of $f$ with
itself.  Thus, $f^2(x) = f(f(x))$, $f^3(x)=f(f(f(x)))$, etc.; in this
context, it is \emph{not} the $k$th power of $f$.

\begin{definition}
A $k$-\emph{cycle} is the orbit $\{a, f(a), \dots, f^{k-1}(a)\}$ of a
$k$-periodic point $a$.
\end{definition}

\begin{definition} \label{d.cycle-stability}
We say that a $k$-periodic point $a$ is \emph{stable},
\emph{asymptotically stable}, or \emph{unstable} if $a$ is a stable,
asymptotically stable, or unstable equilibrium point of
\eqref{e.kth-iterate}, respectively.
\end{definition}

\begin{remark}
It can be shown that all the $k$-periodic points in a $k$-cycle have
the same stability property.  Thus, we can say ``$k$-cycle'' instead
of ``$k$-periodic point'' in Definition~\ref{d.cycle-stability}.
\end{remark}

\begin{remark}
It may happen that an orbit starting at some initial point $b$ becomes
periodic only after a finite number of steps.  Such orbits are called
\emph{eventually periodic}.
\end{remark}

Since the discussion of periodic points can be reduced to the
discussion of equilibrium points for the map \eqref{e.kth-iterate},
everything we have said about stability of equilibrium points
naturally translates into results on the stability of cycles.  For
future reference, we state the following elementary stability
criterion for $2$-cycles.

\begin{theorem} \label{t.two-cycle}
Suppose that the difference equation \eqref{e.simple2} has a $2$-cycle
$\{a,b\}$ and that $f$ is continuously differentiable.  Then
\begin{enumerate}
\item the $2$-cycle is asymptotically stable if $\lvert f'(a) \,
f'(b)\rvert<1$,
\item the $2$-cycle is unstable if $\lvert f'(a) \,
f'(b)\rvert>1$.
\end{enumerate}
\end{theorem}

The proof is a direct application of Theorem~\ref{t.stability} to
\eqref{e.kth-iterate}.  The details are left as an exercise to the
reader. 


\section{The logistic equation}

Difference equations are frequently used to model phenomena of
population growth, spread of disease, or market dynamics as in the
example of the previous section.  One of the simplest such population
models is the \emph{logistic difference equation}, which we introduce
in the following.

\begin{problem}
Describe the growth of a population of a single species in a habitat
of constrained resources.
\end{problem}

Let us discuss a set of modeling assumptions that allow us to obtain a
fully specified, yet simple quantitative description.  As for the
market model of Section~\ref{s.stability}, we assume discrete periods
of equal length with regard to physical time, labeled with integer
$n$.  The size of the population in period $n$ is denoted $p_n$.

When resources are plenty, the number of offspring per model period
will be a fraction $b$ of the size of the population in existence.
Likewise, the number of deaths per model period will be another
fraction $d$ of the population size.  Thus,
\begin{equation}
  p_{n+1} = (1+b-d) \, p_n \equiv \mu \, p_n \,.
  \label{e.unconstrained}
\end{equation}
Since this equation describes good, unconstrained conditions, we
expect that $b>d$ and therefore $\mu>1$.  The solution to
\eqref{e.unconstrained}, 
\begin{equation}
  p_n = \mu^n \, p_0 \,,
\end{equation}
then grows without bounds as $n \to \infty$.

A model for the growth of a population in a habitat of constrained
resources must therefore contain a term that limits the growth when
the population size becomes large.  There are many effects that may be
taken into account and models can become very complicated.  However,
with the following simplifying assumptions, we may obtain a very
simple, yet powerful model.  To facilitate the discussion, we assume
that the habitat consists of $s$ discrete sites.  (This assumption is
not essential; we would reach the same conclusion in a corresponding
continuum set-up.)  We further assume that
\begin{enumerate}
\item In each model period, there is a loss of population, for example
due to fights or starvation, proportional to the probability that two
or more individuals occupy the same site.
\item Individuals move through the habitat at random, i.e.\ each
individual is equally likely to occupy each site.
\item The habitat is large and sparsely populated.  More precisely, we
assume that $1 \ll p \ll \sqrt s$. \label{a3}
\end{enumerate}
Thus, it is our task to compute the probability of two individuals
occupying the same site, and to simplify the resulting expression
according to assumption~\eqref{a3}.

We recognize that the problem is precisely the ``birthday problem'' we
discussed earlier.  As in the solution of the birthday problem, we
first compute the probability that no two individuals occupy the same
site, 
\begin{align}
  P(\text{no two on same site})
  & = \frac{s \cdot (s-1) \cdot \dots \cdot (s-p+1)}{s^p} 
      \notag \\
  & = \frac{s!}{s^p \, (s-p)!} 
      \notag \\
  & \sim \frac{\sqrt{2 \pi s} \, s^s \, \mathrm{e}^{-s}}%
              {s^p \, \sqrt{2 \pi (s-p)} \, (s-p)^{s-p} \, 
               \mathrm{e}^{-(s-p)}} 
      \notag \\
  & = \sqrt{\frac{s}{s-p}} \, 
      \exp \biggl( (p-s) \, \ln \frac{s-p}{s} - p \biggr)
\end{align}
We now approximate this expression up to first order terms in the
small quantity $p/s$.  We use two-term Taylor approximations of the
functions
\begin{equation}
  \ln (1+x) \approx x - \frac12 \, x^2 \,, \quad
  \mathrm{e}^x \approx 1 + x \,, \quad
  \frac1{1-x} \approx 1 + x \,, \quad \text{and} \quad
  \sqrt{1+x} \approx 1 + \frac12 \, x \,.
\end{equation}
Then
\begin{equation}
  \ln \frac{s-p}{s} 
  = \ln \biggl( 1 - \frac{p}s \biggr)
  \approx - \frac{p}s - \frac12 \, \frac{p^2}{s^2} \,,
\end{equation}
so that, since $1 \gg p^2/s \gg p^3/s^2$, 
\begin{equation}
  \exp \biggl( (p-s) \, \ln \frac{s-p}{s} - p \biggr)
  \approx \exp \biggl( 
                 (s-p) \, 
                 \biggl( 
                   \frac{p}s + \frac12 \, \frac{p^2}{s^2}
                 \biggr) - p 
               \biggr)
  \approx \exp \biggl( -\frac{p^2}{2s} \biggr)
  \approx 1 - \frac{p^2}{2s} \,.
\end{equation}
Further,
\begin{equation}
  \sqrt{\frac{s}{s-p}} 
  = \sqrt{\frac1{1-p/s}}
  \approx \sqrt{1+\frac{p}s}
  \approx 1 + \frac12 \, \frac{p}{s} \,.
\end{equation}
Altogether, noting that $1 \gg p^2/s \gg p/s \gg p^3/s^2$, 
\begin{equation}
  P(\text{no two on same site})
  \approx \biggl(
            1 + \frac12 \, \frac{p}{s}
          \biggr) \cdot
          \biggl(
            1 - \frac{p^2}{2s} 
          \biggr)
  \approx 1 - \frac{p^2}{2s} 
\end{equation}
and therefore
\begin{equation}
  P(\text{at least two on same site})
  = 1 - P(\text{no two on same site})
  \approx \frac{p^2}{2s} \,.
\end{equation}
We conclude the following.

\begin{theorem} \label{t.sparse} The probability that at least two of
$p$ randomly placed individuals share the same site in a large
sparsely populated habitat is asymptotically proportional to $p^2$.
\end{theorem}

Theorem~\ref{t.sparse} suggests that we should augment the
unconstrained growth equation \eqref{e.unconstrained} by a loss term
proportional to $p_n^2$.  The resulting difference equation is the
\emph{logistic difference equation}
\begin{equation}
  p_{n+1} = \mu \, p_n - \nu \, p_n^2
  \label{e.logistic-unscaled} 
\end{equation}
with positive parameters $\mu$ and $\nu$.  The qualitative behavior is
essentially governed by a single parameter.  Indeed, setting $p_n =
\mu/\nu \, x_n$, we obtain the logistic equation in its standard form
\begin{equation}
  x_{n+1} = \mu \, x_n \, (1 - x_n) \,.
  \label{e.logistic}
\end{equation}

\begin{remark}
In this model we allow the number of individuals to be any positive
real number.  When the number of individuals is large, this is a
useful simplification.  However, when the number of individuals is
small, the discreteness of $p$ becomes significant.  Moreover, in real
population dynamics, random fluctuations can have dramatic impact on
small populations.
\end{remark}

\begin{remark} \label{r.time}
If we reconsider the derivation leading to Theorem~\ref{t.sparse}
carefully, we see that the constant $\nu$ must be small and inversely
proportional to the artificial parameter $s$.  For purposes of this
discussion, let $\nu=1/s$.  Then a nontrivial equilibrium point of the
logistic difference equation \eqref{e.logistic-unscaled} is $p^*=s \,
(\mu-1)$.  (A more detailed discussion of equilibrium points is given
in the next section.)  Since the equilibrium point should not depend
on the artificial parameter $s$, we must therefore require that $\mu$
be very close to $1$ as $s$ becomes large.  This corresponds to
choosing units of time where one iteration of the logistic map
corresponds to a very small amount of real time.
\end{remark}

\begin{remark}
The statement of Remark~\ref{r.time} can also be re-interpreted as
follows.  Choose the arbitrary normalization $p^*=1$, so that
$\mu=1+1/s$.  Then \eqref{e.logistic-unscaled} coincides with the
so-called explicit Euler scheme \cite{Euler} for the numerical
solution of the \emph{logistic differential equation}
\begin{equation}
  \frac{\mathrm{d} p}{\mathrm{d} t} = p - p^2 \,,
\end{equation}
commonly written in the form
\begin{equation}
  p_{n+1} = p_n + h \, (p_n - p_n^2) \,,
\end{equation}
where $h=1/s$ is the so-called step size.  In other words, our
derivation of the model is valid precisely in the parameter regime
where the logistic difference equation is a good approximation to the
logistic differential equation.
\end{remark}

\section{Bifurcations}

In the following, we study the qualitative behavior of the first order
difference equations of the form \eqref{e.simple-iteration} as a
function of a parameter.  By \emph{qualitative behavior} we mean the
number of equilibrium points, of $k$-cycles, and their stability.  We
speak of a \emph{bifurcation} when a small change in the parameter
causes a change in the qualitative behavior of the equation.

We illustrate the concept using the logistic equation in the form
\eqref{e.logistic} where $\mu$ is the parameter we will vary.  

\paragraph{Equilibrium points.}

We first determine the equilibrium points of the logistic difference
equation.  Setting
\begin{equation}
  f(x) = \mu \, x \, (1 - x) \,,
\end{equation}
the fixed point equation $x=f(x)$ has two solutions
\begin{equation}
  x^*_1 = 0 \quad \text{and} \quad x^*_2 = \frac{\mu-1}\mu \,.
\end{equation}
Since $f'(x) = \mu - 2 \, \mu \, x$, we have, in particular,
\begin{equation}
  f'(x^*_1) = \mu \quad \text{and} \quad 
  f'(x^*_2) = 2 - \mu \,.
\end{equation}
By Theorem~\ref{t.stability},
\begin{itemize}
\item $x^*_1$ is an asymptotically stable equilibrium point for
$\lvert \mu \rvert < 1$ and an unstable equilibrium point for $\lvert
\mu \rvert > 1$.

\item $x^*_2$ is an asymptotically stable equilibrium point for
$1<\mu<3$ and an unstable equilibrium point for $\mu < 1$ or $\mu >
3$. 
\end{itemize}
The points $\mu=1$ and $\mu=3$ are bifurcation points since the
stability of the equilibrium points is changing as the parameter $\mu$
varies across.

\paragraph{$2$-cycles.}

It can be shown (homework exercise!) that the logistic map has a
$2$-cycle whenever $\mu>3$ with orbit $\{a,b\}$ where
\begin{gather}
  a = \frac{1+\mu - \sqrt{(\mu-3)(\mu+1)}}{2 \, \mu} \,, \\
  b = \frac{1+\mu + \sqrt{(\mu-3)(\mu+1)}}{2 \, \mu} \,. \\
\end{gather}
Moreover, by Theorem~\ref{t.two-cycle}, this $2$-cycle is
asymptotically stable for $3 < \mu < 1 + \sqrt 6$ and is unstable for
$\mu > 1 + \sqrt 6$.

Thus, not only does the equilibrium point $x^*_2$ lose stability at
$\mu=3$, but a stable $2$-cycle emerges.  There is a third bifurcation
point at $\mu = 1 + \sqrt 6$ where the $2$-cycle loses stability and,
as can be shown by lengthy calculation, a stable $4$-cycle emerges.

The behavior analyzed so far is summarized in Figure~\ref{f.1}.  This
graph is called a \emph{bifurcation diagram}.  The horizontal axis
shows the parameter while the vertical axis shows the location of
equilibrium points and cycles.  In standard terminology, the
bifurcation at $\mu=1$ is called a \emph{transcritical bifurcation}
and the bifurcation at $\mu=3$ a \emph{pitchfork bifurcation}.

\begin{figure}[htb]
\centering
  \includegraphics[width=15cm]{bifurcation.pdf}
\caption{Partial bifurcation diagram for the logistic map.  Shown are
the branches which we have explicitly computed.  Unstable branches are
dotted, stable branches indicated by solid lines.  The equilibrium
points are in black, the $2$-cycle in gray.  At the point where the
$2$-cycle becomes unstable, a stable $4$-cycle emerges.  This
$4$-cycle and its further bifurcations are \emph{not} shown.}  
\label{f.1}
\end{figure}

\paragraph{Higher cycles and chaos.}

The analysis of the bifurcation structure beyond the $2$-cycle is
rather complicated; \cite{W} has some more information and references
to the research literature.  In particular, there are no closed form
expressions for most branches of the bifurcation diagram.

On the other hand, it is comparatively easy to find the stable cycles
by numerical computation.  For each value of $\mu$, choose an
arbitrary starting point and iterate the logistic map a large number
of times.  Since an asymptotically stable cycle is attracting, the
value of $x_n$ for $n$ large will now be very close to the stable
cycle.  Keep iterating and plot the next many iterates---they will now
approximately trace out the stable cycle in the bifurcation diagram.
Figure~\ref{f.2} shows the result of such a computation (with the
additional twist that the graph is shaded according to the probability
of a certain pixel being visited during the iteration).

\begin{figure}[htb]
\centering
  \includegraphics[angle=270,width=15cm]{LogisticMap_BifurcationDiagram.pdf}
\caption{Numerically computed bifurcation diagram for the logistic
map.  The computation can only show the attracting sets.  The
graph is taken from \cite{logistic} where the method of computation is
explained in detail.  Their parameter $r$ is the same as $\mu$ in
\eqref{e.logistic}.} 
\label{f.2}
\end{figure}

Note that for each initial value you will approach exactly one stable
equilibrium point or cycle.  For the logistic map it is known that
only one starting point  will suffice.  In general, additional analysis
or systematic search is required to verify that all possible
attracting sets have been found.

Through a combination of analysis and numerical computation, the
following picture has been established \cite{S00,W,logistic}.

\begin{itemize}
\item As $\mu$ is increased, each $2^k$-cycle has a bifurcation point
where it becomes unstable and bifurcates into a stable
$2^{k+1}$-cycle.  

The first such bifurcation point is the break-up of the stable
equilibrium point or, equivalently, $2^0$-cycle into the $2^1$-cycle
at $\mu=3$ as shown in Figure~\ref{f.1}.  All other such bifurcation
points are qualitatively similar.

\item If $\ell_k$ denote the interval on the $\mu$-axis in which the
$2^k$ cycle is stable, then
\begin{equation}
  \lim_{k \to \infty} \frac{\ell_k}{\ell_{k+1}}
  = \delta = 4.6692016\dots \,.
\end{equation}
This number is called the \emph{Feigenbaum constant} and is universal
across large classes of maps.  In particular, the fact that this limit
exists implies that this cascade of bifurcations, often referred to as
a \emph{period doubling cascade}, takes place in a finite interval on
the $\mu$-axis.  For the logistic map, the cascade ends at $\mu
\approx 3.57$, as can be seen in Figure~\ref{f.2}.

\item For $3.57 \lesssim \mu < 4$, cycles of all periods exist.  For
most values of $\mu$, all cycles are unstable and we speak of
\emph{chaos}.  However, in the interior of this interval there are
\emph{islands of stability} in which stable $3, 5, \dots$-periodic
cycles exist, each with their own period doubling cascade.

\item Beyond $\mu=4$, solutions diverge for almost all initial
values. 

\end{itemize}

We close this section with a short remark on chaos.  Mathematical
definitions of ``chaos'' are somewhat problem dependent, but usually
involve the following features.

\begin{itemize}
\item \emph{Sensitive dependence on initial conditions.}  More
precisely, if the system is started from any two nearby initial
points, then the distance between the resulting solutions grows
exponentially in time, or with the number of iterations, until it is
on the order of the system size.

\item \emph{Dense orbits.} Every orbit of the system comes arbitrarily
close to every point in the domain of definition, or at least in a
``massive'' subset of the domain of definition, called a \emph{strange
attractor}. 

\end{itemize}

The notion of chaos is important for a number of reasons.  First, it
is a mathematically challenging task to describe and quantify the
behavior of chaotic systems in meaningful ways.  A lot of new and
important mathematics has emerged from the pursuit of such questions.
Second, it shows that very simple systems can have very complicated
behavior---the logistic difference equation being a case in point.
Third, the presence of chaotic dynamics implies severe practical
limits to the predictive power of a model.  Since orbits diverge and
mix at an exponential rate, doubling the accuracy of the initial
values (which are typically measurements) will not double the time
interval over which an accurate predication can be achieved, but only
increase this time interval by a fixed amount.  This severely
restricts the ability to achieve accurate long-range predictions.


\section{Generating functions revisited}

Let us now consider difference equations of order greater than one.
Earlier this semester, we discussed the method of generating functions
for the solution of certain difference equations.  In the following,
we review the method and state the theory from a different viewpoint.

In this section, we assume that the difference equation is
\emph{linear} and \emph{homogeneous} with \emph{constant
coefficients}.  For simplicity, we will only consider \emph{second
order} equations, although everything we do will generalizes to higher
order in a natural way.  Such an equation can be written in the form
\begin{equation}
  \label{e.o2-linear}
  a_{n+1} + c_1 \, a_n + c_2 \, a_{n-1} = 0 \,;
\end{equation}
for the time being, we leave the initial conditions unspecified.  Now
define the generating function
\begin{equation}
  \phi(x) = \sum_{i=0}^\infty a_i \, x^i \,.
\end{equation}
Then, by expanding the left hand expression and ordering terms by
powers in $x$, we find that
\begin{equation}
  (1 + c_1 \, x + c_2 \, x^2) \, \phi(x)
  = a_0 + (a_1 + c_1 \, a_0) \, x \,,
\end{equation}
or
\begin{equation}
  \phi(x) = \frac{a_0 + (a_1 + c_1 \, a_0) \, x}%
                 {1 + c_1 \, x + c_2 \, x^2} \,.
  \label{e.phi1}
\end{equation}
Thus, we need to find the explicit power series expansion of
\eqref{e.phi1} whose coefficients are the solution to the difference
equation \eqref{e.o2-linear}.  Let us recall that, in the ring of
formal power series (cf.\ \cite{I}, p.\ 22),
\begin{equation}
  \frac1{1-x} = 1 + x + x^2 + x^3 + \dots 
  \label{e.geometric}
\end{equation}
and
\begin{equation}
  \frac1{(1-x)^2} = 1 + 2 \, x + 3 \, x^2 + 4 \, x^3 + \dots \,.
  \label{e.series2}
\end{equation}
Let us denote the roots of the polynomial $1 + c_1 \, x + c_2 \, x^2$
by $\mu_1$ and $\mu_2$.  (In the field of complex numbers, these roots
always exist.  In the following, we will consider examples with real
roots only; the general discussion, however, applies equally to the
complex case.)  We must distinguish two cases.

\paragraph{Two distinct roots.}

This is the case we encountered previously, e.g.\ when determining the
general form of the Fibonacci sequence.  In this case, we know that
$\phi$ has a partial fraction decomposition of the form
\begin{align}
  \phi(x) & = \frac{A'}{\mu_1 - x} + \frac{B'}{\mu_2 - x} 
              \notag \\
          & = \frac{A}{1 - \dfrac{x}{\mu_1}} + 
              \frac{B}{1 - \dfrac{x}{\mu_2}}
              \notag \\
          & = A \, \biggl(
                     1 + \frac{x}{\mu_1} + \frac{x^2}{\mu_1^2} +
                     \dots 
                   \biggr) +
              B \, \biggl(
                     1 + \frac{x}{\mu_2} + \frac{x^2}{\mu_2^2} +
                     \dots 
                   \biggr) \,,
  \label{e.partial-fraction}
\end{align}
where $A$ and $B$ are some constants which depend on $a_0$ and $a_1$
and which are determined by the partial fraction decomposition, and
where $A'=A \, \mu_1$ and $B'=B \, \mu_2$.  We now take the point of
view that, although we could work out the values of the constants
through the computation above, we simply note that this computation
shows that the general form of the solution to the difference equation
\eqref{e.o2-linear} is
\begin{equation}
  a_n = A \, \mu_1^{-n} + B \, \mu_2^{-n} \,.
  \label{e.sol-two-roots}
\end{equation}
The constants $A$ and $B$ can then be determined by the condition that
the initial conditions are replicated by \eqref{e.sol-two-roots}.

\paragraph{Double Roots.}

When the two roots $\mu_1 = \mu_2 \equiv \mu $ coincide---we speak of
a \emph{double root}---a partial fraction decomposition of the form
\eqref{e.partial-fraction} is not possible.  Instead, we refer to
\eqref{e.series2} and write
\begin{align}
  \phi(x) & = \frac{A' \, (\mu-x) + B'}{(\mu-x)^2}
              \notag \\
          & = \frac{A}{1 - \dfrac{x}\mu} + 
              \frac{B}{\left( 1 - \dfrac{x}\mu \right)^2}
              \notag \\
          & = A \, \biggl(
                     1 + \frac{x}{\mu} + \frac{x^2}{\mu^2} +
                     \dots 
                   \biggr) +
              B \, \biggl(
                     1 + 2 \, \frac{x}{\mu} + 3 \, \frac{x^2}{\mu^2} +
                     \dots 
                   \biggr) \,.
\end{align}
As in the case of two distinct roots, $A$ and $B$ are two constants
which depend on $a_0$ and $a_1$.  In principle, they could be
determined by explicitly carrying all constants through the
computation above.  Instead, however, we only note that the general
form of the solution to the difference equation \eqref{e.o2-linear} is
\begin{equation}
  a_n = A \, \mu^{-n} + B \, (n+1) \, \mu^{-n} \,.
  \label{e.sol-double-root}
\end{equation}
As before, the constants $A$ and $B$ can be determined by the
condition that the initial conditions are replicated by
\eqref{e.sol-double-root}.

We have now effectively proved a theorem on the general solution of
\eqref{e.o2-linear}.  We can rewrite the final result in a form more
commonly encountered by making the following modifications to the
expressions above.  First, notice that we can make the prefactor of
the second term in \eqref{e.sol-double-root} proportional to $n$
rather than $n+1$.  The ``missing'' bit can then be grouped together
with the first term, forcing a redefinition of the constant $A$.
Second, $\mu$ is a root of $1+ c_1 \, x + c_2 \, x^2$ if and only if
$\lambda = 1/\mu$ is a root of $x^2 + c_1 \, x + c_2$.  We can
therefore express the solution in terms of roots of the latter
polynomial.  Finally, we can rescale the difference equation
\eqref{e.o2-linear} by an arbitrary constant without changing any of
the above procedure.  With these modifications, we restate our result
as follows.

\begin{definition}
The \emph{characteristic polynomial} of the difference equation
\begin{equation}
  \label{e.o2-linear-2}
  c_0 \, a_{n+1} + c_1 \, a_n + c_2 \, a_{n-1} = 0 
\end{equation}
is the polynomial
\begin{equation}
  p(x) = c_0 \, x^2 + c_1 \, x + c_2 \,.
  \label{e.cp}
\end{equation}
\end{definition}

\begin{theorem} \label{t.o2} Let $\lambda_1$ and $\lambda_2$ denote
the roots of the characteristic polynomial \eqref{e.cp}.  Then any
solution to \eqref{e.o2-linear-2} is of the form
\begin{enumerate}
\item when $\lambda_1 \neq \lambda_2$, then $a_n = A \, \lambda_1^n + B
\, \lambda_2^n$,
\item when $\lambda_1 = \lambda_2 \equiv \lambda$, then $a_n = A \,
\lambda^n + B \, n \, \lambda^n$,
\end{enumerate}
for arbitrary constants $A$ and $B$.
\end{theorem}

\begin{remark}
Strictly speaking, our derivation breaks down when one of the roots is
zero.  For a true second order equation, however, $c_2 \neq 0$ and
this cannot happen.  Otherwise we are back to the case of first order
equations discussed in Section~\ref{s.first-order}.
\end{remark}

\begin{remark}
\emph{Linear Algebra} provides an alternative way to prove
Theorem~\ref{t.o2}.  The difference equation can be written as a
system of two first-order difference equations.  It turns out that the
characteristic polynomial of the difference equation equals the
characteristic polynomial of the resulting matrix.  When the roots are
distinct, the system can be solved by diagonalization; in the case of
a double root, the Jordan normal form is used to solve the equation.
\end{remark}

We now give two typical application examples of the theory.  The first
example also illustrates that Theorem~\ref{t.o2} is particularly
useful when the difference equation is augmented not by initial
conditions, but by some other conditions on the solution.  In this
case, a straightforward application of the method of generating
functions would not be possible.

\begin{problem}[Gambler's ruin]
A gambler is repeatedly playing a game in which he is to gain $1$ Euro
with probability $q$ and lose $1$ Euro with probability $1-q$.  Assume
that the gambler has $n$ Euros to play.  When he is down to zero, we
say that the gambler is ruined and cannot continue gambling.  He will
also quit when he has reached a pre-set target of $N$ Euros.  What is
the probability that the gambler is ruined?
\end{problem}

Let $p_n$ denote the probability of being ruined when owning $n$
Euros.  If $n=0$, the gambler is ruined for sure, thus $p_0=1$.  When
$n=N$, the gambler quits without being ruined, so $p_N=0$.  For
$0<n<N$, the gambler keeps playing.  If he wins, he will be ruined
with probability $p_{n+1}$ in the next iteration of the game; if he
loses, he will be ruined with probability $p_{n-1}$.  Since the first
alternative occurs with probability $q$ and the second alternative
occurs with probability $1-q$, we deduce that
\begin{equation}
  p_n = q \, p_{n+1} + (1-q) \, p_{n-1} \,.
  \label{e.gambler}
\end{equation}
This is clearly a difference equation of the form
\eqref{e.o2-linear-2} with characteristic polynomial
\begin{equation}
  p(x) = q \, x^2 - x + (1-q)
\end{equation}
whose roots are 
\begin{equation}
  \lambda_{1,2} 
  = \frac{1 \pm \sqrt{1 - 4 \, q \, (1-q)}}{2 \, q} 
  = \frac{1 \pm \sqrt{(1-2 \,q)^2}}{2 \, q} 
\end{equation}
so that 
\begin{equation}
  \lambda_1 = 1 
  \qquad \text{and} \qquad
  \lambda_2 = \frac{1-q}{q} \,.
  \label{e.roots}
\end{equation}
We can now apply Theorem~\ref{t.o2}.  There are two distinct cases.

\paragraph{Fair game.} 

When the probability of losing and of winning the game is both
$\tfrac12$, we speak of a \emph{fair game}.  In this case, both roots
of the characteristic polynomial equal $\lambda=1$, so that the
general solution to \eqref{e.gambler} is
\begin{equation}
  p_n = A + B \, n \,.
\end{equation}
Since $p_0 = 1$, we must have $A=1$.  Then $p_N=0$ gives $1+B\,N=0$ so
that $B=-1/N$ and therefore
\begin{equation}
  p_n = 1 - \frac{n}{N} \,.
\end{equation}
In particular, if the gambler owns a finite amount of money $n$, but
does not quit gambling, we can represent this by taking the limit $N
\to \infty$.  In this case, the gambler will be ruined with
probability $1$.

\paragraph{Biased game.} 

When $q \neq \tfrac12$, the two roots in \eqref{e.roots} are distinct.
In this case,
\begin{equation}
  p_n = A + B \, \biggl( \frac{1-q}{q} \biggr)^n \,.
  \label{e.pnbiased}
\end{equation}
Since $p_0 = 1$, we must have $A + B = 1$.  Moreover, since $p_N=0$,
we find that
\begin{equation}
  A + B \, \biggl( \frac{1-q}{q} \biggr)^N = 0 \,.
\end{equation}
Solving for $A$ and $B$, we obtain
\begin{equation}
  A = - \frac{\biggl( \dfrac{1-q}{q} \biggr)^N}%
             {1 - \biggl(  \dfrac{1-q}{q} \biggr)^N} 
  \qquad \text{and} \qquad
  B = \frac1{1 - \biggl( \dfrac{1-q}{q} \biggr)^N} \,.
\end{equation}
Inserting this result back into \eqref{e.pnbiased}, we finally obtain
\begin{equation}
  p_n = \frac{\biggl( \dfrac{1-q}{q} \biggr)^n - 
              \biggl( \dfrac{1-q}{q} \biggr)^N}{1 - 
              \biggl( \dfrac{1-q}{q} \biggr)^N} \,.
\end{equation}
Let us again ask the question what happens when the gambler does not
quit gambling, i.e.\ take the limit $N \to \infty$ with $n$ fixed.
If $q < \tfrac12$, then $(1-q)/q>1$ so that
\begin{equation}
  p_n (N) \to 1 \quad \text{as } N \to \infty \,,
\end{equation}
i.e., the gambler is ruined for sure; if $q > \tfrac12$, then
$(1-q)/q<1$ so that
\begin{equation}
  p_n (N) \to \biggl( \frac{1-q}{q} \biggr)^n
  \quad \text{as } N \to \infty \,.
\end{equation}
We see that in a game which is biased in favor of the gambler, there
is non-zero probability of ``survival'' even if he will never quit.

We summarize these results as follows.

\begin{theorem}[Gambler's theorem]
If you don't quit playing a game that's not in your favor, you will be
ruined for sure.
\end{theorem}

We close this section with an example from \cite{S00} which
illustrates interpretation and possible qualitative behavior of a
second order linear difference equation or, equivalently, of a
two-variable first order system of difference equations.

\begin{problem}[Romeo and Juliet]
Model the onset of a love/hate relationship between Romeo and Juliet,
taking into account 
\begin{enumerate}
\item the change of feeling for each other while apart, and
\item the change of feeling for each other in response to the other's
emotions. 
\end{enumerate}
\end{problem}

Again, we construct a discrete time model with $n$ denoting units of
time, say days.  Romeo's feelings for Juliet on day $n$ are quantified
by the variable $R_n$, positive values representing love and negative
values representing hate.  Similarly, Juliet's feelings for Romeo are
quantified by the variable $J_n$.  Before Romeo and Juliet meet on day
$n=0$, their feelings for each other are neutral, so $J_0=R_0=0$.  We
will first formulate a simple model for the dynamics of their
relationship and then investigate whether zero is a stable equilibrium
point, or whether a small perturbation may evolve into something
better (or worse).

When Romeo and Juliet are on their own, $R_n$ and $J_n$ evolve
independent of each other.  Let us therefore focus on $R_n$ only.  We
assume a simple linear model,
\begin{gather}
  R_{n+1} = a_R \, R_n 
\end{gather}
where $a_R$ is a positive parameter.  If $a_R$ was negative, Romeo's
feeling toward Julia would alternate from love to hate on consecutive
days, which is not what we would expect.  When $0<a_R<1$, then Romeo's
intrinsic love for Juliet fades out, when $a_R>1$, it increases with
time.

In the presence of Julia, Romeo also reacts to her feelings toward
him.  This feedback is modeled by a second term in Romeo's equation,
\begin{gather}
  R_{n+1} = a_R \, R_n + p_R \, J_n \,,
\end{gather}
where the feedback coefficient $p_R$ can be either positive or
negative.  The different possible combinations of $a_R$ and $p_R$
give rise to four romantic styles:
\begin{center}
\begin{tabular}{c|cc}
 & $p_R < 0 \displaystyle\vphantom\int$ & $p_R>0$ \\
\hline
$0 < a_R < 1 \displaystyle\vphantom\int$ & hermit & cautious lover \\ 
$a_R>1\displaystyle\vphantom\int $       & likes to tease, but not to
please & eager beaver  
\end{tabular}
\end{center}

The equation for Juliet's feeling toward Romeo shall be symmetric to
Romeo's equation, so that the combined system of difference equations
reads
\begin{subequations}
\begin{gather}
  R_{n+1} = a_R \, R_n + p_R \, J_n \,, \label{e.r} \\
  J_{n+1} = a_J \, J_n + p_J \, R_n \,. 
\end{gather}
\end{subequations}

We can now investigate the outcome of their love affair by studying
these difference equations using Theorem~\ref{t.o2}.  The crucial
observation is that a system of two difference equations can always be
formulated as a single second order difference equation, and vice
versa.  Instead of proving this general statement, we proceed by
example.  Shifting the index of the first equation and multiplying the
second by $p_R$, we obtain
\begin{subequations}
\begin{gather}
  R_{n+2} = a_R \, R_{n+1} + p_R \, J_{n+1} \,, \\
  p_R \, J_{n+1} = p_R \, a_J \, J_n + p_R \, p_J \, R_n \,,
\end{gather}
\end{subequations}
so that
\begin{gather}
  R_{n+2} 
  = a_R \, R_{n+1} + p_R \, a_J \, J_n + p_R \, p_J \, R_n \,
\end{gather}
The remaining $J_n$ can be eliminated by using the unshifted version
of \eqref{e.r}, and we finally obtain
\begin{gather}
  R_{n+2} - (a_R + a_J) \, R_{n+1} 
  + (a_J \, a_R - p_J \, p_R) \, R_n = 0 \,.
  \label{e.r2}
\end{gather}
The characteristic polynomial of this difference equation is
\begin{equation}
  p(x) = x^2 - (a_R + a_J) \, x + a_J \, a_R - p_J \, p_R
\end{equation}
which has roots
\begin{align}
  \lambda_{1,2}
  & = \frac{a_J + a_R \pm 
            \sqrt{(a_R+a_J)^2 - 4 \, (a_J \, a_R - p_J \, p_R)}}2
      \notag \\
  & = \frac{a_J + a_R \pm \sqrt{(a_R-a_J)^2 + 4 \, p_J \, p_R}}2 \,.
  \label{e.rjroots}
\end{align}
A comprehensive analysis of the solution is rather involved and best
done in the linear algebra framework (as this would allow to make
simultaneous statements about Romeo's and Juliet's feelings).  However,
we can make a few basic observations from \eqref{e.rjroots}:
\begin{itemize}
\item When $a_R$ and $a_J$ are similar in magnitude and $p_J$ and
$p_R$ have opposite signs (i.e.\ exactly one of the two counter-reacts
to the other's feelings), then the expression under the square root
may be negative.  In this case, the difference equation can still be
solved, but the roots are complex numbers which manifests itself in
oscillatory behavior of the solutions.  The relationship then gets
into a love-hate cycle which may grow or decay in amplitude; see the
right hand graph in Figure~\ref{f.romeo-juliet}.

\item A relationship can only develop if at least one of the roots is
greater than one, so that the general solution from Theorem~\ref{t.o2}
has at least one growing component.

\item We cannot infer from \eqref{e.r2} that if Romeo's love for
Juliet is growing that Juliet's love for Romeo is also growing.  We
need to explicitly look at the corresponding Juliet solution.

\item Initialization of this model is very important.  In a parameter
regime where love can grow, hate can also grow.  As the left graph in
Figure~\ref{f.romeo-juliet} shows, initial configurations which are
close can result in completely different outcomes.

\item This model can only represent the onset of an affair, as there
are no limiting terms which would keep initially growing solutions
bounded as $n\to\infty$.  This situation is very much like for the
simple, unconstrained model \eqref{e.unconstrained} as opposed to the
full logistic equation for the dynamics of a one-species population.

\end{itemize}

\begin{figure}[tb]
\centering
\includegraphics[width=15cm]{romeo-juliet.pdf}
\caption{The Romeo-Juliet system with neutral self-induced growth
  coefficients $a_R=a_J=1$ and positive coupling coefficients
  $p_R=p_J=0.1$ (left) vs.\ mixed coupling coefficients $p_R=-0.1$,
  $p_J=0.1$ (right).  The shaded region in the left hand graph
  indicates initial configurations which evolve into growing mutual
  hate.  On the right, the relationship evolves into a love-hate cycle
  independent of the initial conditions.}
\label{f.romeo-juliet}
\end{figure}

\section{Inhomogeneous linear equations}

The theory developed in the previous section may seem rather special.
However, there are two important principles that show that the study
of linear homogeneous problems is relevant in a much more general
context.

The first principle concerns the general solution to inhomogeneous
linear equations of the form \eqref{e.inhomogeneous}, and will be
formulated below.  The next principle concerns the study of
equilibrium points to nonlinear equations and will be discussed in
Section~\ref{s.nonlinear} below.

To formulate the statement in succinct terms, we need to introduce the
concept of linear independence.  Let $x_n^{(1)}, \dots, x_n^{(j)}$
denote $j$ sequences (for example, $j$ solutions to a difference
equation with different starting values).  Then we say that
$x_n^{(1)}, \dots, x_n^{(j)}$ are \emph{linearly independent} or
simply \emph{independent} provided that
\begin{equation}
  \alpha_1 \, x_n^{(1)} + \dots + \alpha_j \, x_n^{(j)} = 0
\end{equation}
for every $n \in \N$ implies that $\alpha_1 = \dots = \alpha_j = 0$.
Without discussing the geometric content of this statement (this can
be found in any book on Linear Algebra), we simply note that in the
case of $j=2$ sequences, the condition for independence reduces to the
simple statement that one sequence is not a multiple of the other.

Recall the inhomogeneous linear equation \eqref{e.inhomogeneous},
\begin{gather}
  x_{n+1} = a_0 (n) \, x_n + a_1 (n) \, x_{n-1} + \dots + 
            a_{k-1}(n) \, x_{n-k+1} + b (n) \,.
  \label{e.i}
\end{gather}
We say that the equation 
\begin{gather}
  x_{n+1} = a_0 (n) \, x_n + a_1 (n) \, x_{n-1} + \dots + 
            a_{k-1}(n) \, x_{n-k+1} \,.
  \label{e.h}
\end{gather}
is the \emph{associated homogeneous} equation.  Then the following is
true.

\begin{theorem}[Superposition Principle]
Let $p_n$ solve the inhomogeneous equation \eqref{e.i} for some
arbitrary starting values.  Further, let $y_n^{(1)}, \dots, y_n^{(k)}$
solve the associated homogeneous equation \eqref{e.h}.  Set
\begin{equation}
  x_n = p_n + \alpha_1 \, y_n^{(1)} + \dots + \alpha_n \, y_n^{(k)} 
  \label{e.general-solution}
\end{equation}
for some $\alpha_1, \dots, \alpha_k$.  Then
\begin{enumerate}
\item $x_n$ is a solution to the inhomogeneous equation \eqref{e.i}.
\item Vice versa, any solution to \eqref{e.i} can be written in the
form \eqref{e.general} provided that $y_n^{(1)}, \dots, y_n^{(k)}$ are
independent.
\end{enumerate}
\end{theorem}

Part (i) is easily proved by direct substitution of
\eqref{e.general-solution} into \eqref{e.i}.  We leave this as an
exercise to the reader.  The proof of part (ii) involves nontrivial
Linear Algebra and shall therefore be omitted.

When the coefficients $a_1, \dots, a_k$ depend on $n$, this theorem,
though true and important for the development of the theory, is of not
much practical help.  In fact, it can be argued that such equations
are really nonlinear ``in disguise'' and that we therefore should not
expect to be able to develop a universal solution procedure.

If, on the other hand, the coefficients $a_1, \dots, a_k$ are
independent of $n$, then Theorem~\ref{t.o2} (respectively its
generalization to order $k>2$) provides a method for solving the
homogeneous problem in a systematic way.  So one only needs to guess a
single inhomogeneous solution, which is often possible, to
automatically obtain \emph{all} solution to the inhomogeneous
equation.

\begin{example}[Model for GDP growth]
An idealized model for the evolution of the gross domestic product
(GDP) of a country assumes that the GDP in year $n$ is given as the
sum of three components,
\begin{equation}
  y_n = c_n + g_n + i_n \,,
\end{equation}
where $c_n$ denotes the consumer spending, assumed to be proportional
to the current economic climate as quantified by the GDP, $g_n$
denotes the government spending, assumed constant, and $i_n$ the
induced investment, assumed to be proportional to the growth rate of
the GDP during the previous year.  Thus, we can write
\begin{equation}
  y_n = \alpha \, y_n + \beta \, (y_{n-1} - y_{n-2}) + \gamma \,,
\end{equation}
for some nonnegative parameters $\alpha$, $\beta$, and $\gamma$ with
$\alpha<1$.   
\end{example}

The superposition principle for this model takes a very simple form.
A particular solution is given by the equilibrium solution, where 
\begin{equation}
  y^* = \alpha \, y^* + \beta \, (y^* - y^*) + \gamma \,,
  \label{e.gdp}
\end{equation}
so that
\begin{equation}
  y^* = \frac\gamma{1-\alpha} \,.
\end{equation}
Two independent solutions to the homogeneous problem can be obtained
via the method of generating functions, e.g.\ in the form of
Theorem~\ref{t.o2}.  We first illustrate this for the special case
where $\alpha=\tfrac12$ and $\beta=\gamma=1$.  Then $y^*=2$ is the
equilibrium solution, and the characteristic polynomial for the
homogeneous equation is given by
\begin{equation}
  p(x) = \tfrac12 \, x^2 - x + 1
\end{equation}
with roots
\begin{equation}
  \lambda_{1,2} = 1 \pm \i = \sqrt2 \, \mathrm{e}^{\i \pi/4} \,.
\end{equation}
The general solution is therefore
\begin{equation}
  y_n = 2 + c_1 \, \lambda_1^n + c_2 \, \lambda_2^n \,.
\end{equation}
Since $\lambda_1=\overline \lambda_2$, where $\overline \lambda_2$
denotes the complex conjugate of $\lambda_2$, the solution is
real-valued if and only if $c_1 = \overline c_2$.  Setting $c_1 = a +
\i b$ and recalling that $\mathrm{e}^{\i\theta} = \cos \theta +\i \,
\sin \theta$, we then find that
\begin{equation}
  y_n = 2 + 2 \Re (c_1 \, \lambda_1^n)
      = 2 + 2^{n/2+1} \, 
        \bigl(a \, \cos (n \pi/4) - b \, \sin (n \pi/4) \bigr) \,,
\end{equation}
where $a$ and $b$ are determined by the initial conditions.  We see
that the solution is oscillating about the equilibrium point with
growing amplitude.  Thus, the given parameter values describe an
unstable economy.

More generally, the roots of the characteristic polynomial to the
homogeneous part of \eqref{e.gdp} are given by 
\begin{equation}
  \lambda_{1,2} 
  = \frac{\beta \pm \sqrt{\beta^2 - 4(1-\alpha)\beta}}%
         {2(1-\alpha)} \,. 
\end{equation}
In the oscillatory case, when $4(1-\alpha)>\beta$, so that
$\lambda_{1,2}$ are a pair of complex conjugate roots, we compute for
$i=1,2$ that
\begin{equation}
  \lvert \lambda_i \rvert^2
  = \lambda_1 \, \lambda_2 = \frac\beta{1-\alpha} \,,
\end{equation}
so that the amplitude of the oscillations is growing whenever
$\beta>1-\alpha$, remains constant for $\beta=1-\alpha$, and decreases
for $\beta<1-\alpha$.  We conclude that the economy is stable if and
only if $\beta \leq 1-\alpha$.

The nonoscillatory case can be analyzed similarly; this shall be left
as an exercise to the reader.


\section{Equilibrium points for nonlinear equation of higher order}
\label{s.nonlinear}

For nonlinear equation of higher order, explicit solutions are rare.
Thus, we would like to extend the study of equilibrium points and
cycles to this situation, so that at least some qualitative properties
can be deduced in the absence of fully explicit solutions.

In analogy with the setting of Section~\ref{s.stability}, particularly
equation \eqref{e.general2}, we consider constant coefficient
nonlinear second order equations in the form
\begin{equation}
  x_{n+2} = f(x_{n+1}, x_n) \,.
  \label{e.nl-order2}
\end{equation}
An equilibrium point $x^*$ is then characterized by
\begin{equation}
  x^* = f(x^*, x^*) \,.
  \label{e.epoint2}
\end{equation}
We are again interested in the question of stability of such
equilibrium points.  To provide an intuitive argument, which could
easily be made rigorous, we write
\begin{equation}
  x_n = x^* + \delta_n 
  \label{e.perturb}
\end{equation}
where we think of $\delta_n$ as an initially small perturbation to an
equilibrium point; we want to study whether the perturbation
$\delta_n$ grows or decreases.  Plugging \eqref{e.perturb} into
\eqref{e.nl-order2}, we can obtain an approximate simpler expression
on the right hand side by using a tangent plane approximation---the
first two terms in a multivariate Taylor expansion---at the point
$x^*$, whence
\begin{equation}
  x^* + \delta_{n+2} 
  \approx f(x^*, x^*)
          + \frac{\partial f(x^*, x^*)}{\partial x} \, \delta_{n+1} 
          + \frac{\partial f(x^*, x^*)}{\partial y} \, \delta_{n} \,,
  \label{e.linearized}
\end{equation}
where $\partial f/\partial x$ denotes the derivative of $f$ with
respect to its first argument with the second argument held constant,
and $\partial f/\partial y$ denotes the derivative of $f$ with respect
to its second argument with the first argument held constant.  Due to
\eqref{e.epoint2}, we obtain an approximate equation for $\delta_n$,
called the \emph{linearization of \eqref{e.nl-order2} about the
equilibrium point} $x^*$,
\begin{equation}
  \delta_{n+2} 
  = \frac{\partial f(x^*, x^*)}{\partial x} \, \delta_{n+1} 
  + \frac{\partial f(x^*, x^*)}{\partial y} \, \delta_{n} \,.
\end{equation}
This equation is a linear, second order, constant coefficient,
homogeneous difference equation.  Thus, the growth of the perturbation
will depend on the absolute value of the roots of its characteristic
polynomial.  This motivates the following analog of
Theorem~\ref{t.stability}.

\begin{theorem}
Let $x^*$ be an equilibrium point of the difference equation
\eqref{e.nl-order2} with continuously differentiable right hand side
$f$.  Let $\lambda_{1}$ and $\lambda_{2}$ denote the roots of the
characteristic polynomial of its linearization about the equilibrium
point.
\begin{enumerate}
\item If $\lvert \lambda_1 \rvert < 1$ and $\lvert \lambda_2 \rvert
  <1$, then $x^*$ is an asymptotically stable equilibrium point.
\item If $\lvert \lambda_1 \rvert > 1$ or $\lvert \lambda_2 \rvert
  >1$, then $x^*$ is an unstable equilibrium point.
\end{enumerate}
\end{theorem}

As for Theorem~\ref{t.stability}, there are cases where the criterion
is inconclusive, namely when $\lvert \lambda_{i} \rvert \leq 1$ with
identity for at least one of the roots.  

A complete proof is more difficult, but essentially similar to the
proof of Theorem~\ref{t.stability}.


\begin{thebibliography}{9}

\bibitem{E96}
  S.N. Elaydi, ``An introduction to difference equations'', second
  edition, Springer-Verlag, New York, 1996.

\bibitem{I}
  O.A. Ivanov, 
  ``Easy as $\pi$'', Springer-Verlag, 1998.

\bibitem{S00}
  S. Strogatz,
  ``Nonlinear Dynamics and Chaos'',
  Perseus Books, 2000.

\bibitem{W}
  E.W. Weisstein,
  \textit{Logistic Map},
  From MathWorld, Wolfram Research. \\
  \url{http://mathworld.wolfram.com/LogisticMap.html}

\bibitem{Euler}
  ``Euler method.'' Wikipedia, The Free Encyclopedia. 8 Dec 2007,
  17:14 UTC. \\ 
  \url{http://en.wikipedia.org/w/index.php?title=Euler_method&oldid=176593257}
  
\bibitem{logistic}
  ``Logistic map.'' Wikipedia, The Free Encyclopedia. 30 Nov 2006,
  15:45 UTC. \\
  \url{http://en.wikipedia.org/w/index.php?title=Logistic_map&oldid=91166548}
\end{thebibliography}

\end{document}
