| /* integration/tests.c |
| * |
| * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or (at |
| * your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, but |
| * WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| * General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, write to the Free Software |
| * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. |
| */ |
| |
| #include <config.h> |
| #include <math.h> |
| #include <gsl/gsl_math.h> |
| |
| #include "tests.h" |
| |
| /* These are the test functions from table 4.1 of the QUADPACK book */ |
| |
| /* f1(x) = x^alpha * log(1/x) */ |
| /* integ(f1,x,0,1) = 1/(alpha + 1)^2 */ |
| |
| double f1 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(x,alpha) * log(1/x) ; |
| } |
| |
| /* f2(x) = 4^-alpha / ((x-pi/4)^2 + 16^-alpha) */ |
| /* integ(f2,x,0,1) = arctan((4-pi)4^(alpha-1)) + arctan(pi 4^(alpha-1)) */ |
| |
| double f2 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(4.0,-alpha) / (pow((x-M_PI/4.0),2.0) + pow(16.0,-alpha)) ; |
| } |
| |
| /* f3(x) = cos(2^alpha * sin(x)) */ |
| /* integ(f3,x,0,pi) = pi J_0(2^alpha) */ |
| |
| double f3 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return cos(pow(2.0,alpha) * sin(x)) ; |
| } |
| |
| /* Functions 4, 5 and 6 are duplicates of functions 1, 2 and 3 */ |
| /* .... */ |
| |
| /* f7(x) = |x - 1/3|^alpha */ |
| /* integ(f7,x,0,1) = ((2/3)^(alpha+1) + (1/3)^(alpha+1))/(alpha + 1) */ |
| |
| double f7 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(fabs(x - (1.0/3.0)),alpha) ; |
| } |
| |
| /* f8(x) = |x - pi/4|^alpha */ |
| /* integ(f8,x,0,1) = |
| ((1 - pi/4)^(alpha+1) + (pi/4)^(alpha+1))/(alpha + 1) */ |
| |
| double f8 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(fabs(x - (M_PI/4.0)),alpha) ; |
| } |
| |
| /* f9(x) = sqrt(1 - x^2) / (x + 1 + 2^-alpha) */ |
| /* integ(f9,x,-1,1) = pi/sqrt((1+2^-alpha)^2-1) */ |
| |
| double f9 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return 1 / ((x + 1 + pow(2.0,-alpha)) * sqrt(1-x*x)) ; |
| } |
| |
| /* f10(x) = sin(x)^(alpha - 1) */ |
| /* integ(f10,x,0,pi/2) = 2^(alpha-2) ((Gamma(alpha/2))^2)/Gamma(alpha) */ |
| |
| double f10 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(sin(x), alpha-1) ; |
| } |
| |
| /* f11(x) = log(1/x)^(alpha - 1) */ |
| /* integ(f11,x,0,1) = Gamma(alpha) */ |
| |
| double f11 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(log(1/x), alpha-1) ; |
| } |
| |
| /* f12(x) = exp(20*(x-1)) * sin(2^alpha * x) */ |
| /* integ(f12,x,0,1) = |
| (20 sin(2^alpha) - 2^alpha cos(2^alpha) + 2^alpha exp(-20)) |
| /(400 + 4^alpha) */ |
| |
| double f12 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return exp(20*(x-1)) * sin(pow(2.0,alpha) * x) ; |
| } |
| |
| /* f13(x) = cos(2^alpha * x)/sqrt(x(1 - x)) */ |
| /* integ(f13,x,0,1) = pi cos(2^(alpha-1)) J_0(2^(alpha-1)) */ |
| |
| double f13 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return cos(pow(2.0,alpha)*x)/sqrt(x*(1-x)) ; |
| } |
| |
| double f14 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return exp(-pow(2.0,-alpha)*x)*cos(x)/sqrt(x) ; |
| } |
| |
| double f15 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return x*x * exp(-pow(2.0,-alpha)*x) ; |
| } |
| |
| double f16 (double x, void * params) { |
| double alpha = *(double *) params ; |
| if (x==0 && alpha == 1) return 1 ; /* make the function continuous in x */ |
| if (x==0 && alpha > 1) return 0 ; /* avoid problems with pow(0,1) */ |
| return pow(x,alpha-1)/pow((1+10*x),2.0) ; |
| } |
| |
| double f17 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return pow(2.0,-alpha)/(((x-1)*(x-1)+pow(4.0,-alpha))*(x-2)) ; |
| } |
| |
| /* f454(x) = x^3 log|(x^2-1)(x^2-2)| */ |
| /* integ(f454,x,0,inf) = 61 log(2) + (77/4) log(7) - 27 */ |
| |
| double f454 (double x, void * params) { |
| double x2 = x * x; |
| double x3 = x * x2; |
| params = 0 ; |
| return x3 * log(fabs((x2 - 1.0) * (x2 - 2.0))) ; |
| } |
| |
| /* f455(x) = log(x)/(1+100*x^2) */ |
| /* integ(f455,x,0,inf) = -log(10)/20 */ |
| |
| double f455 (double x, void * params) { |
| params = 0 ; |
| return log(x) / (1.0 + 100.0 * x * x) ; |
| } |
| |
| /* f456(x) = log(x) */ |
| /* integ(f456*sin(10 pi x),x,0,1) = -(gamma + log(10pi) - Ci(10pi))/(10pi) */ |
| |
| double f456 (double x, void * params) { |
| params = 0 ; |
| if (x == 0.0) |
| { |
| return 0; |
| } |
| return log(x) ; |
| } |
| |
| /* f457(x) = 1/sqrt(x) */ |
| /* integ(f457*cos(pi x / 2),x,0,+inf) = 1 */ |
| |
| double f457 (double x, void * params) { |
| params = 0 ; |
| if (x == 0.0) |
| { |
| return 0; |
| } |
| return 1/sqrt(x) ; |
| } |
| |
| /* f458(x) = 1/(1 + log(x)^2)^2 */ |
| /* integ(log(x) f458(x),x,0,1) = (Ci(1) sin(1) + (pi/2 - Si(1)) cos(1))/pi |
| = -0.1892752 */ |
| |
| double f458 (double x, void * params) { |
| params = 0 ; |
| |
| if (x == 0.0) |
| { |
| return 0; |
| } |
| else |
| { |
| double u = log(x); |
| double v = 1 + u * u; |
| |
| return 1.0 / (v * v) ; |
| } |
| } |
| |
| /* f459(x) = 1/(5 x^3 + 6) */ |
| /* integ(f459/(x-0),x,-1,5) = log(125/631)/18 */ |
| |
| double f459 (double x, void * params) { |
| params = 0 ; |
| return 1.0 / (5.0 * x * x * x + 6.0) ; |
| } |
| |
| /* myfn1(x) = exp(-x - x^2) */ |
| /* integ(myfn1,x,-inf,inf) = sqrt(pi) exp(-1/4) */ |
| |
| double myfn1 (double x, void * params) { |
| params = 0; |
| return exp(-x - x*x) ; |
| } |
| |
| /* myfn2(x) = exp(alpha*x) */ |
| /* integ(myfn2,x,-inf,b) = exp(alpha*b)/alpha */ |
| |
| double myfn2 (double x, void * params) { |
| double alpha = *(double *) params ; |
| return exp(alpha*x) ; |
| } |