|
Packit |
67cb25 |
/* cheb/test.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_test.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_ieee_utils.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_chebyshev.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_T0 (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return 1.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_T1 (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return x;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_T2 (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return 2*x*x - 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_sin (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return sin(x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_DP (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return 2.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_P (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return 2.0*x + 3.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_IP1 (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return 30*(x+5.0)/10.0; /* first order approximation to integral over -5,5 */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double f_IP2 (double x, void * p) {
|
|
Packit |
67cb25 |
p = 0;
|
|
Packit |
67cb25 |
return x*x + 3*x - 10.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
test_dim (const size_t n, const double a, const double b,
|
|
Packit |
67cb25 |
gsl_function * F, gsl_function * DF, gsl_function *IF)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tol = 100.0 * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
double x;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_series * cs = gsl_cheb_alloc(n);
|
|
Packit |
67cb25 |
gsl_cheb_series * csd = gsl_cheb_alloc(n);
|
|
Packit |
67cb25 |
gsl_cheb_series * csi = gsl_cheb_alloc(n);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_init(cs, F, a, b);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=a; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(cs, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, GSL_FN_EVAL(F, x), tol, "gsl_cheb_eval, F(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test derivative */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_calc_deriv(csd, cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=a; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(csd, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, GSL_FN_EVAL(DF, x), tol, "gsl_cheb_eval, deriv F(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test integral */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_calc_integ(csi, cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=a; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(csi, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, GSL_FN_EVAL(IF, x), tol, "gsl_cheb_eval, integ F(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_free(csi);
|
|
Packit |
67cb25 |
gsl_cheb_free(csd);
|
|
Packit |
67cb25 |
gsl_cheb_free(cs);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main(void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double tol = 100.0 * GSL_DBL_EPSILON;
|
|
Packit |
67cb25 |
double ftol = 20.0;
|
|
Packit |
67cb25 |
double x;
|
|
Packit |
67cb25 |
size_t i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_series * cs = gsl_cheb_alloc(40);
|
|
Packit |
67cb25 |
gsl_cheb_series * csd = gsl_cheb_alloc(40);
|
|
Packit |
67cb25 |
gsl_cheb_series * csi = gsl_cheb_alloc(40);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_function F_sin, F_T0, F_T1, F_T2, F_DP, F_P, F_IP1, F_IP2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_sin.function = f_sin;
|
|
Packit |
67cb25 |
F_sin.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_T0.function = f_T0;
|
|
Packit |
67cb25 |
F_T0.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_T1.function = f_T1;
|
|
Packit |
67cb25 |
F_T1.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_T2.function = f_T2;
|
|
Packit |
67cb25 |
F_T2.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_P.function = f_P;
|
|
Packit |
67cb25 |
F_P.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_DP.function = f_DP;
|
|
Packit |
67cb25 |
F_DP.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_IP1.function = f_IP1;
|
|
Packit |
67cb25 |
F_IP1.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
F_IP2.function = f_IP2;
|
|
Packit |
67cb25 |
F_IP2.params = 0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_ieee_env_setup();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_init(cs, &F_T0, -1.0, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
size_t expected = 40;
|
|
Packit |
67cb25 |
size_t order = gsl_cheb_order (cs);
|
|
Packit |
67cb25 |
size_t size = gsl_cheb_size (cs);
|
|
Packit |
67cb25 |
double * p = gsl_cheb_coeffs (cs);
|
|
Packit |
67cb25 |
gsl_test(order != expected, "gsl_cheb_order");
|
|
Packit |
67cb25 |
gsl_test(size != expected + 1, "gsl_cheb_size");
|
|
Packit |
67cb25 |
gsl_test(p != cs->c, "gsl_cheb_coeffs");
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i<cs->order; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c_exp = (i == 0) ? 2.0 : 0.0;
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_0(x)", i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_init(cs, &F_T1, -1.0, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i<cs->order; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c_exp = (i == 1) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_1(x)", i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_init(cs, &F_T2, -1.0, 1.0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i<cs->order; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double c_exp = (i == 2) ? 1.0 : 0.0;
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[i], c_exp, tol, "c[%d] for T_2(x)", i);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_init(cs, &F_sin, -M_PI, M_PI);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[0], 0.0, tol, "c[0] for F_sin(x)");
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[1], 5.69230686359506e-01, tol, "c[1] for F_sin(x)");
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[2], 0.0, tol, "c[2] for F_sin(x)");
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[3], -6.66916672405979e-01, tol, "c[3] for F_sin(x)");
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[4], 0.0, tol, "c[4] for F_sin(x)");
|
|
Packit |
67cb25 |
gsl_test_abs (cs->c[5], 1.04282368734237e-01, tol, "c[5] for F_sin(x)");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(cs, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval, sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_err(cs, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_err, sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_err, error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval_n(cs, 25, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, sin(x), tol, "gsl_cheb_eval_n, sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_n_err(cs, 25, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, sin(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-sin(x)) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_n_err, error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test derivative */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_calc_deriv(csd, cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(csd, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval, deriv sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef TEST_DERIVATIVE_ERR
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_err(csd, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, cos(x), tol, "gsl_cheb_eval_err, deriv sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_err, deriv error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval_n(csd, 25, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, cos(x), 1600 * tol, "gsl_cheb_eval_n, deriv sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef TEST_DERIVATIVE_ERR
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_n_err(csd, 25, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, cos(x), 100.0 * tol, "gsl_cheb_eval_n_err, deriv sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-cos(x)) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_n_err, deriv error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test integral */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_calc_integ(csi, cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval(csi, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval, integ sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef TEST_INTEGRAL_ERR
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_err(csi, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_err, integ sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_err, integ error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r = gsl_cheb_eval_n(csi, 25, x);
|
|
Packit |
67cb25 |
gsl_test_abs(r, -(1+cos(x)), tol, "gsl_cheb_eval_n, integ sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#ifdef TEST_INTEGRAL_ERR
|
|
Packit |
67cb25 |
for(x=-M_PI; x
|
|
Packit |
67cb25 |
double r, e;
|
|
Packit |
67cb25 |
gsl_cheb_eval_n_err(csi, 25, x, &r, &e);
|
|
Packit |
67cb25 |
gsl_test_abs(r, -(1+cos(x)), 100.0 * tol, "gsl_cheb_eval_n_err, integ sin(%.3g)", x);
|
|
Packit |
67cb25 |
gsl_test_factor(fabs(r-(-1-cos(x))) + GSL_DBL_EPSILON, e, ftol,
|
|
Packit |
67cb25 |
"gsl_cheb_eval_n_err, integ error sin(%.3g)", x);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
#endif
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_cheb_free(csi);
|
|
Packit |
67cb25 |
gsl_cheb_free(csd);
|
|
Packit |
67cb25 |
gsl_cheb_free(cs);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Test low order cases */
|
|
Packit |
67cb25 |
test_dim (2, -5.0, 5.0, &F_P, &F_DP, &F_IP2);
|
|
Packit |
67cb25 |
test_dim (1, -5.0, 5.0, &F_P, &F_DP, &F_IP1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
exit (gsl_test_summary());
|
|
Packit |
67cb25 |
}
|