Cubic Spline interpolation

Cubic Spline interpolation between each data pair x_i, x_i+1 is done through a third-grade polynomial function, under the condition that the function f(x), as well as its first and second derivatives f(x)¹ and f(x)², have to be continuum in the limits of the interval.

$x_{i} \leq x \leq x_{i + 1}$

$f_{i} (x) = f_{i} (x_{i}) + (x - x_{i}) f_{i}^{1} (x) + \frac{1}{2} (x - x_{i})^{2} f_{i}^{2} (x_{i}) + \frac{1}{6} (x - x_{i})^{3} f_{i}^{3} (x_{i})$

and, if the third-order derivative is written as a function of the second-order derivative,

$f_{i} (x) = f_{i} (x_{i}) + (x - x_{i}) f_{i}^{1} (x) + \frac{1}{2} (x - x_{i})^{2} f_{i}^{2} (x_{i}) + \frac{1}{6 h} (x - x_{i})^{3} (f_{i}^{2} (x_{i + 1}) - f_{i}^{2} (x_{i}))$

where h is:

$x_{i + 1} = x_{i} + h$

and the continuity conditions:

$\begin{matrix} f_{i} (x_{i + 1}) = f_{i + 1} (x_{i + 1}) \\ f_{i}^{1} (x_{i + 1}) = f_{i + 1}^{1} (x_{i + 1}) \\ f_{i}^{2} (x_{i + 1}) = f_{i + 1}^{2} (x_{i + 1}) \end{matrix}$

Library functions (Numerical Recipes):

spline():Go to function call: it takes six arguments; vectors x[], y[] (which contain the experimental points between which the interpolation will take place), the number of data points n, the variables yp1 and ypn (which contain the values of the first derivative in the first point x_i and in the last one x_N, respectively), and the vector y2[] (which returns the second order derivatives in each of the x_i points).
splint():Go to function call: it takes six arguments; vectors xa[], ya[] (which contain the experimental points between which the interpolation will take place), the vector y2a[] (which returns the second order derivatives in each of the x_i points, obtained from spline()), the number of data points n, the variable x (which is the abscissa for which we want to interpolate), and the variable y (used to return the the interpolated value).

The files for this calculation are in folder breitwigner/2_CubicSplines. It uses the Numerical Recipes libraries, spline.c and splint.c.

The program splinesc.c is very similar to lagrange.c and works in a similar way as explained in the Lagrange interpolations page.

Once we have all the inputs, we can call the spline() function (you may want to check the definition again) and then we call splint() inside a loop through the points (See the definition).

	/***************************************************
	 * Calculating the first derivative,
	 * evaluating in E1 and E1 and
	 * calling the functions
	 ***************************************************/
	yp1 = -0,00001;
	ypn = -0,5;
	
	spline(xa, ya, N, yp1, ypn, y2);
	for (i=1; i<=n; i++)
		splint(xa, ya, y2, N, x[i], &y[i]);

Results

And so the following values are obtained:

σ₀ = 70297.063
E_r = 76.100 MeV
γ = 57.986 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 Cubic splines

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

Breit-Wigner curve obtained with Cubic splines. Zoom in

Zoom in of the previous plot.