// ********************************************************************* 
//                      METODOS NUMERICOS                              *
// *********************************************************************
//                  Ejemplo de integración numérica                    *
// *********************************************************************


// Definicion de la funcion a integrar
deff('y=f(x)','y=2 + sin(2*sqrt(x))');

// Definicion de la solucion real
deff('y=If(x)','y= 2*x - sqrt(x) * cos(2*sqrt(x)) + (1/2)*sin(2*sqrt(x))');

// Limites de Integracion
a = 1;
b = 6;

//  Hallando el valor real de la integral
I_real = feval(b,If) - feval(a,If)
pause

//  Hallando el valor de la integral por aproximacion numerica

//  Analicemos el efecto del tamaño del paso de integracion

N = [3 10 20 40 80 160];    // El intervalo [1,6] divideremos en N sub-intervalos
casos = length(N);          // Definiendo el numero de casos que seran estudiados
result_trap = [];
result_simp = [];
result_trap_scilab = [];

// Cargando funciones que contienen los algoritmos
getf('trapecio.sci')       // Cargando funcion que efectua integracion numerica por R. del Trapecio
getf('simpson.sci')       // Cargando funcion que efectua integracion numerica por R. de Simpson

for i = 1:casos
   n = N(i);																// N° de intervalos
   h = (b-a)/n;                   // Tamaño del intervalo (paso de integracion)
   I_trap = trapecio(f,a,b,n);    // Llamando a la funcion que efectua la integracion por Regla del Trapecio
   E_trap = I_real - I_trap;      // Hallando el error
  	result_trap = [result_trap; n h I_trap E_trap];  // Matriz de resultado

   I_simpson = simpson(f,a,b,n);   // Llamando a la funcion que efectua la integracion por Regla del Trapecio
   E_simpson = I_real - I_simpson;  // Hallando el error
	result_simp = [result_simp; n h I_simpson E_simpson];    // Matriz de resultado
end

disp('Resultado de Integracion por Regla del Trapecio')
disp('   N       h            Int_calc      Error')
disp(result_trap)

pause

disp('Resultado de Integracion por Regla de Simpson')
disp('   N       h            Int_calc      Error')
disp(result_simp)

// Graficando el Valor Absoluto del Error Obtenido - Usando escala semilogaritmica
plot2d([result_trap(:,2) result_simp(:,2)],[abs(result_trap(:,4)) abs(result_simp(:,4))],[1,2],leg="Trapecio@Simpson",logflag ="nl");
xtitle('Error en Metodos Numericos de Integracion', 'Paso de Integracion h', 'Valor Absoluto Error')

pause


// Comparando ahora los resultados obtenidos con las funciones propias de Scilab

disp('Usando la funcion Inttrap de Scilab')

// Comparando con la funcion inttrap
for i = 1:casos
   n = N(i);																		 // N° de intervalos
   h = (b-a)/n;                   				 // Tamaño del intervalo (paso de integracion)
   x = a:h:b;                     				 // Creando puntos de la funcion x
   y = f(x);                      				 // Hallando los valores correspondientes de y=f(x)
   I_trap_scilab = inttrap(x, y);   	 // Llamando a la funcion de Scilab
   E_trap_scilab = I_real - I_trap_scilab;      // Hallando el error
  	result_trap_scilab = [result_trap_scilab; n h I_trap_scilab E_trap_scilab];  // Matriz de resultado
end


result_trap_scilab


// Comparancion con la funcion intg
disp('Usando la funcion Intg de Scilab')
[v,err]=intg(a,b,f);             	           // Llamando a la funcion de Scilab    
I_intg_scilab = v
E_intg_scilab = I_real - I_intg_scilab      // Hallando el error
   

// Fin de archivo