\documentclass{article}

\usepackage{amsmath,amssymb,amsthm,graphicx,html}
\numberwithin{equation}{section}

\theoremstyle{plain}
\newtheorem{lemma}{Lemma}
\newtheorem{theorem}[lemma]{Theorem}
\newtheorem{corollary}[lemma]{Corollary}

\theoremstyle{remark}
\newtheorem{remark}{Remark}
\newtheorem{example}{Example}
\newtheorem{project}{Project}
\newtheorem{homework}{Homework}

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

\newcommand*  {\const}{\operatorname{const}}

\newcommand*  {\e}{{\mathrm e}}
\renewcommand*{\d}{{\mathrm d}}

\newcommand*  {\A}  {{\boldsymbol A}}
\newcommand*  {\Q}  {{\boldsymbol Q}}
\renewcommand*{\P}  {{\boldsymbol P}}

\newcommand*  {\R} {{\mathbb R}}
\newcommand*  {\Z} {{\mathbb Z}}

\newcommand{\link}[1]{\htmladdnormallink{\texttt{#1}}{#1}}


\begin{document}

\title{Sequential Data Assimilation by Nudging}
\author{Marcel Oliver}
\date{October 2004}

\maketitle

\begin{abstract}
\emph{Data assimilation} is the process of initializing a forecast
(for example a weather forecast) from incomplete observation.  We
consider a system described by system of differential equations where
not all degrees of freedom are observable, either because this is not
practical, as, for example, in the atmosphere where wind speed and
temperature information above ground level is expensive to measure, or
because these degrees of freedom are fundamentally inaccessible.  We
also assume that the system of differential equations is an accurate
model of the real physics (leaving aside issues of model uncertainty),
and that measured quantities are exact (leaving aside measurement
errors).  Thus, the task is to use the differential equation to
determine the complete state at a fixed time from time series data on
the accessible degrees of freedom.

The simplest possible form of this so-called \emph{sequential data
assimilation} is called \emph{nudging} or \emph{dynamical
relaxation}---we introduce additional forces into the differential
equation that ``nudge'' the solution toward the correct internal state
as time evolves.  These lectures will explore nudging in one of the
simplest possible setting for which the technique is interesting: the
three-dimensional Lorenz equations with one accessible observable.
\end{abstract}

\tableofcontents

\noindent
\emph{Instructions for students:} The homework exercises are marked
``Homework'' and are due one week after the last lecture of this
module.  The questions titled ``Project'' are suggestions for the term
paper, or general food for thought.


\section{A crash course in differential equations}

\subsection{Basic setting and examples}

Let $D\subset \R^n$ be the set of admissible states, which we call the
\emph{phase space}.  For a given \emph{vector field} $f \colon D \to
\R^n$ and \emph{initial value} $y_0 \in D$, we seek a continuously
differentiable function $y \colon [0,T) \to D$ which satisfies
\begin{subequations}
  \label{e.ode}
\begin{gather}
  \dot y(t) = f(y(t)) \,, \label{e.ode.a} \\
  y(0) = y_0 \,.
\end{gather}
\end{subequations}
The dot shall denote derivatives with respect to $t$, which is common
practice whenever the physical meaning of $t$ is time.  The interval
of existence of a solution may be unbounded, i.e.\ $T=\infty$, or the
solution solution may cease to exist after a finite interval of time;
see the remarks on blowup in Section~\ref{s.picard} below.  We will
often drop explicit arguments and write $y=y(t)$.  On the other hand,
when the dependence on the initial value is important, we sometimes
write $y(t;y_0)$.

\begin{remark}
The form of equation \eqref{e.ode} seems special in several respects.
However, the formulation is actually very general, as many more
differential equations can be reduced to this particular case.  For
example:
\begin{enumerate}
\item Our vector field $f$ does not explicitly depend on $t$.  Such
equations are called \emph{autonomous}.  Equations with explicit
$t$-dependence are called \emph{non-autonomous}.  By regarding $t$ as
the $(n+1)$st component of the unknown function which satisfies the
equation $\dot t=1$, a non-autonomous equation can always be written
in autonomous form.  However, non-autonomous equations often warrant
special treatment regardless of this formal trick.

\item Our system is first order---no higher order derivatives appear.
Any higher order system can be written as a first order system by
introducing separate dependent variables for the intermediate
derivatives.  This procedure is illustrated in
Example~\ref{s.vanderpol} below.

\item The system is linear and non-degenerate in $\dot y$.  If it is
not, one can sometimes use the implicit function theorem to locally
write the equation as \eqref{e.ode}.  In general, this is difficult,
and further complications arise if this process degenerates.  Such
issues will not be considered here.

\item All data is specified at the initial time; the equation is an
\emph{initial value problem}.  There are also \emph{boundary value
problems} where data is specified at different points in time.  Such
equations often lead to quite different mathematics and will not be
considered here.
\end{enumerate}
\end{remark}

\begin{example} \label{x.exp}
(Exponential growth or decay.) Take $D=\R$;
\begin{subequations}
\begin{gather}
  \dot y = a \, y \,, \\
  y(0) = y_0 \,,
\end{gather}
\end{subequations}
where $a\in\R$ is a given number.  The solution is known from calculus
and can be checked by differentiation,
\begin{equation}
  y(t) = \e^{at} \, y_0 \,.
\end{equation}
\end{example}

\begin{example}
(Linear autonomous system.) Take $D=\R^n$;
\begin{subequations}
\begin{gather}
  \dot y = A \, y \,, \\
  y(0) = y_0 \,,
\end{gather}
\end{subequations}
where $A$ is a given $n\times n$ matrix.  We can formally write
\begin{equation}
  \label{e.linear}
  y(t) = \e^{At} \, y_0 \,.
\end{equation}
It is not difficult to make sense of this expression by defining the
exponential of a matrix through its power series.  This will be
studied in more detail in a third year lecture course on differential
equations or dynamical systems.  If the diagonalization or the Jordan
normal form of $A$ is known, its matrix exponential can be computed
explicitly.  In this sense, linear autonomous systems are ``boring''
as the solution procedure reduces to a standard problem of linear
algebra.
\end{example}

Points $\bar y \in D$ where the solution $y(t;\bar y)$ does not change
in time are called \emph{equilibrium points} (sometimes also
\emph{stationary points} or \emph{fixed points} or \emph{critical
points}).  We observe that if $y$ does not change, then $\dot y$ must
be zero.  Thus, according to \eqref{e.ode.a}, equilibrium points are
characterized by $f(\bar y) = 0$.

\begin{example}
(Logistic equation.) Take $D=\R$;
\begin{subequations}
\begin{gather}
  \dot y = y - a \, y^2 \,, \\
  y(0) = y_0 \,,
\end{gather}
\end{subequations}
Solving this equation is a standard exercise in calculus.  However,
much can be said about the behavior of the solution without actually
finding $y(t)$ explicitly.

The equilibrium points are clearly $\bar y_1=0$ and $\bar y_2=1/a$.
If $0 < y < 1/a$ then $f(y)>0$ and therefore $y(t)$ must be
increasing.  Outside the interval marked by the two equilibria, $y(t)$
must be decreasing.  This completely characterizes the possible
qualitative behavior of the solutions, as is illustrated in
Figure~\ref{f.logistic}.
\end{example}

\begin{figure}[tb]
\centering
\includegraphics{logistic.pdf}
\caption{Solutions to the logistic equation with $a=1$ and for three
different initial values.  The $y$-axis is the one-dimensional phase
space where the direction of the vector field $f$ is indicated.  In
one dimension, $y(t)$ is necessarily monotonic and the equilibrium
points at $y=0$ and $y=1$ act as topological barriers.}
\label{f.logistic}
\end{figure}

\subsection{Existence and uniqueness; blowup}
\label{s.picard}

Differential equations of the form \eqref{e.ode} can be shown to
possess unique solutions under relatively mild assumptions on the
vector field $f$.  This is the content of the following theorem.

\begin{theorem}[Picard--Lindel\"of] \label{t.picard}
Assume that $f$ is uniformly Lipshitz continuous on $D$, i.e.\ there
exists some $L>0$ such that
\begin{equation}
  \lVert f(x) - f(y) \rVert < L \, \lVert x-y \rVert
\end{equation}
for every $x, y \in D$.  Then there exists a time $T>0$ depending only
on $L$ and $\lVert y_0 \rVert$ such that
\begin{enumerate}
\item there exists a continuously differentiable $y \colon [0,T) \to
  D$ which solves the ordinary differential equation \eqref{e.ode};
\item the solution is unique;
\item $y(t;y_0)$ depends continuously on $y_0$ for every fixed $t \in
  [0,T)$.
\end{enumerate}
\end{theorem}

Here and in the following we use $\lVert v \rVert$ to denote the
Euclidean length of a vector $v \in \R^n$, i.e.
\begin{equation}
  \lVert v \rVert^2 = v_1^2 + \dots + v_n^2 \,.
\end{equation}

The proof of this theorem is based on a fixed point argument and is
thus nonconstructive.  Moreover, the time interval of existence
$[0,T)$ that arises in the proof of the theorem is typically much
shorter than the true maximal interval of existence.  However, a
careful analysis of the proof (not given here) reveals that the time
of existence is uniform in $\lVert y_0 \rVert$.  In other words, so
long as the solution remains in a fixed bounded set, one can
repeatedly apply the argument, thereby extending the interval of
existence.  The only way that the solution can cease to exist is by
``blowing up'', as stated more formally in the following.

\begin{theorem} \label{t.blowup}
Provided that $f$ is defined on $\R^n$ and uniformly Lipshitz
continuous on bounded subsets of $\R^n$, then either $T=\infty$ or
\begin{equation}
  \lVert y(t) \rVert \to \infty \quad \text{as } t \to T \,.
\end{equation}
\end{theorem}

\begin{example} \label{x.blowup}
Consider the differential equation
\begin{subequations}
\begin{gather}
  \dot y = y^2 \,, \\
  y(0) = 1 \,,
\end{gather}
\end{subequations}
Formal separation of variables gives
\begin{equation}
  \label{e.separate}
  \frac{\d y}{y^2} = \d t
\end{equation}
(such use of differentials can be made rigorous, or checked by direct
computation on the final answer), and upon integration
\begin{equation}
  \label{e.integrated}
  \int_{y(0)}^{y(t)} \frac{\d y}{y^2} = \int_0^t \d t \,.
\end{equation}
(The necessity to integrate from $y(0)$ to $y(t)$ on the left may seem
odd, but can be derived as follows: write the left hand side of
\eqref{e.separate} initially as $y^{-2} \, \dot y \, \d t$, integrate
from $0$ to $t$ on both sides, and then change variables from $t$ to
$y$ on the left.)  {From} \eqref{e.integrated}, we obtain
\begin{equation}
  - \frac1y \bigg|_{y(0)}^{y(t)} = t \,.
\end{equation}
Solving for $y(t)$ and inserting the initial condition, we arrive at
\begin{equation}
  y(t) = \frac1{1-t}
\end{equation}
which clearly has blow-up time $T=1$.
\end{example}

\begin{example}
Linear differential equations do not blow up, as solution formula
\eqref{e.linear} shows.
\end{example}

In general, we have to prove that equations do not blow up, as (most
of the time) this means that the model is running out of the
physically relevant regime.  This has to be done for each type of
equation separately.  Theorems~\ref{t.picard} and~\ref{t.blowup}
suggest how to go about it: First, use the existence argument to be
sure to have a solution to work with for at least a short interval of
time, and then try to show that this solution remains in a bounded
subset of $\R^n$, contradicting blowup.  The latter usually requires
some tricks, and we'll look at the issue in detail in
Section~\ref{s.lorenz}.

\begin{homework} \label{h.1}
Perform a computation similar to Example~\ref{x.blowup} that shows
that the solution to the equation
\begin{subequations}
\begin{gather}
  \dot y = \sqrt y \,, \\
  y(0) = 0 \,,
\end{gather}
\end{subequations}
is not unique.  (Note that the function $f(y)=\sqrt y$ does not
satisfy a Lipshitz condition at $y=0$, thus Theorem~\ref{t.picard}
does not apply.)
\end{homework}

\begin{project}
Write up the proof of the Picard--Lindel\"of theorem, which can be
found in most textbooks on differential equations.  Is the requirement
of the Lipshitz condition ``sharp'', or can it be weakened?  Explore
examples of non-uniqueness like Homework~\ref{h.1}, and try to
close the gap between examples and theorem.
\end{project}


\subsection{Nonlinear equations in two dimensions}
\label{s.vanderpol}

To sharpen our intuition, we will consider an equation that is only
slightly more complicated than the equations discussed before.  It is
called \emph{van der Pol equation}, and describes oscillations in an
an electrical circuit involving a cathode ray tube.  Usually written
as the scalar second order equation
\begin{equation}
  \ddot x + x = \mu \, (1-x^2) \, \dot x \,,
\end{equation}
with $\mu$ a positive parameter, it can also be written as a
two-component system of first order equations by defining $y=\dot x$,
so that
\begin{subequations}
\begin{align}
  \dot x &= y \,, \\
  \dot y &= -x + \mu \, (1-x^2) \, y \,,
\end{align}
\end{subequations}
augmented with initial conditions on $x$ and $y$.

\begin{figure}[tb]
\centering
\includegraphics{vanderpol.pdf}
\caption{Phase plane of the van der Pol oscillator.  The thick orbit
is the stable limit cycle, the dashed lines are two other orbits that
converge to the limit cycle from the inside and the outside,
respectively.} 
\label{f.vanderpol}
\end{figure}

The van der Pol system illustrates that even in a two dimensional
phase space there are strong topological obstructions on the dynamics:
For large times, an orbit can only converge to an equilibrium point, a
limit cycle, or diverge.  This result, called the
\emph{Poincar\'e--Bendixon Theorem}, should be intuitively clear; see,
for example, Figure~\ref{f.vanderpol}.  A formal proof, however, is
more involved and not the subject of these lectures.


\subsection{The Lorenz equations and chaos}
\label{s.lorenz}

Beyond the topological constraints of dimensions one and two, very
complicated behavior is possible.  In the following, we will study the
\emph{Lorenz equations}
\begin{subequations} \label{e.lorenz}
\begin{align}
  \dot x &= \sigma \, (y-x) \,, \\
  \dot y &= \rho x - y - x z \,, \\
  \dot z &= x y - \beta z \,, \label{e.lorenz.c}
\end{align}
\end{subequations}
where $\sigma$, $\rho$, and $\beta$ are positive parameters.  The
equations represent the dynamics of the fundamental modes of a
two-dimensional model of fluid convection; see \cite{W}, for example.

Part of a single trajectory in phase space is shown in
Figure~\ref{f.attr}.  This type of qualitative behavior is called
\emph{deterministic chaos}.  The word ``deterministic'' indicates that
each trajectory follows a well posed initial value problem and is
therefore at least in principle computable.  ``Chaos'' refers to the
following.

\begin{figure}[!tb]
\centering
\includegraphics{lorenzattr.pdf}
\caption{The Lorenz attractor.  Solution curve of the Lorenz equation
in phase space with $\sigma=10$, $\rho=35$, $\beta=\tfrac83$,
$x(0)=10$, $y(0)=0$, and $z(0)=10$ for the time interval
$t\in[0,30]$.} 
\label{f.attr}
\end{figure}

\begin{enumerate}
\item The state of the system after some time \emph{depends
sensitively on the initial data} to an extent that makes long term
prediction practically impossible.

\item Trajectories approach a \emph{strange attractor} as $t\to
\infty$.  This is a so-called fractal set with a dimension typically
non-integer and larger than one.  We will not go into the mathematics
of these concepts here.

\item \emph{Intermittency}---solutions have long relatively regular
states, interrupted by intermittent periods of rapid change; see
Figure~\ref{f.yseries}.
\end{enumerate}

\begin{figure}[tb]
\centering
\includegraphics{lorenzy.pdf}
\caption{Time series for $y(t)$ for the solution to the Lorenz
  equations corresponding to the computation of Figure~\ref{f.attr}.}
\label{f.yseries}
\end{figure}

\begin{remark}
The historic significance of the Lorenz equations is summarized in the
following quote from \cite{DM}.
  
\begin{quote}
\small
The significance of these equations, which were discovered by Edward
Lorenz back in the 60s, is that relatively simple systems such as
these could exhibit rather complex (specifically, chaotic)
behavior. The chaotic aspect of this system demonstrates that, despite
being given the initial conditions to any arbitrary degree of
accuracy, one cannot predict a sufficiently advanced state of the
system. In other words, the system has built into itself the property
of amplifying small perturbations until they become so significant
they affect the accuracy of the results.

For some years after Lorenz published his paper, mathematicians and
physicists did not believe such behavior was even possible at
all---and even for those who thought otherwise, conventional wisdom
suggested that an extremely complex model would be required.

Lorenz, a meteorologist, made his discovery by observing weather
phe\-no\-me\-na---in particular, convection of fluids (and to a
weatherman, the ``fluid'' he's most interested in is air). He took
various mathematical models of fluid convection and simplified them
into a system of differential equations that basically led him to the
now-famous Lorenz equations.

This stripped-down system doesn't really model any real-world
situations to a degree that would interest an applied scientist, but
Lorenz wasn't concerned with modeling; rather, he wanted to examine
what he felt were the fundamental, intrinsic properties of such
nonlinear systems. He was more interested in the complexity of
behavior than in how this behavior might translate to real-world
phenomena.

That said, it took years for his work to receive the attention it so
well deserved. Many colleagues and peers initially either (1) refused
to consider his work because, being interested in how mathematics
models the world around them, they gave no thought to abstract models
that ``don't behave anything like real weather''; or (2) did not
believe such behavior was indeed exhibited by these equations, and
attributed the results to computational errors.

Well, nowadays every mathematician and scientist understands and
accepts the idea first put forth by Lorenz, that chaotic behavior can
arise naturally from simple nonlinear systems, and by extension,
real-world phenomena that are modeled by similar nonlinear systems.
\end{quote}
\end{remark}

In the remainder of this section, we will prove that solutions to the
Lorenz equations do not blow up.  In fact, we will show something
stronger, namely that solutions to the Lorenz equations enter a set
which is bounded independent of the magnitude of the initial data
after a finite transient time (which can well depend on the initial
condition).

Writing $Y=(x,y,z)$, we first note that we need to find a bound on the
quantity $\lVert Y \rVert^2 = x^2 + y^2 + z^2$.  The plan is to find a
(scalar!)\ differential equation for this quantity.  However, even for
simple examples it is usually impossible to achieve exact closure;
usually we are forced apply some inequalities to close this equation.
In this case, we will eventually even replace the norm by a different
dynamic quantity to obtain the strongest possible result.

We begin by multiplying the components of the Lorenz system
\eqref{e.lorenz} by $x$, $y$, and $z$, respectively, thereby obtaining
\begin{subequations}
  \label{e.components2}
\begin{align}
  \frac12 \, \frac\d{\d t} x^2
    &= \sigma \, x y - \sigma \, x^2 \,, \\ 
  \frac12 \, \frac\d{\d t} y^2
    &= \rho \, x y - y^2 - x y z \,, \\
  \frac12 \, \frac\d{\d t} z^2
    &= x y z - \beta \, z^2 \,. 
\end{align}
\end{subequations}
Adding the three equations, we obtain
\begin{equation}
  \frac12 \, \frac\d{\d t} \lVert Y \rVert^2
  = - \sigma \, x^2 - y^2 - \beta \, z^2 + (\sigma+\rho) \, x y
  \,. \label{e.1}
\end{equation}
To proceed, note the following.

\begin{lemma} \label{l.youngs}
For arbitrary real numbers $a$ and $b$,
\begin{equation}
  \label{e.youngs}
  \lvert a  b \rvert \leq \tfrac12 \, a^2 + \tfrac12 \, b^2 \,.
\end{equation}
\end{lemma}

\begin{proof}
The binomial formula gives
\begin{equation}
  0 \leq (a \pm b)^2 = a^2 \pm 2 a b + b^2 \,.
\end{equation}
Solving for $\mp ab$ yields \eqref{e.youngs}.
\end{proof}

We use this lemma with $a=(\sigma + \rho)\, x$ and $b = y$ to estimate
the cross term on the right of \eqref{e.1},
\begin{equation}
   \lvert (\sigma+\rho) \, x y \rvert
   \leq \tfrac12 (\sigma + \rho)^2 \, x^2 + \tfrac12 \, y^2
\end{equation}
and therefore
\begin{equation}
  \frac12 \, \frac\d{\d t} \lVert Y \rVert^2
  \leq \bigl( \tfrac12 (\sigma + \rho)^2 - \sigma \bigr) \, x^2
  - \tfrac12 \, y^2 - \beta \, z^2 \,.
  \label{e.2}
\end{equation}
Depending on the values of the parameters, two cases are possible.  If
\begin{equation}
  \tfrac12 (\sigma + \rho)^2 - \sigma < 0 \,,
\end{equation}
then there exists a constant $\gamma>0$ (take the smallest in absolute
value of the coefficients on the right side of equation~\ref{e.2})
such that
\begin{equation}
  \frac12 \, \frac\d{\d t} \lVert Y \rVert^2
  \leq - \gamma \, \lVert Y \rVert^2 \,. 
\end{equation}
This differential inequality has the same simple form as the
differential equation in Example~\ref{x.exp}, and therefore yields an
exponentially decreasing bound on the norm of $Y$.  In other words, no
matter where in the phase space we start out, the solution will
converge to zero.  This is obviously not a very interesting choice of
parameters for the Lorenz model.

In the general case, we can still find a constant $c>0$ (take the
largest of the coefficients on the right side of equation~\ref{e.2})
such that
\begin{equation}
  \frac12 \, \frac\d{\d t} \lVert Y \rVert^2
  \leq c \, \lVert Y \rVert^2 \,. 
\end{equation}
This differential inequality implies an exponentially increasing bound
on the norm of $Y$.  In other words, the norm of $Y$ can grow, but not
so fast as to blow up in finite time.

However, the computer generated Figures~\ref{f.attr}
and~\ref{f.yseries} suggest that not only do trajectories not blow up,
but that also they are restricted to a bounded region in phase space.
Looking at the above argument once again, we see that the troublesome
term is the cross-term proportional to $xy$ on the right side of
\eqref{e.1}.  On the other hand, this same term appears on the right
hand side of the $z$-equation of the Lorenz model as well, and we can
add an appropriate lower order term proportional to $z$ to the norm
which we intend to estimate.  More precisely, we set
\begin{equation}
  V = \tfrac12 \rho \, x^2
    + \tfrac12 \sigma \, y^2
    + \tfrac12 \sigma \, z^2 - 2 \rho \sigma z
  \,.
\end{equation}
The choice of prefactors in front of the quadratic terms is just a
convenience, however it is crucial to find the matching coefficient in
front of $z$.

By adding the appropriate multiples of (\ref{e.components2}a--c) and
\eqref{e.lorenz.c}, we directly obtain
\begin{align}
  \frac{\d V}{\d t}
  & = - \rho \sigma \, x^2 - \sigma \, y^2 - \beta \sigma \, z^2
      + 2 \beta \rho \sigma z \,.
      \notag \\
  & = - \bigl(
          \rho \sigma \, x^2
          + \sigma \, y^2
          + \tfrac12 \beta \sigma \, z^2
        \bigr)
      - \tfrac12 \beta \sigma \, z^2
      + 2 \beta \rho \sigma z
      \notag \\
  & \leq -\gamma \, V 
      - \tfrac12 \beta \sigma \, z^2
      + 2 \rho \sigma \,(\gamma + \beta) \, z \,,
      \label{e.3}
\end{align}
where $\gamma$ is chosen small enough such that the last inequality
holds, for example $\gamma = \min \{2 \sigma, 2, \beta \}$.  To find a
usable bound for the last term on the right of \eqref{e.3}, we apply
Lemma~\ref{l.youngs} with $a=2\rho (\gamma+\beta) \sqrt{\sigma/\beta}$
and $b=\sqrt{\beta \sigma} \, z$.  The resulting differential
inequality reads
\begin{equation}
  \label{e.de}
  \frac{\d V}{\d t} \leq c - \gamma \, V
\end{equation}
with
\begin{equation}
  c = \frac{2 \sigma \rho^2 (\gamma+\beta)^2}\beta \,.
\end{equation}
Thus, if $V$ exceeds $c/\gamma$, then the right hand side of
\eqref{e.de} is negative and $V$ is forced to decrease.

Tying together the remaining loose ends in Homework~\ref{h.2}
and~\ref{h.3}, we have altogether proved the following.

\begin{theorem} \label{t.bound}
There exists a bound $B>0$ depending only on the parameters $\sigma$,
$\rho$, and $\beta$ such that for any solution $Y=(x,y,z)$ to the
Lorenz equations \eqref{e.lorenz} there exists a time $\tau \geq 0$,
which may also depend on the initial data, such that
\begin{equation}
  \lVert Y (t) \rVert \leq B 
\end{equation}
for all $t \geq \tau$.  In particular, there is no finite-time
blow-up. 
\end{theorem}

\begin{remark}
This is just one result on the Lorenz equations which will be
important later.  For more information, see modern textbooks on
differential equations, for example \cite{V}.
\end{remark}

\begin{homework} \label{h.2}
Solve the differential inequality \eqref{e.de} explicitly.
\end{homework}

\begin{homework} \label{h.3}
Show that $V(t)$ is bounded if and only if $\lVert Y (t) \rVert$ is
bounded.
\end{homework}

\begin{project}
Repeat the above argument for the van der Pol equation.  Note that the
first part of the argument, i.e.\ the proof that the norm of the
solution grows at most exponentially, is easy.  The hard part is to
show that trajectories enter some bounded set independent of their
initial value.
\end{project}

\section{Nudging}

\subsection{Direct insertion}

We assume that one of the variables, $y(t)$ say, is observable, and
that we have complete exact knowledge of its value for all $t \geq 0$.
We also assume that the dynamics is accurately represented by the
Lorenz model \eqref{e.lorenz}---we are not discussing modeling errors
here.  However, the remaining variables are considered inaccessible;
in particular, we assume no knowledge about their initial values.  Our
task is to use the available information to reconstruct the correct
``internal'' states $x(t)$ and $z(t)$.

A natural and simple approach is to insert the known time series for
$y$ into the equations for $x$ and $z$, then compute the resulting
evolution with arbitrarily chosen initial data.  I.e., we solve
\begin{subequations}
  \label{e.y-nudging}
\begin{align}
  \dot \xi   &= \sigma \, (y-\xi) \,, \\
  \dot \zeta &= \xi y - \beta \zeta \,.
\end{align}
with arbitrarily chosen
\begin{align}
  \xi (0)   &= \xi_0 \,, \\
  \zeta (0) &= \zeta_0 \,.
\end{align}
\end{subequations}
For convenience, we assume that the initial state is already within
the finite bounded region afforded by Theorem~\ref{t.bound}.  This
assumption is not essential, but simplifies the discussion in several
places.

If we initialized \eqref{e.y-nudging} with the exact (but unknown)
initial data $\xi(0)=x_0$ and $\zeta(0)=z_0$, we would replicate our
reference solution exactly.  For other initial data we hope that our
reconstructions $\xi$ and $\zeta$ will at least converge to the
reference solution.  The numerical experiment shown in
Figure~\ref{f.y-nudging} suggests that this is indeed so.

\begin{figure}[tb]
\centering
\includegraphics{y-nudging.pdf}
\caption{Comparison of the $z$-component of the Lorenz solution
(solid curve) to the solution of the ``nudged'' equations
\eqref{e.y-nudging} with $y$ taken as the observable (dashed curve).
The parameters are as in Figure~\ref{f.attr}, and the nudged equation
is initialized with $\xi_0=\zeta_0=0$.}
\label{f.y-nudging}
\end{figure}

We begin our analysis by setting $u=x-\xi$, $v=y-\eta$, and
$w=z-\zeta$.  (The reconstruction $\eta$ of the $y$-component does not
appear in system \eqref{e.y-nudging}, but is convenient for future
use.)  By direct computation, using the Lorenz
equation~\eqref{e.lorenz} and the ``nudged'' system
\eqref{e.y-nudging}, we find
\begin{subequations}
\begin{align}
  \dot u & = \dot x - \dot \xi
           = \sigma \, (y-x) - \sigma \, (y - \xi)
           = - \sigma u \,, \\
  \dot w & = \dot z - \dot \zeta
           = xy - \beta \, z - (\xi y - \beta \zeta)
           = yu - \beta w \,.
           \label{e.w.b}
\end{align}
\end{subequations}
The first equation does not depend on the second and can be solved
immediately,
\begin{equation}
  u(t) = \e^{-\sigma t} \, u_0 \,.
\end{equation}
Thus, $u(t)\to 0$ exponentially as $t \to \infty$, independent of the
initial data.  We insert this result into \eqref{e.w.b}, and write
\begin{equation}
  \dot w + \beta w = y \, \e^{-\sigma t} \, u_0 \,.
\end{equation}
This equation can be solved after multiplication with a so-called
\emph{integrating factor}: We observe that the left hand side is the
result of a product rule after multiplying both sides with $\exp(\beta
t)$.  Hence,
\begin{equation}
  \frac\d{\d t} \bigl( \e^{\beta t} \, w)
  = \e^{\beta t} \, y \, \e^{-\sigma t} \, u_0 \,. 
\end{equation}
Integration in time gives
\begin{equation}
  \e^{\beta s} \, w(s) \bigg|_0^t
  = u_0 \int_0^t \e^{(\beta-\sigma)s} \, y(s) \, \d s \,.
\end{equation}
The integral on the right can obviously not be solved explicitly since
we do not even have an explicit expression for $y$.  However, we only
need an estimate that will prove convergence.  We thus solve for $w$,
take the absolute value, apply the triangle inequality, and note that
Theorem~\ref{t.bound} guarantees the existence of a bound $\lvert y(t)
\rvert \leq B$ for all time:
\begin{align}
  \lvert w (t) \rvert
  & =    \biggl|
           \e^{-\beta t} \, w_0 + u_0 \, \e^{-\beta t}
           \int_0^t \e^{(\beta-\sigma)s} \, y(s) \, \d s
         \biggr|
         \notag \\
  & \leq \e^{-\beta t} \, \lvert w_0 \rvert
         + \lvert u_0 \vert \, \e^{-\beta t}
         \int_0^t \e^{(\beta-\sigma)s} \, \lvert y(s) \rvert \, \d s
         \notag \\
  & \leq \e^{-\beta t} \, \lvert w_0 \rvert
         + \lvert u_0 \vert \, B \, \e^{-\beta t}
         \int_0^t \e^{(\beta-\sigma)s} \, \d s
         \notag \\
  & =    \e^{-\beta t} \, \lvert w_0 \rvert
         + \lvert u_0 \vert \, B \,
         \frac{\e^{-\sigma t} - \e^{-\beta t}}{\beta-\sigma} 
         \notag \\
  & \to  0 \text{ as } t \to \infty \,.
\end{align}
(Note that in the last integration we have considered the case when
$\sigma\neq\beta$.  If $\sigma=\beta$, an even simpler argument leads
to the same conclusion.)  We conclude the following.

\begin{theorem}
When applying \emph{nudging by direct insertion} to the Lorenz model
with observable $y$,
\begin{equation}
  \lvert x(t) - \xi(t) \rvert \to 0
  \qquad\text{and}\qquad
  \lvert z(t) - \zeta(t) \rvert \to 0
\end{equation}
as $t \to \infty$ for any choice of initial values $\xi_0$ and
$\zeta_0$. 
\end{theorem}

We now proceed to see if it is possible to use $z$ as the only
observable.  The corresponding nudged equations are
\begin{subequations}
  \label{e.z-nudging}
\begin{align}
  \dot \xi  &= \sigma \, (\eta-\xi) \,, \\
  \dot \eta &= \rho \xi - \eta - \xi z \,,
\end{align}
\end{subequations}
again with unknown initial values $\xi_0$ and $\eta_0$.  By direct
computation,
\begin{subequations}
\begin{align}
  \dot u &= \sigma \, (v-u) \,, \\
  \dot v &= \rho u - v - u z \,,
\end{align}
\end{subequations}
so that multiplication of these equations with $u$ and $v$,
respectively, gives
\begin{subequations}
\begin{align}
  \frac12 \, \frac{\d}{\d t} u^2
  &= \sigma uv - \sigma \, u^2 \,, \\
  \frac12 \, \frac{\d}{\d t} v^2
  &= \rho uv - v^2 - u v z \,.
\end{align}
\end{subequations}
To obtain an estimate on the norm of the difference between the nudged
and the reference solution, we add the two equations:
\begin{equation}
  \frac12 \, \frac{\d}{\d t} \bigl( u^2 + v^2 \bigr)
  = - \sigma \, u^2 - v^2 + (\sigma + \rho - z) \, uv \,.
  \label{e.norm-equality}
\end{equation}
If we want to control the cross term, which is the only term that may
allow $u$ and $v$ to increase, the best we can do without additional
information is to apply Lemma~\ref{l.youngs} with $a=(\sigma + \rho -
z) \, u$ and $b=v$.  We find that
\begin{equation}
  (\sigma + \rho - z) \, uv
  \leq \tfrac12 \, (\sigma + \rho - z)^2 \, u^2 + \tfrac12 \, v^2 \,,
\end{equation}
so that
\begin{equation}
  \frac12 \, \frac{\d}{\d t} \bigl( u^2 + v^2 \bigr)
  \leq \bigl(\tfrac12 \, (\sigma + \rho - z)^2 - \sigma \bigr) \, u^2 
       - \tfrac12 \, v^2 \,.
\end{equation}
We have no control on the sign of the first term (unless the
parameters $\sigma$ and $\rho$ are very small, and we can show that
the corresponding bound on $z$ given by Theorem~\ref{t.bound} is also
small).  Thus, if $u$ and $v$ have the same sign,
\eqref{e.norm-equality} seems to indicate that the difference between
nudged and reference trajectory increases rather than decreases, while
for $u$ and $v$ of opposite sign, the difference decreases as desired.
Note that we encountered the same type of difficulty when we tried to
prove boundedness of trajectories in Section~\ref{s.lorenz}; see
equation \eqref{e.1}.  The modification of the argument we proposed
there would not suffice here as it yields only bounds of finite
magnitude, not convergence to zero.  However, it may be possible to
prove a weaker statement, see Project~\ref{p.4}.

The argument above does not prove, strictly speaking, that we cannot
use $z$ as the only observable.  It demonstrates, however, that our
usual method of proving convergence fails.  And indeed, the numerical
experiment shown in Figure~\ref{f.z-nudging} verifies our suspicion.

\begin{figure}[tb]
\centering
\includegraphics{z-nudging.pdf}
\caption{Comparison of the $y$-component of the Lorenz solution
(solid curve) to the solution of the ``nudged'' equations
\eqref{e.y-nudging} with $z$ taken as the observable (dashed curve).
The parameters are as in Figure~\ref{f.attr}, and the nudged equation
is initialized with $\xi_0=0$ and $\eta_0=10$.}
\label{f.z-nudging}
\end{figure}

\begin{homework}
Prove that nudging by direct insertion with $x$ as the observable also
works.
\end{homework}

\begin{project} \label{p.4}
Figure~\ref{f.z-nudging} suggests that although the two curves do not
converge, they remain within a small distance of each other.  Can you
prove that this distance is significantly smaller than the typical
amplitude of the solution, e.g.\ as measured by the bound $B$ from
Theorem~\ref{t.bound}? 
\end{project}


\subsection{Time-discrete data}
\label{s.discrete}

Often, we have to assimilate data which is only available at discrete
times.  There are two ways to proceed.  First, we could find a
function which interpolates the given data points and perform
continuous nudging as described previously using the interpolant.  In
this case, we would need to analyze the interpolation error, a task
similar to analyzing measurement errors.  We will not proceed further
in this direction here.  Second, we could solve the full system in
between measurements, re-initializing with new data when it becomes
available.  We will study this procedure in detail, again in the
context of the Lorenz equations.

We assume that $x(t)$ is observable at discrete, equidistant times
$t_n=nh$, where $n=0,1,\dots$ and $h$ denotes the distance between the
sampling points.  The functions $y$ and $z$ are assumed to be
non-observable.  We proceed as follows.
\begin{enumerate}
\item When $t=0$, set $\xi_0=x(0)$ and initialize $\eta$ and $\zeta$
arbitrarily.

\item Let $\xi$, $\eta$, and $\zeta$ evolve from time $t$ to time
$t+h$ according to the exact Lorenz equations
\begin{subequations} \label{e.lorenz.greek}
\begin{align}
  \dot \xi &= \sigma \, (\eta-\xi) \,, \\
  \dot \eta &= \rho \xi - \eta - \xi \zeta \,, \\
  \dot \zeta &= \xi \eta - \beta \zeta \,.
\end{align}
\end{subequations}

\item At time $t+h$, keep the values of $\eta$ and $\zeta$, but set
$\xi$ to the new observation $x(t+h)$.  Go to step (ii) and repeat
with new initial time $t+h$.
\end{enumerate}

Again, we ask the question if, or under which conditions, do the
reconstructions $\xi$, $\eta$, and $\zeta$ converge to the correct
internal states.  Numerical experimentation indicates that for small
enough step sizes, convergence is achieved
(Figure~\ref{f.discrete}), but if the step size grows too large,
convergence fails (Figure~\ref{f.discrete-large}).

\begin{figure}[tb]
\centering
\includegraphics{discrete-nudging.pdf}
\caption{Discrete nudging.  Shown are the observable $x$ and internal
state $y$ with solid lines, and their reconstructions $\xi$ and $\eta$
with dashed lines.  $\xi$ is clearly distinguished by its
discontinuity through re-initialization at the beginning of each step.
The parameters are as in Figure~\ref{f.attr}, and the internal states
are initialized with $\eta_0=10$ and $\zeta_0=0$.}
\label{f.discrete}
\end{figure}

\begin{figure}[tb]
\centering
\includegraphics{discrete-nudging-large.pdf}
\caption{Non-convergence for large step size $h$.  The setting is the
same as in Figure~\ref{f.discrete}, but the step size is four times
larger.}
\label{f.discrete-large}
\end{figure}

We will now analyze a single step of this procedure and prove that for
$h$ small enough the difference to the reference solution reduces by
at least a constant factor smaller than one.  As in the previous
section, we compute, using \eqref{e.lorenz} and
\eqref{e.lorenz.greek},
\begin{subequations}
\begin{align}
  \dot u &= \sigma \, (v-u) \,, \\
  \dot v &= \rho u - v - xz + \xi \zeta
          = \rho u - v - xw - u \zeta \,, \\
  \dot w &= xy - \xi \eta - \beta w
          = x v - u \eta - \beta w \,.
\end{align}
\end{subequations}
Multiplying these equations with $u$, $v$, and $w$, respectively, we
find
\begin{subequations}
\begin{align}
  \frac12 \, \frac\d{\d t} u^2
  &= \sigma uv - \sigma \, u^2 \,, \label{e.u2} \\
  \frac12 \, \frac\d{\d t} v^2
  &= \rho uv - v^2 - xvw - u v \zeta \,, \label{e.u2.b}\\
  \frac12 \, \frac\d{\d t} w^2
  &= x v w - u w \eta - \beta \, w^2 \label{e.u2.c}\,.
\end{align}
\end{subequations}
Setting $V^2=v^2+w^2$, we obtain by adding \eqref{e.u2.b} and
\eqref{e.u2.c} that
\begin{equation}
  \frac12 \, \frac\d{\d t} V^2
  = (\rho - \zeta) \, uv - v^2 - u w \eta - \beta \, w^2 \,.
  \label{e.V2}
\end{equation}
We estimate the cross term in \eqref{e.u2} using Lemma~\ref{l.youngs}
with $a=u$ and $b=v$.  Similarly, the first cross term in \eqref{e.V2}
is estimated by Lemma~\ref{l.youngs} with $a=(\rho-\zeta) \, u$ and
$b=v$; for the second cross term in \eqref{e.V2} we take $a=\eta
u/\sqrt{\beta}$ and $b=\sqrt{\beta} \, w$.  Altogether, we obtain the
system of differential inequalities
\begin{subequations}
\begin{align}
  \frac12 \, \frac\d{\d t} u^2
  & \leq \frac12 \, \sigma \, v^2 - \frac12 \, \sigma \, u^2 \,, \\
  \frac12 \, \frac\d{\d t} V^2
  & \leq \frac12 \, (\rho - \zeta)^2 \, u^2 - \frac12 \, v^2
       + \frac12 \, \frac{\eta^2}\beta \, u^2
       - \frac12 \, \beta \, w^2 \,.
\end{align}
\end{subequations}
Theorem~\ref{t.bound} provides upper bounds on $\eta$ and $\zeta$, so
that $(\rho-\zeta)^2 \leq (\rho+B)^2$ and $\eta^2 \leq B^2$.  Setting
\begin{equation}
  c = \max \{ \sigma, (\rho+B)^2 + B^2/\beta \}
\end{equation}
and
\begin{equation}
  \varepsilon = \min \{ 1, \sigma, \beta \} 
\end{equation}
we see that $u$ and $V$ satisfy the system of differential
inequalities
\begin{subequations}
  \label{e.dineq}
\begin{align}
  \frac\d{\d t} u^2
  & \leq c \, V^2 - \varepsilon \, u^2 \,, \\
  \frac\d{\d t} V^2
  & \leq c \, u^2 - \varepsilon \, V^2 \,.
\end{align}
\end{subequations}
For simplicity, we assume $c=1$; the general case can be recovered
through an appropriate rescaling of time.  We will first solve this
system as a differential equation, and then discuss how this solution
provides bounds on the behavior of functions satisfying the
differential inequalities \eqref{e.dineq}.  Setting
\begin{equation}
  U = 
  \begin{pmatrix}
    u^2 \\ V^2
  \end{pmatrix}
  \qquad \text{and} \qquad
  A =
  \begin{pmatrix}
    - \varepsilon & 1 \\ 1 & -\varepsilon
  \end{pmatrix} \,,
\end{equation}
we now have to solve the linear system
\begin{equation}
  \dot U = A \, U \,.
  \label{e.U}
\end{equation}
Note that $A$ has the diagonalization $A=SDS^{-1}$ with
\begin{equation}
  D =
  \begin{pmatrix}
    1 - \varepsilon & 0 \\
    0 & -1 - \varepsilon 
  \end{pmatrix}
  \qquad \text{and} \qquad
  S = \frac1{\sqrt2}
  \begin{pmatrix}
    1 & 1 \\ 1 & -1
  \end{pmatrix} = S^{-1} \,,
\end{equation}
so that
\begin{align}
  \e^{At}
  & = S \, \e^{Dt} \, S^{-1}
    = S \,
      \begin{pmatrix}
        \e^{(1-\varepsilon)t} & 0 \\
        0 & \e^{-(1+\varepsilon)t}
      \end{pmatrix} \, S^{-1}
      \notag \\
  & = \frac12 
      \begin{pmatrix}
        1 & 1 \\ 1 & -1
      \end{pmatrix}
      \begin{pmatrix}
        \e^{(1-\varepsilon)t} & \e^{(1-\varepsilon)t} \\
        \e^{-(1+\varepsilon)t} & -\e^{-(1+\varepsilon)t}
      \end{pmatrix}
      \notag \\
  & = \frac12
      \begin{pmatrix}
        \e^{(1-\varepsilon)t} + \e^{-(1+\varepsilon)t} &
        \e^{(1-\varepsilon)t} - \e^{-(1+\varepsilon)t} \\
        \e^{(1-\varepsilon)t} - \e^{-(1+\varepsilon)t} &
        \e^{(1-\varepsilon)t} + \e^{-(1+\varepsilon)t}
      \end{pmatrix}
\end{align}
Since we are initializing $\xi$ exactly, $u(0)=0$, hence $U(0)=
(0,V_0^2)$, so that the solution to \eqref{e.U} is given by
\begin{equation}
  U(t) = \e^{At} \, U(0)
       = \frac12
      \begin{pmatrix}
        \e^{(1-\varepsilon)t} - \e^{-(1+\varepsilon)t} \\
        \e^{(1-\varepsilon)t} + \e^{-(1+\varepsilon)t}
      \end{pmatrix} \, V_0^2 
\end{equation}
and, in particular,
\begin{equation}
  V^2(h) = \tfrac12 \,
           \bigl(
             \e^{(1-\varepsilon)h} + \e^{-(1+\varepsilon)h}
           \bigr) \, V_0^2 \,.
\end{equation}
It is easy to see that the prefactor on the right is bounded by a
constant $\mu<1$ for $h$ not too large (see Homework~\ref{h.7}), so
that, after $n$ iterations of the argument, our error measure has
decreased to
\begin{equation}
  V(nh) \leq \mu^{n/2} \, V_0 \,.
\end{equation}
This bound converges to zero as $n\to \infty$.

It can be shown that the result of this computation truly provides a
bound on $u$ and $V$ so long as $V$ is decreasing (see
Project~\ref{p.7}), even though we had replaced the inequalities in
\eqref{e.dineq} by equalities.  Thus, we have altogether proved the
following.

\begin{theorem} \label{t.discrete}
There exists a maximum step size $H$ so that nudging the Lorenz model
with discrete data in the $x$ variable at times $t_n=nh$ with $h \in
[0,H]$ yields an approximation $(\xi,\eta,\zeta)$ that converges to
$(x,y,z)$ as $n\to\infty$. 
\end{theorem}

\begin{remark}
Our proof gives a sufficient condition for convergence.  Numerical
experiments, on the other hand, indicate that even when $h$ is much
larger than the bounds given above, convergence still occurs
(Figure~\ref{f.errors}).  The zoom shown in Figure~\ref{f.zoom} shows
that this is not due to ``being too generous'' in our estimates when
analyzing a single step: There are clearly steps in which the computed
error does grow.  Convergence in this situation is thus a consequence
of subtle averaging over many steps, and a rigorous analysis requires
completely different techniques.  The case when
Theorem~\ref{t.discrete} truly applies is shown in
Figure~\ref{f.tiny}: the error decreases monotonically from step to
step, although it may increase within a single step.
\end{remark}

\begin{figure}[tb]
\centering
\includegraphics{discrete-nudging-error.pdf}
\caption{Error measures $u^2$ (solid line) and $V^2$ (dashed line) for
the solutions plotted in Figure~\ref{f.discrete}.}
\label{f.errors}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics{discrete-nudging-zoom.pdf}
\caption{Zoom into the first few steps shown in Figure~\ref{f.errors}.
One can clearly see that in each nudging interval the error initially
decreases, but the step size $h$ is too large to guarantee a decrease
in error up to the end of the interval.}
\label{f.zoom}
\end{figure}

\begin{figure}[tb]
\centering
\includegraphics{discrete-nudging-tiny.pdf}
\caption{Discrete nudging with tiny step size.  In this case,
Theorem~\ref{t.discrete} applies and the error measure $V^2$ (dashed
line) decreases across each nudging interval.  The error measure $u^2$
(solid line) is so small that it is magnified by a factor of $200$.}
\label{f.tiny}
\end{figure}

\begin{homework} \label{h.7}
Prove that such $\mu<1$ exists as claimed in the text, and give an
estimate in terms of $\varepsilon$ of the maximum possible value for
$h$.
\end{homework}

\begin{project} \label{p.7}
Give an argument why it was possible to take the solution of
\eqref{e.U} as an upper bound for $u$ and $V$ satisfying the
differential \emph{inequalities} \eqref{e.dineq}.
\end{project}


\subsection{Newtonian relaxation}

In practical applications, nudging is often not done by direct
insertion, but by introducing additional forces into the full system
that try to ``relax'' the system to the desired trajectory.  For
example, when nudging in $x$, we may solve
\begin{subequations} 
\begin{align}
  \dot \xi &= \sigma \, (\eta-\xi) + \frac{x-\xi}\tau \,, \\
  \dot \eta &= \rho \xi - \eta - \xi \zeta \,, \\
  \dot \zeta &= \xi \eta - \beta \zeta \,.
\end{align}
\end{subequations}

\begin{project}
Analyze this procedure following the proofs in the previous two
sections.  By varying $\tau$, one can trade off speed of convergence
versus robustness with regard to noise in the data.  Explore this
claim by computer experiment and/or analysis.
\end{project}

\section{Comments}

Nudging is arguably the simplest form of sequential data assimilation.
It is easy to implement, fast, and, depending on the particular
situation, relatively robust and transparent.  However, it has also
some distinct disadvantages.  First, nudging needs relatively large
amounts of data.  It cannot distribute the reconstruction error
uniformly across the time interval of interest, and is therefore
suboptimal in the usage of the supplied data.  Second, it relies
fundamentally on dissipativity, i.e.\ the nonconservation of energy,
in the equation.  However, one sometimes would like to assimilate data
for systems that are not dissipative, or are only very slightly so.

An alternative to nudging is \emph{variational data assimilation}, in
which some measure for the error, for example the function $U$ of
Section~\ref{s.discrete}, is actively minimized.  Such a procedure is
naturally more expensive and more difficult to implement.  On the
other hand, variational data assimilation performs much better when
data is sparse and the equations are not very dissipative.  It also
involves interesting mathematics, but this story will be told another
time.  A few entry points for further reading, available online, are
\cite{R,S,M}.

\begin{project}
Expand this critique of nudging, and survey methods of variational
data assimilation.
\end{project}



\begin{thebibliography}{99}

\bibitem{R}
  M.S. Roulston,
  \textit{Data assimilation in chaotic systems}, \\
  \link{http://gfd.whoi.edu/proceedings/1998/PDF/rep\%5Froulston.pdf}

\bibitem{S}
  T.W. Schlatter,
  \textit{Four-dimensional data assimilation}, in ``The Analysis of
  Observations with Applications in Atmospheric Science``, \\
  \link{http://www.asp.ucar.edu/colloquium/1992/notes/}
  
\bibitem{V}
  F. Verhulst,
  ``Nonlinear Differential Equations and Dynamical Systems'',
  Springer-Verlag, 1996.

\bibitem{W}
  Eric W. Weisstein, \textit{Lorenz Attractor},
  MathWorld, Wolfram Research, \\
  \link{http://mathworld.wolfram.com/LorenzAttractor.html}

\bibitem{DM}
  \textit{Lorenz equations},
  The Math Forum, 
  Drexel University, \\
  \link{http://mathforum.org/library/drmath/view/61791.html}

\bibitem{M}
  \textit{Understanding data assimilation: how models create their
  initial conditions},
  University Corporation for Atmospheric Research, \\
  \link{http://meted.ucar.edu/nwp/pcu1/ic6/}

\end{thebibliography}

\end{document}

