Rheolef  7.2
an efficient C++ finite element environment
jacobi_roots.icc
Go to the documentation of this file.
1 // zeta = roots of the jacobi polynomial(alpha,beta) of order R
21 #include "rheolef/compiler.h"
22 #include "rheolef/compiler_eigen.h"
23 
24 namespace rheolef {
25 template <class Size, class T, class OutputIterator>
26 void jacobi_roots (Size R, T alpha, T beta, OutputIterator zeta) {
27  if (R == 0) return;
28  if (R == 1) { zeta [0] = (beta-alpha)/(alpha+beta+2); return; }
29  using namespace Eigen;
30  Matrix<T,Dynamic,1> diag(R), subdiag(R-1);
31  diag[0] = (beta-alpha)/(alpha+beta+2);
32  for (size_t r = 1; r < R; r++) {
33  diag[r] = (sqr(beta) - sqr(alpha))/((alpha+beta+2*r)*(alpha+beta+2*r+2));
34  }
35  subdiag[0] = 2/(alpha+beta+2)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
36  for (size_t r = 2; r < R; r++) {
37  subdiag[r-1] = 2/(alpha+beta+2*r)
38  *sqrt(r*(alpha+r)*(beta+r)*(alpha+beta+r)/((alpha+beta+2*r-1)*(alpha+beta+2*r+1)));
39  }
40  SelfAdjointEigenSolver<Matrix<T,Dynamic,1> > es;
41  es.computeFromTridiagonal (diag, subdiag, EigenvaluesOnly);
42  for (size_t i = 0; i < R; ++i) zeta[i] = es.eigenvalues()(i);
43  if (alpha == beta && R%2 == 1) zeta[R/2] = 0.;
44 }
45 
46 } // namespace rheolef
Expr1::float_type T
Definition: field_expr.h:230
Float alpha[pmax+1][pmax+1]
Definition: bdf.icc:28
This file is part of Rheolef.
csr< T, M > diag(const vec< T, M > &d)
Definition: csr.cc:56
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
Float beta[][pmax+1]