Numerical Recipes applied
to
Thermodynamics

Non-linear least squares fitting: Levemberg-Marquardts

Given a group of M pairs of experimental data (x1, y1), (x2, y2), ... , (xM, yM), having to fulfill certain function, and which depend on the variable x as well as on a group of N parameters aj

y(x)=y(x,a1,a2,...,aN)

If the dependency of the aj parameters is non-linear, the y(x) function is expanded by means of a Taylor series of first order around the point aj0 (the initial values of the coefficients)

y(x)=y0(x)+j=1Ny0(x)ajδaj

where

δaj=aj-aj0

and yj represents the value of the function for aj.

The least squares method requires to minimize the expression:

χ2=i=1M1σi2yi-y0(xi)-j=1Ny0(xi)ajδaj2

with respect to each one of the coefficients, i.e.:

χ2δak=0        k=1,...,N

which leads to the following system of equations:

i=1Mj=1N1σi2yo(xi)ajyo(xi)akδaj=i=1M1σi2yi-y0(xi)yo(xi)ak        k=1,...,N

which can be written in matrix notation:

α δa=β

where:

αjk=i=1M1σi2y0(xi)ajy0(xi)ak

βk=i=1M1σi2yi-y0(xi)y0(xi)ak

The method of Marquardt exchanges matrix α in the previous equation by matrix α':

α' δa=β

α'jj=αjj(1+λ)

α'ij=αij        (ij)

The files for this calculation are in folder breitwigner/3_Marquardt. It uses the Numerical Recipes libraries mentioned above, plus covsrt.c and gausj.c.

The program marquardt.c reads the experimental data from data.txt and writes the results in an output.txt file. The format of this file is given in three columns, where the first represents the energy, the second, the cross-section measured at that energy, and the third, the uncertainties. The values are taken from the experimental data of neutron scattering cross-sections provided in the main page.

In this case we have 3 parameters, σ, Er and γ. We have to provide initial values for them to start the iteration, and they will be adjusted by the mrqmin() function. The program fits all the data points by default, with ia[i] = 1.

	/***************************************************
	 * Initial values for the parameter vectors
	 * For example:
	 *    a[1] = 70000.0;
	 *    a[2] = 75.0;
	 *    a[3] = 60.0;
	 ***************************************************/
	for (i=1; i<=N; i++) {
		printf("Enter initial value of parameter %d: ", i);
		scanf("%lf", &a[i]);
	}

	// Fit all
	for (i=1; i<=N; i++)
		ia[i] = 1;

Once we have all the inputs, we can call the mrqmin() function (you may want to check the definition again).

	/***************************************************
	 * ITERATION LOOP
	 ***************************************************/
	char answer;
	do {
		printf("\nHow many iterations?: ");
		scanf("%d", &k);

		for(i=1; i<=k; i++) {
			iter++;
			ochisq = chisq;

			mrqmin(x, y, sig, M, a, ia, N, covar, alpha, &chisq, funcs, &alamda);

			if (chisq > ochisq)
				n=0;
			else if (fabs(ochisq - chisq) < 0.1)
				n++;

			if(n >= 4)
				break;
		}

		if(n < 4) {
			printf("\nIterate again? (Y/N)?: ");
			scanf("%c", &answer);
		}
	} while (answer == 's' || answer == 'S');

Finally, we can use the user defined function func() to recalculate the experimental points (see func() definition ).

	printf(     "\nData points recalculated from the fit:\n");
	fprintf(fp, "\nData points recalculated from the fit:\n");
	for (i=1; i<=M; i++) {

		(*funcs)(x[i], a, &ymod[i], dyda, N);
	
		sum = 0.0;
		for (j=1; j<=N; j++)
			sum += pow(dyda[j]*sqrt(covar[j][j]), 2);

		printf(     "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
		fprintf(fp, "\nSigma[%d] = %3.3lf\t\tdSigma[%d] = %3.3lf", i, ymod[i], i, sqrt(sum));
	}

Results

The following values with their errors are obtained for the parameters of interest:

Breit-Wigner curve obtained with the method of least-squares