Loading...
Searching...
No Matches
jacobi_roots.icc
Go to the documentation of this file.
// zeta = roots of the jacobi polynomial(alpha,beta) of order R
#include "rheolef/compiler.h"
#include "rheolef/compiler_eigen.h"
namespace rheolef {
template <class Size, class T, class OutputIterator>
void jacobi_roots (Size R, T alpha, T beta, OutputIterator zeta) {
if (R == 0) return;
if (R == 1) { zeta [0] = (beta-alpha)/(alpha+beta+2); return; }
using namespace Eigen;
Matrix<T,Dynamic,1> diag(R), subdiag(R-1);
diag[0] = (beta-alpha)/(alpha+beta+2);
for (size_t r = 1; r < R; r++) {
diag[r] = (sqr(beta) - sqr(alpha))/((alpha+beta+2*r)*(alpha+beta+2*r+2));
}
subdiag[0] = 2/(alpha+beta+2)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
for (size_t r = 2; r < R; r++) {
subdiag[r-1] = 2/(alpha+beta+2*r)
*sqrt(r*(alpha+r)*(beta+r)*(alpha+beta+r)/((alpha+beta+2*r-1)*(alpha+beta+2*r+1)));
}
SelfAdjointEigenSolver<Matrix<T,Dynamic,1> > es;
es.computeFromTridiagonal (diag, subdiag, EigenvaluesOnly);
for (size_t i = 0; i < R; ++i) zeta[i] = es.eigenvalues()(i);
if (alpha == beta && R%2 == 1) zeta[R/2] = 0.;
}
} // namespace rheolef
T
Expr1::float_type T
Definition field_expr.h:230
Eigen
Definition field_eigen.h:35
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::jacobi_roots
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
Definition jacobi_roots.icc:26
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56