MÉTODO DE SIMPSON DE 1/3

Este método sirve para calcular las integrales   por ejemplo de tipo exponencial elevada al cuadrado como la que se trabajo en clase. Recordemos que el calcular una integral  nos significa otra cosa mas que el área que se encuentra debajo de la curva.

Algo muy importante que debemos tener en cuenta es que para este método recuerden que "n" debe ser PAR.. se preguntaran ¿Por Qué?... esto por que para generar una parábola  es necesario contar con 2 puntos.

Se podría calcular  si "n" fuera impar pero se utilizaría otro método el del trapecio y para obtener dicha área se realizaría  la suma del método de simpson 1/3 y el método del trapecio.

Un ejemplo del código:


#include<iostream.h>
#include<conio.h>
#include<stdio.h>
#include<math.h>
main()
{
clrscr();
long double coef[20],a,b,h, x[20],y[20], I, div, promedio;
long double deriv[20];
int exp[20],n;// grado, n;
#ifndef BOOLEAN_H
#define BOOLEAN_H
      //typedef int bool;
      //const int true  = 1;
      //const int false = 0;
#endif
 // int numero, cont = 0;
     int cont=0;
   double grado;
   bool continuar;
   do
   {
      continuar = false;
      cin.clear();
      if(cont > 0) cin.ignore(1024, '\n');
cout<<"........Programa que hace  un tercio de simson........."<<endl<<endl;
cout<<"Dame el grado del polinomio a integrar:"<<endl;
cin>>grado;
cont++;
      if(cin.fail() && cin.rdstate())
      {
      cout << "NO ES UN NUMERO" << endl;
      continuar = true;
      }
}while (continuar);
cout<<endl<<endl;
//*******************************************1
for(int i=0;i<=grado; i++)
{

     cout<<"Dame el coeficiente de x^"<<i<<endl;
     cin>>coef[i];
     deriv[i]=coef[i];
     exp[i]=i;
//deriv[i]=coef[i];
//exp[i]=i;

}
//*******************************************1
cout<<"La funcion quedara de la siguiente forma:"<<endl<<endl;
cout<<"\tf(x)= ";

//*******************************************2
for(int c=0; c<=grado;c++)
{
if(exp[c]!=0)
{
cout<<"+";
}
if(exp[c]!=0)
{
cout<<"("<<coef[c]<<"x^"<<exp[c]<<")";
}
if(exp[c]==0)
{
cout<<"("<<coef[c]<<")";
}
}
//*******************************************2
cout<<endl;
    int  cont2 = 0;
   // double a;
    bool continuar2;
    do
    {
    continuar2= false;
    cin.clear();

    if(cont2 > 0) cin.ignore(1024, '\n');

cout<<"Dame el punto 'a' donde se va a evaluar la funcion"<<endl;
cin>>a;
    cont2++;

    if(cin.fail() && cin.rdstate())
    {
    cout << "NO ES UN NUMERO" << endl;
    continuar2 = true;
    }
}while (continuar2);
//++++++++++
    int  cont3 = 0;
  //  double C;
    bool continuar3;
    do
    {
    continuar3= false;
    cin.clear();

    if(cont3 > 0) cin.ignore(1024, '\n');

cout<<"Dame el punto 'b' donde se va a evaluar la funcion"<<endl;
cin>>b;
    cont3++;

    if(cin.fail() && cin.rdstate())
    {
    cout << "NO ES UN NUMERO" << endl;
    continuar3 = true;
    }
}while (continuar3);
cout<<endl<<"En cuantos segmentos se va a dividir"<<endl;
scanf("%d",&n);
while (pow(-1,n)<0)  //Esto es por si es impar
{
printf("\n -------Lo SABIAMOs  siempre hay alguien q no sabe q es par-------- ");
scanf("%d",&n);
//while 1
}

h=((b-a)/(n));
//*******************************************3

x[0]=a;

cout<<"Valores:"<<endl<<"\t"<<"x0="<<x[0]<<endl<<endl;
for(int d=1; d<=n; d++)
{

x[d]=x[d-1]+h;

cout<<"\tx"<<d<<"="<<x[d]<<endl<<endl;

}
//*******************************************3


//*******************************************4


for(int e=0; e<=n; e++)
{
long double fx=0;

for(int f=0; f<=grado; f++)
 {

if(f!=0)
   {
fx=(coef[f]*(pow(x[e],f)))+fx;
   }
if(f==0)
   {
fx=coef[f]+fx;
   }


 }
y[e]=fx;
cout<<endl<<"El valor de:"<<endl<<endl<<"\t f(x"<<e<<")="<<y[e]<<endl<<endl;

}
//*******************************************4

long double impar,par;
//******************************************5
long double auximpar=0;
long double auxpar=0;

for(int g=1; g<=(n-1); g+=2)
{

auximpar=y[g]+auximpar;
//cout<<endl<<auximpar<<endl;

}
impar=auximpar*4;
//cout<<endl<<impar<<endl;

//******************************************5

//******************************************6

for(int j=2; j<=(n-2); j+=2)
{

auxpar=y[j]+auxpar;
//cout<<endl<<auxpar<<endl;

}
par=auxpar*2;
//cout<<endl<<par<<endl;
//******************************************6

div=((y[0]+par+impar+y[n])/(3*n));

I=((b-a)*div)    ;

cout<<"El resultado de la integral es: "<<I<<endl;





//derivar*******************************************************************

//int ;

//auxgrado=grado;
//****************************************7
for(int l=grado; l>(grado-4);l--)
{


for(int k=1; k<=l; k++)
{

deriv[k-1]=(deriv[k]*k);



}


}
//*****************************************7




cout<<endl<<endl<<"La cuarta derivada es:"<<endl;
//****************************************8
cout<<"fIV(x)=";
for(int m=0; m<=(grado-4);m++)
{


if(exp[m]!=0)
{
cout<<"+";
}
if(exp[m]!=0)
{
cout<<"("<<deriv[m]<<"x^"<<exp[m]<<")";
}
if(exp[m]==0)
{
cout<<"("<<deriv[m]<<")";
}

}
//***************************************8

promedio=((b-a)/2);


//*******************************************9



long double fxaux=0;

for(int o=0; o<=(grado-4); o++)
 {

if(o!=0)
   {
fxaux=(deriv[o]*(pow(promedio,o)))+fxaux;
   }
if(o==0)
   {
fxaux=deriv[o]+fxaux;
   }


 }
//y[e]=fx;

long double numerador, dominador, conciente, erroraprox;
numerador=(pow((b-a),5));
dominador=((pow(n,4))*180);
conciente=((numerador)/(dominador))*(-1);
erroraprox=(conciente)*(fxaux);

cout<<"El valor del error aproximado es: "<<endl;
cout<<erroraprox;

//*******************************************


getch();
return 0;
}
 Imagen:

En este enlace podemos encontrar una explicacion detallada de como se resolveria paso a paso y la expliacion del método bien especificada