Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
jacobi_roots.icc
Go to the documentation of this file.
1
20// zeta = roots of the jacobi polynomial(alpha,beta) of order R
21#include "rheolef/compiler.h"
22#include "rheolef/compiler_eigen.h"
23
24namespace rheolef {
25template <class Size, class T, class OutputIterator>
26void 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
This file is part of Rheolef.
void jacobi_roots(Size R, T alpha, T beta, OutputIterator zeta)
csr< T, M > diag(const vec< T, M > &d)
Definition csr.cc:56