\documentclass{article}

\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{float} % an [H] forces a float (such as algorithm) to appear HERE!
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}
 
\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Matlab,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={(*@}{@*)}          % if you want to add a comment within your code
}

\begin{document}

\begin{center}
  \Huge{Implicit ODE Methods\\
  Backward Euler and stiff equations}\\
  \Large{MATH1090: Directed Study in Differential Equations}
  \normalsize{http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1090\_2020/implicit/implicit.pdf}
\end{center}

\hrule

\begin{center}
  \includegraphics[scale=0.25]{switchback.png}
  \vskip 0.1in
  {\it{A zigzag path is a sign of trouble!}}
\end{center}
\begin{tcolorbox}[colback=red!5,title=Implicit methods and stiff equations]
{\it{A stiff differential equation can cause an ODE solver to produce oscillatory errors unless a
very small step size is used.  Using an implicit method instead may kill the oscillations.
}}
\end{tcolorbox}

\section{A surprisingly difficult ODE}

Consider the following differential equation, which we will give the nickname {\it{``stiff''}}:
\begin{align*}
  &y' = 50 ( \cos(t) - y )\\
  &y(0.0) = 0.0
\end{align*}
for which the exact solution is
\begin{displaymath}
  y(t) = 50 \, \frac{\sin(t) + 50 \cos(t) - 50 e^{-50t}}{2501}
\end{displaymath}
and suppose we want estimate the solution over the interval $[0,1]$ using the Euler method.
\vskip 0.1in
We could write the solution as a single program, that allows us to choose the number
of steps $n$, and returns the resulting set of $t, y$ values for tabulating or plotting:
\begin{lstlisting}[caption = Forward Euler method for stiff equation.]
function [ t, y ] = stiff_euler ( n )

  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, 1 );

  a = 0.0;
  b = 1.0;
  dt = ( b - a ) / ( n - 1 );

  t(1) = a;
  y(1) = 0.0;
%
%  We can call a function for dydt, or write the formula.
%
  for i = 1 : n
    t(i+1) = t(i) + dt;
    y(i+1) = y(i) + dt * stiff_prime ( t(i), y(i) );
%   y(i+1) = y(i) + dt * 50.0 * ( cos ( t(i) ) - y(i) );
    
  end

  return
end
\end{lstlisting}
For convenience, we might write the derivative as a function:
\begin{lstlisting}[caption = stiff\_prime.m evaluates the right hand side.]
function dydt = stiff_prime ( t, y )
  dydt = 50.0 * ( cos ( t ) - y );
  return
end
\end{lstlisting}
and the exact solution:
\begin{lstlisting}[caption = stiff\_exact.m evaluates the exact solution.]
function value = stiff_exact ( t )
  value = 50.0 * ( sin ( t ) + 50.0 * cos(t) - 50.0 * exp ( - 50.0 * t ) ) / 2501.0;
  return
end
\end{lstlisting}
\vskip 0.1in
A value of $n=27$ might seem reasonable to get from $t=0.0$ to $t=1.0$.  We can compute the solution,
and plot it versus the true solution (on a finer grid) as follows:
\begin{lstlisting}[caption = stiff\_euler\_test.m computes a forward Euler solution.]
function stiff_euler_test ( n )
  [ t1, y1 ] = stiff_euler ( n );

  t2 = linspace ( 0.0, 1.0, 101 );
  y2 = stiff_exact ( t2 );

  plot ( t1, y1, 'ro-', t2, y2, 'b-' )
  return
end
\end{lstlisting}
The resulting approximate zigzag curve (red) looks nothing like the true smooth solution:
\begin{center}
  \includegraphics[scale=0.60]{stiff_euler.png}
  \vskip 0.1in
  {\it{Explicit Euler solution with 27 steps.}}
