% SIMIODE-TeacherStudent-Template.tex and simiode.cls for LaTeX2e
% Version 22 November 2016.
% 
% This template is to be used with the  class file "simiode.cls". 
%  
% The latest version of this class and template files 
% can be found at
% 
%       <http://www.SIMIODE.org> in the Author section.
%
% A "Quick Start Guide" appears below.
%
% This file must be processed with version LaTeX2e or higher.
% Appreciation to Mark Frantz, IUPUI; John Thoo, Yuba College; and Eric Sullivan, Carroll College.
%
\documentclass{simiode}
%
% With this version of the tex file and the associated SIMIODE.cls file you can produce BOTH a Student Version and a Teacher Version of your modeling scenario or technique narrative.
%
% Student Version contains the title STUDENT VERSION along with the NAME OF SCENARIO, with no abstract, no keywords, no tags, no author identifiers, and no COMMENTS section. The Student Version contains only the STATEMENT of the modeling scenario only with appropriate figures, data, instructions, and perhaps a bibliography appropriate for student resourcing, but nothing in the bibliography which would "give away" the nature of the modeling scenario, etc.
%
%The Teacher Version contains the material in the Student Version, namely the title NAME OF SCENARIO, now headed with   TEACHER VERSION, but also full author  identification, abstract, key words, tag words, and, following the STATEMENT, a section of COMMENTS which can be discussion of presentation  ideas, "solutions," techniques, etc.
%
% Authors can prepare a full Teacher Version using the notions immediately below  to create a pdf for review, but in review please eliminate all information about the author as we use a double-blind referee system. Such information will be restored when accepted and published in SIMIODE at www.simiode.org.
%
%***********   Toggle the following line to turn the teacher version on (1) or off (0). ***********
%
\def\TeacherEdition{0} % Toggle this line's value.
\newcommand{\teacher}[1]{\ifnum\TeacherEdition=1 #1 \fi}
%
%
% The \teacher{ } command wraps any text that will appear only in the teacher version of
% the document. See this file for examples. This has been done for you in the document below.   
% Rename the pdf after you compile each final version.
%
\usepackage{amsmath,amssymb}
\usepackage{color}
\usepackage{graphicx}
%\usepackage{hyperref}
\usepackage{url}
%
\renewcommand{\baselinestretch}{1.3}
%
% NOTE regarding macros %%
%
% Any macros defined for your paper should be contained in 
% the top matter. Likewise, any environment definitions such 
% as \newtheorem or \newenvironment should also go in the 
% top matter. 
% 
% Do not \input your macros as separate files within your 
% final version because more files make it harder for the 
% editors and users to keep track of and/or modify and 
% append your paper.
%
% Top matter is above and the beginning of the document below.
%
\begin{document}

\newpage
%
% You only need to change the NAME OF SCENARIO in the next line and delete the footnote. 
% See above for instructions on toggling from Teacher Version to Student Version.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Top matter is above and 
% the beginning of the document below.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{\ifnum\TeacherEdition=1 TEACHER VERSION \else STUDENT VERSION \fi\\Introduction to Numerical Methods}


\markboth{Introduction to Numerical Methods}{Introduction to Numerical Methods}






\author{Brian Winkel,
SIMIODE\\
Chardon OH 44024  USA\\
BrianWinkel@simiode.org}

\makeStitlePDFLaTex


\begin{abstract}
We develop elementary approaches to numerically solving first order differential equations with Euler's Method and Improved Euler's Method.  We ask students to develop these geometrically and to compute numeric solutions and compare them to the analytic solutions for a number of different step sizes.  
\end{abstract}


\section*{INTRODUCTION}

\subsection*{Introduction}

We  try to solve differential equations analytically (i.e. obtain a closed form solution) for the most part by attempting to maneuver the equation into a position where we can ``see" an antiderivative emerge and then jump on the opportunity by integrating both sides of the differential equation. This solution strategy, while difficult at times, offers us a nice formula with structure we can examine and tells us where our parameters are involved. For example the analytic solution for the exponential growth equation (with $k >0$),

\begin{equation}\label{eq:straight}
y'(t) = k y(t)\, , \quad y(0) = y_0 \,, 
\end{equation}
\noindent is simply $y(t) = y_0 e^{k\, t}$. We see that $k$ is the exponent and if we make $k$ positive and bigger then $y(t)$ will grow faster. It is all right in front of  us when we can analytically solve. Our emphasis is on modeling and we will delegate analytic solutions, for the most part, to {\it Mathematica\/}'s {\tt DSolve} command or other software commands. Of course, other comparable software packages have their own differential equation solving commands. However, we will find out that even a powerful computer algebra system will not always be able to pull a solution for $y(t)$ out of the proverbial hat. So, it may not be possible to solve the differential equation analytically (or in closed form) or it might be just beyond the software's  integration skills or our own by-hand strategies and methods. Indeed, there exist differential equations on which mathematicians have spent considerable effort proving  that there exist no analytic solutions.  For example consider 

\begin{equation}\label{eq:normal}
\displaystyle{y'(t) = e^{- t^2}\, ,\qquad  y(0) = y_0} \,\, . 
\end{equation}

\noindent Perhaps in your calculus class a wise instructor challenged students to integrate the function $\displaystyle{e^{-t^2}}$ which is necessary to solve (\ref{eq:normal}). Integrating this function is something statisticians often need to do, for if you plot the function $x(t) = \displaystyle{e^{-t^2}}$ you will see the familiar bell-shaped curve which describes many phenomena in the world and is a simplified form of the {\it normal distribution}  used in  probability and statistics courses. However, there are proofs that no anlaytic or closed form solution exists for (\ref{eq:normal}).

\subsection*{What if there is no analytic solution?}
So what do we do in  case we cannot find an analytic solution to a differential equation,
 
\begin{equation}\label{eq:generic}
\displaystyle{y'(t) = f(t, y(t))\, ,\qquad  y(0) = y_0} \,\, , 
\end{equation}
\noindent where $f(t, y(t))$ is perhaps some ``wild'' function of $t$ and $y(t)$?
We revert to how we got the differential equation in the first place, i.e. the derivative $y'(t)$ represents some rate of change,

\begin{equation}\label{eq:RateChange}
\displaystyle{y'(t) = \lim_{h \to 0} \frac{y(t+h) - y(t)}{h}} \, .
\end{equation}

Recall that for very small values of $h$ the ratio $\displaystyle{\frac{y(t+h) - y(t)}{h}}$ is approximately the derivative, $y'(t)$. So perhaps if we used this approximation and kept $h$ small we could render an approximation, maybe even a good approximation to the solution for (\ref{eq:RateChange}).  Let us  do just that by attempting to solve the  {\it difference\/}  equation~(\ref{eq:RateChange})  approximation to the actual {\it differential\/} equation~(\ref{eq:generic}). 

\begin{equation}\label{eq:RateChange}
\displaystyle{\frac{y(t+h) - y(t)}{h} = f(t, y(t))\, , \quad y(0) = y_0} \, .
\end{equation}

We will attempt to develop a scheme  and apply it to a differential equation whose solution we actually know to see just how good our scheme really is and how changing the size of $h$, often referred to as the {\it step size\/}, alters our accuracy.  We consider this simple differential equation:
\begin{equation}\label{eq:simple}
\displaystyle{y'(t) = .05*y(t)\, ,\qquad  y(0) = 100} \, . 
\end{equation}

\noindent One interpretation for this model is growth on a continuously compounding bank account paying annual interest of 5\% with a \$100 initial deposit.  We know the solution to this differential equation, by separating the variables and integrating, is  $y(t) = 100\, e^{0.05 t}$.

\subsection*{Developing an approximation scheme for solving numerically}

We  now turn to developing an approximation scheme for (\ref{eq:RateChange}).  Scientists will often develop algorithms or test their theories in known situations, minimizing the effects of extraneous variables. This is  known as solving a {\it toy problem\/} in mathematics.  Figure~\ref{fig:ModelingOneAEulerMethod} might  give you some ideas.

\begin{figure}[h]
\centering
\includegraphics[width=8cm]{ModelingOneAEulerMethod}

%\vspace*{-1.7 true cm}

\parbox{12.9 true cm}{\caption{\label{fig:ModelingOneAEulerMethod} Possible estimate plot for solution of differential equation~(\ref{eq:straight}). The dashed curve is the analytic solution to (\ref{eq:straight}) while the linear piecewise curve is the graphical summary of  Euler's Method, treated below, with estimated points indicated and connections drawn.}}
\vspace*{-.5 true cm}
\end{figure}

\begin{itemize}
\item [a)] Starting at point $(0,100)$ for \$100 at time $t=0$ begin numerically estimating the solution to the difference equation~(\ref{eq:interest}). $h$ is known as the {\it step size\/} of our iterations.

\begin{equation}\label{eq:interest}
\displaystyle{\frac{y(t+h) - y(t)}{h}  = .05*y(t)\, ,\qquad  y(0) = 100} \,\, . 
\end{equation}

Perhaps the formulation in (\ref{eq:Eulerinterest}) where we solve for $y(t+h)$  will help. If we know a value or an approximation for the function $y(t)$, in this case at $t = 0$, can we
 use (\ref{eq:Eulerinterest}) to estimate $y(t)$ at $t = 5$? Here  $h = 5$. We solve (\ref{eq:interest}) for $y(t+h)$ to find:


\begin{equation}\label{eq:Eulerinterest}
y(t+h)  = y(t) + h*0.05*y(t)\, .
\end{equation}

Can we use our knowledge of slope and line to estimate $y(5)$?  Indeed, this first iteration is akin to moving along the tangent line at $(0,0)$ with slope $y'(0)$ until  the estimated value of $y(x)$ at  $x = 5$ is obtained. We shall call this estimate $\hat{y}(5)$ (pronounced ``y hat at .5") to denote an  estimate. Write out the exact arithmetic necessary  to obtain $\hat{y}(5)$. Compare this to  the known solution in our case, i.e.  $y(5) = 100 e^{.05*5} = \$128.40$. Identify the relevant points on the plot in Figure~\ref{fig:ModelingOneAEulerMethod} for this new found estimate and for the actual value of our solution to the differential equation (\ref{eq:straight}). Write out your thoughts  about the comparison, including your observations of what is in Figure~\ref{fig:ModelingOneAEulerMethod}.

\item [b)] Now we have an estimate, $\hat{y}(5)$, of the solution of  (\ref{eq:interest}) at $t = 5$. Enter your results in Table~\ref{table:Euler}.  Then use this new found value to estimate the solution to  (\ref{eq:interest}) at $t = 10$. Write out explicitly the arithmetic statement that will get you from $\hat{y}(5)$ to $\hat{y}(10)$ using $h = 5$, being sure to identify each term in the situation in terms of $y$, $t$ and $h$. Also identify where the estimate is  on the plot in Figure~\ref{fig:ModelingOneAEulerMethod}.

\item [c)] Use your estimate, $\hat{y}(10)$, to render an estimate to the solution of the differential equation~(\ref{eq:interest}) at $t = 15$, i.e. compute $\hat{y}(15)$. Compare your estimate, using $h = 5$,  to the actual value of the solution of the differential equation.

\item [d)] Now do several more computations to fill out Table~\ref{table:Euler}. Also make your own plot of the actual solution and the estimates you have computed. For the former use {\it Mathematica\/}'s {\tt Plot} command and label the plot {\tt SolutionPlot} and for the latter use {\tt ListPlot} of your data set and label the data plot {\tt DataPlot}. Finally, {\tt Show[SolutionPlot, DataPlot]} to see your estimates in comparison with the analytic solution. Comments are in order. What can you say about what is happening and why it is happening?

 
\begin{table} [h!]

\begin{center}

\begin{tabular}{|c|p{1.8 true cm}|p{1.8 true cm}|p{1.8 true cm}|p{1.8 true cm}|p{1.8 true cm}|p{1.8 true cm}|} \hline
$t$ & \multicolumn{1}{|c|}{0}& \multicolumn{1}{|c|}{5}& \multicolumn{1}{|c|}{10}& \multicolumn{1}{|c|}{15}& \multicolumn{1}{|c|}{20}& \multicolumn{1}{|c|}{25} \\ \hline
$\hat{y}(t)$ &\multicolumn{1}{|c|}{\$100.00} & & & & & \\ \hline
${y}(t)$ &\multicolumn{1}{|c|}{\$100.00} & & & & & \\ \hline
\end{tabular}
\parbox{13 true cm}{\caption{\label{table:Euler} Your estimates for  the  solution of differential equation~(\ref{eq:Eulerinterest}) for several values of $t$.}}
\end{center}

\end{table}


\end{itemize}

\vspace*{-.6 true cm}
We can summarize what we have done with (\ref{eq:RateChange}) in one general formula as shown in (\ref{eq:EulerForm})

\begin{equation}\label{eq:EulerForm}
\displaystyle{\underbrace{\hat{y}(t+h)}_{\rm new} =  \underbrace{\hat{y}(t)}_{\rm old} + \underbrace{h*f(t, \hat{y}(t)) }_{\rm change} }\, .
\end{equation}

We simply find a new estimate, $\hat{y}(t+h)$ by adding the change $f(t, \hat{y}(t))*h$ to the old estimate  $\hat{y}(t)$. This method is an old one and is known as {\it Euler's Method\/}, after Leonhard Euler \cite{EulerWiki2010}, the great (no make that awesome!) Swiss mathematician of the 18$^{\rm th}$ century. 

\begin{itemize}
\item [d)] Here is some {\it Mathematica\/} code that will help you play with this situation.  We use $h = 0.1$ to gain more accuracy and we iterate for 300 times to solve our differential equation~(\ref{eq:Eulerinterest}) over the interval $[0, 300*0.1 = 30]$:


\begin{verbatim}
h = 0.1; 
data = {{0, 100}}; 
Do[AppendTo[data, 
      {data[[i - 1, 1]] + h, data[[i - 1, 2]] + .05 data[[i - 1, 2]]*h}], 
{i, 2, 300}]
\end{verbatim}

An explanation of this sequence of commands in {\it Mathematica\/} is in order.  
\begin{itemize}
\item [$\bullet$] {\tt h = 0.1} is our step size in our estimating process.
\item  [$\bullet$] {\tt data} is where we  store our successive estimates - two coordinates for each data point, 
$(t_i, \hat{y}(t_i))$ for the $i^{\rm th}$ estimate.
\item [$\bullet$] We {\tt Append} our new estimate to our data set.
\item [$\bullet$] {\tt data[[i - 1, 1]]} is the first coordinate of the $(i-1)^{\rm th}$  data point (we start at $i=2$) while {\tt data[[i - 1, 2]]} is the  second coordinate of the $(i-1)^{\rm th}$ data point and {\tt  .05 data[[i - 1, 2]]*h} is  the change at $t_i$.
\item [$\bullet$] We do this for $i$ from 2 to 300  with step size $h = 0.1$, which will get us to $t = 0.1*300 = 30$ years.
\end{itemize}

Of course we can change $h$, but then we need to change the number of iterations  to be sure we go as far along as we really want.

Use this code and your plotting expertise to estimate the value $\hat{y}(30)$. Try values as big as $h = 30$ and as small as $h = 0.001$. Hint: Change  {\tt \{i, 2, 30\}} to {\tt \{i, 2, 30/h\}} to automatically account for just how many estimates are necessary to compute to get to $t=30$ with this given value of $h$. Render a plot of these estimates for a reasonable set of  values of $h$ and then for each $h$  selected plot the analytic solution and the estimate data on the same axes. Comment on what these plots show, say for $h = 0.001$ and $h = 10$.

\item [e)] Summarize what we have just learned in part (d). Try to think beyond this one example to how this method might apply to other equations. Visualize the geometry of the estimating process of the straight line segments used in the estimating and of the error when we arrive at  the time we want, in this case $t = 30$, in terms of the step size $h$ we elect.
\item [f)] We take some time to get a handle on just how good this process is at estimating solutions to differential equations.  Consider the following differential equations with initial conditions:
\begin{enumerate}
\item [i)] $y'(t) = 0.03 y(t)\,, \quad y(0) = 10$.
\item [ii)] $y'(t) = 0.003 y(t)^2\,, \quad y(0) = 2$.
\item [iii)] $y'(t) = 0.05 y(t) + 3 \,, \quad y(0) = 10$.
\item [iv)] $y'(t) = 5 - t +.02 y(t)\, ,  \quad y(0) = 1$.
\end{enumerate}

