gauss_jacobi.icc
Go to the documentation of this file.
#include "rheolef/compiler.h"
#include "rheolef/gamma.h"
#include "rheolef/jacobi.h"
#include "rheolef/jacobi_roots.h"
#include <iterator>
namespace rheolef {
template <class Size, class OutputIterator1, class OutputIterator2>
void
gauss_jacobi (Size R,
typename std::iterator_traits<OutputIterator1>::value_type alpha,
typename std::iterator_traits<OutputIterator1>::value_type beta,
OutputIterator1 zeta, OutputIterator2 omega)
{
typedef typename std::iterator_traits<OutputIterator1>::value_type T;
T num = pow(T(2.), alpha+beta+3)/sqr(alpha+beta+T(1.*R)+1);
if (alpha == floor(alpha) && beta == floor(beta))
for (Size k = 1; k <= size_t(static_cast<int>(beta)); k++)
num *= (T(1.*R)+T(1.*k))/(alpha+T(1.*R)+T(1.*k));
else
num *= (my_gamma(alpha+T(1.*R)+1)/my_gamma(alpha+beta+T(1.*R)+1))
*(my_gamma(beta+T(1.*R)+1)/my_gamma(T(1.*R)+1));
jacobi_roots (R, alpha, beta, zeta);
jacobi<T> P (R-1, alpha+1, beta+1);
for (Size r = 0; r < R; r++)
omega[r] = num/((1-sqr(zeta[r]))*sqr(P(zeta[r])));
}
} // namespace rheolef
rheolef::my_gamma
T my_gamma(const T &x)
Definition: gamma.icc:25
rheolef::gauss_jacobi
void gauss_jacobi(Size R, typename std::iterator_traits< OutputIterator1 >::value_type alpha, typename std::iterator_traits< OutputIterator1 >::value_type beta, OutputIterator1 zeta, OutputIterator2 omega)
Definition: gauss_jacobi.icc:29
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::jacobi
Definition: jacobi.icc:24
rheolef::jacobi_roots
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
Definition: jacobi_roots.icc:26
rk::beta
Float beta[][pmax+1]
Definition: runge_kutta_semiimplicit.icc:60
T
Expr1::float_type T
Definition: field_expr.h:218