% 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{epsfig}
%\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.
%
\title{\ifnum\TeacherEdition=1 TEACHER VERSION  \else STUDENT VERSION \fi\\Numerical Methods of Various Orders: A Comparison}





\markboth{Numerical Methods of Various Orders}{Numerical Methods of Various Orders} 


\author{Swarn Singh\\
Department of Mathematics\\ 
Sri Venkateswara College, University of Delhi\\
New Delhi  110021 INDIA\\
sswarn2005@yahoo.co.in, ssingh@svc.ac.in}

\makeStitlePDFLaTex

\begin{abstract}
It is  not always possible to solve a differential equation analytically.This material makes an effort to teach the basics of numerical methods for first order differential equations by following  graphical and numerical approaches. Here we also discuss  the order of accuracy of the methods and compare their CPU times.
\end{abstract}


\section*{INTRODUCTION}

\subsection*{BACKGROUND}  %%%%%%%%%%%%%%%%%%%
When it is not possible to find analytic solutions for a differential equation, numerical methods are often used.  When an analytic solution of  a differential equation is possible, but it is time-consuming, then the approximation obtained with a numerical solution could be appropriate. To solve a given differential equation, we replace the derivatives by finite differences, based on the limit definition of the derivative. The methods we offer provide  algebraic equations for implementation on computers. 


\subsection*{METHODS AND IMPLEMENTATION}\label{sec:METHODS and IMPLEMENTATION}
The steps involved in finite difference methods are:
\begin{itemize}
\item  discretize the continuous domain into a grid,
\item approximate the derivatives by finite differences, based on the limit definition of the derivative,
\item substitute the approximations in the differential equation to obtain algebraic equations, and
\item solve the resulting system of finite difference equations.
\end{itemize}
Before  proceeding further, let us first clarify what  we mean by accuracy and order of accuracy. We first consider absolute error:

\centerline{Absolute error = $| {\rm analytic\,\, solution} - {\rm numerical\,\, solution}|$.\\} 

A numerical method is said to be {\it r-th order (or order r)\/} if there is a number $C$ independent of $h$ such that the error $E(h)$ behaves like 
$$
E(h)\approx Ch^r
$$


\noindent Hence, the larger the order of accuracy, the faster the error reduces as $h$  decreases. 
One can also find the order of the method using the formula
$$
r=\frac{\ln(E(h_2)/E(h_1))}{\ln(h_2/h_1)}\, ,
$$
\noindent where $h_1$ and $h_2$ are two small step sizes for numerical solutions.

Consider the initial value problem 
$$
y'=f(x,y),  y(0)=y_0\, .
$$
To approximate the solution of the initial value problem we use the following methods based on finite differences.

