#include #include #include #include #include using namespace std; template void recurs(list& x, list& f, const functor& fun, list::iterator xit, list::iterator fit, double precision, int Max_level, int level){ // Recursive algorithm inserts points into the mesh if the function values // differ more than precision // Recursion can not go deeper than Max_level list::iterator fip = fit; fip++; list::iterator xip = xit; xip++; if (fabs(*(fit)-*(fip))Max_level){ return; } else{ double x_new = 0.5*((*xit)+ (*xip)); list::iterator xin = x.insert(xip, x_new); list::iterator fin = f.insert(fip, fun(x_new)); recurs(x, f, fun, xit, fit, precision, Max_level, level+1); recurs(x, f, fun, xin, fin, precision, Max_level, level+1); } } template void mesh(list& x, list& f, const functor& fun, double a, double b, int Nmin=9, int Max_level=20, double precision=1e-3) { // This function constructs an adaptive mesh for function f. // The mesh is recursively subdivided until the difference between the function values is less than precision // or there were more than Max_level subdivisions - recursion steps. // This function creates a list for mesh and a list of function values. for (int i=0; i::iterator xit = x.begin(); list::iterator fit = f.begin(); vector::iterator> xiter(x.size()); vector::iterator> fiter(x.size()); for (int i=0; i& x, const list& f) {// given two lists : x and f, it computes integral using simpson's rule. list::const_iterator f0, f1, f2; list::const_iterator i0, i1, i2; f0 = f.begin(); i0 = x.begin(); double sum=0; for (; i0!= x.end();){ double a = *i0; i1 = i0; double b = *(++i1); if (i1==x.end()) break; i2 = i1; double c = *(++i2); if (i2==x.end()) break; double fa = *f0; f1 = f0; double fb = *(++f1); f2 = f1; double fc = *(++f2); double dh1 = b-a; double dh2 = c-b; sum += (dh1+dh2)/(6*dh1*dh2)*((3*b-2*a-c)*dh2*fa + (c-a)*(c-a)*fb + (2*c-3*b+a)*dh1*fc); i0++; i0++; f0++; f0++; } if (x.size()%2==0){ i1 = x.end(); f1 = f.end(); double b = *(--i1); double fb = *(--f1); i0 = i1; f0 = f1; double a = *(--i0); double fa = *(--f1); sum += 0.5*(b-a)*(fa+fb); } return sum; } double F0(double x) { return (1-2*sqrt(fabs(x)))*exp(-x*x);} double F1(double x) { return cos(x);} double F3(double x) { return 1/sqrt(x); } int main() { double IF0_exact = -0.3391897770124197*2; list x, f; mesh(x, f, F0, -5., 5.); double integral = simps(x,f); cout<<"Total number of points used "<::const_iterator i=x.begin(); i!=x.end(); i++){ fdat1<<*i<<" "<::const_iterator i=x.begin(); i!=x.end(); i++){ fdat2<<*i<<" "<::const_iterator i=x.begin(); i!=x.end(); i++){ fdat3<<*i<<" "<