In each case solve the differential equation analytically (use {\it Mathematica\/}'s {\tt DSolve} command) and compute specifically the value of $y(10)$. Then, using $h =  0.5$, 0.25, 0.2, 0.1, 0.05, and 0.01 as step sizes, compute estimates, $\hat{y}(10)$   using Euler's Method and compare  the  answer you got with the analytic solution; a plot of the analytic solution along with the data points for each increment $h$ would be in order. Just how good an approximation Euler's Method  is can be assessed by examining how good it is at predicting a specific value of our function, say, $\hat{y}(10)$ for $y(10)$.  Fill in  Table~\ref{table:results} and come to some conclusions as to the accuracy in Euler's Method.

\begin{table}[h]
\begin{center}


\begin{tabular}{|p{2 true cm}|p{2 true cm}|p{2 true cm}|p{2 true cm}|} \hline \hline
\multicolumn{1}{|c|}{$h$} & \multicolumn{1}{|c|}{$y(10)$} & \multicolumn{1}{|c|}{$\hat{y}(10)$} & \multicolumn{1}{|c|}{$y(10)-\hat{y}(10)$}\\ \hline\hline
\multicolumn{1}{|c|}{0.50}& & & \\ \hline
\multicolumn{1}{|c|}{0.25}& & & \\ \hline
\multicolumn{1}{|c|}{0.20}& & & \\ \hline
\multicolumn{1}{|c|}{0.10}& & & \\ \hline
\multicolumn{1}{|c|}{0.05}& & & \\ \hline
\multicolumn{1}{|c|}{0.01}& & & \\ \hline
\end{tabular}

\parbox{12.4 true cm}{\caption{\label{table:results} A summary of errors for the numerical solution using Euler's Method for one of the four ((i), (ii), (iii), or (iv)) differential equations for various step sizes $h$ in computing $\hat{y}(10)$.}}
\end{center}

\vspace*{-.5 true cm}

\end{table}

For each differential equation, (i), (ii),  (iii), and  (iv), plot the error between the estimate and the actual value of the solution at $t=10$ vs.~step size $h$. What can we infer from these plots? State this as a mathematical conjecture.
  
\item [g)] Now use this new  method to  approximate the solution of the differential equation (\ref{eq:complex}) over the time interval $t \in [0,4]$:

\begin{equation}\label{eq:complex}
\displaystyle{y'(t) = 0.2\, y(t) \frac{\left(10 +5 \sin(25 t)-y(t)\right)}{100 +5\sin(25 t)}\,, \quad y(0) = 50}\, .
\end{equation}
and plot your solution. Of course, you might try to {\tt DSolve}, but we do not hold out much  hope that an analytic solution for (\ref{eq:complex}) can be found. We might then consider a numerical solution with {\it Mathematica\/}'s {\tt NDSolve}, a very powerful, advance numerical solution technique which gives a comparable function (only as a list of points and interpolations).  Compare your results from Euler's Method with {\it Mathematica\/}'s {\tt NDSolve} solution.  Also, consider what this differential equation might model.

\item [h)] Show that a concise formula for Euler's Method with step size $h$ for approximating the solution of the differential equation $y'(t) = f(t, y(t))$ with initial condition $y(t_0) = y_0$  is given by (\ref{eq:EM}):
\vspace*{-.1 true cm}
\begin{equation}\label{eq:EM}
y(t_{k+1}) = y(t_k) + h* f(t_k, y(t_{k+1})) \, ,
\end{equation}
\noindent where $t_k = t_0 + k*h$, $k = 0,1,2, \ldots, $.

One way to do this is to convince yourself, with pictures perhaps, that 

\vspace*{-.3 true cm}
$$y(t_{1}) = y(t_0) + h* f(t_0, y(t_{0})) \, , $$

\vspace*{-.3 true cm}
\noindent and then that 

\vspace*{-.3 true cm}
$$y(t_{2}) = y(t_1) + h* f(t_1, y(t_{1})) \, , $$

\vspace*{-.2 true cm}
\noindent as starters.
\vspace*{-.2 true cm}
\end{itemize}

What we are doing here is repeatedly taking tangent line predictions at each iteration from $t = t_k$ to $t = t_{k+1}$ and, apparently, getting worse and worse results at each step.

\subsection*{Improving on Euler's Method}

Let us return to  Figure~\ref{fig:ModelingOneAEulerMethod} 
and see if we can improve on our approximation strategy. In Figure~\ref{fig:ModelingOneAEulerMethod} we see  the actual solution to the differential equation~(\ref{eq:simple})  (curved line) 

\begin{equation}\label{eq:simple}
\displaystyle{y'(t) = 0.05*y(t)\, ,\qquad  y(0) = 100} \,\, , 
\end{equation}

\noindent and the stepwise approximations from $t=0$ to $t=5$ and thence to $t=10$.

Notice that when we go from the point $(5, \hat{y}(5))$ to $(10, \hat{y}(10))$ our accuracy in our estimation gets worse. However, if we were to connect the initial point $(0,100)$ to the estimate point $(10,\hat{y}(10))$ to form a line segment and then take the estimate at $t = 5$ to be the value on this line segment at $t=5$ we can see that we would have a better estimate of $y(5)$. Call this new estimate $\overline{y}(5)$. In effect, we are using the estimate $\hat{y}(10)$ to tell our estimation process to get ready to  correct for the sharper turn between $(5, \hat{y}(5))$ and $(10, \hat{y}(10))$.




\begin{figure}[h]
\centering
\includegraphics[width=8cm]{ModelingOneAFirstImprovedEulersMethod}

%\vspace*{-1.7 true cm}

\parbox{11.6 true cm}{\caption{\label{fig:ModelingOneAFirstImprovedEulersMethod} Improved Euler's Method estimate at $t=5$  for the solution of differential equation~(\ref{eq:simple}). The dashed curve is the analytic solution to (\ref{eq:simple}); the linear piecewise curve is the graphical summary of  Euler's Method showing estimated points, $(5, \hat{y}(5))$ and $(10, \hat{y}(10))$,    with  connections drawn and solution $(5, \hat{y}(5))$; and the dotted line is the intermediate step in computing the Improved Euler's Method and solution $(5,\overline{y}(5))$.}}
\vspace*{-.6 true cm}
\end{figure}
 
\begin{itemize}
\item [i)] Use the point slope method for determining the equation of a linear function $y = g(t)$  connecting $(0,100)$ to the estimate point $(10,\hat{y}(10))$ and  then evaluate $g(5)$ as an improved estimate for our differential equation (\ref{eq:simple}) at $t = 5$.  Call this improved estimate, $\overline{y}(t)$. Compare $y(5)$, $\hat{y}(5)$, and $\overline{y}(5)$. Is it really improved?

Here is another, less complicated way to view the situation and arrive at a formula for the Improved Euler's Method.  Let us  average the two slopes (1) $m_1$ the slope of the line segment from $(0, y(0))$ to our point at which we use Euler's Method for a first estimate, $(5, \hat{y}(5))$, and (2) $m_2$, the slope of the line segment from $(5, \hat{y}(5))$ to the second Euler's Method estimate  $(10, \hat{y}(10))$ (since the horizontal step size in each segment is the same, namely $h = 5$).  Call this average slope $M_1$. If we then go from our initial point $(0,y(0)$ toward the point $(10, \hat{y}(10))$  along a line of this new slope, $M_1$, we will arrive at our Improved Euler's Method point of $(5, \overline{y}(5))$.  Here is the formula for this Improved Euler's Method and we ask you to explain each step of the expressions:
\vspace*{-.3 true cm}


\begin{equation}\label{eq:IEM}
\overline{y}(5) = y(0) + h \frac{1}{2}(m_1+m_2) \, , \,\,\,  {\rm here} \,\, h=5\,  .
\end{equation}

\noindent Hint: We suggest you compute $m_1$ and $m_2$ first.

\item [j)] Build a concise formula for this Improved Euler's Method with step size $h$ for approximating the solution of the differential equation $y'(t) = f(t, y(t))$ with initial condition $y(t_0) = y_0$ and use your results to fill in (\ref{eq:IIEM}):

\begin{equation}\label{eq:IIEM}
y(t_{k+1}) = y(t_k) + \rule{3 true cm}{.01 cm} \, ,
\end{equation}
\noindent where $t_k = t_0 + k*h$, $k = 0,1,2, \ldots , $.  Figure~\ref{fig:ModelingOneAGeneralImprovedEulerMethod} might help, as it depicts the general step from estimate $\overline{y}(t_k)$ to the next estimate $\overline{y}(t_{k+1})$.
 


\begin{figure}[h]
\centering
\includegraphics[width=8cm]{ModelingOneAGeneralImprovedEulerMethod}

%\vspace*{-1.7 true cm}

\parbox{11.6 true cm}{\caption{\label{fig:ModelingOneAGeneralImprovedEulerMethod} General diagram to help develop iteration formula for Improved Euler's Method estimate of solution at $(t_{k+1},\overline{y}(t_{k+1})$ from that at $t_k$ for solution of differential equation~(\ref{eq:simple}). Also shown is the Euler's Method estimate for solution at $(t_{k+1},\hat{y}(t_{k+1}))$.}}
\vspace*{-.7 true cm}
\end{figure}

\item [k)] Summarize what you have just learned in part (i) - (j). Try to think beyond this one example to how this method might apply to other differential equations. Visualize the geometry of the estimating process via the straight line segments and  the corrections used in the estimating process, as well as  the error when we change  the selected step size $h$.

 Let us now take some time to get a handle on just how good this process is at estimating solutions to differential equations.  Consider the following differential equations with initial conditions:
\begin{enumerate}
\item [i)] $y'(t) = 0.03 y(t)\,, \quad y(0) = 10$.
\item [ii)] $y'(t) = 0.003 y(t)^2\,, \quad y(0) = 2$.
\item [iii)] $y'(t) = 0.05 y(t) + 3 \,, \quad y(0) = 10$.
\item [iv)] $y'(t) = 5 - t +0.02 y(t)\, ,  \quad y(0) = 1$.
\end{enumerate}

