f(x,y) = (1-x)² + 100(y-x²)²
Least-squares fit.
Suppose an experiment measuring the activity of a radioactive substance as function of time gives the following result,
# t_i y_i e_i 0.47 5.49 0.26 1.41 4.08 0.12 2.36 3.54 0.27 3.30 2.61 0.10 4.24 2.09 0.15 5.18 1.91 0.11 6.13 1.55 0.13 7.07 1.47 0.07 8.01 1.45 0.15 8.95 1.25 0.09or
double t[]= {0.47,1.41,2.36,3.30,4.24,5.18,6.13,7.07,8.01,8.95}; double y[]= {5.49,4.08,3.54,2.61,2.09,1.91,1.55,1.47,1.45,1.25}; double e[]= {0.26,0.12,0.27,0.10,0.15,0.11,0.13,0.07,0.15,0.09}; int n = sizeof(t)/sizeof(t[0]);where the time ti, activity yi, and the experimental error σi are given in some rescaled units.
Do the least-squares fit of the function
f(t)=A*exp(-t/T)+B(where the initial amount A, the lifetime T, and the background noise B are the fitting parameters) to the data and estimate the lifetime, T, of the substance.
Plot the experimental data and the fitting function.
Hints:
It is probably easier to use the algorithm without derivatives
gsl_multimin_fminimizer_nmsimplex2
.
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.
It can be programmed like this,
struct experimental_data {int n; double *t,*y,*e;}; double function_to_minimize (const gsl_vector *x, void *params) { double A = gsl_vector_get(x,0); double T = gsl_vector_get(x,1); double B = gsl_vector_get(x,2); struct experimental_data *p = (struct experimental_data*) params; int n = p->n; double *t = p->t; double *y = p->y; double *e = p->e; double sum=0; #define f(t) A*exp(-(t)/T) + B for(int i=0;i<n;i++) sum+ = pow( f(t[i]) - y[i] )/e[i] ,2); return sum; }
(Optional) Do the same least-squares fit using one of the algorithms with derivatives. It might be more effective, that is, it might arrive at the optimum with fewer evaluations of the function-to-minimize.