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


\begin{center}
\Large \bf
Computational Methods and Data Analysis 2002\\
\large 
\end{center}
\vspace{0.4cm}
\noindent 
{\bf Ex.~2:  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 Newton's equation of motion read:
\begin{equation}
m {d^2x \over dt^2}= mg - kv^2
\end{equation}

First we examine only the behaviour of the velocity. In this case 
we have a first-order differential equation for  the 
rate of change of velocity  
\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. If you divide both sides by $m$ 
and call v=y(t) one sees that this is of the form 
\begin{equation}
{dy \over dt}= f(y,t)
\end{equation}
with $f(y,t)=g-k/m v^2$

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 a simple analytical form for $v_{lim}$ obtained by setting 
$dv/dt=0$ to compare to the numerical solution. It has also
an analytical solution for $v(t)$ which can be used to check the numerical 
error. The analytical solution is 
\begin{equation}
v(t)= \sqrt{mg \over k} {exp(2 \sqrt{kg \over m}t -1) \over 
exp(2 \sqrt{kg \over m}t +1)}
\end{equation}

where we have used the initial values $v(t=0)=0$. This solution is easily 
obtained by defining a variable $w=\sqrt{k/(mg)} v$ and by rewriting
$1/(1-w^2)=1/2(1/(1-w)+1/(1+w))$. You can try and check it.

\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.
In this case we call $y_1(t)=v$, $y_2(t)=x$ and we get a system of two 
coupled differental equations.
\begin{eqnarray}
{dy_1 \over dt}&=& f(y,t)\\
{dy_2 \over dt}&=&y_1
\end{eqnarray}

It is up to you to decide whether to solve first the problem for a single 
equations (for $v$) and then the coupled system of two equations 
(for $v$ and $x$) or write directly the general solution for a system of 
equations. Use $x(t=0)=0$ as initial condition. 


\end{itemize}

In writing the program, keep the function to integrate the equations in a 
separate file.

\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. You have to implement 
at least one of the two.

\begin{description}

\item[b)]  
Use the  Runge Kutta method with adaptive step size 
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. Optional: generalize to a system of two coupled 
equations 
to calculate velocity and position as a function of time.

\item[c)] Same as in b) but using the Runge Kutta Fehlberg method. 
\end{description}
\end{document}


\end{document}