@cindex basis splines, B-splines@cindex splines, basisThis chapter describes functions for the computation of smoothingbasis splines (B-splines). The header file @file{gsl_bspline.h}contains prototypes for the bspline functions and related declarations.@menu* Overview of B-splines::* Initializing the B-splines solver::* Constructing the knots vector::* Evaluation of B-spline basis functions::* Evaluation of B-spline basis function derivatives::* Example programs for B-splines::* References and Further Reading::@end menu@node Overview of B-splines@section Overview@cindex basis splines, overviewB-splines are commonly used as basis functions to fit smoothingcurves to large data sets. To do this, the abscissa axis isbroken up into some number of intervals, where the endpointsof each interval are called @dfn{breakpoints}. These breakpointsare then converted to @dfn{knots} by imposing various continuityand smoothness conditions at each interface. Given a nondecreasingknot vector@c{$t = \{t_0, t_1, \dots, t_{n+k-1}\}$}@math{t = @{t_0, t_1, @dots{}, t_@{n+k-1@}@}},the @math{n} basis splines of order @math{k} are defined by@tex\beforedisplay$$B_{i,1}(x) = \left\{\matrix{1, & t_i \le x < t_{i+1}\cr                            0, & else}\right.$$$$B_{i,k}(x) = \left[(x - t_i)/(t_{i+k-1} - t_i)\right] B_{i,k-1}(x) +             \left[(t_{i+k} - x)/(t_{i+k} - t_{i+1})\right] B_{i+1,k-1}(x)$$\afterdisplay@end tex@ifinfo@exampleB_(i,1)(x) = (1, t_i <= x < t_(i+1)             (0, elseB_(i,k)(x) = [(x - t_i)/(t_(i+k-1) - t_i)] B_(i,k-1)(x) + [(t_(i+k) - x)/(t_(i+k) - t_(i+1))] B_(i+1,k-1)(x)@end example@end ifinfofor @math{i = 0, @dots{}, n-1}. The common case of cubic B-splinesis given by @math{k = 4}. The above recurrence relation can beevaluated in a numerically stable way by the de Boor algorithm.If we define appropriate knots on an interval @math{[a,b]} thenthe B-spline basis functions form a complete set on that interval.Therefore we can expand a smoothing function as@tex\beforedisplay$$f(x) = \sum_{i=0}^{n-1} c_i B_{i,k}(x)$$\afterdisplay@end tex@ifinfo@examplef(x) = \sum_i c_i B_(i,k)(x)@end example@end ifinfogiven enough @math{(x_j, f(x_j))} data pairs. The @math{c_i} canbe readily obtained from a least-squares fit.@node Initializing the B-splines solver@section Initializing the B-splines solver@cindex basis splines, initializingUsing B-splines requires a gsl_bspline_workspace:@deftypefun {gsl_bspline_workspace *} gsl_bspline_alloc (const size_t @var{k}, const size_t @var{nbreak})This function allocates a workspace for computing B-splines of order@var{k}. The number of breakpoints is given by @var{nbreak}. Thisleads to @math{n = nbreak + k - 2} basis functions. Cubic B-splinesare specified by @math{k = 4}. The size of the workspace is@math{O(5k + nbreak)}.@end deftypefun@deftypefun void gsl_bspline_free (gsl_bspline_workspace * @var{w})This function frees the memory associated with the workspace @var{w}.@end deftypefunEvaluation of B-spline basis function derivatives additionally requiresa @code{gsl_bspline_deriv_workspace}:@deftypefun {gsl_bspline_deriv_workspace *} gsl_bspline_deriv_alloc (const size_t @var{k})This function allocates a workspace for computing the derivatives of aB-spline basis function of order @var{k}.  The size of the workspaceis @math{O(2k^2)}.@end deftypefun@deftypefun void gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * @var{w})This function frees the memory associated with the derivativeworkspace @var{w}.@end deftypefun@node Constructing the knots vector@section Constructing the knots vector@cindex knots@deftypefun int gsl_bspline_knots (const gsl_vector * @var{breakpts}, gsl_bspline_workspace * @var{w})This function computes the knots associated with the given breakpointsand stores them internally in @code{w->knots}.@end deftypefun@deftypefun int gsl_bspline_knots_uniform (const double a, const double b, gsl_bspline_workspace * @var{w})This function assumes uniformly spaced breakpoints on @math{[a,b]}and constructs the corresponding knot vector using the previouslyspecified @var{nbreak} parameter. The knots are stored in@code{w->knots}.@end deftypefun@node Evaluation of B-spline basis functions@section Evaluation of B-splines@cindex basis splines, evaluation@deftypefun int gsl_bspline_eval (const double @var{x}, gsl_vector * @var{B}, gsl_bspline_workspace * @var{w})This function evaluates all B-spline basis functions at the position@var{x} and stores them in @var{B}, so that the @math{i}th elementof @var{B} is @math{B_i(x)}. @var{B} must be of length@math{n = nbreak + k - 2}. This value may also be obtained by calling@code{gsl_bspline_ncoeffs}.It is far more efficient to compute all of the basis functions atonce than to compute them individually, due to the nature of thedefining recurrence relation.@end deftypefun@deftypefun int gsl_bspline_eval_nonzero (const double @var{x}, gsl_vector * @var{Bk}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w})This function evaluates all potentially nonzero B-spline basisfunctions at the position @var{x} and stores them in @var{Bk}, sothat the @math{i}th element of @var{Bk} is @math{B_(istart+i)(x)}.The last element of @var{Bk} is @math{B_(iend)(x)}.  @var{Bk} must beof length @math{k}.  Only returning nonzero basis functions allows usersto more cheaply perform tasks requiring linear combinations of the basisfunctions, e.g.  when evaluating an interpolated function.@end deftypefun@deftypefun size_t gsl_bspline_ncoeffs (gsl_bspline_workspace * @var{w})This function returns the number of B-spline coefficients given by@math{n = nbreak + k - 2}.@end deftypefun@node Evaluation of B-spline basis function derivatives@section Evaluation of B-spline derivatives@cindex basis splines, derivatives@deftypefun int gsl_bspline_deriv_eval (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, gsl_bspline_workspace * @var{w}, gsl_bspline_deriv_workspace * @var{dw})This function evaluates all B-spline basis function derivatives of orders@math{0} through @math{nderiv} (inclusive) at the position @var{x}and stores them in @var{dB}.  The @math{(i,j)}th element of @var{dB}is @math{d^j/dx^j B_i(x)}.  @var{dB} must be of size @math{n = nbreak +k - 2} by @math{nderiv + 1}.  The value @math{n} may also be obtainedby calling @code{gsl_bspline_ncoeffs}.  Note that function evaluationsare included as the @math{0}th order derivatives in @var{dB}.It is far more efficient to compute all of the basis functions derivativesat once than to compute them individually, due to the nature of thedefining recurrence relation.@end deftypefun@deftypefun int gsl_bspline_deriv_eval_nonzero (const double @var{x}, const size_t @var{nderiv}, gsl_matrix * @var{dB}, size_t * @var{istart}, size_t * @var{iend}, gsl_bspline_workspace * @var{w}, gsl_bspline_deriv_workspace * @var{dw})This function evaluates all potentially nonzero B-spline basis functionderivatives of orders @math{0} through @math{nderiv} (inclusive) atthe position @var{x} and stores them in @var{dB}.  The @math{(i,j)}thelement of @var{dB} is @math{d^j/dx^j B_(istart+i)(x)}.  The last rowof @var{dB} contains @math{d^j/dx^j B_(iend)(x)}.  @var{dB} must beof size @math{k} by at least @math{nderiv + 1}.  Note that functionevaluations are included as the @math{0}th order derivatives in @var{dB}.Only returning nonzero basis functions allows users to more cheaplyperform tasks requiring linear combinations of the basis functions, e.g.when evaluating an interpolated function.@end deftypefun@node Example programs for B-splines@section Example programs for B-splines@cindex basis splines, examplesThe following program computes a linear least squares fit to data usingcubic B-spline basis functions with uniform breakpoints. The data isgenerated from the curve @math{y(x) = \cos{(x)} \exp{(-x/10)}} on@math{[0, 15]} with gaussian noise added.@example@verbatiminclude examples/bspline.c@end exampleThe output can be plotted with @sc{gnu} @code{graph}.@example$ ./a.out > bspline.datchisq/dof = 1.118217e+00, Rsq = 0.989771$ graph -T ps -X x -Y y -x 0 15 -y -1 1.3 < bspline.dat > bspline.ps@end example@iftex@sp 1@center @image{bspline,3.4in}@end iftex@node References and Further Reading@section References and Further ReadingFurther information on the algorithms described in this section can befound in the following book,@itemize @asisC. de Boor, @cite{A Practical Guide to Splines} (1978), Springer-Verlag,ISBN 0-387-90356-9.@end itemize@noindentA large collection of B-spline routines is available in the@sc{pppack} library available at @uref{http://www.netlib.org/pppack},which is also part of @sc{slatec}.