29 typename std::iterator_traits<OutputIterator1>::value_type alpha,
30 typename std::iterator_traits<OutputIterator1>::value_type beta,
31 OutputIterator1 zeta, OutputIterator2 omega)
33 typedef typename std::iterator_traits<OutputIterator1>::value_type
T;
34 check_macro (R >= 2,
"gauss_lobatto_jacobi: node number " << R <<
" may be >= 2");
35 T num =
pow(
T(2), alpha+beta+1)/
T(1.*R-1);
36 T w0 =
pow(
T(2), alpha+beta+1)*(beta+1);
37 T wf =
pow(
T(2), alpha+beta+1)*(alpha+1);
38 if (alpha == floor(alpha) && beta == floor(beta)) {
39 num *=
T(1.*R)/((alpha+R)*(beta+R));
40 w0 *= 1/(
T(R-1.)*(alpha+R));
41 wf *= 1/(
T(R-1.)*(beta+R));
42 for (Size k = 1; k <= size_t(static_cast<int>(beta)); k++) {
43 num *=
T(1.*R+k)/(alpha+
T(1.*R+k));
44 w0 *=
T(sqr(
int(k)))/(
T(1.*R-1+k)*(alpha+R+k));
46 for (Size k = 1; k <= size_t(static_cast<int>(alpha)); k++) {
47 wf *=
T(sqr(
int(k)))/(
T(1.*R-1+k)*(beta+
T(1.*R+k)));
61 for (Size r = 1; r < R-1; ++r)
62 omega[r] = num/sqr(P1(zeta[r]));
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void gauss_lobatto_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)