Rheolef  7.1
an efficient C++ finite element environment
jacobi.icc
Go to the documentation of this file.
1 #include <cstdlib>
22 namespace rheolef {
23 template <class T>
24 class jacobi {
25  public:
26  jacobi (size_t R1, T alpha1, T beta1)
27  : R(R1), alpha(alpha1), beta(beta1) {}
28  T operator() (T x) const {
29  if (R == 0) return 1;
30  T P_prev = 1;
31  T P = ((alpha+beta+2)*x + alpha - beta)/2;
32  for (size_t r = 1; r < R; r++) {
33  T P_save = P;
34  T a = (beta*beta - alpha*alpha)/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+2)));
35  T b = 2*(alpha+T(1.*r))*(beta+T(1.*r))/((alpha+beta+T(2.*r))*(alpha+beta+T(2.*r+1)));
36  T c = 2*T(r+1.)*(alpha+beta+T(r+1.))/((alpha+beta+T(2.*r+1))*(alpha+beta+T(2.*r+2)));
37  P = ((x-a)*P - b*P_prev)/c;
38  P_prev = P_save;
39  }
40  return P;
41  }
42  protected:
43  size_t R;
45 };
46 } // namespace rheolef
jacobi(size_t R1, T alpha1, T beta1)
Definition: jacobi.icc:26
T operator()(T x) const
Definition: jacobi.icc:28
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.