help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Iryna Feuerstein
Subject: Re: [Help-gsl] How to calculate integral of a b-spline?
Date: Wed, 08 Feb 2012 12:22:15 +0100
User-agent: Mozilla/5.0 (Windows NT 6.1; WOW64; rv:10.0) Gecko/20120129 Thunderbird/10.0

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



reply via email to

[Prev in Thread] Current Thread [Next in Thread]