Hello Rhys,
thank you very much for the example. It's the thing, I was looking for,
and I can use this approach in my program now.
Best wishes,
Iryna.
Am 07.02.2012 23:51, schrieb 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;
}