#include #define LOGGING #ifdef LOGGING #define Logg(message) std::clog << message << std::endl; #else // LOGGING #define Logg(message) #endif /* NO_ARG_CHECK */ template T integrate_trapez_simple(functor& F, double a, double b, int N) {// Trapezoid rule for fixed number of function evaluations N at points a, a+(b-a)/(N-1),... double dh = (b-a)/(N-1); T sum = 0.5*F(a); double x = a+dh; for (int i=1; i T integrate4(functor& F, double a, double b, int Np) { // Forth order integration routine with equidistant mesh static double coeff[4] = {17./48., 59./48.,43./48.,49./48.}; double dh = (b-a)/(Np-1), x = a; T sum = 0; for (int i=0; i<4; i++, x+=dh) sum += coeff[i]*F(x); for (int i=4; i T integrate_trapez(functor& F, double a, double b, double precision=1e-5, int max_recursion=22) { T olds = 0, s = 0.5*(F(a)+F(b))*(b-a); // lowest order approximation int m=1; // number of subdivisions for (int i=0; i 4 && fabs(s-olds) <= precision*fabs(olds)){ //Avoid spurious early convergence. Logg("Trapezoid: The integral step used =" << 0.5*dh << " with " << m << " divisions"); return s; } olds=s; } std::cerr<<"Too many steps in routine integrate_trapez"< T integrate_simps(functor& F, double a, double b, double precision=1e-5, int max_recursion=22) { T oldsm = 0, sm = 0; // simpson approximation T olds = 0, s = 0.5*(F(a)+F(b))*(b-a); // trapezoid approximation, lowest order approximation int m=1; // number of subdivisions for (int i=0; i 4 && fabs(sm-oldsm) <= precision*fabs(oldsm)){ //Avoid spurious early convergence. Logg("Simpson: The integral step used ="<<0.5*dh<<" with "<(F0, -5, 5, 200); double v1 = integrate_trapez(F0, -5, 5,1e-4); double v2 = integrate_simps(F0, -5, 5,1e-4); cout<(F1, 0.0, 51*M_PI/2., 200); double p1 = integrate_trapez(F1, 0.0, 51*M_PI/2., 1e-4); double p2 = integrate_simps(F1, 0.0, 51*M_PI/2., 1e-4); double p3 = integrate4(F1, 0.0, 51*M_PI/2., 200); cout<