\subsection*{1.~Euler's method}
Divide the solution domain $[a,b]$ into $n$ subintervals.
The idea behind Euler's method is to use the tangent line to the solution at the independent variables $x_{i+1} $ as intervals in $[a,b]$.
To find the coordinates of the points in our numerical solution, we use the initial condition as our starting point and generate the rest of the solution by using the iterative formulae,

\begin{eqnarray}\label{eq:EulerMeth}
x_{i+1}&=&x_i+h\\ \nonumber
y_{i+1}&=&y_i+h*f(x_i,y_i)\, .
\end{eqnarray}

\noindent Terminate this process when we have reached the right hand end of the desired interval.

Let us work on a simple numerical method example completely by hand first. 

\begin{equation}\label{eq:first}
y'= x + 2y, \quad y(0) = 0
\end{equation}
 to find the value for the approximate numerical solution for (\ref{eq:first}) at $x = 1$.
 
 Incidentally, the analytic solution of  (\ref{eq:first}) can be found using integrating factors or a computer algebra system to be
 \begin{equation}\label{eq:analytic}
 y(x) = \frac{1}{4} \left(-2 x+e^{2
   x}-1\right)\, ,
 \end{equation}
\noindent and we can compare both analytic and numerical solutions.

We now apply Euler's method to estimate  a solution in $[0,1]$.  Here $f(x, y) = x + 2y$ and the initial condition tells us the values of the coordinates of our starting point:
$x_0 = 0, y_0 = 0$.

We  use  Euler's method formulae, (\ref{eq:EulerMeth}), to generate values for $x_1$ and $y_1$ using steps of size $h = 0.25$.

For the first iteration, when $i=0$ we have
$
x_1 = x_0 + h
$
or
$
x_1 = x_0 + 0.25 =0.25
$
so,
$
y_1 = y_0 + h f(x_0, y_0)
=0
$. Thus,  the numerical solution at $x_1 = 0.25$ goes through the point $(0.25, 0)$. When you get to your estimate for $y(0.25)$ how does it compare to the analytic solution from (\ref{eq:analytic}), i.e.  $y(0.25) = 0.03718$?


Now for the second iteration
$
x_2 = x_1 + h
$
and 
$
x_2 = 0.25+ 0.25 =0.5
$
so,
$
y_2 = y_1+ h f(x_1, y_1)
=0+0.25*(x_1+2*y_1)=0.25*0.25=0.0625
$.
So the numerical solution at $x_2 = 0.5$ goes next through the point  $(0.5, 0.0625)$. This compares with the analytic solution from (\ref{eq:analytic}), i.e.  $y(0.5) = 0.17957$?

We now move on to get the third point in the solution $(x_3, y_3)$ and go further to our final point $(x_4 = 1, y_4 =\,  ?)$. When you get to your estimate, $y_4$,  for $y(1)$ how does it compare to the analytic solution from (\ref{eq:analytic}), i.e.  $y(1) = 1.09726$?
We compare our numerical estimate approach  with the analytic solution graphically  in Figure~\ref{fig:plot} and note the increasing errors in our approximation as $x$ increases.

\begin{figure}[h!]
\begin{center}
\includegraphics[width=9cm]{plot.jpg}
\caption{\label{fig:plot} Comparison of numerical estimate approach (straight line segments with iteration points in large dots) with the analytic solution (solid curved line).}
\end{center}
\end{figure}


\subsubsection*{Activities}

\begin{enumerate}
\item Summarize the results of all the calculations in a tabular form, as in Table 1.
\item Plot the graph of the solution. 
\item Compare the solution obtained with the analytic solution.
\end{enumerate}


\begin{table}

\label{ex5:tab1} 
\begin{center}
\begin{tabular}{|c|p{1.5 cm}|p{1.5 cm}|}\hline
$i$&\hspace*{.35 true cm} $x_i$& \hspace*{.35 true cm} $y_i$\\ \hline
0&&\\ \hline
1&&\\ \hline
2&&\\ \hline
3&&\\ \hline
4&&\\ \hline
\end{tabular}
\caption{Table of Numerical Solution using Euler's method.}

\end{center}
\end {table}

%%%%%%%%%%%
\subsection*{2.~Improved Euler's Second Order method}
Use the formula from Euler's method to obtain a first approximation to the solution $y_{i+1}$ as
$$
y^*_{i+1}=y_i+h*f(x_i,y_i)\, .
$$
Now improve this approximation  once more by applying Euler's method and using the average of $f(x_i,y_i)$ and $f(x_{i+1},y^*_{i+1})$ as
$$
y_{i+1}=y_i+\frac{1}{2}h*[f(x_i,y_i)+f(x_{i+1},y^*_{i+1})],
$$

\subsubsection*{Activities}

\begin{enumerate}
\setcounter{enumi}{3}
\item Solve the example given in the previous section using improved Euler's method and compile the solution in a tabular form.
\item Plot the graph of the solution.
\item Compare the solution obtained with the analytic solution and compute the absolute error between the numerical and analytic solution.
\item Obtain the order of accuracy of the method.
\end{enumerate}


\subsection*{3.~Runge Kutta Third Order method}

The Runge Kutta Third Order (RK3) method is as follows. 

\begin{eqnarray*}
y_{i+1}&=&y_i+\frac{1}{6}(k_1+4*k_2+k_3)\, \quad {\rm where}\\
k_1&=&h*f(x_i,y_i),\\
k_2&=&h*f(x_i+\frac{h}{2},y_i+\frac{k_1}{2}),\,\,\, {\rm and}\\
k_3&=&h*f(x_i+h,y_i-k_1+2*k_2)\, .
\end{eqnarray*}

\subsubsection*{Activities}

\begin{enumerate}
\setcounter{enumi}{7}
\item Solve the example given in the section 1 using RK3 method and compile the solution in a tabular form.
\item Plot the graph of the solution.
\item Compare the solution obtained with the analytic solution and compute the absolute error.
\item Obtain the order of accuracy of the method.
\end{enumerate}

\subsection*{4.~Runge Kutta fourth order method}

The Runge Kutta fourth order (RK4) method is as follows.

\begin{eqnarray*}
y_{i+1}&=&y_i+\frac{1}{6}(k_1+2*k_2+2*k_3+k_4)\, \quad {\rm where}\\
k_1&=&h*f(x_i,y_i),\\
k_2&=&h*f(x_i+\frac{h}{2},y_i+\frac{k_1}{2}),\\
k_3&=&h*f(x_i+\frac{h}{2},y_i+\frac{k_2}{2}),\,\,\, {\rm and} \\
k_4&=&h*f(x_i+h,y_+k_3)\, .
\end{eqnarray*}


\subsubsection*{Activities}

\begin{enumerate}
\setcounter{enumi}{11}
\item  Solve the example using RK4 method.
\item  Plot the table and graph of the solution.
\item  Compare the solution obtained with the analytic solution and compute the absolute error between the numerical and the analytic solution.
\item  Obtain the order of accuracy of RK4 method.
\end{enumerate}

\subsection*{COMPARISON OF THE METHODS }\label{sec:Comparison}

\subsubsection*{Activities}

\begin{enumerate}
\setcounter{enumi}{15}
\item Solve the example given in Section 1 with all the four methods and calculate the absolute errors obtained between the analytic and numerical solutions. Make a single table of these.
\item Compare the accuracy of all the methods and check the order of accuracy. Describe the importance of order of accuracy of the methods.
\item For all the methods plot the graphs of error versus iteration. 
\item  Solve the example by using software packages Matlab, Mathematica, Maple, etc.  or see (\ref{eq:first}).
\item  Calculate CPU time taken by each method to get the desired accuracy and complete  Table 2. 
\item  Compare the accuracy and the CPU time taken for all the methods.
\end{enumerate}

\subsubsection*{Final Activity}

Consider the non-linear IVP
$$
y'=(\sin^2 y+\cos x)y, \quad y(0)=1\, .
$$
Implement the above methods and calculate absolute error and CPU time.


\begin{table}

\label{ex5:tab1} 
\centering
\begin{center}
\begin{tabular}{|c|p{1.5 cm}|p{1.5 cm}|p{1.5 cm}|p{2 cm}|}\hline
$i$&\hspace*{.35 true cm} $x_i$& \hspace*{.35 true cm} $y_i$ &\hspace*{.15 true cm} Error & \hspace*{.05 true cm} CPU Time\\ \hline
0&&&& \\ \hline
1&&&& \\ \hline
2&&&& \\ \hline
3&&&& \\ \hline
4&&&& \\ \hline
\end{tabular}
\caption{Table of Numerical Solution and CPU time.}

\end{center}
\end{table} 





\subsection*{PREREQUISITES}
Before starting this activity, it may be helpful (but not required) for students to  know how to
\begin{itemize}
\item  find analytic solutions of ordinary differential equations and
\item use software packages such as Matlab/Mathematica/Maple to solve first order ODE's
\end{itemize}

\subsection*{OBJECTIVES}
Upon successful completion of this activity, students will be able to
\begin{itemize}
\item  know about the importance of learning higher order methods with respect to fewer number of iterations involved to obtain the desired accuracy and
\item know about the costs involved in applying higher order methods in terms of CPU time.
\end{itemize}



\begin{thebibliography}{98}


\bibitem{Bradie}Bradie B.  2007.  {\it  A Friendly Introduction to Numerical Analysis\/}.  Delhi INDIA: Pearson Education.

\bibitem{Jain}Jain M.K., S.R.K.~Iyengear, and R.K.~Jain.  1996. {\it Numerical Methods for Scientific and Engineering Computation\/}. Delhi INDIA: New Age International Pvt. Ltd.

\end{thebibliography}




\teacher{\section*{NOTES FOR TEACHERS}
This activity is designed for a course on numerical analysis or differential equations. The students may  have done a course in differential equations, but the nature of the material here does not require deep knowledge of differential equations. It can be used as a group assignment. Moreover examples where analytical solutions are unknown can be considered to solve numerically.
}
\end{document}