\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}

\DeclareMathOperator*{\argmin}{arg\,min}

\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{RK:\\
  Runge-Kutta ODE Solvers}\\
  \Large{MATH1902: Numerical Solution of Differential Equations}
  \normalsize{http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1902\_2020/rk/rk.pdf}
\end{center}
\hrule
\vskip 0.1in
\begin{center}
  \includegraphics[scale=1.0]{runge_kutta.png}
  \vskip 0.1in
  {\it{A fourth order Runge Kutta step involves several initial test steps.}}
\end{center}

\begin{tcolorbox}[colback=red!5,title=ODE's by a Runge-Kutta method]
{\it{If the Euler method requires too many steps, we can select a more accurate
solver from the Runge-Kutta family.}}
\end{tcolorbox}

\section{How accurate is the Euler method?}

We are interested in approximately solving an ordinary differential equation 
with an initial condition:
\begin{align*}
y' &= f(t,y) \quad \textit{(a given derivative function we will also call ``dydt()'')}\\
y(t_0) &= y_0 \quad \textit{(a given initial value)}
\end{align*}
and so far we have used the forward Euler method.  Given the pair of values
$(t_0,y_0)$ and a stepsize which we call $dt$ or $h$, the method produces
the approximation $(t_1,y_1)$ for the next time, where 
\begin{align*}
  t_1 &= t_0 + dt;\\
  y_1 &= y_0 + dt * dydt(t_0,y_0);
\end{align*}
Now where did the Euler estimate come from?  Let's suppose the exact solution $y(t)$
of the ODE is ``reasonably smooth'' - in this case, assume it has at least two
continuous derivatives over the time interval we are interested in. 
Then by Taylor's theorem we have
\begin{displaymath}
  y(t_1) = y(t_0) + (t_1 - t_0) \, y'(t_0) + \frac{1}{2} \, (t_1 - t_0)^2 \, y"(\xi)
\end{displaymath}
for some point $ t_0 \leq \xi \leq t_1$.
\vskip 0.1in 
If we compare the Euler method and the Taylor result, we can evaluate the absolute
value of the error in our first Euler step, known as the ``local error'':
\begin{displaymath}
  \text{Euler local error} = | y(t_1) - y_1 | = \frac{1}{2} \, dt^2 \, | y"(\xi) |
\end{displaymath}
This tells us that if we cut the stepsize in half, the error bound on our first
Euler step would be multiplied by a fourth.  This suggests why a smaller step
size helps.  For the Euler method, we say the the local error behaves like 
$dt^2$, which we can symbolize by
\begin{displaymath}
  \text{Euler local error} \propto dt^2
\end{displaymath}
or we say this local error is $O(dt^2)$, pronounced ``order of $dt$ squared''.
\vskip 0.1in
However, we are interested in approximating the solution from $t_0$ all the way to
a fixed final time. say $T$, and then estimating the ``global error'', that is,
an upper bound on the absolute value of the difference between our list of approximate
values and the exact solution at those times.  If we know the local error behavior,
we can actually make a statement about the global error.
\vskip 0.1in
We can see how this might be true in the case of the Euler method.
If we cut the stepsize in half, then we double 
the number of steps we have to take to reach $T$.  Roughly speaking, we take
twice as many steps, each having one quarter the local error, which turns out to give
about half the global error at the final time.  
\vskip 0.1in
{\it{Assuming the exact solution is twice continuously differentiable over the interval
$[t_0,T]$, then the local error for the forward Euler method is $O(dt^2)$ and the
global error is $O(dt)$.}}
\vskip 0.1in
In other words, cutting the stepsize in half will divide the very first error by 4,
but overall, this will only tend to reduce the error everywhere by half.  Unfortunately, 
if we have a difficult ODE to solve, and our first attempt
at solving it with the Euler method gives us a very large global error, then we may have
to reduce the stepsize many times before we see a reasonable approximation.  We saw
this problem already in the predator/prey ODE.
\vskip 0.1in
Thus, if we encounter an ODE for which the Euler method does not seem to be giving
us a good approximation, and we suspect that we would need a very small time step
to proceed, we should be interested in using other methods which have a smaller
local error behavior, because they should get us good results faster.
\vskip 0.1in
As it turns out, the Euler method can be regarded as the first member of 
a family known as Runge-Kutta methods.  When we think of it this way, the 
Euler method can be renamed as {\it{rk1}}, where the ``1'' indicates that it
produces a global error characterized as $O(dt^1)$.  
We will now look at 
the next Runge-Kutta method, known as {\it{rk2}}, for which the global error 
is $O(dt^2)$.  We expect this method to do better on our predator-prey problem.

\section{{\it{rk2.m}}, an ODE solver with global error $O(dt^2)$}

The Euler method worked by evaluating the ODE right hand side at the starting point,
and then behaving as this slope was a good estimate of the behavior of the curve
over the interval $[t,t+dt]$.  If we want an improved solution at $t+dt$, we need
a better estimate for the slope over the interval.  The Runge-Kutta method 
{\it{rk2}} does this by trying to average the slopes at the two endpoints.
\vskip 0.1in
In order to evaluate the slope at the far endpoint,
we need a value for $y$ there.  We don't have such a value, but we can use the
Euler method to estimate it.  So in other words, the {\it{rk2}} method starts out
by using the Euler method to estimate the next solution.  It uses the estimated
next solution to evaluate the derivative at the far end.  Then it averages the
near and far derivatives to get a slope estimate over the interval.  And it uses
that improved estimate to take the step.
\vskip 0.1in
Here is a sketch of the process.  We assume we have computed a solution estimate
$(t_n,y_n)$, which we think of as the ``near'' data.  We want to make an estimate
$(t_{n+1},y_{n+1})$ at the ``far'' end of our current interval:
\begin{flalign*}
  f_{near} &= y'(t_n,y_n)           \quad \quad \quad\quad \quad \quad \quad \quad {\textit{Slope at near end}}\\
  y_{far} &= y_n + dt \, y'(t_n,y_n) \quad \quad \quad\quad {\textit{Estimate far solution using Euler step}}\\
  f_{far} &= y'(t_{n+1},y_{far})    \quad \quad \quad\quad \quad \quad {\textit{Estimate slope at far end}}\\
  y_{n+1} &= y_n + dt \, (f_{near}+f_{far})/2  \quad {\textit{Forward step using improved slope estimate}}\\
\end{flalign*}
\vskip 0.1in
Here is a MATLAB version of this algorithm:
\begin{lstlisting}[caption = MATLAB code for rk2 method.]
function [ t, y ] = rk2 ( dydt, tspan, y0, n )

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

  t(1,1) = tspan(1);
  y(1,:) = y0(:);
 
  for i = 1 : n
    f1 = dydt ( t(i,1),      y(i,:)           );
    f2 = dydt ( t(i,1) + dt, y(i,:) + dt * f1' );
    t(i+1,1) = t(i,1) + dt;
    y(i+1,:) = y(i,:) + dt * ( f1' + f2' ) / 2.0;
  end

  return
end
\end{lstlisting}
{\bf{Try to create a copy of {\it{rk2.m}} yourself.}}  Make sure you think about what is
going on in each line.  In particular, notice that the solution vector {\tt{y()}}
is written as a 2-dimensional vector.  What is stored along a row, and what is stored
along a column?
\vskip 0.1in
Because the derivative vectors {\tt{f1}} and {\tt{f2}}  are column vectors, we
have to transpose with an apostrophe sign when we want to form expressions with
the row vector information in {\tt{y}}.  MATLAB makes a distinction between row
and column vectors, and the MATLAB ODE interface expects column vectors for
derivatives and row vectors for individual solutions.  Don't worry too much 
about this for now, but you should start becoming aware of how MATLAB vectors
work, and be able to identify the difference between row and column vectors.

\section{{\it{predator\_rk2.m}} uses rk2 to solve the predator ODE}

Because the {\it{rk2.m}} code has the same input and output format as {\it{euler.m}}, we
can rerun the predator/prey ODE with this new solver.  We can make a copy of
{\it{predator\_euler.m}} and call it {\it{predator\_rk2.m}}, and replace the
word {\it{euler}} by {\it{rk2}}.  In particular, our call to the ODE solver becomes:
\begin{lstlisting}[caption = Calling the rk2 ODE solver.]
[ t, y ] = rk2 ( dydt, tspan, y0, n );
\end{lstlisting}
\vskip 0.1in
For the {\it{euler()}} solver, we had to use as many as $n=100,000$ steps to get a reasonable
solution to the predator prey problem.  Because the global error in the {\it{rk2()}} solver 
behaves quadratically, we could hope to get equivalent accuracy with somewhere between
100 and 1,000 steps.  {\it{Can you explain why a method of order 2 might be as accurate
as a method of order 1, while using roughly the square root of the number of steps?}}  
\vskip 0.1in
{\bf{Investigate this claim!}} Use the {\it{rk2}} solver on the predator/prey problem 
with steps $n = 100, 200, 300, ...$ until you get a solution that you believe is 
accurate.  
\vskip 0.1in
To choose a satisfactory value of $n$, consider the evidence of the three types of plots:
\begin{packed_enumerate}
  \item{the plots of y(t) and y'(t) should each have repeated peaks of almost the same height;}
  \item{the phase plot of y(t) versus y'(t) should look like a smooth closed curve;}
  \item{the convervation plot should show very litle change.}
\end{packed_enumerate}
The conservation plot should give you the best report on whether the solution is accurate.
\vskip 0.1in
You should be able to conclude that {\it{rk2()}} got you a good answer with much less work
than the Euler method aka {\it{rk1()}}.

\section{{\it{rk4.m}}, an ODE solver with global error $O(dt^4)$}

If we are interested in even higher accuracy, there is a sequence of Runge-Kutta 
ODE solvers of increasing order.  One of the popular methods is {\it{rk4}}, whose
global error tends to behave like $O(dt^4)$.  We will use this method in a few
examples; the user input/output interface is the same as we have seen for 
{\it{euler/rk1}} and {\it{rk2}}.  On the other hand, the actual MATLAB code is 
complicated and more difficult to explain, so we will simply display it, and remark
that it uses four samples of the right hand side function in order to make a very
accurate approximation of the next solution.
\begin{lstlisting}[caption = The rk4 ODE solver.]
function [ t, y ] = rk4 ( dydt, tspan, y0, n )

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

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

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

  for i = 1 : n

    f1 = dydt ( t(i,1),            y(i,:) );
    f2 = dydt ( t(i,1) + dt / 2.0, y(i,:) + dt * f1' / 2.0 );
    f3 = dydt ( t(i,1) + dt / 2.0, y(i,:) + dt * f2' / 2.0 );
    f4 = dydt ( t(i,1) + dt,       y(i,:) + dt * f3' );

    t(i+1,1) = t(i,1) + dt;
    y(i+1,:) = y(i,:) + dt * ( f1' + 2.0 * f2' + 2.0 * f3' + f4' ) / 6.0;

  end

  return
end
\end{lstlisting}
{\bf{Unless you are very patient and careful, I suggest you copy rk4.m from the web 
site, rather than trying to type it in!}}
\vskip 0.1in
As we saw in the {\it{rk2.m}} source code, the right hand side values are column
vectors, and have to be transposed with an apostrophe before being used in computations
involving the row vectors containing solution values.

\section{{\it{expsin\_solution.m}}, a new test ODE}

So far we have judged the improved accuracy of Runge Kutta methods by looking
at plots.  Now it's time to get seriously numerical.  To do so, we'll consider
a simple test function called {\it{expsin}}, whose formula is
\begin{displaymath}
  y(t) = e^{\sin(t)}
\end{displaymath}
\vskip 0.1in
{\bf{Write a corresponding MATLAB function {\it{expsin\_solution.m}} with the 
structure}}:
\begin{lstlisting}[caption = Pseudocode for expsin function.]
function value = expsin_solution ( t )
  value = ?;
  return
end
\end{lstlisting}
Plot this function over the interval $0 \leq t \leq 10$, and observe that it
has a simple periodic behavior.  

\section{ {\it{expsin\_deriv.m}}, the derivative function}

Determine the formula for $y'(t)$.
\vskip 0.1in
{\bf{Write a corresponding MATLAB function {\it{expsin\_deriv.m}} with the structure:}}
\begin{lstlisting}[caption = Pseudocode for expsin ODE right hand side.]
function value = expsin_deriv ( t, y )
  value = ?;
  return
end
\end{lstlisting}

\section{ {\it{expsin\_errors.m}}, solve the expsin ODE 3 ways}

You will need to write a MATLAB function {\it{expsin\_errors.m}} which solves the 
{\it{expsin}} problem.  You are going to solve it using each of the methods
{\it{rk1}}, {\it{rk2}}, {\it{rk4}}, and for each method you are going to use several values of the
number of steps, {\tt{n=5, 10, 20, 40, 80}}.  It turns out you could do this by writing
a program that does all 15 computations, or a program that does one particular
computation, which you have to modify for each new computation.
\vskip 0.1in
Solve the problem over the interval $0 \leq t \leq 10$, using an initial value 
of {\tt{y0=1.0}}.
\vskip 0.1in
Each time that you compute the approximate solution, you must also use
{\it{expsin\_solution(t)}} to compute the exact solution at the same points,
and use {\it{rms.m}} to compute the RMS norm of the difference between the 
approximate and exact solutions.
\vskip 0.1in
\begin{lstlisting}[caption = An output of expsin\_errors.m for rk1]
function expsin_errors ( n )
  dydt = ?
  t0 = ?
  tstop = ?
  tspan = ?
  y0 = ?
  
  [ t, y1 ] = rk1 ( ? );

  y = ?

  e = rms ( y - y1 );

  fprintf ( 1, '  n = %d, ||y-y1|| = %g\n', n, e );

  return
end
\end{lstlisting}
\vskip 0.1in
{\bf{Write the MATLAB code necessary to solve the {\it{expsin}} ODE with each of
the three Runge-Kutta methods, and the several values of {\tt{n}}.}}
\vskip 0.1in
Thus, you will be filling out the entries in the following table,
where {\tt{y1}}, {\tt{y2}} and {\tt{y4}} are the approximate solutions returns by 
{\it{rk1}}, {\it{rk2}} and {\it{rk4}},
and {\tt{y}} is the exact solution vector.
\begin{verbatim}
         rk1         rk2         rk4
         rms(y-y1)   rms(y-y2)   rms(y-y4)   
n =  5:  .........   .........   .........
n = 10:  .........   .........   .........
n = 20:  .........   .........   .........
n = 40:  .........   .........   .........
n = 80:  .........   .........   .........
\end{verbatim}
I expect that you will observe several things in this table:
\begin{packed_itemize}
  \item{for a given {\tt{n}}, the RMS error decreases reading from
  left to right (increasing Runge-Kutta order);}
  \item{for a given Runge Kutta solver, the RMS error decreases as we read from
  top to bottom (increase {\tt{n}};}
  \item{for a given Runge Kutta solver, doubling {\tt{n}} divides the error roughly
  by 2, 4, or 16;}
\end{packed_itemize}

\section { Homework \#4}

Submit a copy of your completed  table to {\tt{trenchea@pitt.edu}} 
 
\end{document}
