Ojo! esta tarea hace un montón de cosas que __no tienen que ver con integración__

integracion.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define ERROR_ACEPTABLE (0.00001)

double segundointegrando(double x){
  return (1+(x*x))*log((x*x))*log((x*x))*2*sqrt((x*x))/((1-(x*x)*(x*x)*(x*x)));
}

double segundaintegral(int N){
  double sumaimpar = 0.0;
  double sumapar = 0.0;

  double a = 0;
  double b = 1;

  /* modificamos N con el proposito de que sea par */
  if(N%2 != 0) N = N+1;

  double h = (b-a)/N;

  double xi;
  int i;
  for(i=1, xi=a+h; i<=N-1; i++, xi+=h){
    double fi = segundointegrando(xi);

    if(i%2 == 0){ /* i es par */
      sumapar += fi;
    } else { /* i es impar */
      sumaimpar += fi;
    }
  }

  /* se aplica la formula de simpson 1/3 finalmente
     usando que analiticamente f(a) y f(b) es 0) */
  double valorfinal = (0 + 4.0*sumaimpar + 2.0*sumapar + 0)*h/3.0;

  return valorfinal;
}


double integrandogamma01(double x){
  double x10 = x*x*x*x*x*x*x*x*x*x;
  return pow(tan(x10),(0.1-1.0))*exp(-tan(x10))*(1+tan(x10)*tan(x10))*10*pow((x10),(9.0/10));
}

double integrandogamma05(double x){
  double x2 = x*x;
  return pow(tan(x2),(0.5-1.0))*exp(-tan(x2))*(1+tan(x2)*tan(x2))*2*pow((x2),(1.0/2));
}

double integrandogamma095(double x){
  double x2 = x*x;
  return pow(tan(x2),(0.95-1.0))*exp(-tan(x2))*(1+tan(x2)*tan(x2))*2*pow((x2),(1.0/2));
}

double integrandogamma3(double x){
  return pow(tan(x),(3.0-1.0))*exp(-tan(x))*(1+tan(x)*tan(x));
}

double gamma01(int N){
  double sumaimpar = 0.0;
  double sumapar = 0.0;

  double a = 0;
  double b = pow(M_PI/2.0, 0.1);

  /* modificamos N con el proposito de que sea par */
  if(N%2 != 0) N = N+1;

  double h = (b-a)/N;

  double xi;
  int i;
  for(i=1, xi=a+h; i<=N-1; i++, xi+=h){
    double fi = integrandogamma01(xi);

    if(i%2 == 0){ /* i es par */
      sumapar += fi;
    } else { /* i es impar */
      sumaimpar += fi;
    }
  }

  /* se aplica la formula de simpson 1/3 finalmente
     usando que analiticamente f(a)=10 y f(b)=0) */
  double valorfinal = (10.0 + 4.0*sumaimpar + 2.0*sumapar + 0)*h/3.0;

  return valorfinal;
}

double gamma05(int N){
  double sumaimpar = 0.0;
  double sumapar = 0.0;

  double a = 0;
  double b = pow(M_PI/2.0, 0.5);

  /* modificamos N con el proposito de que sea par */
  if(N%2 != 0) N = N+1;

  double h = (b-a)/N;

  double xi;
  int i;
  for(i=1, xi=a+h; i<=N-1; i++, xi+=h){
    double fi = integrandogamma05(xi);

    if(i%2 == 0){ /* i es par */
      sumapar += fi;
    } else { /* i es impar */
      sumaimpar += fi;
    }
  }

  /* se aplica la formula de simpson 1/3 finalmente
     usando que analiticamente f(a)=2 y f(b)=0) */
  double valorfinal = (2.0 + 4.0*sumaimpar + 2.0*sumapar + 0)*h/3.0;

  return valorfinal;
}

