Numerical Recipes applied
to
Thermodynamics

Method of the Lagrange interpolators

Given a group of N pairs of experimental data (x1, y1), (x2, y2), ... , (xN, yN), we want to obtain the values (x, y) interpolated from the experimental values using the Lagrange polynomial of grade N-1:

y=(x-x2)(x-x3)(x-xN)(x1-x2)(x1-x3)(x1-xN)y1+(x-x1)(x-x3)(x-xN)(x2-x1)(x2-x3)(x2-xN)y2+(x-x1)(x-x2)(x-xN-1)(xN-x1)(xN-x2)(xN-xN-1)yN

y=k=1Nj=1jkN(x-xj)(xk-xj)yk

The formula was first published by Waring (1779), rediscovered by Euler in 1783, and published by Lagrange in 1795 (Jeffreys and Jeffreys 1988).

The files for this calculation are in folder breitwigner/1_Lagrange. It uses the Numerical Recipes libraries and polint.c.

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

Then, a call is made to the polint.c routine to make the interpolation in an interval fabricated in execution time. To generate the interval, the program asks for a minimum and maximum value. Ideally, those values should be around the peak of the resonance, so that we can have more points in that region where the curve changes the most. Looking at the figure in the main page, the peak is around E = 75 MeV, so good limit values could be Emin = 70.1 MeV and Emax =79.9 MeV.

To choose the step size, we need to keep in mind that when constructing interpolating polynomials, there is a trade-off between having a better fit and having a smooth well-behaved fitting function. The more data points that are used in the interpolation, the higher the degree of the resulting polynomial, and therefore the greater oscillation it will exhibit between the data points. Therefore, a high-degree interpolation may be a poor predictor of the function between points, although the accuracy at the data points will be "perfect". In this case, taking into account the limit values previously choosen, a good step size could be 0.5.

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

	/***************************************************
	 * Creating the interpolation interval
	 ***************************************************/
	printf(" Lower limit Emin of the interpolation interval: ");
	scanf("%lf", &jmin);	
	printf(" Upper limit Emax of the interpolation interval: ");
	scanf("%lf", &jmax);
	printf(" Step size: ");
	scanf("%lf", &dx);
	
	fprintf(fp,"\n Lower limit Emin of the interpolation interval: Emin = %3.3lf", jmin);
	fprintf(fp,"\n Upper limit Emax of the interpolation interval: Emax = %3.3lf", jmax);
	fprintf(fp,"\n Step size: dE = %lf", dx);
		
	n  = (int)( (jmax - jmin) / dx) + 2;	
	x  = dvector(1,n);
	y  = dvector(1,n);
	dy = dvector(1,n);
	getchar();
	
	for (i=1, j=jmin; i<=n-1, j<jmax; i++, j=j+dx)
		x[i] = j;
	x[n] = jmax;
	/***************************************************
	 * Calling the interpolation function
	 ***************************************************/
	for (i=1; i<=n; i++)
		polint(xa, ya, N, x[i], &y[i], &dy[i]);

Once we have obtained the interpolated values of the cross-section with their errors from the generated energies, we can calculate the typical properties of the Breit-Wigner resonance, and their errors:

(See figure in the problem's page)

	m      = maximum (y, n);
	gamma  = 2.0*x[m] / sqrt( (y[m]/ya[1]) - 1.0);
	dgamma = ( gamma / 2*(1.0 - (ya[1]/y[m])) ) * sqrt(pow(0.1/ya[1],2) + pow(dy[m]/y[m], 2));
	So     = y[m]*pow(gamma,2) / 4.0;
	dSo    = So*sqrt(pow(dgamma/gamma, 2) + pow(2.0*dy[m] / y[m],2));

Finally, the results of the calculation are dumped into an output.txt file.

Results

As stated previously, looking at the figure, the maximum seemed to be around E = 75 MeV. With the calculation, a value of Er = 74.60 MeV with an error of 0.04 MeV is obtained.

For gamma we obtain a value of γ = 56.9 MeV with an error of 0.2 MeV.

The region of the curve where there is more variation is the peak, that's why we take more points there. Repeat the calculation adding points for the rest of the curve.

Breit-Wigner curve obtained with Lagrange interpolators

The peak is the region of higher variation of the curve.

Breit-Wigner curve obtained with Lagrange interpolators. Zoom in

Zoom in of the previous plot.