cheb_approx.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2020, 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_tl
25 */
26 #ifndef O2SCL_CHEB_APPROX_H
27 #define O2SCL_CHEB_APPROX_H
28 
29 #include <cmath>
30 #include <o2scl/funct.h>
31 #include <o2scl/err_hnd.h>
32 
33 #include <boost/numeric/ublas/vector.hpp>
34 
35 #ifndef DOXYGEN_NO_O2NS
36 namespace o2scl {
37 #endif
38 
39  /** \brief Chebyshev approximation (GSL)
40 
41  Approximate a function on a finite interval using a Chebyshev series:
42  \f[
43  f(x) = \sum_n c_n T_n(x)
44  \f]
45  where \f$ T_n(x)=\cos(n \arccos x) \f$
46 
47  See also the \ref ex_cheb_sect .
48  */
49  template<class fp_t=double> class cheb_approx_tl {
50 
51  public:
52 
54 
55  protected:
56 
57  /// Coefficients
59  /// Order of the approximation
60  size_t order;
61  /// Lower end of the interval
62  fp_t a;
63  /// Upper end of the interval
64  fp_t b;
65  /// Single precision order
66  size_t order_sp;
67  /// Function evaluated at Chebyshev points
69  /// True if init has been called
71 
72  /// Desc
73  fp_t pi;
74 
75  public:
76 
77  cheb_approx_tl() {
78  init_called=false;
79  pi=boost::math::constants::pi<fp_t>();
80  }
81 
82  /// Copy constructor
84  order=gc.order;
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;
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, fp_t a1, fp_t 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  fp_t two=2;
136  fp_t half=1/two;
137  fp_t bma=(b-a)/two;
138  fp_t bpa=(b+a)/two;
139  fp_t fac=two/((fp_t)(order+1));
140 
141  for(k=0;k<=order;k++) {
142  fp_t y=cos(pi*(k+half)/(order+1));
143  f[k]=func(y*bma+bpa);
144  }
145 
146  for(j=0;j<=order;j++) {
147  fp_t sum=0.0;
148  for(k=0;k<=order;k++) {
149  sum+=f[k]*cos(pi*j*(k+half)/((fp_t)(order+1)));
150  }
151  c[j]=fac*sum;
152  }
153 
154  init_called=true;
155 
156  return;
157  }
158 
159  /// Create an approximation from a vector of coefficients
160  template<class vec_t> void init(fp_t a1, fp_t b1,
161  size_t ord, vec_t &v) {
162  order=ord;
163  order_sp=order;
164  a=a1;
165  b=b1;
166  c.resize(order+1);
167  for(size_t i=0;i<order+1;i++) c[i]=v[i];
168 
169  init_called=true;
170 
171  return;
172  }
173 
174  /// Create an approximation from a vector of function values
175  template<class vec_t> void init_func_values(fp_t a1, fp_t b1,
176  size_t ord, vec_t &fval) {
177  size_t k, j;
178 
179  if(a>=b) {
180  b=a1;
181  a=b1;
182  } else {
183  a=a1;
184  b=b1;
185  }
186  order=ord;
187  order_sp=ord;
188  c.resize(order+1);
189  f.resize(order+1);
190 
191  fp_t two=2;
192  fp_t half=1/two;
193  fp_t bma=(b-a)/two;
194  fp_t bpa=(b+a)/two;
195  fp_t fac=two/((fp_t)(order+1));
196 
197  for(j=0;j<=order;j++) {
198  fp_t sum=0.0;
199  for(k=0;k<=order; k++) {
200  sum+=fval[k]*cos(pi*j*(k+half)/((fp_t)(order+1)));
201  }
202  c[j]=fac*sum;
203  }
204 
205  init_called=true;
206 
207  return;
208  }
209  //@}
210 
211  /// \name Evaulation methods
212  //@{
213 
214  /** \brief Evaluate the approximation
215  */
216  fp_t eval(fp_t x) const {
217 
218  if (init_called==false) {
219  O2SCL_ERR("Series not initialized in cheb_approx::eval()",
221  return 0.0;
222  }
223 
224  size_t i;
225  fp_t d1 = 0.0;
226  fp_t d2 = 0.0;
227 
228  fp_t y = (2*x-a-b)/(b-a);
229  fp_t y2 = 2*y;
230 
231  for (i=order; i >= 1; i--) {
232  fp_t temp = d1;
233  d1 = y2*d1-d2+c[i];
234  d2 = temp;
235  }
236 
237  return y*d1-d2+c[0]/2;
238  }
239 
240  /** \brief Evaluate the approximation
241  */
242  fp_t operator()(fp_t x) const {
243  return eval(x);
244  }
245 
246  /** \brief Evaluate the approximation to a specified order
247  */
248  fp_t eval_n(size_t n, fp_t x) const {
249  size_t i;
250  fp_t d1 = 0.0;
251  fp_t d2 = 0.0;
252 
253  size_t eval_order;
254  if (n<order) eval_order=n;
255  else eval_order=order;
256 
257  fp_t y = (2*x-a-b)/(b-a);
258  fp_t y2 = 2*y;
259 
260  for (i = eval_order; i >= 1; i--) {
261  fp_t temp = d1;
262  d1 = y2*d1-d2+c[i];
263  d2 = temp;
264  }
265 
266  return y*d1-d2+c[0]/2;
267  }
268 
269  /** \brief Evaluate the approximation and give the uncertainty
270  */
271  void eval_err(fp_t x, fp_t &result, fp_t &abserr) {
272 
273  size_t i;
274  fp_t d1 = 0.0;
275  fp_t d2 = 0.0;
276 
277  fp_t y = (2*x-a-b)/(b-a);
278  fp_t y2 = 2*y;
279 
280  fp_t absc = 0.0;
281 
282  for (i = order; i >= 1; i--) {
283  fp_t temp = d1;
284  d1 = y2*d1-d2+c[i];
285  d2 = temp;
286  }
287 
288  result = y*d1-d2+c[0]/2;
289 
290  /* Estimate cumulative numerical error */
291 
292  for (i = 0; i <= order; i++) {
293  absc += o2scl::o2abs(c[i]);
294  }
295 
296  /* Combine truncation error and numerical error */
297 
298  fp_t dbl_eps=std::numeric_limits<fp_t>::epsilon();
299 
300  abserr = o2scl::o2abs(c[order])+absc*dbl_eps;
301 
302  return;
303  }
304 
305  /** \brief Evaluate the approximation to a specified order
306  and give the uncertainty
307  */
308  void eval_n_err(size_t n, fp_t x, fp_t &result, fp_t &abserr) {
309  size_t i;
310  fp_t d1 = 0.0;
311  fp_t d2 = 0.0;
312 
313  fp_t y = (2*x-a-b)/(b-a);
314  fp_t y2 = 2*y;
315 
316  fp_t absc = 0.0;
317 
318  size_t eval_order;
319  if (n<order) eval_order=n;
320  else eval_order=order;
321 
322  for (i = eval_order; i >= 1; i--) {
323  fp_t temp = d1;
324  d1 = y2*d1-d2+c[i];
325  d2 = temp;
326  }
327 
328  result = y*d1-d2+c[0]/2;
329 
330  /* Estimate cumulative numerical error */
331 
332  for (i = 0; i <= eval_order; i++) {
333  absc += o2scl::o2abs(c[i]);
334  }
335 
336  fp_t dbl_eps=std::numeric_limits<fp_t>::epsilon();
337 
338  /* Combine truncation error and numerical error */
339  abserr = o2scl::o2abs(c[eval_order])+absc*dbl_eps;
340 
341  return;
342  }
343  //@}
344 
345  /// \name Maniupulating coefficients and endpoints
346  //@{
347  /** \brief Get a coefficient
348 
349  Legal values of the argument are 0 to \c order (inclusive)
350  */
351  fp_t get_coefficient(size_t ix) const {
352  if (ix<order+1) {
353  return c[ix];
354  }
355  O2SCL_ERR
356  ("Requested invalid coefficient in cheb_approx::get_coefficient()",
358  return 0.0;
359  }
360 
361  /** \brief Set a coefficient
362 
363  Legal values of the argument are 0 to \c order (inclusive)
364  */
365  void set_coefficient(size_t ix, fp_t co) {
366  if (ix<order+1) {
367  c[ix]=co;
368  return;
369  }
370  O2SCL_ERR
371  ("Requested invalid coefficient in cheb_approx::get_coefficient()",
373  return;
374  }
375 
376  /// Return the endpoints of the approximation
377  void get_endpoints(fp_t &la, fp_t &lb) {
378  la=a;
379  lb=b;
380  return;
381  }
382 
383  /** \brief Get the coefficients
384  */
385  template<class vec_t> void get_coefficients(size_t n, vec_t &v) const {
386  for(size_t i=0;i<order+1 && i<n;i++) {
387  v[i]=c[i];
388  }
389  return;
390  }
391 
392  /** \brief Set the coefficients
393  */
394  template<class vec_t> void set_coefficients(size_t n, const vec_t &v) {
395  for(size_t i=0;i<order+1 && i<n;i++) {
396  c[i]=v[i];
397  }
398  return;
399  }
400  //@}
401 
402  /// \name Derivatives and integrals
403  //@{
404  /// Make \c gc an approximation to the derivative
405  void deriv(cheb_approx_tl &gc) const {
406 
407  size_t n=order+1;
408 
409  const fp_t con=2/(b-a);
410  size_t i;
411 
412  gc.init_called=true;
413  gc.order=order;
414  gc.order_sp=order;
415  gc.a=a;
416  gc.b=b;
417  gc.c.resize(n);
418  gc.f.resize(n);
419 
420  gc.c[n-1]=0.0;
421 
422  if (n > 1) {
423 
424  gc.c[n-2]=2*(n-1)*c[n-1];
425 
426  for(i=n;i>=3;i--) {
427  gc.c[i-3]=gc.c[i-1]+2*(i-2)*c[i-2];
428  }
429 
430  for(i=0;i<n;i++) {
431  gc.c[i]*=con;
432  }
433  }
434 
435  return;
436  }
437 
438  /// Make \c gc an approximation to the integral
439  void integ(cheb_approx_tl &gc) const {
440 
441  size_t n=order+1;
442 
443  const fp_t con=(b-a)/4;
444 
445  gc.init_called=true;
446  gc.order=order;
447  gc.order_sp=order;
448  gc.a=a;
449  gc.b=b;
450  gc.c.resize(n);
451  gc.f.resize(n);
452 
453  if (n == 1) {
454 
455  gc.c[0]=0.0;
456 
457  } else if (n == 2) {
458 
459  gc.c[1]=con*c[0];
460  gc.c[0]=2*gc.c[1];
461 
462  } else {
463 
464  fp_t sum=0;
465  fp_t fac=1;
466 
467  for(size_t i=1;i<=n-2;i++) {
468  gc.c[i]=con*(c[i-1]-c[i+1])/((fp_t)i);
469  sum += fac*gc.c[i];
470  fac=-fac;
471  }
472  gc.c[n-1]=con*c[n-2]/(n-1);
473  sum += fac*gc.c[n-1];
474  gc.c[0]=2*sum;
475  }
476 
477  return;
478  }
479  //@}
480 
481  };
482 
483  /** \brief Double-precision version of \ref o2scl::cheb_approx_tl
484  */
486 
487 #ifndef DOXYGEN_NO_O2NS
488 }
489 #endif
490 
491 #endif
o2scl::cheb_approx_tl::b
fp_t b
Upper end of the interval.
Definition: cheb_approx.h:64
o2scl::cheb_approx_tl::set_coefficient
void set_coefficient(size_t ix, fp_t co)
Set a coefficient.
Definition: cheb_approx.h:365
boost::numeric::ublas::vector
The default vector type from uBlas.
Definition: vector.h:3735
o2scl::cheb_approx_tl::eval
fp_t eval(fp_t x) const
Evaluate the approximation.
Definition: cheb_approx.h:216
o2scl::cheb_approx_tl::cheb_approx_tl
cheb_approx_tl(const cheb_approx_tl &gc)
Copy constructor.
Definition: cheb_approx.h:83
o2scl::cheb_approx_tl::get_coefficients
void get_coefficients(size_t n, vec_t &v) const
Get the coefficients.
Definition: cheb_approx.h:385
o2scl::cheb_approx_tl::deriv
void deriv(cheb_approx_tl &gc) const
Make gc an approximation to the derivative.
Definition: cheb_approx.h:405
o2scl
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
Definition: anneal.h:42
o2scl::cheb_approx_tl::operator()
fp_t operator()(fp_t x) const
Evaluate the approximation.
Definition: cheb_approx.h:242
o2scl::cheb_approx_tl::eval_n
fp_t eval_n(size_t n, fp_t x) const
Evaluate the approximation to a specified order.
Definition: cheb_approx.h:248
o2scl::cheb_approx_tl::eval_n_err
void eval_n_err(size_t n, fp_t x, fp_t &result, fp_t &abserr)
Evaluate the approximation to a specified order and give the uncertainty.
Definition: cheb_approx.h:308
o2scl::cheb_approx_tl::get_coefficient
fp_t get_coefficient(size_t ix) const
Get a coefficient.
Definition: cheb_approx.h:351
o2scl::cheb_approx_tl::operator=
cheb_approx_tl & operator=(const cheb_approx_tl &gc)
Copy constructor.
Definition: cheb_approx.h:94
o2scl::cheb_approx_tl::set_coefficients
void set_coefficients(size_t n, const vec_t &v)
Set the coefficients.
Definition: cheb_approx.h:394
o2scl::exc_einval
@ exc_einval
invalid argument supplied by user
Definition: err_hnd.h:59
o2scl::cheb_approx_tl::pi
fp_t pi
Desc.
Definition: cheb_approx.h:73
o2scl::cheb_approx_tl::init
void init(fp_t a1, fp_t b1, size_t ord, vec_t &v)
Create an approximation from a vector of coefficients.
Definition: cheb_approx.h:160
o2scl::cheb_approx_tl::c
ubvector c
Coefficients.
Definition: cheb_approx.h:58
o2scl::cheb_approx_tl::order_sp
size_t order_sp
Single precision order.
Definition: cheb_approx.h:66
o2scl::cheb_approx_tl::init
void init(func_t &func, size_t ord, fp_t a1, fp_t b1)
Initialize a Chebyshev approximation of the function func over the interval from a1 to b1.
Definition: cheb_approx.h:120
o2scl::cheb_approx_tl::a
fp_t a
Lower end of the interval.
Definition: cheb_approx.h:62
O2SCL_ERR
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
o2scl::cheb_approx_tl::init_called
bool init_called
True if init has been called.
Definition: cheb_approx.h:70
o2scl::o2abs
float o2abs(const float x)
Absolute value for single precision numbers.
o2scl::cheb_approx_tl::integ
void integ(cheb_approx_tl &gc) const
Make gc an approximation to the integral.
Definition: cheb_approx.h:439
o2scl::cheb_approx
cheb_approx_tl< double > cheb_approx
Double-precision version of o2scl::cheb_approx_tl.
Definition: cheb_approx.h:485
o2scl::cheb_approx_tl::order
size_t order
Order of the approximation.
Definition: cheb_approx.h:60
o2scl::cheb_approx_tl
Chebyshev approximation (GSL)
Definition: cheb_approx.h:49
o2scl::cheb_approx_tl::get_endpoints
void get_endpoints(fp_t &la, fp_t &lb)
Return the endpoints of the approximation.
Definition: cheb_approx.h:377
o2scl::cheb_approx_tl::eval_err
void eval_err(fp_t x, fp_t &result, fp_t &abserr)
Evaluate the approximation and give the uncertainty.
Definition: cheb_approx.h:271
o2scl::cheb_approx_tl::init_func_values
void init_func_values(fp_t a1, fp_t b1, size_t ord, vec_t &fval)
Create an approximation from a vector of function values.
Definition: cheb_approx.h:175
o2scl::cheb_approx_tl::f
ubvector f
Function evaluated at Chebyshev points.
Definition: cheb_approx.h:68

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