Use o botão abaixo para reportar erros ou dar sugestões.

Computação Científica - Um Livro Colaborativo

3.2 Problemas


Exemplo 3.2.1. Implemente um programa que aproxime o valor da integral de f(x) no intervalo [a,b] usando o método de Simpson composto, dado por ∫ abf(x)dx ≈ h 3 f(x1) + 2 ∑ i=1n−1f(x 2i+1) + 4 ∑ i=1nf(x 2i) + f(x2n+1) = h 3 f(x1) + 4f(x2) + 2f(x3) + 4f(x4) + 2f(x5) + 4f(x6) + ⋯ + 2f(x2n−1) + 4f(x2n) + f(x2n+1)

onde h = b − a 2n xi = a + (i − 1)h,i = 1,2,⋯,2n + 1.

#include <stdio.h>  
#include <math.h>  
 
double f(double x)  
{  
  return exp(-x);  
}  
 
double Simpson(double a,double b, int n)  
{  
  double x,h,soma;  
  int i;  
  x=a;  
  h=(b-a)/(2.0*n);  
  soma=f(x);  
  for (i=1;i<=n;i++)  
  {  
    x+=h;  
    soma+=4*f(x);  
    x+=h;  
    soma+=2*f(x);  
  }  
  soma-=f(x);  
  soma*=h/3;  
  return soma;  
}  
 
main (void)  
{  
  printf("Integral=%f\n",Simpson(0,1,20));  
}

Nós ainda não estudamos estruturas de vetores para armazenar a malha x1 = a, x2 = a + h, x3 = a + 2h, ⋯. Por isso, no código acima, calculamos dentro do for usando a relação de recorrência xi = xi−1 + h. Observe que calculamos f(x1) = f(a) antes do for. No for, nós corremos i = 1,...,n e a cada iteração calculamos 4f(xi) e 2f(xi + h). Isso significa que, na última iteração, calculamos 4f(x2n) e 2f(x2n+1) = 2f(b). Mas f(b) deveria ser somado uma vez só, por isso, precisamos diminuir f(b) depois de encerrado o for.

No código da solução do exemplo 3.2.1, nós calculamos uma aproximação para ∫ abf(x)dx, mas não se preocupamos com a precisão do cálculo. Vamos introduzir o seguinte critério de parada: se o valor absoluto da diferença entre as aproximações com n e n + 1 pontos na malha for menor que 10−8, então apresenta o resultado dado pela malha com n + 1 pontos. Façamos isso introduzindo uma nova função, como segue abaixo:

#include <stdio.h>  
#include <math.h>  
 
double f(double x)  
{  
  return exp(-x);  
}  
 
double Simpson(double a,double b, int n)  
{  
  double x,h,soma;  
  int i;  
  x=a;  
  h=(b-a)/(2.0*n);  
  soma=f(x);  
  for (i=1;i<=n;i++)  
  {  
    x+=h;  
    soma+=4*f(x);  
    x+=h;  
    soma+=2*f(x);  
  }  
  soma-=f(x);  
  soma*=h/3;  
  return soma;  
}  
 
double Integral(double a, double b,double tolerancia)  
{  
  int n=1;  
  double I_antigo=0,I_atual=2*tolerancia;  
  while ((fabs(I_antigo-I_atual)>=tolerancia))  
  {  
    I_antigo=Simpson(a,b,n);  
    I_atual=Simpson(a,b,n+1);  
    n++;  
  }  
  return I_atual;  
}  
 
main (void)  
{  
printf("Integral=%.12f\n",Integral(0,1,1e-8));  
}

Nessa versão do código, os valores I_antigo=0 e I_atual=2*tolerancia são para que o while execute a primeira vez. Observe que, se I_antigo=I_atual=0, o programa não entra no while. Uma alternativa para implementar o critério de parada é inserir uma variável de controle que, quando a tolerância é atingida, o programa executa mais algumas iterações. Isso dá uma segurança maior para o resultado obtido, evitando que, por exemplo, I_atual e I_antigo fiquem p?oximos, mas longe do valor da integral. Veja a próxima versão do código:

#include <stdio.h>  
#include <math.h>  
 
double f(double x)  
{  
  return exp(-x);  
}  
 
double Simpson(double a,double b, int n)  
{  
  double x,h,soma;  
  int i;  
  x=a;  
  h=(b-a)/(2.0*n);  
  soma=f(x);  
  for (i=1;i<=n;i++)  
  {  
    x+=h;  
    soma+=4*f(x);  
    x+=h;  
    soma+=2*f(x);  
   }  
  soma-=f(x);  
  soma*=h/3;  
  return soma;  
}  
 
double Integral(double a, double b,double tolerancia)  
{  
  int n=1,controle=2;  
  double I_antigo,I_atual;  
  while (controle)  
  {  
    I_antigo=Simpson(a,b,n);  
    I_atual=Simpson(a,b,n+1);  
    if((fabs(I_antigo-I_atual)<tolerancia)) controle--;  
    else controle=2;  
    n++;  
  }  
  return I_atual;  
}  
 
main (void)  
{  
  printf("Integral=%.12f\n",Integral(0,1,1e-8));  
}

É importante validar a rotina antes de considerar o problema resolvido. Para isso, aproxime a integral de várias funções que conhecemos o valor exato da integral e compare.

Creative Commons License Este texto é disponibilizado nos termos da licença Creative Commons Atribuição-CompartilhaIgual 3.0 Não Adaptada (CC-BY-SA 3.0). Página gerada em 19/8/2020 às 17:33:16.

Informe erros ou edite você mesmo!