\documentstyle[12pt]{article}
\textwidth 16cm
\textheight 23cm
\topmargin 0cm
\evensidemargin 0cm
\oddsidemargin 0cm
\parindent=0pt
\parskip=10pt
\begin{document}


\begin{center}
\Large \bf
Numerieke methoden in de natuurwetenschappen\\
\large 
\vspace{0.2cm}
18-Nov-2000\\
\end{center}
\vspace{0.4cm}
\noindent 
{\bf Ex.~5:  Solution of initial value differential equations}\\
\vspace{0.4cm}

\noindent 
In this exercise you will solve the same problem with three different 
algorithms. In the final report you should 
comment about the advantages/disadvantages of each method.

Consider a sphere falling from rest in presence of air resistance. 
The rate of change of velocity is given by 
\begin{equation}
m {dv \over dt}= mg - kv^2
\end{equation}
where $m=10^{-2}$~kg is the mass of the object,
$v$ the velocity, $g=9,8$~ms$^{-2}$ is the acceleration due to gravity
and $k=10^{-4}$~kg/m is the drag coefficient. You know what to expect, 
the sphere will be first accelerated by gravity but with increasing velocity
the friction term will eventually lead to a motion with constant velocity 
$v_{lim}$. 
This simple problem has  been chosen for the following reasons:
\begin{itemize}
\item it has an analytical solution for $v(t)$ to check the numerical error. 
\item there are different regimes of the motion (linear, rapidly varying, 
constant) and it will be clear that algorithms with adjustable time step are particularly suited in such situations.
\item one can calculate both velocity and position of the falling body 
by reducing the second order equation of motion to a system of two coupled 
first order differential equations.
\end{itemize}

In writing the program, keep all subroutines for integration in a separate 
module which does use any other module. All quantities needed are passed as 
arguments of the call and as external functions. In this way, you can directly 
use these subroutines also for other problems in the future. Make also use of 
an input file to read initial conditions, values of friction coefficient, time 
range of integration, number of equations to be integrated  and whatever 
constitutes an input variable.  


\begin{description}
\item[a)]  Use the  Euler and Runge Kutta method of order 2 and 4 
to calculate the velocity for $0 < t < 10$~s distributed over a uniform grid
with spacing $\delta t$. In order to find an appropriate value for the time 
step $\delta t$ check the results against the expected behaviour and value of  
$v_{lim}$.  
\end{description}

As already mentioned the behaviour of the 
sphere goes through different regimes and integration over a uniform grid in 
time is obviously not the most efficient choice. The purpose of the following 
two exercises is to apply other algorithms where the time step is chosen
at every step in order to keep the precision constant.

\begin{description}

\item[b)]  
Use the  Runge Kutta method with adaptive step size (step doubling)
to find the velocity of the sphere for $0 < t < 10$~s
with tolerance 10$^{-3}$ and 10$^{-4}$. Compare the solutions and the
steps used in the calculation. Generalize to a system of two coupled equations 
to calculate velocity and position as a function of time.

\item[c)] Calculate the solution also by the Runge Kutta Fehlberg method 
(also called embedded Runge Kutta method) as in point b).
The method is described in Numerical Recipes in Fortran, ch.16.2. An on-line 
version of this book can be found in our website Telescience.
\end{description}
\end{document}


\end{document}