Simpson.cpp

#include <iostream>
#include <cmath>
#include <fstream>

// Función para integrar 
double funcion(double x){
    return 100*x*x*cos(20*x);
}
// Calculo de la integral aproximada bajo simpson
double simpson(double a,double b,int N){
    double h=(b-a)/N;
    double suma=0.0;
    for(int i=0;i<N;++i){
        double xi = a+i*h;
        if(i>0 && i<N-1){
            if(i % 2){
            suma+=4*funcion(xi);
            }
            else suma+=2*funcion(xi);
        }
        else suma+=funcion(xi);
    }
    return h/3*suma;
}

int main(){
std::ofstream datafile ("resultados_simpson.dat");
for (int N=1;N<=10000;++N){
    double integral = simpson(0.0,1,N);
    double integral_exacta = 4.7459;
    datafile<<N<<" "<<integral<<" "<<integral_exacta<<std::endl;
}

datafile.close();

// Script gnuplot integral aprox
std::ofstream scriptFile1("grafico_integral_vs_simpson.gp");
scriptFile1<<"set term png\n";
scriptFile1<<"set output 'grafico_integral_vs_simpson.png'\n";
scriptFile1<<"set xlabel 'N'\n";
scriptFile1<<"set ylabel 'Integral'\n";
scriptFile1<<"set logscale x\n";
scriptFile1<<"plot 'resultados_simpson.dat' u 1:2 w l title 'Integral Aproximada', '' u 1:3 w l title 'Valor Exacto'\n";
scriptFile1.close();

std::ofstream dataFile2("datos_simpson.dat");
int N_plot =10;
double h_plot=(1.0-0.0)/N_plot;
for(double x=0.0;x<=1.0;x+=h_plot){
    dataFile2<<x<<" "<<funcion(x)<<std::endl;
    dataFile2<<x<<" "<<0.0<<std::endl;
}

dataFile2.close();

std::ofstream scriptFile2("grafico_funcion_rect.gp");
scriptFile2<<"set term png\n";
scriptFile2<<"set output 'grafico_funcion_rect.png'\n";
scriptFile2<<"set xlabel 'N'\n";
scriptFile2<<"set ylabel 'Integral'\n";
scriptFile2<<"set logscale x\n";
scriptFile2<<"plot 'datos_cuadratura_rectangular.dat' w l title 'Funcion', '' w boxes title 'Rectangulo'\n";
scriptFile2.close();

// Ejecutar Gnuplot
system("gnuplot grafico_integral_vs_simpson.gp");
//system("gnuplot grafico_integral_exacta.gp");
system("gnuplot grafico_funcion_rect.gp");

return 0;
}