\end{center}
It doesn't just seem like inaccuracy, but a wild variation in the direction.  This equation
is just one example of a stiff differential equation.  Such equations include a feature that
gives a typical ODE solver make a staggering, zigzag approximation to the solution.  Even if
we don't know the true solution, we can usually recognize that something is wrong.  In general,
if we reduce the stepsize {\it{dt}} or equivalently increase $n$, eventually the solution
will look reasonable.  However, this stepsize reduction may be quite severe, even though the
shape of the solution curve doesn't give us any warning about how difficult it will be for
\vskip 0.1in
The ODE solvers we have looked at so far have all been of the {\bf{explicit}} or {\bf{forward}} 
type.  That is, when we write the formula for the next solution estimate {\tt{y(i+1}}, we put all the 
ODE derivative information on the right hand side, evaluated at the previous time {\tt{t(i)}} 
and previous solution {\tt{y(i)}}.  And that's what the Euler method can be thought of as
drawing the tangent line at the current solution, and following it forward in time a small distance.
\vskip 0.1in
A different approach to solving ODE's is known as the family of {\bf{implicit}} or {\bf{backward}} methods.  To make an
implicit version of the Euler method, we start out by writing the Euler update equation again,
except that we evaluate the right hand side of the ODE at the ``future'' step $i+1$.  In other
words, for our stiff example, we replace
\begin{lstlisting}[caption = Explicit Euler update step.]
  y(i+1) = y(i) + dt * 50.0 * cos ( t(i) ) - dt * 50 * y(i)
\end{lstlisting}
by
\begin{lstlisting}[caption = Implicit Euler update step.]
y(i+1) = y(i) + dt * 50.0 * cos ( t(i+1) ) - dt * 50 * y(i+1)
\end{lstlisting}
Once we rewrite the update step, we have to rearrange it so that it tells us
how to solve for {\tt{y(i+1)}}.  For this case, it's not too hard to work out what
the implicit Euler update state should be:
\begin{verbatim}
y(i+1)                      =   y(i) + dt * 50.0 * cos ( t(i+1) ) - dt * 50 * y(i+1) 
y(i+1) + dt * 50   * y(i+1) =   y(i) + dt * 50.0 * cos ( t(i+1) )
(  1.0 + dt * 50 ) * y(i+1) =   y(i) + dt * 50.0 * cos ( t(i+1) )
y(i+1)                      = ( y(i) + dt * 50.0 * cos ( t(i+1) ) ) / ( 1.0 + 50.0 * dt )  
\end{verbatim}

\begin{lstlisting}[caption = stiff\_euler\_backward.m computes a backward Euler solution.]
function stiff_euler_backward_test ( n )
  [ t1, y1 ] = stiff_euler_backward ( n );

  t2 = linspace ( 0.0, 1.0, 101 );
  y2 = stiff_exact ( t2 );

  plot ( t1, y1, 'ro-', t2, y2, 'b-' )
  return
end
\end{lstlisting}
Using the same number of steps, our backward Euler solver does a much better job of 
using the ODE to approximate the true curve:
\begin{center}
  \includegraphics[scale=0.60]{stiff_euler_backward.png}
  \vskip 0.1in
  {\it{Implicit Euler solution with 27 steps.}}
\end{center}
Just as there are many explicit ODE solvers, there are many implicit ODE solvers.  The backward 
Euler method has only first order accuracy, so if we think an ODE is stiff, and we want high
accuracy, we might look for an implicit method of higher order.

\section{The Midpoint Method}

The midpoint method for solving an ODE is an implicit method.  When describing programs,
we write {\tt{t(i)}} and {\tt{y(i)}}, but mathematically, we write $t_i$ and $y_i$, and
we think of the derivative or ODE right hand side as a function $y'(t,y)$.  Using this
notation, we can mathematically describe the midpoint method as seeking values that
satisfy the following equation:
\begin{displaymath}
  y_{i+1} = y_i + dt * \frac{y'(t_i,y_i) + y'(t_{i+1},y_{i+1})}{2}
\end{displaymath}
In a way, this looks like an explicit Euler method, except that instead of using
$y'(t_i,y_i)$ in the update, we use an average of the derivative at the current
and future points.  Of course, since we don't actually know the future point,
this is really an implicit method.  If we are lucky, for a particular problem
we can rewrite this update formula so that it is simple to see how to determine
$y_{i+1}$.
\vskip 0.1in
For our stiff example, we start with the definition, and then solve for $y_{i+1}$:
\begin{align*}
  y_{i+1} &= y_i + dt * \frac{y'(t_i,y_i) + y'(t_{i+1},y_{i+1})}{2}\\
  y_{i+1} &= y_i + dt * \frac{50*(\cos(t_i)-y_i) + 50*(\cos(t_{i+1})-y_{i+1})}{2}\\
  y_{i+1} + \frac{dt}{2} (50*y_{i+1}) &= y_i + dt * \frac{50*(\cos(t_i)-y_i) + 50*\cos(t_{i+1})}{2}\\
  ( 1.0 + 25*dt ) * y_{i+1} &= y_i + dt * \frac{50*(\cos(t_i)-y_i) + 50*\cos(t_{i+1})}{2}\\
  y_{i+1} &= \frac{ y_i + 25 * dt * ( \cos(t_i)-y_i + \cos(t_{i+1})) }{ 1.0 + 25*dt }
\end{align*}
\vskip 0.1in
It's actually tricky to write this update formula correctly in MATLAB, getting all the
parentheses in the right place.  So I will write out the solver routine in full, and
show the update formula spread out over several lines, so that I can force the
open and close parentheses to line up clearly.
\begin{lstlisting}[caption = stiff\_midpoint.m solves the stiff problem.]
function [ t, y ] = stiff_midpoint ( n )

  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, 1 );

  a = 0.0;
  b = 1.0;
  dt = ( b - a ) / n;

  t(1) = a;
  y(1) = 0.0;

  for i = 1 : n
    t(i+1) = t(i) + dt;
    y(i+1) = ... 
      ( ...
        y(i) + 25.0 * dt * ( ...
                             cos ( t(i) ) - y(i) + cos ( t(i+1) ) ...
                           ) ...
      ) ...
      / ( 1.0 + 25.0 * dt );
  end

  return