double gamma095(int N){
  double sumaimpar = 0.0;
  double sumapar = 0.0;

  double a = 0;
  double b = pow(M_PI/2.0, 0.5);

  /* modificamos N con el proposito de que sea par */
  if(N%2 != 0) N = N+1;

  double h = (b-a)/N;

  double xi;
  int i;
  for(i=1, xi=a+h; i<=N-1; i++, xi+=h){
    double fi = integrandogamma095(xi);

    if(i%2 == 0){ /* i es par */
      sumapar += fi;
    } else { /* i es impar */
      sumaimpar += fi;
    }
  }

  /* se aplica la formula de simpson 1/3 finalmente
     usando que analiticamente f(a)=10 y f(b)=0) */
  double valorfinal = (0.0 + 4.0*sumaimpar + 2.0*sumapar + 0)*h/3.0;

  return valorfinal;
}

double gamma3(int N){
  double sumaimpar = 0.0;
  double sumapar = 0.0;

  double a = 0;
  double b = M_PI/2.0;

  /* modificamos N con el proposito de que sea par */
  if(N%2 != 0) N = N+1;

  double h = (b-a)/N;

  double xi;
  int i;
  for(i=1, xi=a+h; i<=N-1; i++, xi+=h){
    double fi = integrandogamma3(xi);

    if(i%2 == 0){ /* i es par */
      sumapar += fi;
    } else { /* i es impar */
      sumaimpar += fi;
    }
  }

  /* se aplica la formula de simpson 1/3 finalmente
     usando que analiticamente f(a)=10 y f(b)=0) */
  double valorfinal = (0.0 + 4.0*sumaimpar + 2.0*sumapar + 0)*h/3.0;

  return valorfinal;
}

void datosgamma01(){
  FILE *fp = fopen("datosgamma01.dat","w");
  
  double cero = 0.0;
  fprintf(fp, "%10.10f\t%10.10f\n", cero, 10.0);

  int N = 200;
  double a = 0.0;
  double b = pow(M_PI, 0.1);
  double h = (b-a)/N;
  double cont = h;
  
  while(cont < b){
    fprintf(fp,"%10.10f\t%10.10f\n",cont, integrandogamma01(cont));
    cont += h;
  }

  fprintf(fp, "%10.10f\t%10.10f\n", b, 0.0);

  fclose(fp);
}

void datosgamma05(){
  FILE *fp = fopen("datosgamma05.dat","w");
  
  double cero = 0.0;
  fprintf(fp, "%10.10f\t%10.10f\n", cero, 2.0);

  int N = 200;
  double a = 0.0;
  double b = sqrt(M_PI/2.0);
  double h = 1.0/N;
  double cont = h;
  
  while(cont < b){
    fprintf(fp,"%10.10f\t%10.10f\n",cont, integrandogamma05(cont));
    cont += h;
  }

  fprintf(fp, "%10.10f\t%10.10f\n", b, 0.0);

  fclose(fp);
}

void datosgamma095(){
  FILE *fp = fopen("datosgamma095.dat","w");
  
  double cero = 0.0;
  fprintf(fp, "%10.10f\t%10.10f\n", cero, 0.0);

  int N = 200;
  double a = 0.0;
  double b = sqrt(M_PI/2);
  double h = (b-a)/N;
  double cont = h;
  
  while(cont < b){
    fprintf(fp,"%10.10f\t%10.10f\n",cont, integrandogamma095(cont));
    cont += h;
  }

  fprintf(fp, "%10.10f\t%10.10f\n", b, 0.0);

  fclose(fp);
}

void datosgamma3(){
  FILE *fp = fopen("datosgamma3.dat","w");
  
  double cero = 0.0;
  fprintf(fp, "%10.10f\t%10.10f\n", cero, 0.0);

  int N = 200;

  double a = 0.0;
  double b = M_PI/2.0;

  double h = (b-a)/N;
  double cont = h;
  
  while(cont < b){
    fprintf(fp,"%10.10f\t%10.10f\n",cont, integrandogamma3(cont));
    cont += h;
  }

  fprintf(fp, "%10.10f\t%10.10f\n", b, 0.0);

  fclose(fp);
}

void datossegundaintegral(){
  FILE *fp = fopen("datossegundaintegral.dat","w");
  
  double cero = 0.0;
  fprintf(fp, "%10.10f\t%10.10f\n", cero, 0.0);

  int N = 200;

  double a = 0.0;
  double b = 1.0;

  double h = (b-a)/N;
  double cont = h;
  
  while(cont < b){
    fprintf(fp,"%10.10f\t%10.10f\n",cont, segundointegrando(cont));
    cont += h;
  }

  fprintf(fp, "%10.10f\t%10.10f\n", b, 0.0);

  fclose(fp);
}

