Computação Científica - Um Livro Colaborativo
9.3 Problemas
Faça um código para calcular , salve seus valores em arquivo e visualize a função no scilab.
Vamos integrar usando a regra de Boole
onde
Vamos resolver o problema por partes. Primeiro, vamos implementar uma rotina que integra uma função pelo método de Boole.
#include <math.h>
double f(double t)
{
return t*t*t*t*exp(-50*(t-2)*(t-2));
}
void main(void)
{
int N=10,i;
double a=0,b=4,h=(b-a)/(4*N),Int=0;
Int=14./45.*h*f(a);
for(i=0;i<N;i++)
{
Int+=64./45.*h*f(a+h);
Int+=24./45.*h*f(a+2*h);
Int+=64./45.*h*f(a+3*h);
Int+=28./45.*h*f(a+4*h);
a+=4*h;
}
Int-=14./45.*h*f(a);
printf("%f\n",Int);
}
Agora, vamos usar uma estrutura com ,, e .
#include <math.h>
struct INT
{
double a;
double b;
int N;
};
double f(double t)
{
return t*t*t*t*exp(-50*(t-2)*(t-2));
}
double Boole(struct INT *I)
{
int i,N=I->N;
double a=I->a,b=I->b,h=(b-a)/(4*N),Int=0;
Int=14./45.*f(a);
for(i=0;i<N;i++)
{
Int+=64./45.*f(a+h);
Int+=24./45.*f(a+2*h);
Int+=64./45.*f(a+3*h);
Int+=28./45.*f(a+4*h);
a+=4*h;
}
Int-=14./45.*f(a);
return h*Int;
}
void main(void)
{
struct INT I;
I.a=0;
I.b=4;
I.N=10;
printf("%f\n",Boole(&I));
}
Agora vamos fazer a primeira versão do código que resolve o problema acima. Vamos calcular para onze pontos no intervalo .
#include <math.h>
struct INT
{
double a;
double b;
int N;
};
double f(double t)
{
return t*t*t*t*exp(-50*(t-2)*(t-2));
}
double Boole(struct INT *I)
{
int i,N=I->N;
double a=I->a,b=I->b,h=(b-a)/(4*N),Int=0;
Int=14./45.*f(a);
for(i=0;i<N;i++)
{
Int+=64./45.*f(a+h);
Int+=24./45.*f(a+2*h);
Int+=64./45.*f(a+3*h);
Int+=28./45.*f(a+4*h);
a+=4*h;
}
Int-=14./45.*f(a);
return h*Int;
}
void main(void)
{
struct INT I;
I.a=0;
I.N=10;
int i,M=10;
for(i=0;i<=M;i++)
{
I.b=4.*i/M;
printf("%f\n",Boole(&I));
}
}
Veja, até agora calculamos as integrais sem preocupação com a convergência. Então, vamos implementar um código que calcule as integrais até que o erro fica menor que uma tolerância relativa de .
#include <math.h>
struct INT
{
double a;
double b;
int N;
};
double f(double t)
{
return t*t*t*t*exp(-50*(t-2)*(t-2));
}
double Boole(struct INT *I)
{
int i,N=I->N;
double a=I->a,b=I->b,h=(b-a)/(4*N),Int=0;
Int=14./45.*f(a);
for(i=0;i<N;i++)
{
Int+=64./45.*f(a+h);
Int+=24./45.*f(a+2*h);
Int+=64./45.*f(a+3*h);
Int+=28./45.*f(a+4*h);
a+=4*h;
}
Int-=14./45.*f(a);
return h*Int;
}
double Integral(struct INT *I)
{
double I1,I2,aux;
I1=Boole(I);
do
{
(I->N)*=2;
I2=Boole(I);
aux=fabs(I2-I1);
I1=I2;
}while ((aux/fabs(I2))>1e-8);
return I2;
}
double F(double x)
{
struct INT I;
I.a=0;
I.N=10;
I.b=x;
return Integral(&I);
}
void main(void)
{
int i,M=1000;
double res;
FILE *arq;
arq=fopen("saida.dat","w+b");
if (arq==NULL) printf("Não foi possível abrir o arquivo");
else
{
printf("arquivo aberto com sucesso");
for(i=0;i<=M;i++)
{
res=F(4.*i/M);
fwrite(&res,sizeof(double),1,arq);
printf("%f\n",res);
}
fclose(arq);
}
}
Agora vamos imprimir o resultado em arquivo e abrir no scilab:
#include <math.h>
struct INT
{
double a;
double b;
int N;
};
double f(double t)
{
return t*t*t*t*exp(-50*(t-2)*(t-2));
}
double Boole(struct INT *I)
{
int i,N=I->N;
double a=I->a,b=I->b,h=(b-a)/(4*N),Int=0;
Int=14./45.*f(a);
for(i=0;i<N;i++)
{
Int+=64./45.*f(a+h);
Int+=24./45.*f(a+2*h);
Int+=64./45.*f(a+3*h);
Int+=28./45.*f(a+4*h);
a+=4*h;
}
Int-=14./45.*f(a);
return h*Int;
}
double Integral(struct INT *I)
{
double I1,I2,aux;
I1=Boole(I);
do
{
(I->N)*=2;
I2=Boole(I);
aux=fabs(I2-I1);
I1=I2;
}while ((aux/fabs(I2))>1e-8);
return I2;
}
double F(double x)
{
struct INT I;
I.a=0;
I.N=10;
I.b=x;
return Integral(&I);
}
void main(void)
{
int i,M=1000;
double res;
FILE *arq;
arq=fopen("saida.dat","w+b");
if (arq==NULL) printf("Não foi possível abrir o arquivo");
else
{
printf("arquivo aberto com sucesso");
for(i=0;i<=M;i++)
{
res=F(4.*i/M);
fwrite(&res,sizeof(double),1,arq);
printf("%f\n",res);
}
fclose(arq);
}
}
Abrimos no scilab com os seguintes comandos:
x=linspace(0,4,M);
arq=mopen(’home/fulano/pasta1/pasta2/saida.dat’,’rb’)
Snum=mget(M,’d’,arq)
mclose(arq)
plot(x,Snum)
Comentários sobre os códigos:
- O ponto na linha res=F(4.*i/M); é importante. Substituir a linha por res=F(4*i/M); não produz o resultado esperado.
- As linhas return Integral(&I); e I1=Boole(I); passam um ponteiro para a mesma estrutura INT, mas a primera usa &I e a segunda I. Procure entender essa diferença na notação.
Observe que a rotina de integração vai aumentando o número de pontos até a convergência. Uma alternativa mais eficiente é implementar um esquema auto-adaptativo com refinamento local, isto é, adicionar pontos apenas nas regiões onde a integral ainda não convergiu.
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.