Theory.

1. Find the files gsl_multimin.h and gsl_multiroots.h in your box; Hint: gsl-config --cflags
2. Suppose your main funtion writes the string "out" to stdout and the string "err" to stderr. What (and why) will be the content of files "a" and "b" after the following commands,

   ./main >a 2>b
   ./main 1>a 2>b
   ./main >a | tee b
   ./main >a | tee b 1>/dev/null
   ./main 2>&1 >a | tee b
   ./main >a 2>&1 | tee b
   ./main >a 2>/dev/null | tee b
   ./main 2>b 1>2
   

   Hint: redirection (computing).

Practice.

1. Multidimensional minimization with Nelder-Mead algorithm gsl_multimin_fminimizer_nmsimplex2.

   • Suppose that an experiment on measuring the activity of a radioactive substance as function of time gives the following result

     # ti,  yi,   σi
     0.23  4.64  0.42
     1.29  3.37  0.37
     2.35  3.00  0.34
     3.41  2.55  0.31
     4.47  2.29  0.29
     5.53  1.66  0.27
     6.59  1.58  0.26
     7.65  1.69  0.25
     8.71  1.37  0.24
     9.77  1.45  0.24
     

     or

     double t[] = {0.23,1.29,2.35,3.41,4.47,5.53,6.59,7.65,8.71,9.77};
     double y[] = {4.64,3.38,3.01,2.55,2.29,1.67,1.59,1.69,1.38,1.46};
     double e[] = {0.42,0.37,0.34,0.31,0.29,0.27,0.26,0.25,0.24,0.24};
     int N = sizeof(t)/sizeof(t[0]);
     

     where the time t_i, activity y_i, and the experimental error σ_i are given in some rescaled units.

   • Fit the function f(t)=A*exp(-t/T)+B where A (initial amount), T (lifetime), and B (background noise) are the fitting parameters, to the data and determine the lifetime, T, of the substance.

   • The (master) function to minimize is F(A,T,B)=∑i(f(ti)-yi)²/σi². That is, you have a three-dimensional minimization problem in the space of the parameters A, T, and B.

   • Plot the experimental data and the fitting function.

   • Hints:

     The master function to minimimize can be like this:

     struct myparams {int N; double *t,*y,*e;};
     
     double master(const gsl_vector *x, void *params){
     	struct myparams *p = (struct myparams*)params;
     	int     N = p->N;
     	double *t = p->t;
     	double *y = p->y;
     	double *e = p->e;
     	double  A = gsl_vector_get(x,0);
     	double  T = gsl_vector_get(x,1);
     	double  B = gsl_vector_get(x,2);
     	double sum=0;
     	for(int i=0;i<N;i++) sum+ = pow( ( A*exp(-t[i]/T)+B - y[i] )/e[i] ,2);
     	return sum;
     }