26 #ifndef O2SCL_CHEB_APPROX_H 27 #define O2SCL_CHEB_APPROX_H 30 #include <gsl/gsl_chebyshev.h> 31 #include <o2scl/funct.h> 32 #include <o2scl/err_hnd.h> 34 #include <boost/numeric/ublas/vector.hpp> 36 #ifndef DOXYGEN_NO_O2NS 98 if (
this==&gc)
return *
this;
119 template<
class func_t>
120 void init(func_t &func,
size_t ord,
double a1,
double b1) {
135 double bma=0.5*(b-
a);
136 double bpa=0.5*(b+
a);
137 double fac=2.0/(order+1.0);
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);
144 for(j=0;j<=
order;j++) {
146 for(k=0;k<=
order; k++) {
147 sum+=f[k]*cos(M_PI*j*(k+0.5)/(order+1));
158 template<
class vec_t>
void init(
double a1,
double b1,
159 size_t ord, vec_t &v) {
165 for(
size_t i=0;i<order+1;i++) c[i]=v[i];
174 size_t ord, vec_t &fval) {
189 double bma=0.5*(b-
a);
190 double bpa=0.5*(b+
a);
191 double fac=2.0/(order+1.0);
193 for(j=0;j<=
order;j++) {
195 for(k=0;k<=
order; k++) {
196 sum+=fval[k]*cos(M_PI*j*(k+0.5)/(order+1));
214 if (init_called==
false) {
215 O2SCL_ERR(
"Series not initialized in cheb_approx::eval()",
224 double y = (2.0*x - a -
b) / (b - a);
227 for (i=order; i >= 1; i--) {
229 d1 = y2*d1 - d2 + c[i];
233 return y*d1 - d2 + 0.5*c[0];
244 double eval_n(
size_t n,
double x)
const {
250 if (n<order) eval_order=n;
251 else eval_order=
order;
253 double y = (2.0*x - a -
b) / (b - a);
256 for (i = eval_order; i >= 1; i--) {
258 d1 = y2*d1 - d2 + c[i];
262 return y*d1 - d2 + 0.5*c[0];
267 void eval_err(
double x,
double &result,
double &abserr) {
273 double y = (2.*x - a -
b) / (b - a);
278 for (i = order; i >= 1; i--) {
280 d1 = y2*d1 - d2 + c[i];
284 result = y*d1 - d2 + 0.5*c[0];
288 for (i = 0; i <=
order; i++) {
294 double dbl_eps=std::numeric_limits<double>::epsilon();
296 abserr = fabs (c[order]) + absc*dbl_eps;
304 void eval_n_err(
size_t n,
double x,
double &result,
double &abserr) {
309 double y = (2.*x - a -
b) / (b - a);
315 if (n<order) eval_order=n;
316 else eval_order=
order;
318 for (i = eval_order; i >= 1; i--) {
320 d1 = y2*d1 - d2 + c[i];
324 result = y*d1 - d2 + 0.5*c[0];
328 for (i = 0; i <= eval_order; i++) {
332 double dbl_eps=std::numeric_limits<double>::epsilon();
335 abserr = fabs(c[eval_order])+absc*dbl_eps;
352 (
"Requested invalid coefficient in cheb_approx::get_coefficient()",
367 (
"Requested invalid coefficient in cheb_approx::get_coefficient()",
382 for(
size_t i=0;i<order+1 && i<n;i++) {
391 for(
size_t i=0;i<order+1 && i<n;i++) {
405 const double con=2.0/(b-
a);
420 gc.
c[n-2]=2.0*(n-1.0)*c[n-1];
423 gc.
c[i-3]=gc.
c[i-1]+2.0*(i-2.0)*c[i-2];
439 const double con=0.25*(b-
a);
463 for(
size_t i=1; i<=n-2; i++) {
464 gc.
c[i]=con*(c[i-1] - c[i+1])/((
double)i);
468 gc.
c[n-1]=con*c[n-2]/(n-1.0);
469 sum += fac*gc.
c[n-1];
479 #ifndef DOXYGEN_NO_O2NS size_t order_sp
Single precision order.
The main O<span style='position: relative; top: 0.3em; font-size: 0.8em'>2</span>scl O$_2$scl names...
void eval_n_err(size_t n, double x, double &result, double &abserr)
Evaluate the approximation to a specified order and give the uncertainty.
double a
Lower end of the interval.
void integ(cheb_approx &gc) const
Make gc an approximation to the integral.
double operator()(double x) const
Evaluate the approximation.
invalid argument supplied by user
void get_endpoints(double &la, double &lb)
Return the endpoints of the approximation.
void get_coefficients(size_t n, vec_t &v) const
Get the coefficients.
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...
double eval_n(size_t n, double x) const
Evaluate the approximation to a specified order.
double b
Upper end of the interval.
bool init_called
True if init has been called.
size_t order
Order of the approximation.
cheb_approx(const cheb_approx &gc)
Copy constructor.
void set_coefficient(size_t ix, double co)
Set a coefficient.
double get_coefficient(size_t ix) const
Get a coefficient.
void eval_err(double x, double &result, double &abserr)
Evaluate the approximation and give the uncertainty.
void set_coefficients(size_t n, const vec_t &v)
Set the coefficients.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
void init(double a1, double b1, size_t ord, vec_t &v)
Create an approximation from a vector of coefficients.
Chebyshev approximation (GSL)
double eval(double x) const
Evaluate the approximation.
cheb_approx & operator=(const cheb_approx &gc)
Copy constructor.
void init_func_values(double a1, double b1, size_t ord, vec_t &fval)
Create an approximation from a vector of function values.
ubvector f
Function evaluated at Chebyshev points.
void deriv(cheb_approx &gc) const
Make gc an approximation to the derivative.