/* Metodos de cuadraturas: metodo de los rectangulos */

kill(all)$

a:0;b:1;
x[0]:a;
Nu:20;
x[n]:=x[0]+n*(b-a)/Nu;

define(f(x),x^2);
rec:sum(f((x[k]+x[k+1])/2),k,0,Nu-1)*((b-a)/Nu); 
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(rec-exac));

kill(all)$
a:0;b:1;
x[0]:a;
Nu:20;
x[n]:=x[0]+n*(b-a)/Nu;
define(f(x),exp(x));
rec:sum(f((x[k]+x[k+1])/2),k,0,Nu-1)*((b-a)/Nu)$
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(rec-exac));

/* Metodos de cuadraturas: metodo de los trapecios */

kill(all)$
a:0;b:1;
x[0]:a;
Nu:20;
x[n]:=x[0]+n*(b-a)/Nu;
define(f(x),x^2);
tra: (f(a)+f(b)+ 2*sum(f(x[k]),k,1,Nu-1))*((b-a)/(2*Nu))$
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(tra-exac));

kill(all)$
a:0;b:1;
x[0]:a;
Nu:100;
x[n]:=x[0]+n*(b-a)/Nu;
define(f(x),exp(x));
tra: (f(a)+f(b)+ 2*sum(f(x[k]),k,1,Nu-1))*((b-a)/(2*Nu))$
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(tra-exac));

/* Metodos de cuadraturas: metodo de simpson */

kill(all)$
a:0;b:1;x[0]:a;Nu:10;
x[n]:=x[0]+n*(b-a)/(2*Nu);

define(f(x),x^2);
simp: (f(a)+f(b)+ 4*sum(f(x[2*k-1]),k,1,Nu)+2*sum(f(x[2*k]),k,1,Nu-1))*((b-a)/(6*Nu))$
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(simp-exac));

kill(all)$
a:0;b:1;
x[0]:a;
Nu:50;
x[n]:=x[0]+n*(b-a)/(2*Nu);
define(f(x),exp(x));
simp: (f(a)+f(b)+ 4*sum(f(x[2*k-1]),k,1,Nu)+2*sum(f(x[2*k]),k,1,Nu-1))*((b-a)/(6*Nu))$
float(%);

exac:float(integrate(f(x),x,a,b));
float(abs(simp-exac));

quad_qag (x^(2), x, 0, 1, 3);

quad_qag (x^(1/2)*log(1/x), x, 0, 1, 3, 'epsrel=5d-8);

quad_qag (x^(1/2)*log(1/x), x, 0, 1, 3);

intG:quad_qag (x^(1/2)*(1-x)^(3/2), x, 0, 1, 3);

kill(f)$
a:0;b:1;
x[0]:a;
Nu:100;
x[n]:=x[0]+n*(b-a)/(2*Nu);
define(f(x),x^(1/2)*(1-x)^(3/2) );
simp: (f(a)+f(b)+ 4*sum(f(x[2*k-1]),k,1,Nu)+2*sum(f(x[2*k]),k,1,Nu-1))*((b-a)/(6*Nu))$
float(%);

intG[1]-%; /*intG es el valor del paquete quad_qag*/

ex:integrate(x^(1/2)*(1-x)^(3/2), x, 0, 1);

float(ex-intG[1]);

float(ex-simp);

integrate(exp(-x^2),x,0,1);

exerr:float(%);

intEr:quad_qag (exp(-x^2), x, 0, 1, 3);

intEr[1]-exerr;

kill(f)$
a:0;b:1;
x[0]:a;
Nu:100;
x[n]:=x[0]+n*(b-a)/(2*Nu);
define(f(x),exp(-x^2) );
simper: (f(a)+f(b)+ 4*sum(f(x[2*k-1]),k,1,Nu)+2*sum(f(x[2*k]),k,1,Nu-1))*((b-a)/(6*Nu))$
float(%);

%-exerr;

taylor(exp(-x^2),x,0,10);

integrate(%,x,0,1);float(%);

%-exerr;

n:100;
float(sum(1/(n+k),k,1,n));

integrate(1/(1+x),x,0,1);
float(%);

n:500;
su2:float(sum(3*n/(n^2+k^2),k,1,n));

exsu2:integrate(3/(1+x^2),x,0,1);

float(%);

float(exsu2-su2);