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