Hex  1.0
Hydrogen-electron collision solver
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
chebyshev.h
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
2  * *
3  * / / / / __ \ \ / / *
4  * / /__ / / / _ \ \ \/ / *
5  * / ___ / | |/_/ / /\ \ *
6  * / / / / \_\ / / \ \ *
7  * *
8  * Jakub Benda (c) 2014 *
9  * Charles University in Prague *
10  * *
11 \* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
12 
13 #ifndef HEX_CHEBYSHEV
14 #define HEX_CHEBYSHEV
15 
16 #include <cmath>
17 #include <iostream>
18 #include <numeric>
19 #include <vector>
20 
21 #include <fftw3.h>
22 
23 #include "arrays.h"
24 #include "misc.h"
25 
38 template <typename Tin, typename Tout> class Chebyshev
39 {
40 public:
41 
42  Chebyshev () : N(0), C(), xt(0.), m(0.) {}
43  Chebyshev (Chebyshev const & cb) : N(cb.N), C(cb.C), xt(cb.xt), m(cb.m) {}
44 
53  template <class Functor> Chebyshev (Functor const & f, int n, Tin a, Tin b)
54  {
55  generate(f, n, a, b);
56  }
57 
65  Chebyshev (ArrayView<Tout> const & array, Tin a, Tin b)
66  {
67  N = array.size();
68  xt = 0.5 * (b + a);
69  m = 0.5 * (b - a);
70 
71  // copy coefficients
72  C = array;
73  }
74 
88  template <class Functor>
89  void generate (Functor const & f, int n, Tin a, Tin b);
90 
94  Tout operator() (Tin const & x) const
95  {
96  Tout ret = 0.5 * C[0];
97  const double xp = scale (x);
98 
99  for (int k = 1; k < N; k++)
100  {
101  double Tk_x = cos(k * acos(xp));
102  ret += C[k] * Tk_x;
103  }
104  return ret;
105  }
106 
113  Tout clenshaw (Tin const & x, int const & m) const
114  {
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; // due to linearity of 'scale'
118 
119  for (int j = m - 1; j >= 1; j--)
120  {
121  d_j = two_x * d_jp1 - d_jp2 + C[j];
122 
123  d_jp2 = d_jp1;
124  d_jp1 = d_j;
125  }
126 
127  d_j = one_x * d_jp1 - d_jp2 + 0.5 * C[0];
128 
129  return d_j;
130  }
131 
147  int tail (double eps) const
148  {
149  double sum = std::abs(0.5 * C[0]);
150  double abs_Ck;
151 
152  for (int k = 1; k < N; k++)
153  {
154  abs_Ck = std::abs(C[k]);
155  sum += abs_Ck;
156 
157  if (abs_Ck <= eps * sum)
158  return k + 1;
159  }
160 
161  return N;
162  }
163 
170  Tout approx (double x, double eps, int * n = 0) const
171  {
172  Tout ret = 0.5 * C[0];
173  double xp = scale(x);
174 
175  int k;
176  for (k = 1; k < N; k++)
177  {
178  double Tk_x = cos(k * acos(xp));
179  Tout delta = C[k] * Tk_x;
180  ret += delta;
181 
182  if (std::abs(delta) <= eps * std::abs(ret))
183  {
184  k++;
185  break;
186  }
187  }
188 
189  if (n != 0)
190  *n = k;
191 
192  return ret;
193  }
194 
198  typedef enum {
202  } Integration;
203 
224  {
225  Chebyshev ret;
226  ret.xt = xt;
227  ret.m = m;
228  ret.N = N;
229 
230  ret.C.resize(ret.N);
231  ret.C[0] = 0;
232  ret.C[N-1] = ret.m * C[N-2] / (2.*(N-2.));
233 
234  for (int i = 1; i < N - 1; i++)
235  ret.C[i] = ret.m * (C[i-1] - C[i+1]) / (2.*i);
236 
237  // transform coefficients to get a definite integral, if requested
238  switch (itype)
239  {
240  case Integ_Indef:
241  {
242  // already done
243  break;
244  }
245  case Integ_Low:
246  {
247  // compute C[0] as an alternating sum of other coefficients
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];
252  break;
253  }
254  case Integ_High:
255  {
256  // compute C[0] as a sum of other coefficients and negate other coefficients
257  for (int i = 1; i < N; i++)
258  {
259  ret.C[0] += ret.C[i];
260  ret.C[i] = -ret.C[i];
261  }
262  break;
263  }
264  }
265 
266  // take into account the normalization (F = c_0/2 + ...)
267  ret.C[0] *= 2.;
268 
269  return ret;
270  }
271 
281  static Tin root (int N, int k, Tin x1 = 0., Tin x2 = 1.)
282  {
283  return x1 + 0.5 * (1. + cos(M_PI * (k + 0.5) / N)) * (x2 - x1);
284  }
285 
291  std::string str() const
292  {
293  std::ostringstream out;
294 
295  out << "[";
296  if (N > 0)
297  {
298  for (int i = 0; i < N - 1; i++)
299  out << C[i] << ", ";
300  out << C[N-1];
301  }
302  out << "]" << std::endl;
303 
304  return out.str();
305  }
306 
308  NumberArray<Tout> const & coeffs() const
309  {
310  return C;
311  }
312 
318  inline static double node (int k, int N)
319  {
320  return cos(M_PI * (k + 0.5) / N);
321  }
322 
323 private:
324 
326  inline double scale(Tin x) const {
327  return (x - xt) / m;
328  }
329 
331  inline Tin unscale(double x) const {
332  return (xt + m*x);
333  }
334 
336  int N;
337 
340 
342  Tin xt;
343 
345  Tin m;
346 };
347 
348 //
349 // template member specializations
350 //
351 
357 template<> template <class Functor>
358 void Chebyshev<double,double>::generate (Functor const & f, int n, double a, double b)
359 {
360  N = n;
361  xt = 0.5 * (b + a);
362  m = 0.5 * (b - a);
363  C.resize(N);
364 
365  // input array
366  rArray fvals(N);
367 
368  // create the FFTW plan
369  fftw_plan plan = fftw_plan_r2r_1d (N, &fvals[0], &C[0], FFTW_REDFT10, 0);
370 
371  // evaluate nodes and function
372  double pi_over_N = M_PI / N;
373  for (int k = 0; k < N; k++)
374  {
375  double xk = cos(pi_over_N * (k + 0.5));
376  fvals[k] = f(unscale(xk));
377  }
378 
379  // execute the transform
380  fftw_execute(plan);
381  fftw_destroy_plan(plan);
382 
383  // normalize
384  double scal = 1./N;
385  for (double & c : C)
386  c *= scal;
387 }
388 
395 template<> template <class Functor>
396 void Chebyshev<double,Complex>::generate (Functor const & f, int n, double a, double b)
397 {
398  N = n;
399  xt = 0.5 * (b + a);
400  m = 0.5 * (b - a);
401  C.resize(N);
402 
403  // input array
404  cArray fvals(4*N), ftraf(4*N);
405 
406  // create the FFTW plan
407  fftw_plan plan = fftw_plan_dft_1d (
408  4*N,
409  reinterpret_cast<fftw_complex*>(&fvals[0]),
410  reinterpret_cast<fftw_complex*>(&ftraf[0]),
411  FFTW_FORWARD,
412  0
413  );
414 
415  // evaluate nodes and function
416  double pi_over_N = M_PI / N;
417  for (int k = 0; k < N; k++)
418  {
419  double xk = cos(pi_over_N * (k + 0.5));
420  fvals[2*k+1] = fvals[4*N-(2*k+1)] = f(unscale(xk));
421  }
422 
423  // execute the transform
424  fftw_execute(plan);
425  fftw_destroy_plan(plan);
426 
427  // copy normalized coefficients
428  double scal = 1./N;
429  for (int i = 0; i < N; i++)
430  C[i] = ftraf[i] * scal;
431 }
432 
433 #endif
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