int main(){
  datosgamma01();
  datosgamma05();
  datosgamma095();
  datosgamma3();
  datossegundaintegral();

  int Ninicial = 10;
  int N = Ninicial;

  FILE* f1 = fopen("integral1.dat","w");
  FILE* f2 = fopen("integral2.dat","w");
  FILE* f3 = fopen("integral3.dat","w");
  FILE* f4 = fopen("integral4.dat","w");
  FILE* f5 = fopen("integral5.dat","w");

  printf("El valor de gamma(0.1) usando Simpson 1/3 con 10000 puntos es de %10.10f\n", gamma01(10000));
  printf("El valor de gamma(0.5) usando Simpson 1/3 con 10000 puntos es de %10.10f\n", gamma05(10000));
  printf("El valor de gamma(0.95) usando Simpson 1/3 con 10000 es de %10.10f\n", gamma095(10000));
  printf("El valor de gamma(3) usando Simpson 1/3 con 10000 es de %10.10f\n", gamma3(10000));
  printf("\nEl valor de la segunda integral usando 10000 puntos es de %10.10f\n", segundaintegral(10000));

  printf("\nPara el menor N (nro. de intevarlos) con error %g se requieren, por integral:\n", ERROR_ACEPTABLE);
  
  double error = fabs(gamma01(Ninicial) - gamma01(Ninicial*2))/fabs(gamma01(Ninicial*2));
  while (error > ERROR_ACEPTABLE){ /* error = 10^-5 */
    error = fabs(gamma01(N) - gamma01(N*2))/fabs(gamma01(N*2));

    fprintf(f1, "%10.10f\t%10.10f\n", gamma01(N), 1.0/N);

    N = N + 100;
  }

  printf("Cuando N=%i, el error de gamma(0.1) es de %10.10f\n", N, error);

  Ninicial = 10;
  N = Ninicial;

  error = fabs(gamma05(Ninicial) - gamma05(Ninicial*2))/fabs(gamma05(Ninicial*2));
  while(error > ERROR_ACEPTABLE){
    error = fabs(gamma05(N) - gamma05(N*2))/fabs(gamma05(N*2));
    fprintf(f2, "%10.10f\t%10.10f\n", gamma05(N), 1.0/N);
    N = N + 10;
  }

  printf("Cuando N=%i, el error de gamma(0.5) es de %10.10f\n", N, error);

  Ninicial = 10;
  N = Ninicial;

  error = fabs(gamma095(Ninicial) - gamma095(Ninicial*2))/fabs(gamma095(Ninicial*2));
  while(error > ERROR_ACEPTABLE){
    error = fabs(gamma095(N) - gamma095(N*2))/fabs(gamma095(N*2));
    fprintf(f3, "%10.10f\t%10.10f\n", gamma095(N), 1.0/N);
    N = N + 10;
  }

  printf("Cuando N=%i, el error de gamma(0.95) es de %10.10f\n", N, error);

  Ninicial = 10;
  N = Ninicial;

  error = fabs(gamma3(Ninicial) - gamma3(Ninicial*2))/fabs(gamma3(Ninicial*2));
  while(error > ERROR_ACEPTABLE){
    error = fabs(gamma3(N) - gamma3(N*2))/fabs(gamma3(N*2));
    fprintf(f4, "%10.10f\t%10.10f\n", gamma3(N), 1.0/N);
    N = N + 10;
  }

  printf("Cuando N=%i, el error de gamma(3) es de %10.10f\n", N, error);

  Ninicial = 10;
  N = Ninicial;

  error = fabs(segundaintegral(Ninicial) - segundaintegral(Ninicial*2))/fabs(segundaintegral(Ninicial*2));
  while(error > ERROR_ACEPTABLE){
    error = fabs(segundaintegral(N) - segundaintegral(N*2))/fabs(segundaintegral(N*2));
    fprintf(f5, "%10.10f\t%10.10f\n", segundaintegral(N), 1.0/N);
    N = N + 10;
  }

  printf("Cuando N=%i, el error de la segunda integral es de %10.10f\n", N, error);

  return 0;
}