[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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;
}