Rheolef
7.1
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
Eigen
Definition:
field_eigen.h:35
bdf::alpha
Float alpha[pmax+1][pmax+1]
Definition:
bdf.icc:28
rheolef::diag
csr< T, M > diag(const vec< T, M > &d)
Definition:
csr.cc:56
rheolef
This file is part of Rheolef.
Definition:
compiler_eigen.h:37
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