cheb_approx.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 /** \file cheb_approx.h
24  \brief File for definition of \ref o2scl::cheb_approx
25 */
26 #ifndef O2SCL_CHEB_APPROX_H
27 #define O2SCL_CHEB_APPROX_H
28 
29 #include <cmath>
30 #include <gsl/gsl_chebyshev.h>
31 #include <o2scl/funct.h>
32 #include <o2scl/err_hnd.h>
33 
34 #include <boost/numeric/ublas/vector.hpp>
35 
36 #ifndef DOXYGEN_NO_O2NS
37 namespace o2scl {
38 #endif
39 
40  /** \brief Chebyshev approximation (GSL)
41 
42  Approximate a function on a finite interval using a Chebyshev series:
43  \f[
44  f(x) = \sum_n c_n T_n(x)
45  \f]
46  where \f$ T_n(x)=\cos(n \arccos x) \f$
47 
48  See also the \ref ex_cheb_sect .
49 
50  \future Move non-template functions to .cpp file.
51  */
52  class cheb_approx {
53  //: public funct {
54 
55  public:
56 
58 
59  protected:
60 
61  /// Coefficients
62  ubvector c;
63  /// Order of the approximation
64  size_t order;
65  /// Lower end of the interval
66  double a;
67  /// Upper end of the interval
68  double b;
69  /// Single precision order
70  size_t order_sp;
71  /// Function evaluated at Chebyshev points
72  ubvector f;
73  /// True if init has been called
75 
76  public:
77 
78  cheb_approx() {
79  init_called=false;
80  }
81 
82  /// Copy constructor
83  cheb_approx(const cheb_approx &gc) {
84  order=gc.order;
85  init_called=gc.init_called;
86  c=gc.c;
87  f=gc.f;
88  a=gc.a;
89  b=gc.b;
90  order_sp=gc.order_sp;
91  }
92 
93  /// Copy constructor
95 
96  // Check for self-assignment so that we don't
97  // reallocate the vectors and lose the coefficients
98  if (this==&gc) return *this;
99 
100  order=gc.order;
101  init_called=gc.init_called;
102  c=gc.c;
103  f=gc.f;
104  a=gc.a;
105  b=gc.b;
106  order_sp=gc.order_sp;
107 
108  return *this;
109  }
110 
111  /// \name Initialization methods
112  //@{
113  /** \brief Initialize a Chebyshev approximation of the function
114  \c func over the interval from \c a1 to \c b1
115 
116  The interval must be specified so that \f$ a < b \f$ , so
117  a and b are swapped if this is not the case.
118  */
119  template<class func_t>
120  void init(func_t &func, size_t ord, double a1, double b1) {
121  size_t k, j;
122 
123  if(a1>=b1) {
124  b=a1;
125  a=b1;
126  } else {
127  a=a1;
128  b=b1;
129  }
130  order=ord;
131  order_sp=ord;
132  c.resize(order+1);
133  f.resize(order+1);
134 
135  double bma=0.5*(b-a);
136  double bpa=0.5*(b+a);
137  double fac=2.0/(order+1.0);
138 
139  for(k=0;k<=order;k++) {
140  double y=cos(M_PI*(k+0.5)/(order+1));
141  f[k]=func(y*bma+bpa);
142  }
143 
144  for(j=0;j<=order;j++) {
145  double sum=0.0;
146  for(k=0;k<=order; k++) {
147  sum+=f[k]*cos(M_PI*j*(k+0.5)/(order+1));
148  }
149  c[j]=fac*sum;
150  }
151 
152  init_called=true;
153 
154  return;
155  }
156 
157  /// Create an approximation from a vector of coefficients
158  template<class vec_t> void init(double a1, double b1,
159  size_t ord, vec_t &v) {
160  order=ord;
161  order_sp=order;
162  a=a1;
163  b=b1;
164  c.resize(order+1);
165  for(size_t i=0;i<order+1;i++) c[i]=v[i];
166 
167  init_called=true;
168 
169  return;
170  }
171 
172  /// Create an approximation from a vector of function values
173  template<class vec_t> void init_func_values(double a1, double b1,
174  size_t ord, vec_t &fval) {
175  size_t k, j;
176 
177  if(a>=b) {
178  b=a1;
179  a=b1;
180  } else {
181  a=a1;
182  b=b1;
183  }
184  order=ord;
185  order_sp=ord;
186  c.resize(order+1);
187  f.resize(order+1);
188 
189  double bma=0.5*(b-a);
190  double bpa=0.5*(b+a);
191  double fac=2.0/(order+1.0);
192 
193  for(j=0;j<=order;j++) {
194  double sum=0.0;
195  for(k=0;k<=order; k++) {
196  sum+=fval[k]*cos(M_PI*j*(k+0.5)/(order+1));
197  }
198  c[j]=fac*sum;
199  }
200 
201  init_called=true;
202 
203  return;
204  }
205  //@}
206 
207  /// \name Evaulation methods
208  //@{
209 
210  /** \brief Evaluate the approximation
211  */
212  double eval(double x) const {
213 
214  if (init_called==false) {
215  O2SCL_ERR("Series not initialized in cheb_approx::eval()",
217  return 0.0;
218  }
219 
220  size_t i;
221  double d1 = 0.0;
222  double d2 = 0.0;
223 
224  double y = (2.0*x - a - b) / (b - a);
225  double y2 = 2.0*y;
226 
227  for (i=order; i >= 1; i--) {
228  double temp = d1;
229  d1 = y2*d1 - d2 + c[i];
230  d2 = temp;
231  }
232 
233  return y*d1 - d2 + 0.5*c[0];
234  }
235 
236  /** \brief Evaluate the approximation
237  */
238  double operator()(double x) const {
239  return eval(x);
240  }
241 
242  /** \brief Evaluate the approximation to a specified order
243  */
244  double eval_n(size_t n, double x) const {
245  size_t i;
246  double d1 = 0.0;
247  double d2 = 0.0;
248 
249  size_t eval_order;
250  if (n<order) eval_order=n;
251  else eval_order=order;
252 
253  double y = (2.0*x - a - b) / (b - a);
254  double y2 = 2.0*y;
255 
256  for (i = eval_order; i >= 1; i--) {
257  double temp = d1;
258  d1 = y2*d1 - d2 + c[i];
259  d2 = temp;
260  }
261 
262  return y*d1 - d2 + 0.5*c[0];
263  }
264 
265  /** \brief Evaluate the approximation and give the uncertainty
266  */
267  void eval_err(double x, double &result, double &abserr) {
268 
269  size_t i;
270  double d1 = 0.0;
271  double d2 = 0.0;
272 
273  double y = (2.*x - a - b) / (b - a);
274  double y2 = 2.0*y;
275 
276  double absc = 0.0;
277 
278  for (i = order; i >= 1; i--) {
279  double temp = d1;
280  d1 = y2*d1 - d2 + c[i];
281  d2 = temp;
282  }
283 
284  result = y*d1 - d2 + 0.5*c[0];
285 
286  /* Estimate cumulative numerical error */
287 
288  for (i = 0; i <= order; i++) {
289  absc += fabs(c[i]);
290  }
291 
292  /* Combine truncation error and numerical error */
293 
294  double dbl_eps=std::numeric_limits<double>::epsilon();
295 
296  abserr = fabs (c[order]) + absc*dbl_eps;
297 
298  return;
299  }
300 
301  /** \brief Evaluate the approximation to a specified order
302  and give the uncertainty
303  */
304  void eval_n_err(size_t n, double x, double &result, double &abserr) {
305  size_t i;
306  double d1 = 0.0;
307  double d2 = 0.0;
308 
309  double y = (2.*x - a - b) / (b - a);
310  double y2 = 2.0*y;
311 
312  double absc = 0.0;
313 
314  size_t eval_order;
315  if (n<order) eval_order=n;
316  else eval_order=order;
317 
318  for (i = eval_order; i >= 1; i--) {
319  double temp = d1;
320  d1 = y2*d1 - d2 + c[i];
321  d2 = temp;
322  }
323 
324  result = y*d1 - d2 + 0.5*c[0];
325 
326  /* Estimate cumulative numerical error */
327 
328  for (i = 0; i <= eval_order; i++) {
329  absc += fabs(c[i]);
330  }
331 
332  double dbl_eps=std::numeric_limits<double>::epsilon();
333 
334  /* Combine truncation error and numerical error */
335  abserr = fabs(c[eval_order])+absc*dbl_eps;
336 
337  return;
338  }
339  //@}
340 
341  /// \name Maniupulating coefficients and endpoints
342  //@{
343  /** \brief Get a coefficient
344 
345  Legal values of the argument are 0 to \c order (inclusive)
346  */
347  double get_coefficient(size_t ix) const {
348  if (ix<order+1) {
349  return c[ix];
350  }
351  O2SCL_ERR
352  ("Requested invalid coefficient in cheb_approx::get_coefficient()",
354  return 0.0;
355  }
356 
357  /** \brief Set a coefficient
358 
359  Legal values of the argument are 0 to \c order (inclusive)
360  */
361  void set_coefficient(size_t ix, double co) {
362  if (ix<order+1) {
363  c[ix]=co;
364  return;
365  }
366  O2SCL_ERR
367  ("Requested invalid coefficient in cheb_approx::get_coefficient()",
369  return;
370  }
371 
372  /// Return the endpoints of the approximation
373  void get_endpoints(double &la, double &lb) {
374  la=a;
375  lb=b;
376  return;
377  }
378 
379  /** \brief Get the coefficients
380  */
381  template<class vec_t> void get_coefficients(size_t n, vec_t &v) const {
382  for(size_t i=0;i<order+1 && i<n;i++) {
383  v[i]=c[i];
384  }
385  return;
386  }
387 
388  /** \brief Set the coefficients
389  */
390  template<class vec_t> void set_coefficients(size_t n, const vec_t &v) {
391  for(size_t i=0;i<order+1 && i<n;i++) {
392  c[i]=v[i];
393  }
394  return;
395  }
396  //@}
397 
398  /// \name Derivatives and integrals
399  //@{
400  /// Make \c gc an approximation to the derivative
401  void deriv(cheb_approx &gc) const {
402 
403  size_t n=order+1;
404 
405  const double con=2.0/(b-a);
406  size_t i;
407 
408  gc.init_called=true;
409  gc.order=order;
410  gc.order_sp=order;
411  gc.a=a;
412  gc.b=b;
413  gc.c.resize(n);
414  gc.f.resize(n);
415 
416  gc.c[n-1]=0.0;
417 
418  if (n > 1) {
419 
420  gc.c[n-2]=2.0*(n-1.0)*c[n-1];
421 
422  for(i=n;i>=3;i--) {
423  gc.c[i-3]=gc.c[i-1]+2.0*(i-2.0)*c[i-2];
424  }
425 
426  for(i=0;i<n;i++) {
427  gc.c[i]*=con;
428  }
429  }
430 
431  return;
432  }
433 
434  /// Make \c gc an approximation to the integral
435  void integ(cheb_approx &gc) const {
436 
437  size_t n=order+1;
438 
439  const double con=0.25*(b-a);
440 
441  gc.init_called=true;
442  gc.order=order;
443  gc.order_sp=order;
444  gc.a=a;
445  gc.b=b;
446  gc.c.resize(n);
447  gc.f.resize(n);
448 
449  if (n == 1) {
450 
451  gc.c[0]=0.0;
452 
453  } else if (n == 2) {
454 
455  gc.c[1]=con*c[0];
456  gc.c[0]=2.0*gc.c[1];
457 
458  } else {
459 
460  double sum=0.0;
461  double fac=1.0;
462 
463  for(size_t i=1; i<=n-2; i++) {
464  gc.c[i]=con*(c[i-1] - c[i+1])/((double)i);
465  sum += fac*gc.c[i];
466  fac=-fac;
467  }
468  gc.c[n-1]=con*c[n-2]/(n-1.0);
469  sum += fac*gc.c[n-1];
470  gc.c[0]=2.0*sum;
471  }
472 
473  return;
474  }
475  //@}
476 
477  };
478 
479 #ifndef DOXYGEN_NO_O2NS
480 }
481 #endif
482 
483 #endif
size_t order_sp
Single precision order.
Definition: cheb_approx.h:70
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
void eval_n_err(size_t n, double x, double &result, double &abserr)
Evaluate the approximation to a specified order and give the uncertainty.
Definition: cheb_approx.h:304
double a
Lower end of the interval.
Definition: cheb_approx.h:66
void integ(cheb_approx &gc) const
Make gc an approximation to the integral.
Definition: cheb_approx.h:435
double operator()(double x) const
Evaluate the approximation.
Definition: cheb_approx.h:238
invalid argument supplied by user
Definition: err_hnd.h:59
void get_endpoints(double &la, double &lb)
Return the endpoints of the approximation.
Definition: cheb_approx.h:373
void get_coefficients(size_t n, vec_t &v) const
Get the coefficients.
Definition: cheb_approx.h:381
void init(func_t &func, size_t ord, double a1, double b1)
Initialize a Chebyshev approximation of the function func over the interval from a1 to b1...
Definition: cheb_approx.h:120
double eval_n(size_t n, double x) const
Evaluate the approximation to a specified order.
Definition: cheb_approx.h:244
double b
Upper end of the interval.
Definition: cheb_approx.h:68
bool init_called
True if init has been called.
Definition: cheb_approx.h:74
size_t order
Order of the approximation.
Definition: cheb_approx.h:64
cheb_approx(const cheb_approx &gc)
Copy constructor.
Definition: cheb_approx.h:83
ubvector c
Coefficients.
Definition: cheb_approx.h:62
void set_coefficient(size_t ix, double co)
Set a coefficient.
Definition: cheb_approx.h:361
double get_coefficient(size_t ix) const
Get a coefficient.
Definition: cheb_approx.h:347
void eval_err(double x, double &result, double &abserr)
Evaluate the approximation and give the uncertainty.
Definition: cheb_approx.h:267
void set_coefficients(size_t n, const vec_t &v)
Set the coefficients.
Definition: cheb_approx.h:390
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
void init(double a1, double b1, size_t ord, vec_t &v)
Create an approximation from a vector of coefficients.
Definition: cheb_approx.h:158
Chebyshev approximation (GSL)
Definition: cheb_approx.h:52
double eval(double x) const
Evaluate the approximation.
Definition: cheb_approx.h:212
cheb_approx & operator=(const cheb_approx &gc)
Copy constructor.
Definition: cheb_approx.h:94
void init_func_values(double a1, double b1, size_t ord, vec_t &fval)
Create an approximation from a vector of function values.
Definition: cheb_approx.h:173
ubvector f
Function evaluated at Chebyshev points.
Definition: cheb_approx.h:72
void deriv(cheb_approx &gc) const
Make gc an approximation to the derivative.
Definition: cheb_approx.h:401

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).