\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{Euler:\\
  A basic ODE solver}\\
  \Large{MATH1902: Numerical Solution of Differential Equations}
  \normalsize{http://people.sc.fsu.edu/$\sim$jburkardt/classes/math1902\_2020/euler/euler.pdf}
\end{center}

\hrule

\begin{center}
  \includegraphics[scale=1.50]{euler_sketch.png}
  \vskip 0.1in
  {\it{The Euler method steps off in the right direction, and gradually drifts away.}}
\end{center}
\begin{tcolorbox}[colback=red!5,title=Euler]
{\it{Writing a simple, general Euler code to solve a differential equation will give you an idea
of the kind of mathematical problems we are interested in, and the computational methods we use
to solve them.
}}
\end{tcolorbox}

\section{Introduction}

This is an introduction for all of us.  You would like to know what we are working on.  We would like
to know how well prepared you are for the computational efforts in this directed study class.
\vskip 0.1in
We ask you to read this introductory material, which outlines a method for approximating
the solution of an ordinary differential equation (ODE).  This method can be expressed 
in the MATLAB programming language, and so we ask you to try to write {\it{euler.m}}, 
a file of MATLAB commands that carry out the procedure for a general ODE.  Then we
consider a specific ODE, which describes the shape of a curve by giving a starting value
and the direction of motion at all subsequent times, and we write a MATLAB file
{\it{humps\_deriv.m}} which expresses the direction information.  Now we are ready to solve the 
problem, but we need one more file {\it{euler\_humps.m}}, which specifies the initial value,
tells the euler method which problem to solve, and makes a plot of the results.
\vskip 0.1in
We encourage you to try to follow the directions below that guide you in writing these
programs and using them.  Please let us know if you can't figure out to form a particular
MATLAB expression, or can't understand an error message.  
\vskip 0.1in
We hope that you will be able to write the three programs, run them together, and
get a nice plot of the results.  You will be asked to submit your plot file as
evidence of your work.

\section{A differential equation}

Consider the following differential equation, which can be thought of as defining
the slope of a curve $y=f(t)$ at every point.
\begin{displaymath}
\frac{dy}{dt} = \frac{- 2 ( t - 0.3 ) }{ ( ( t - 0.3 )^2 + 0.01 )^2}
            + \frac{- 2 ( t - 0.9 ) }{ ( ( t - 0.9 )^2 + 0.04 )^2}
\end{displaymath}
If we also know the value of the function at a particular time, then theoretically
we may have enough information to recover the formula for the original function.
However, often this process is not possible.  Nonetheless, we can use this information
to construct an approximate picture of the solution, by estimating the value of the
function at a sequence of points.  This is the meaning of the numerical solution
of a differential equation.
\vskip 0.1in
For this example, we will assume the initial condition $y(0) = 5.1765$, and we will
want to construct a plot of $y(t)$ over the interval $0 \leq t \leq 2$.

\section{{\it{humps\_deriv.m}}: the derivative function}

In MATLAB, a derivative function has a typical structure.  If we follow this structure, then we can
use the same derivative function with many different ODE solvers that MATLAB provides.  Even though
we aren't ready to do that yet, we will follow the pattern, and this will pay off in the future.
\vskip 0.1in
Write a MATLAB function named {\it{humps\_deriv.m}}, with the following format, filling in the details:
\begin{lstlisting}[caption = Outline of the derivative function humps\_deriv.m.]
function dydt = humps_deriv ( t, y )
  dydt = (the formula above);
  return
end
\end{lstlisting}
\vskip 0.1in
In MATLAB, the symbols {\bf{*}}, {\bf{/}}, and {\bf{\^{}}} are used to indicate
multiplication, division, and exponentiation.
You will need to be careful to use parentheses to correctly define your fractions.
In general, you should terminate each MATLAB command with a semicolon; otherwise
MATLAB will print out the value computed by that command, which is usually not
what you want.
If your text line is very long, you can break it into two by using three periods,
as in:
\begin{lstlisting}[caption = You can use a triple period to break a long line.]
  x = 1 + 2 + 3 + 4 + 5 + 6;
  
  x = 1 + 2 + 3 ...
    + 4 + 5 + 6;
\end{lstlisting}

\section{{\it{euler.m}}: approximate ODE solver}

The forward Euler method for solving an ODE is very simple:
\begin{packed_itemize}
  \item{Assume your position is (t,y);}
  \item{Assume you have chosen a stepsize dt;}
  \item{Compute the current slope dy/dt;}
  \item{Update your position to t+dt, y+dy/dt*dt;}
  \item{Repeat as often as desired.}
\end{packed_itemize}
\vskip 0.1in
Now we will write a corresponding MATLAB function, named {\it{euler.m}}.  Note that this
function is ``general'', that is, it is not written for a specific problem, but allows the
user to specify, through input variables, the details of the particular problem to be solved:
\vskip 0.1in
The function will have the rough form:
\begin{lstlisting}[caption = Outline of the general ODE solver euler.m.]
function [ t, y ] = euler ( dydt, tspan, y0, n )
  (initialization)
  for i = 1 : n
    compute next t and y
  end
  return
end
\end{lstlisting}

Here, the input variables are:
\begin{packed_itemize}
  \item{{\tt{dydt}} is the name of a derivative function.}
  \item{{\tt{tspan}} is a vector of length 2.  {\tt{t(1)}} is the start, and {\tt{t(2)}} the final time.}
  \item{{\tt{y0}} contains the value of $y$ at the initial time.}
  \item{{\tt{n}} is the number of equal steps to take from {\tt{t(1)}} to {\tt{t(2)}}.}
\end{packed_itemize}
and the output is:
\begin{packed_itemize}
  \item{{\tt{t}} a vector of length {\tt{n+1}}, containing the initial time, and the {\tt{n}} later times.}
  \item{{\tt{y}} is a vector of length {\tt{n+1}}, containing the initial y0, and the {\tt{n}} later estimates.}
\end{packed_itemize}
\vskip 0.1in

\section{Filling in the details}

The initialization part of the code will look something like this
\begin{lstlisting}[caption = The initialization commands for euler.m]
  t = ?;
  y = ?;
  t(1) = ?;
  y(1) = ?;
  dt = ?;
\end{lstlisting}
Specifically, we need to do the following things:
\begin{packed_itemize}
  \item{create vectors t and y, to store the results.  We can use the MATLAB command
\begin{verbatim}
    vector = zeros ( rows, 1 );
\end{verbatim}
  We are going to take n steps, but we also want to include the initial values, 
  so both vectors should have n+1 rows;}
  \item{Set the first entries of t (contained in tspan(1)) and y (the value y0).}
  \item{Set the time step dt.  What is the formula for dt so that n steps of size
  dt will take us from tspan(1) to tspan(2)?}
\end{packed_itemize}
\vskip 0.1in
Inside the {\tt{for}} loop, we need to evaluate the next values of the variables
t(i+1) and y(i+1).
\begin{lstlisting}
  t(i+1) = t(i) + dt;
  y(i+1) = y(i) + dt * dydt ( t(i), y(i) );
\end{lstlisting}
\vskip 0.1in
Now the euler code should be ready to go.

\section{Use the solver and the derivative functions to make an estimate}

Let's try a quick test run of your codes to see if things are correct.
Don't end this command with a semicolon.  We want to see the value of the output:
You can issue the following command:
\begin{lstlisting}[caption = Try 10 steps from 0 to 2.]
[ t, y ] = euler ( @humps_deriv, [0.0,2.0], 5.1765, 10 )
\end{lstlisting}

\begin{lstlisting}[caption = Output from a 10 step run.]
t =

   0.00000
   0.20000
   0.40000
   0.60000
   0.80000
   1.00000
   1.20000
   1.40000
   1.60000
   1.80000
   2.00000

y =

   5.1765e+00
   1.7675e+01
   1.1867e+02
   2.1050e+01
   1.6150e+01
   2.9192e+01
   1.2072e+01
   4.4357e+00
   1.7620e+00
   5.8524e-01
  -3.0498e-02
\end{lstlisting}
Now, to see a picture of your results, type:
\begin{lstlisting}[caption = Plotting your data.]
plot ( t, y );
\end{lstlisting}

\section{{\it{euler\_humps.m}}: Collecting commands into a script}

Computing is an experimental process.  In a MATLAB interactive session, you may need to type 6 or 10 commands in a row before you get a result.  If you want to see the result again, you have to remember what those commands were, and type them all over again correctly.  If you want to repeat the experiment, but with one parameter changed, you again may have to type in your commands again.  This happens especially when the result of one command depends on the output of previous commands.
\vskip 0.1in
So while interactive experimentation is very helpful in feeling your way towards a solution, I strongly recommend that, once you know what you want to do, you create a {\it{script file}} that lists your commands in order.  That way, you can easily execute the entire experiment with a single command.  If you want to change a small part of the experiment, you use an editor to modify just that one item, and run it again.  If there are always one or two variables that should be specified interactifely, these can simply be input to your script.
\vskip 0.1in
Here's how it would work for the example we just ran.  First, we choose a name for the script file.  Here, I would choose the name {\it{euler\_humps.m}}.  As it turns out, the number of steps to be taken really ought to be chosen at the last minute, so we will make that an input variable:
\begin{lstlisting}[caption = Script version 1.]
function euler_humps ( n )
  return
end
\end{lstlisting}
\vskip 0.1in
Now we stick in the commands that we used before, except that, instead of using the value
10, we now use the input value {\tt{n}}:
\begin{lstlisting}[caption = Script version 2.]
function euler_humps ( n )

  [ t, y ] = euler ( @humps_deriv, [0.0,2.0], 5.1765, n );
  plot ( t, y );
  
  return
end
\end{lstlisting}
At this point, you can run your script with 20 steps.  This should immediately produce
a plot, and it should be a little smoother than the previous one:
\begin{lstlisting}[caption = Run the script.]
euler_humps ( 20 );
\end{lstlisting}
\vskip 0.1in
Now our script includes some quantities that show up only as numbers.  I prefer to use
named variables, that suggest what these quantities mean.  So, even though it's more work,
I would rewrite the script as:
\begin{lstlisting}[caption = Script version 3.]
function euler_humps ( n )

  dydt = @ humps_deriv;
  tstart = 0.0;
  tstop = 2.0;
  tspan = [ tstart, tstop ];
  y0 = 5.1765;
  
  [ t, y ] = euler ( dydt, tspan, y0, n );

  plot ( t, y );
  
  return
end
\end{lstlisting}
and you can see that my call to {\tt{euler()}} uses the same names as were used in the text of the {\tt{euler()}} function.
\vskip 0.1in
Finally, it's possible to improve the look of the plot.  MATLAB plotting commands can be a bit
mysterious, so here I will simply list them, and note that the {\tt{print()}} command actually
allows us to save a copy of the plot as a {\bf{png}} file.
\begin{lstlisting}[caption = Script version 4.]
function euler_humps ( n )

  dydt = @ humps_deriv;
  tstart = 0.0;
  tstop = 2.0;
  tspan = [ tstart, tstop ];
  y0 = 5.1765;
  
  [ t, y ] = euler ( dydt, tspan, y0, n );

  plot ( t, y, 'linewidth', 3 );
  grid ( 'on' );
  xlabel ( '<-- T -->' );
  ylabel ( '<-- Y(T) -->' );
  title ( 'Euler approximation to humps ODE' );
  print ( '-dpng', 'euler_humps.png' );
  
  return
end
\end{lstlisting}
Try this improved version of the script, using {\tt{n=40}} as input.  You should see the plot on the screen, but you should also find a PNG plot file {\it{euler\_humps.png}} in your current directory.  Remember the {\tt{print()}} command, since we will often ask you to verify a homework problem by sending us a PNG plot of the results!

\section{Homework \#1}

If you followed the commands in the previous section, then you have created a PNG plot file
{\it{euler\_humps.png}} for the case {\tt{n=40}}.  As a verification, send this file
as an attachment to {\tt{trenchea@pitt.edu}}.

\end{document}