In each case solve the differential equation analytically (use {\it Mathematica\/}'s {\tt DSolve} command) and compute specifically the value of $y(10)$. Then using $h =   0.5$, 0.25, 0.2, 0.1, 0.05, and 0.01 as step sizes, compute estimates, $\overline{y}(t)$   using Improved Euler's Method. Then compare  the  answer you got with the analytic solution; a plot of the analytic solution along with the data points for each increment $h$ would be in order.  How well the  Improved Euler's Method  performs can be assessed by examining how well  it  predicts a specific value of our function, say, $y(10)$.  Fill in Table~\ref{table:ImprovedEuler} and come to some conclusions as to the accuracy in te Improved Euler's Method.

\begin{table}[h]
\begin{center}


\begin{tabular}{|p{2 true cm}|p{2 true cm}|p{2 true cm}|p{2 true cm}|} \hline \hline
\multicolumn{1}{|c|}{$h$} & \multicolumn{1}{|c|}{$y(10)$} & \multicolumn{1}{|c|}{$\hat{y}(10)$} & \multicolumn{1}{|c|}{$y(10)-\hat{y}(10)$}\\ \hline\hline
\multicolumn{1}{|c|}{0.50}& & & \\ \hline
\multicolumn{1}{|c|}{0.25}& & & \\ \hline
\multicolumn{1}{|c|}{0.20}& & & \\ \hline
\multicolumn{1}{|c|}{0.10}& & & \\ \hline
\multicolumn{1}{|c|}{0.05}& & & \\ \hline
\multicolumn{1}{|c|}{0.01}& & & \\ \hline
\end{tabular}


\parbox{12.4 true cm}{\caption{\label{table:ImprovedEuler} A summary of errors for the numerical solution using the Improved Euler's Method for one of the four ((i), (ii), (iii), or (iv)) differential equations for various step sizes $h$ in computing $\hat{y}(10)$.}}
\end{center}

\vspace*{-.5 true cm}

\end{table}

For your differential equations (i), (ii),  (iii), or  (iv) plot the error between the estimate and the actual value of the solution at $t=10$ vs. step size $h$. What can you infer from these plots? State this as a mathematical conjecture.
  
\item [l)] Compare your results and your conclusions about the Euler's Method (Table~\ref{table:results}) and Improved Euler's Method (Table~\ref{table:ImprovedEuler}). Offer a general mathematical summary  and compare Euler's Method and  the Improved Euler's Method.

\item [m)] Now use the Improved Euler's Method to solve the following differential equation (\ref{eq:complex1}):
\begin{equation}\label{eq:complex1}
\displaystyle{y'(t) = .2 y(t) \frac{\left(10 +5 \sin(25 t)-y(t)\right)}{100 +5\sin(25 t)}\,, \quad y(0) = 50}\, .
\end{equation} and plot your solution,  comparing it with the solution you got using Euler's Method in (g). 

\end{itemize}

\subsection*{Advanced Numerical Methods}

There exist  some very sophisticated algorithms for estimating solutions for differential equations. Indeed, some corporations and individuals  make very good profits from producing computer code to approximate solutions to these  differential equations, often with hundreds, sometimes thousands, of differential equations in the systems they seek to analyze.

We could pursue additional methods here. For example the Runge-Kutta approach is quite good.  There also exist varying step size algorithms, which essentially take large steps when the derivative, $y'(t)$, is small and smaller step sizes when the derivative, $y'(t)$, is large. Why would one want to vary step sizes in this way?

{\it Mathematica\/} has a very powerful and accurate numerical solution technique for solving differential equations numerically, which is invoked with the command   {\tt NDSolve} as opposed to {\tt DSolve}, which gives  the analytic solution if one exists -- well, if it exists in the strategies {\it Mathematica\/} employs to integrate!

 {\tt NDSolve} has lots of variations in itself, but its basis is a sophisticated algorithm which renders an {\it Interpolating Function\/}, i.e. a set of points with interpolation between these points. In general, {\tt NDSolve}'s capabilities are so good that once we have picked off the numerical solution to a differential equation from {\tt NDSolve} output, we can do just about anything we want with this interpolating function, e.g., differentiate, integrate, plot, use in other equations, etc. What we cannot see from the output of {\tt NDSolve} is a nice closed form solution with parameters visible (which {\tt DSolve} does offer); {\tt NDSolve} renders the solution as a set of numbers and rules for interpolating, but it is a solution that we can use with confidence. 

\begin{itemize}
\item [n)] 
This activity is to show that {\it Mathematica\/}'s {\tt NDSolve} command is very good in returning an interpolation function. In fact, it is so good that you can differentiate and integrate the resulting interpolating function very successfully.

Compare plots of the two solution strategies, {\tt DSolve} and {\tt NDSolve}, for the differential equation~(\ref{eq:compare}).
\vspace*{-.3 true cm}


\begin{equation}\label{eq:compare}
y'(t) = 12 - 5 \sin(3t) + y(t)\, , \quad y(0) = 3.
\end{equation}
Here is the {\it Mathematica\/} code to solve and pick out the solution of (\ref{eq:compare}) using {\tt DSolve}

\begin{verbatim}
yd[t_ ] = y[t] /. DSolve[{y'[t] == 12 - 5 Sin[3 t] + y[t], 
                y[0] == 3}, y[t], t][[1]]
\end{verbatim}



\noindent Here is the {\it Mathematica\/} code to solve and pick out the solution of (\ref{eq:compare}) using {\tt NDSolve}


\begin{verbatim}
ynd[t_ ] = y[t] /. NDSolve[{y'[t] == 12 - 5 Sin[3 t] 
                 + y[t], y[0] == 3}, y[t], {t, 0, 3}][[1]]
\end{verbatim}
\end{itemize}

\section*{Conclusion}


We now have more differential equation solving capability with an understanding of the difference between {\tt DSolve} and {\tt NDSolve} and some idea when to use which command. More importantly, we have a better understanding of how computers numerically solve differential equations which they cannot solve analytically. 

\begin{thebibliography}{1}

\bibitem{EulerWiki2010} Wikipedia. 2010. {\it Leonhard Euler\/}.  {\tt http://en.wikipedia.org/wiki/Leonhard\_Euler}. Accessed 7 January 2019.

\end{thebibliography}

\teacher{
\section*{NOTES FOR THE TEACHER}

We have  found that giving students some insight into just what a numerical method is and how an analytic and a numerical solution compare is very helpful in giving students a perspective on what is desirable for an adequate solution in various modeling situations. Often, the algebraic form in an analytic solution is helpful in seeing the phenomena and the role of the key parameters while we are trying to model, e.g., when we get to second order differential equations we will see the steady-state solution and the transient solution portions of a solution which will be very helpful in understanding both the form of the solution and the meaning of the solution.

There are many variations one could ask students to perform just in Euler and Improved Euler methods alone, e.g., effects of concavity  of the solution function on under- or overestimating numerical solutions. Moreover, one could push further into Runge-Kutta and other advanced methods, but we have found that on our journey of teaching differential equations through modeling and use of technology that it is enough to stop here in this first pass through such a course to give students a feeling for what a numerical method can offer and its value in modeling and its limitations. Moreover, we seek to get students to see the difference between the analytic solution command and the powerful numerical solution command of our computer algebra system, in this case {\tt DSolve} and {\tt NDSolve}, respectively, in {\it Mathematica\/}.

In practice we hope to have students develop a sense of when they can expect to obtain analytic solutions and when they will only have numerical solutions. In the latter situation we show them the exactness of {\tt NDSolve} and the ability to manipulation the numerical solutions offered by {\it Mathematica\/}, e.g., differentiate and integrate such solutions, and thus their real power to model effectively.
}

\end{document}