end
\end{lstlisting}
Now by calling the function {\it{stiff\_midpoint\_test(n)}}, we can compute a new estimate and plot it:
\begin{center}
  \includegraphics[scale=0.60]{stiff_midpoint.png}
  \vskip 0.1in
  {\it{Midpoint solution with 27 steps.}}
\end{center}
We can guess from the plot that the midpoint method has better accuracy than the
backward Euler method.  In fact, while the backward Euler method
has accuracy $O(dt)$, the midpoint method's accuracy is $O(dt^2)$.  This means that, for
a given stepsize, we are likely to get better results with the midpoint rule, and if we
cut the stepsize in half, the backward Euler error is divided by 2, but the midpoint error
is divided by 4.  
The midpoint method has some other features that make it very attractive for research computing,
and so we will be discussing it more in the future.

\section{Report: Convergence Comparison}

Since we know the exact solution of the {\it{stiff}} problem, we can determine the error we make
at each step.  We summarize the error using a norm.  A good norm to
use is the RMS norm, which is just the $\ell^2$ norm adjusted for the vector size.  
\vskip 0.1in
We want to verify that the forward Euler and backward Euler methods have an accuracy that is 
$O(dt)$, whereas the midpoint method's accuracy is $O(dt^2)$.   To do this, we will solve the {\it{stiff}}
problem repeatedly, increasing the number of steps from $n=20$ to 40, 80, 160, 320.  Each
time we will compute the RMS norm of the error, and make a table.  Next week we will talk
about how to analyze the results of the table.
\vskip 0.1in
You can modify the codes {\it{stiff\_euler\_test.m}}, {\it{stiff\_euler\_backward\_test.m}}, and 
{\it{stiff\_midpoint\_test.m}} to solve the equation, returning values {\tt{t1}} and {\tt{y1}}.  
Evaluate the exact solution {\tt{y2}} at the same set of {\tt{t1}} values, and then compute {\tt{e}}, the RMS norm
of {\tt{y2-y1}}.  For each value of {\tt{n}} and each ODE solver, record the value of {\tt{e}}.
\begin{verbatim}
                    Backward
 n     Euler          Euler         Midpoint
---  ----------     ----------     ----------
 20  ..........     ..........     ..........
 40  ..........     ..........     ..........
 80  ..........     ..........     ..........
160  ..........     ..........     ..........
320  ..........     ..........     ..........
\end{verbatim}
Bring your completed table to our next meeting, at 2:00pm, Thursday, 13 February, room Thackery 624.
\end{document}
