38 template <
typename Tin,
typename Tout>
class Chebyshev
53 template <
class Functor>
Chebyshev (Functor
const & f,
int n, Tin a, Tin b)
88 template <
class Functor>
89 void generate (Functor
const & f,
int n, Tin a, Tin b);
96 Tout ret = 0.5 * C[0];
97 const double xp = scale (x);
99 for (
int k = 1; k < N; k++)
101 double Tk_x =
cos(k * acos(xp));
115 Tout d_j = 0, d_jp1 = 0, d_jp2 = 0;
116 const double one_x = scale(x);
117 const double two_x = 2 * one_x;
119 for (
int j = m - 1; j >= 1; j--)
121 d_j = two_x * d_jp1 - d_jp2 + C[j];
127 d_j = one_x * d_jp1 - d_jp2 + 0.5 * C[0];
152 for (
int k = 1; k < N; k++)
157 if (abs_Ck <= eps * sum)
170 Tout
approx (
double x,
double eps,
int * n = 0)
const
172 Tout ret = 0.5 * C[0];
173 double xp = scale(x);
176 for (k = 1; k < N; k++)
178 double Tk_x =
cos(k * acos(xp));
179 Tout delta = C[k] * Tk_x;
232 ret.C[N-1] = ret.m * C[N-2] / (2.*(N-2.));
234 for (
int i = 1; i < N - 1; i++)
235 ret.C[i] = ret.m * (C[i-1] - C[i+1]) / (2.*i);
248 for (
int i = 1; i < N; i += 2)
249 ret.C[0] += ret.C[i];
250 for (
int i = 2; i < N; i += 2)
251 ret.C[0] -= ret.C[i];
257 for (
int i = 1; i < N; i++)
259 ret.C[0] += ret.C[i];
260 ret.C[i] = -ret.C[i];
281 static Tin
root (
int N,
int k, Tin x1 = 0., Tin x2 = 1.)
283 return x1 + 0.5 * (1. +
cos(M_PI * (k + 0.5) / N)) * (x2 - x1);
293 std::ostringstream out;
298 for (
int i = 0; i < N - 1; i++)
302 out <<
"]" << std::endl;
318 inline static double node (
int k,
int N)
320 return cos(M_PI * (k + 0.5) / N);
326 inline double scale(Tin x)
const {
331 inline Tin unscale(
double x)
const {
357 template<>
template <
class Functor>
369 fftw_plan plan = fftw_plan_r2r_1d (N, &fvals[0], &C[0], FFTW_REDFT10, 0);
372 double pi_over_N = M_PI / N;
373 for (
int k = 0; k < N; k++)
375 double xk =
cos(pi_over_N * (k + 0.5));
376 fvals[k] = f(unscale(xk));
381 fftw_destroy_plan(plan);
395 template<>
template <
class Functor>
404 cArray fvals(4*N), ftraf(4*N);
407 fftw_plan plan = fftw_plan_dft_1d (
409 reinterpret_cast<fftw_complex*>(&fvals[0]),
410 reinterpret_cast<fftw_complex*>(&ftraf[0]),
416 double pi_over_N = M_PI / N;
417 for (
int k = 0; k < N; k++)
419 double xk =
cos(pi_over_N * (k + 0.5));
420 fvals[2*k+1] = fvals[4*N-(2*k+1)] = f(unscale(xk));
425 fftw_destroy_plan(plan);
429 for (
int i = 0; i < N; i++)
430 C[i] = ftraf[i] * scal;
Chebyshev(Chebyshev const &cb)
Definition: chebyshev.h:43
Definition: chebyshev.h:200
NumberArray< T > cos(NumberArray< T > const &A)
Return per-element cosine.
Definition: arrays.h:1747
virtual size_t resize(size_t n)
Resize array.
Definition: arrays.h:685
Definition: chebyshev.h:201
Tout approx(double x, double eps, int *n=0) const
Evaluate expansion.
Definition: chebyshev.h:170
int tail(double eps) const
Get significant number of coefficients.
Definition: chebyshev.h:147
static double node(int k, int N)
Chebyshev node from (0,1).
Definition: chebyshev.h:318
Tout operator()(Tin const &x) const
Definition: chebyshev.h:94
void generate(Functor const &f, int n, Tin a, Tin b)
Compute the transformation.
size_t size() const
Length of the array (number of elements).
Definition: arrays.h:276
T sum(const ArrayView< T > v)
Sum elements in array.
Definition: arrays.h:1770
Chebyshev()
Definition: chebyshev.h:42
Integration
Definition: chebyshev.h:198
std::string str() const
Convert to string.
Definition: chebyshev.h:291
Chebyshev integrate(Integration itype=Integ_Indef) const
Get expansion of integral of the approximated function.
Definition: chebyshev.h:223
Chebyshev(ArrayView< Tout > const &array, Tin a, Tin b)
Constructor from reference to array.
Definition: chebyshev.h:65
Chebyshev(Functor const &f, int n, Tin a, Tin b)
Constructor.
Definition: chebyshev.h:53
Chebyshev approximation.
Definition: chebyshev.h:38
NumberArray< Tout > const & coeffs() const
Return reference to the coefficient array.
Definition: chebyshev.h:308
rArray abs(const cArrayView u)
Definition: arrays.cpp:19
Tout clenshaw(Tin const &x, int const &m) const
Evaluate expansion.
Definition: chebyshev.h:113
static Tin root(int N, int k, Tin x1=0., Tin x2=1.)
Chebyshev node in interval.
Definition: chebyshev.h:281
Definition: chebyshev.h:199