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

[Rhys Ulerich]

> 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;
}