\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{Midpoint and ODE45\\
  Two more ODE Solvers}\\
  \Large{MATH1902: Numerical Solution of Differential Equations}
  \normalsize{http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1902\_2020/midpoint/midpoint.pdf}
\end{center}

\hrule

\begin{center}
  \includegraphics[scale=0.25]{midpoint_method.png}
  \vskip 0.1in
  {\it{The midpoint method approximates the slope in the middle.}}
\end{center}
\begin{tcolorbox}[colback=red!5,title=More ODE Solvers]
{\it{
Our current choices for ODE solvers include the fixed stepsize versions of
the Euler method, Runge Kutta solvers rk1, rk2, rk4, the implicit Euler method,
and some simple adaptive stepsize solvers based on rk1 and rk2.  We introduce
the midpoint method, a more powerful alternative to the backward Euler method,
and ode45(), an ODE solver provided with MATLAB, which has high accuracy,
powerful error estimation and stepsize adaptation.
}}
\end{tcolorbox}

\section{The Midpoint Method ``improves'' Backward Euler}

We introduced the backward Euler method as a technique that could handle stiff
differential equations.  The special characteristic of the method is that it is
implicit, so that each step requires setting up and solving an implicit equation.
This extra work was easy to set up using MATLAB's {\it{fsolve()}} function.
\vskip 0.1in
Thus the backward Euler method was able to get a good approximation to the solution
to our stiff ODE example, while the usual Euler method produced a ridiculously
zigzagging result.   However, both these methods are only first order - that is,
the error decreases like $dt$.  We also saw, in the discussion on conservation,
that for some ODE's, it may be important to use a solver that can conserve a
quantity such as energy or mass.  The backward Euler method can only do this
for a conservation quantity that is a linear function of the variables, whereas
many physical systems have a conservation law that is a quadratic function.
\vskip 0.1in
The midpoint method is an enhanced version of the backward Euler method.
It also involves solving an implicit equation, but it has second order accuracy
(the error now descreases like $dt^2$, and it can capture quadratic conservation
laws.  We have already seen the midpoint method in action in the previous
class, but now we will explain how it is defined, and try some more examples.
\vskip 0.1in
Mathematically, the midpoint method for solving $y'(t)=f(t,y)$ using a
fixed stepsize $dt$ produces a sequence of approximate values, whose $i+1$-th
member should 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}
To see where this came from, rewrite it as
\begin{displaymath}
  \frac{dy}{dt} \approx \frac{y_{i+1} - y_i}{dt} = \frac{y'(t_i,y_i) + y'(t_{i+1},y_{i+1})}{2}
\end{displaymath}
We can interpret this as the statement that $y'(t_{i+1/2})$ should be well approximated
by the average of the derivatives at $t_i$ and $t_{i+1}$.  In fact, this approximation
is of one higher degree of accuracy than we had for the forward or backward Euler methods.
\vskip 0.1in
To figure out the value of $y_{i+1}$, we have to solve this implicit equation,
so in general we will need to use something like MATLAB's {\it{fsolve()}} to crack 
the implicit equation that results.  However, if the ODE is linear, then we can
actually work out the solution for this special case.  Let's start by looking
at such a simpler example.

\section{{\it{stiff\_midpoint\_explicit()}}: The Midpoint Method for a Linear ODE}

Our stiff ODE example has the form:
\begin{displaymath}
  y' = 50 ( \cos(t) - y )
\end{displaymath}
with initial condition:
\begin{displaymath}
  y(0) = 0
\end{displaymath}
\vskip 0.1in
Suppose we wish to apply the midpoint method to this problem, using a stepsize $dt$.
Then we have
\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 = midpoint\_explicit.m solves the stiff problem.]
function [ t, y ] = midpoint_explicit ( n )

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

  t0 = 0.0;
  tstop = 1.0;
  dt = ( tstop - t0 ) / n;

  t(1) = t0;
  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}
