\documentclass[12pt]{article}
\textwidth 16cm
\textheight 23cm
\topmargin 0cm
\evensidemargin 0cm
\oddsidemargin 0cm
\usepackage{epsfig}
\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.~3: General nonlinear least squares fitting}\\
\vspace{0.4cm}

In many experiments it happens that a given line of resonant absorption splits 
under the effect of a given perturbation, such as a magnetic or electric field. In this 
exercise the purpose is to extract the value of the splitting by applying a 
general  nonlinear least squares fitting to the data. In figure we present a 
case rather simple but illustrative of a whole class of problems.
\begin{figure}[h]
\epsfig{file=figex3.eps}
%\caption{caption}
\end{figure}
The solid line is the function $f_0^{exp}$ due to resonant absorption at 
frequency $\omega=2.0$ measured in 
absence of external fields,  the dotted line the function $f^{exp}$ due to a 
splitting of the 
line into two contributions of half intensity and an amount $\omega_s$ 
apart (the subscript $s$ stays for splitting). 
Physical reasons  tell us the line  $f^{exp}$ could be fitted by a 
function $f^{calc}$ which can be expressed in terms of the unperturbed 
lineshape $f_0^{exp}$ as:
\begin{equation}
f^{calc}(\omega_i,\omega_s)= {1 \over 2}
[ f_0^{exp}(\omega_i + {\omega_s \over 2}) + 
  f_0^{exp}(\omega_i - {\omega_s \over 2})]
\end{equation}

The interesting point of this exercise is that you do not need to know the  
explicit form giving the shape of the absorption line but you can use the 
unperturbed data $f_0^{exp}$ and formula (1) as fitting function to 
extract the fitting parameter $\omega_s$ 
by the data representing $f^{exp}(\omega_i)$.

The file splitting.dat contains 4000 values of $\omega$, $f_0^{exp}$ 
and $f^{exp}$ in 
the first, second and third column respectively, that is to say 
$\omega_i$, $f_0^{exp}(\omega_i)$ 
and $f^{exp}(\omega_i)$ for $i=1,..., N$ with $N=4000$. 
These values will be used
to find the value of $\omega_s$ of the splitting which best describes the data
with two digit precision. You can read this file with matlab and store the 
values in vectors with the following commands:

\begin{tabular}{lll}
load && 'splitting.dat'\\
omega& =& splitting(:,1)\\ 
f0   & =& splitting(:,2)\\
fexp & =& splitting(:,3)\\
\end{tabular}

Consider that you have data in units of 0.001 and express your tentative 
values of $\omega_s$ and of $\delta\omega_s$ in multiples of this value. 
Besides, since you have to shift $f_0^{exp}$ to evaluate 
$f^{calc}(\omega_i,\omega_s)$ consider you have less than 4000 points available 
for the fit, namely 4000-($\omega_s$/0.001) points. 

In order to facilitate comparison with the formula's given in the notes 
distributed in the lecture, consider that you have only one fitting parameter 
so that the vector {\bf p}
has only one element $p_1=\omega_s$. Consequently, the matrix $\it{M}$
has only one element so that you have to solve just one equation iteratively.

The values $\omega_i$ take the place of $x_i$, the function $f^{exp}(\omega_i)$ takes the place of $y_i(x_i)$ and 
$f^{calc}(\omega_i,\omega_s)$ that of $y(x_i, {\bf p})$, $i$ running over the 
available experimental points.
First you choose an initial guess for $\omega_s^0$ and for its increment
$\delta\omega_s^0$. With this you can calculate the numerical derivative 
 by forward difference so that you can calculate the only element $E_1$ 
of the vector {\bf E} and $M_{11}$
of the matrix $\it{M}$. Their ratio gives the 
new value of the increment $\delta\omega_s$. Add this increment 
(or better $\sim$half of it in order to avoid large jumps of the value of 
the parameter) to the initial $\omega_s^0$ and 
repeat until convergence is reached. If you have chosen a reasonable first 
guess this should happen in a few iterations. If you do not reach convergance within about 10 iterations stop and try another initial value of $\omega_s$.

To keep the iterative procedure as simple and readable as possible it is 
useful to define a separate function which calculates 
the function $f^{calc}(\omega_i,\omega_s)$
and its derivative with respect to $\omega_s$


\end{document}