help-gsl
[Top][All Lists]

## Re: [Help-gsl] How to calculate integral of a b-spline?

 From: Rhys Ulerich Subject: Re: [Help-gsl] How to calculate integral of a b-spline? Date: Tue, 7 Feb 2012 16:51:59 -0600

> I need to calculate an integral of the following type:
>
> \int f(x)*B(i,k)(x) dx,
>
> where B(i,k) is a basis spline of the degree k  with a De Boor point i.

Perhaps the following routine could be a good starting point.  It has
been thoroughly tested.  It computes \int B(i,k)(x) dx using the
lowest cost Gauss-Legendre integration rule that should exactly
recover the analytical results.  The error handling (SUZERAIN_ERROR,
etc) very much resembles that found in the GSL.

As you know the support of B(i,k)(x), you also know the support of
f(x)*B(i,k)(x).  The trick would be either (a) knowing something about
f(x) so you could pick the correct Gauss-Legendre integration order or
(b) substituting a more general integration process.

Best of luck,
Rhys

/**
* Compute the coefficients \f$\gamma_{i} \f$ for <code>0 <= i < w->n</code> *
* such that \f$\vec{\gamma}\cdot\vec{\beta} = \int \sum_{i} \beta_{i} * B_{i}^{(\mbox{nderiv})}(x) \, dx\f$.
*
* @param[in]  nderiv The derivative to integrate.
* @param[out] coeffs Real-valued coefficients \f$\gamma_{i} \f$.
* @param[in]  inc Stride between elements of \c x
* @param[in]  dB Temporary storage to use of size <tt>w->k</tt> by
*             no less than <tt>nderiv + 1</tt>.
* @param[in]  w Workspace to use (which sets the integration bounds).
* @param[in]  dw Workspace to use.
*
* @return ::SUZERAIN_SUCCESS on success.  On error calls suzerain_error() and
*      returns one of #suzerain_error_status.
*/
int
suzerain_bspline_integration_coefficients(
const size_t nderiv,
double * coeffs,
size_t inc,
gsl_matrix *dB,
gsl_bspline_workspace *w,
gsl_bspline_deriv_workspace *dw);

int
suzerain_bspline_integration_coefficients(
const size_t nderiv,
double * coeffs,
size_t inc,
gsl_matrix *dB,
gsl_bspline_workspace *w,
gsl_bspline_deriv_workspace *dw)
{
/* Obtain an appropriate order Gauss-Legendre integration rule */
gsl_integration_glfixed_table * const tbl
= gsl_integration_glfixed_table_alloc((w->k - nderiv + 1)/2);
if (tbl == NULL) {
SUZERAIN_ERROR("failed to obtain Gauss-Legendre rule from GSL",
SUZERAIN_ESANITY);
}

/* Zero integration coefficient values */
for (size_t i = 0; i < w->n; ++i) {
coeffs[i * inc] = 0.0;
}

/* Accumulate the breakpoint-by-breakpoint contributions to coeffs */
double xj = 0, wj = 0;
for (size_t i = 0; i < (w->nbreak - 1); ++i) {

/* Determine i-th breakpoint interval */
const double a = gsl_bspline_breakpoint(i,   w);
const double b = gsl_bspline_breakpoint(i+1, w);

for (size_t j = 0; j < tbl->n; ++j) {

/* Get j-th Gauss point xj and weight wj */
gsl_integration_glfixed_point(a, b, j, &xj, &wj, tbl);

/* Evaluate basis functions at point xj */
size_t kstart, kend;
gsl_bspline_deriv_eval_nonzero(xj, nderiv,
dB, &kstart, &kend, w, dw);

/* Accumulate weighted basis evaluations into coeffs */
for (size_t k = kstart; k <= kend; ++k) {
coeffs[k * inc] += wj * gsl_matrix_get(dB,
k - kstart, nderiv);
}
}
}

/* Free integration rule resources */
gsl_integration_glfixed_table_free(tbl);

return SUZERAIN_SUCCESS;
}