This solver is actually part of the file {\it{stiff\_midpoint\_explicit.m}}, which uses 
it to to compute and plot a solution.  We can compare the result from backward Euler result to
our new midpoint calculation: 
\begin{center}
  \includegraphics[scale=0.40]{stiff_backward_euler_explicit.png}
  \includegraphics[scale=0.40]{stiff_midpoint_explicit.png}
  \vskip 0.1in
  {\it{Explicit backward euler and midpoint solutions with $n=25$ 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{ {\it{midpoint()}}: A general midpoint solver using {\it{fsolve()}} }

Except when the right hand side is linear in the variable, most ODE's don't
allow us to convert the midpoint method into an explicit formula.  So, just
as we did for the backward Euler method, we are interested in creating a
MATLAB function that can apply the midpoint method to a general problem, 
while looking ``simple'' to the user (same interface as for {\it{euler}},
{\it{rk2()}} and so on), while taking care of the hard work of solving the
implicit equation by calling the MATLAB function {\it{fsolve()}}.
\vskip 0.1in
\begin{lstlisting}[caption = Signature for the midpoint ODE solver.]
function [ t, y ] = midpoint ( dydt, tspan, y0, n )
\end{lstlisting}
\vskip 0.1in
The midpoint function will look similar to our other ODE solvers, except that 
when we need to approximate the solution at the next time, we have to call a
hidden function which solves the implicit midpoint equation:
\begin{lstlisting}[caption = Text of midpoint.m]
function [ t, y ] = midpoint ( f, tspan, y0, n )

%  Ask fsolve() not to print messages.

  options = optimoptions ( 'fsolve', 'Display', 'off' );

  m = length ( y0 );
  t = zeros ( n + 1, 1 );
  y = zeros ( n + 1, m );

  dt = ( tspan(2) - tspan(1) ) / n;

  t(1,1) = tspan(1);
  y(1,:) = y0(:);

  for i = 1 : n

    to = t(i,1);
    yo = y(i,:);
    
    tp = to + dt; 
    yp = yo + dt * f ( to, yo )';  %  Initial guess is forward Euler estimate.

    yp = fsolve ( @(yp)midpoint_residual(f,to,yo,tp,yp), yp, options );

    t(i+1,1) = tp;
    y(i+1,:) = yp;

  end

  return
end
\end{lstlisting}
\vskip 0.1in
Meanwhile {\it{midpoint\_residual()}} defines a function whose value would be
zero if $yp$ is the correct midpoint prediction.  It is simply this:
\begin{lstlisting}[caption = Code to evaluate the midpoint residual.]
function value = backward_euler_residual ( f, to, yo, tp, yp )
  value = yp - yo - ( tp - to ) * ( f ( t0, y0) + f ( tp, yp ) )' / 2.0;
  return
end
\end{lstlisting}
You don't need to understand why we have to write things this way, but
you should be able to see that if {\it{fsolve()}} finds a value $yp$ for which
{\it{midpoint\_residual()}} returns a value of 0, then we have found our desired
result.

\section{{\it{stiff\_midpoint()}} repeats our computation}

Let's use our general midpoint method to solve the stiff ODE again, and make sure
that the explicit and implicit results are identical for this problem.
\vskip 0.1in
The important thing to note is that, by using fsolve() to handle the implicit equation,
we now have a midpoint ODE solver that can be used in the exact same way as our
familiar codes like {\it{euler()}} and {\it{rk2()}}:
\begin{lstlisting}
function stiff_midpoint ( n )

  dydt = @ stiff_deriv;
  t0 = 0.0;
  tstop = 1.0;
  tspan = [ t0, tstop ];
  y0 = 0.0;
%
%  Same format as if we were calling euler() or rk2() or other of our solvers...
%
  [ t1, y1 ] = midpoint ( dydt, tspan, y0, n );

  return
end
\end{lstlisting}
\begin{center}
  \includegraphics[scale=0.40]{stiff_midpoint_explicit.png}
  \includegraphics[scale=0.40]{stiff_midpoint.png}
  \vskip 0.1in
  {\it{Explicit and implicit midpoint solutions with 25 steps.}}
\end{center}

\section{{\it{quadex\_midpoint()}} handles the QUADEX ODE}

The homework for the class on stiff ODE's asked you to apply the backward Euler method
to the QUADEX ODE.  This is a difficult problem, and we know that small errors
can cause the calculation to go bad.  Even with 2,500 steps, the backward Euler
method can't hang onto the solution for the whole test range.
\vskip 0.1in
The midpoint method is also designed to handle stiff problems, but it has the
extra advantage of greater accuracy.  Using only $n=200$ steps over the same
interval, we get a much more accurate picture of this tricky problem.
\begin{center}
  \includegraphics[scale=0.50]{quadex_backward_euler.png}
  \includegraphics[scale=0.50]{quadex_midpoint.png}
  \vskip 0.1in
  {\it{QUADEX solution by backward Euler (2500 steps) and midpoint (200 steps).}}
\end{center}

\section{The Degree of Accuracy}

We have claimed that the midpoint method has degree 2, or order of accuracy 2,
meaning that the global error should behave like $dt^2$.  It is more clear to 
say that if we double $n$ or equivalently, cut $dt$ in half, the error should
go down by a factor of $1/4$.  The forward and backward Euler methods have
degree 1, meaning halving the stepsize only halves the error.
\vskip 0.1in
To get a feeling for what this looks like, and how to investigate the order of
accuracy of any ODE solver, we will consider the EXP ODE:
\begin{align*}
  y'   &= y;    &\text{ODE}\\
  y(0) &= 1;    &\text{Initial condition}\\
  y(t) &= e^t;  &\text{Exact solution}\\
\end{align*}
We will solve this problem over the interval $0 \leq t \leq 5$, and report the
RMS error in our solution, using $n=10, 20, 40, 80$.  You can see that the error 
goes down by roughly a factor of 2 for the Euler codes, but by a factor of 4
for the midpoint rule:
\begin{verbatim}
Method              n=10      n=20      n=40      n=80
---------------     -----     -----     -----     -----
euler()             33.15     20.21     11.44     5.14
backward_euler()   301.01     52.64     18.16     7.72
midpoint()           6.07      1.27      0.29     0.07
\end{verbatim}
This chart is evidence that our accuracy claim is correct.

\section{The Midpoint Method and Research Questions}

The midpoint method for solving an ODE is of special interest to Professor Trenchea.
Here are a few points of interest.
\begin{packed_itemize}
  \item{People disagree about what method should be called ``the midpoint method''.
  Professor Trenchea very firmly believes in the method we have used here.  Textbooks and
  other researchers have used this same name for methods that look similar, differ
  in small or large details;}
  \item{Because the midpoint method has a good conservation and accuracy properties, 
  it is preferred by researchers in areas such as chemistry, physics, magnetism, and 
  even glacial modeling; Professor Trenchea is trying to demonstrate its advantages
  to other scientists who are using inferior methods;}
  \item{Professor Trenchea is very interested in using an adaptive stepsize scheme
  (more powerful than the ones we looked at) with the midpoint method.  This becomes
  a difficult research problem when the ODE switches rapidly between smooth and
  varying behavior.  You can ask Professor Trenchea about the van der Pol ODE, or
  look it up on Wikipedia!}
\end{packed_itemize}

\begin{center}
  \includegraphics[scale=0.40]{vanderpol_ode45_plot.png}
  \vskip 0.1in
  {\it{An ``easy'' version of the van der Pol ODE.}}
\end{center}

\section{ {\it{ode45()}}: MATLAB's ODE Solver}

The user interface to {\it{ode45()}} is almost the same as 
I have been using for our example solvers like {\it{euler()}} and {\it{rk2()}}.
In particular, calling {\it{ode45()}} for a solution looks almost the same
as we are used to:
\begin{lstlisting}
  [ t, y ] = ode45 ( f, tspan, y0 );
\end{lstlisting}
where {\tt{f}} is the name of the user derivative function, {\tt{tspan}} is
a vector of length 2 containing the initial and final times, and {\tt{y0}} is
an $ny$-vector containing the values of the dependent variables at time {\tt{t0}}.
\vskip 0.1in
Our previous ODE solvers included a final argument {\tt{n}}, the number of steps
to take.  But {\it{ode45()}} does not accept that value; this is because it uses
an adapative stepsize algorithm, and so the number of steps to be taken cannot
be determined in advance.  But by allowing the stepsize to vary, {\tt{ode45()}}
aims to produce a solution that is as accurate as desired, as efficiently as possible.
After the solution is returned, you can determine the number of steps that were
taken by a command like this:
\begin{lstlisting}
  n = length ( t ) - 1;
\end{lstlisting}
Here, we are subtracting 1 because we want $n$ to count the number of steps, and so
we don't wish to include the initial value as a step.
\vskip 0.1in
Another thing to keep in mind is how {\tt{ode45()}} handles a problem involving
a system of $ny$ equations.  First, it is important that the derivative routine
returns a {\it{column}} vector; {\it{ode45()}} will refuse to proceed if a row
vector is returned instead.  Thus, for our pendulum ODE problem, the derivative
function looks like this:   
\begin{lstlisting}
function dydt = pendulum_deriv ( t, y )

  [ g, l, m, t0, y0 ] = pendulum_parameters ( );

  u = y(1);
  v = y(2);

  dudt = v;
  dvdt = - ( g / l ) * u;

  dydt = [ dudt; dvdt ];   %  <-- the semicolon makes dydt a COLUMN vector.

  return
end
\end{lstlisting}
Secondly, note that the output of {\it{ode45()}} will be the $nt+1$-vector {\tt{t}}
and the $(nt+1) \times ny$-array {\tt{y}}.  The pendulum solution values at time $i$ are
{\tt{y(i,1)}} and {\tt{y(i,2)}}, or in general {\tt{y(i,1:ny)}}.  To plot the angular 
deflection of the pendulum problem over time, we want to plot {\tt{y(1:nt+1,1)}} or 
{\tt{y(:,1)}}.
\vskip 0.1in
Be sure you know the shape of your data!
MATLAB can make a mess of things if you give it a row vector when it wants a column vector;
sometimes it prints an error message, and sometimes it thinks you are doing a very fancy
linear algebra computation and computes a result you didn't want!
\vskip 0.1in
The {\it{ode45()}} code adapts the stepsize by estimating the error it has made
in the most recent step.  It does this by comparing the predictions of Runge Kutta
methods of degrees 4 and 5, which is where the ``45'' comes from in the name.
Note that, if you cut the stepsize in half, the error in a solver of degree 5 
is divided by 32.  This high  accuracy means  the solver usually does
not have to take very small steps to get good results.
\vskip 0.1in
The {\it{ode45()}} code uses a sophisticated error measurement.  After computing
the difference of two solution estimates to get an error estimate, it decides whether
the error is small or large using a relative scale based on the size of the solution. 
Our simple error estimation ideas would fail such a test.
\vskip 0.1in
Because {\it{ode45()}} uses high degree Runge Kutta solvers, it is expecting that
the derivative function and solution curve will be smooth, with at least 5 continuous
and bounded derivatives.  Some ODE's do not have such smooth solutions, and will
cause the algorithm to be inefficient or to fail.  MATLAB provides alternatives
that use lower degree solvers, such as {\it{ode23()}} for such cases.
\vskip 0.1in
Because {\it{ode45()}} uses explicit Runge Kutta solvers, it may perform poorly
on stiff problems.  When a problem seems to be stiff, MATLAB provides several
alternatives that can work better, including {\it{ode23s()}}.

\section{Homework \#9}

Set up the pendulum and predator-prey problems using {\it{ode45()}}, using the same
data that we had for the conservation class.  For both problems, count $n$, the number
of steps taken, compute and plot the conserved quantities $h()$, and print the
last value $h(n+1)$.  (For the pendulum problem, $h()$ is computed by the
function {\it{pendulum\_energy()}}.)

Fill in the chart to compare the {\it{ode45()}} results
against {\it{midpoint()}}.  What does this chart suggest about the {\it{ode45()}}
efficiency and conservation?

\begin{verbatim}
Problem      initial    midpoint     midpoint     ode45            ode45
               h(1)      h(n+1)        n          h(n+1)             n
            
pendulum       5.37       5.37        400         ..........     ..........

predator       3.22       3.25        400         ..........     ..........
\end{verbatim}

Submit your completed table via email to {\bf{trenchea@pitt.edu}}

\end{document}
