Incrementos Adaptativos (Fehlberg)

Prévia do material em texto

Me´todos Computacionais da F´ısica
Trabalho 2
Fabio Rasera Figueiredo
244430
8 de abril de 2016
1 Introduc¸a˜o ao problema
Neste trabalho se analisa o controle de erro e incrementos adaptativos nos
me´todos de Runge-Kutta. Aqui buscaremos avaliar a precisa˜o dos me´todos
quando na˜o se possui uma soluc¸a˜o exata para comparac¸a˜o. No entanto,
utilizaremos uma equac¸a˜o diferencial com soluc¸a˜o exata conhecida para que
possamos verificar graficamente que as ta´ticas utilizadas sa˜o eficientes. A
equac¸a˜o diferencial a ser resolvida sera´:
x′ = x− t2 + 1 (1)
Em que a soluc¸a˜o exata e´:
x(t) = (t+ 1)2 − 0.5et (2)
2 Controle de erro e incrementos adaptativos
Um me´todo simples de controle consiste em estimar o erro corrente ao avaliar
a diferenc¸a entre a soluc¸a˜o calculada com um passo h e a soluc¸a˜o calculada
com 2 passos h/2. Se o erro calculado encontra-se acima do erro ma´ximo
tolera´vel que fora estipulado, o valor do incremento h e´ modificado para um
menor. Assim se prossegue ate´ que a soluc¸a˜o esteja conforme o erro aceita´vel
pelo usua´rio.
Pode-se, tambe´m, aumentar o incremento quando o erro esta´ abaixo
do limite aceita´vel, assim o programa calcula a soluc¸a˜o em passos maiores
quando poss´ıvel e pode vir a poupar tempo de processamento. Para tal
se utiliza o mesmo processo descrito anteriormente, mas, ale´m da opc¸a˜o de
diminuir o passo ao calcular-se um erro acima do aceita´vel se acrescenta
a opc¸a˜o de aumentar o passo ao detectar-se que o erro esta´ menor que o
aceita´vel. Desta forma o incremento se autoregula, adaptando-se a` necessi-
dade corrente.
1
O incremento novo sera´ calculado pela seguinte expressa˜o:
hnovo = hS
(
�
�c
) 1
m
(3)
Onde m = n+1, sendo n a ordem do me´todo utilizado, � e´ o erro ma´ximo
permitido pelo usua´rio, �c e´ o erro corrente e S sera´ um paraˆmetro para
definir o tamanho do aumento ou diminuic¸a˜o de h. Aqui utilizo S1 = 0.9
para que a diminuic¸a˜o do incremento seja um pouco menor do que a exigida
para que o passo esteja dentro da precisa˜o requerida e S2 = 2.0 para que o
incremento aumente o dobro do exigido.
Esta estrate´gia foi implementada na resoluc¸a˜o da Equac¸a˜o 1 utilizando
os me´todos RK2 e RK4. O co´digo segue abaixo.
//AUTOR: Fabio Rasera
#include <stdio.h>
#include <math.h>
double RK4(double tn, double hn, double xn);
double RK2(double tm, double hm, double xm);
double f(double t, double x);
double sol(double t);
int main(){
double ERRO, ERRO1, x0_rk2, x0_rk4, x0, h, x, x_rk2, x_rk4, x1_rk2, x1_rk4, t ,k1, k2, k3, k4, e;
//Arquivos de saı´da
FILE* pt;
pt=fopen("RK2.txt","w");
FILE* pt1;
pt1=fopen("RK4.txt","w");
e=pow(10,-6); //Para^metro de tolera^ncia para o erro corrente
//Dados iniciais
x0=sol(0); x0_rk2=x0;
t=0.0;
h=0.1;
//---RK2---//
do{
//passo h
x_rk2=RK2(t,h,x0_rk2);//RK2
//passo h/2
x1_rk2=RK2(t,(h/2.0),x0_rk2);
x1_rk2=RK2((t+(h/2.0)), (h/2.0), x1_rk2);
ERRO=fabs(x1_rk2-x_rk2);
if(ERRO>=e) //Testa se o erro e´ inaceita´vel
{
h=h*0.9*pow(e/ERRO , 1.0/3.0);
}
else //Caso o erro seja aceita´vel, utiliza o valor calculado e segue para o pro´ximo ponto
{
x0_rk2=x_rk2;
t+=h;
fprintf(pt, "%12e %12e %12e %.12e %.12e\n",t,x_rk2,h,ERRO,(fabs(sol(t)-x_rk2)) / fabs(sol(t)) );
h=h*2*pow(e/ERRO , 1.0/3.0); //Aumento de incremento
}
}while(t<=10);
fclose(pt);
//Dados iniciais
x0_rk4=x0;
2
h=0.1;
t=0.0;
//---RK4---//
do{
//Passo h
x_rk4=RK4(t,h,x0_rk4);
//Passo h/2
x1_rk4=RK4(t,(h/2.0),x0_rk4);
x1_rk4=RK4((t+(h/2.0)), (h/2.0), x1_rk4);
ERRO1=fabs(x1_rk4-x_rk4);
if(ERRO1>=e)//Testa se o erro e´ inaceita´vel
{
h=h*0.9*fabs(pow(e/ERRO1 , 1.0/5.0));
}
else //Caso o erro seja aceita´vel, utiliza o valor calculado e segue para o pro´ximo ponto
{
x0_rk4=x_rk4;
t+=h;
fprintf(pt1,"%12e %12e %12e %.12e %.12e\n",t,x_rk4,h,ERRO1,(fabs(sol(t)-x_rk4)) / fabs(sol(t)) );
h=h*2*fabs(pow(e/ERRO1 , 1.0/5.0)); //Aumento de incremento
}
}while(t<=10);
fclose(pt1);
return(0);
}
//Func¸~ao derivada
double f(double t, double x){
return(x-(t*t)+1);
}
//Soluc¸~ao analı´tica
double sol(double t){
return((t+1)*(t+1)-0.5*exp(t));
}
//Me´todo RK2
double RK2(double tm, double hm, double xm){
double k1=f(tm,xm);
double k2=f(tm+0.5*hm, xm+0.5*hm*k1);
double xy= xm+k2*hm;
return(xy);
}
//Me´todo RK4
double RK4(double tn, double hn, double xn){
double k1=f(tn,xn);
double k2=f(tn+0.5*hn, xn+0.5*k1*hn);
double k3=f(tn+0.5*hn, xn+0.5*k2*hn);
double k4=f(tn+hn, xn+k3*hn);
double xf=xn+(k1 + 2*k2 + 2*k3 + k4)*hn/6.0;
return(xf);
}
A seguir encontram-se os gra´ficos dos pontos calculados para cada me´todo
plotados sobre a soluc¸a˜o anal´ıtica.
3
 0
 1
 2
 3
 4
 5
 6
 7
 0 0.5 1 1.5 2 2.5 3
x
(t
)
t
RK2
Condições iniciais: x(0)=0.5, h0=0.1, Erro máx=10
−6
RK2 num. c/ inc. adapt.
RK2 num s/ inc. adapt.
Exato
Figura 1: Soluc¸a˜o atrave´s do meto´do RK2 utilizando incremento adaptativo.
 0
 1
 2
 3
 4
 5
 6
 7
 0 0.5 1 1.5 2 2.5 3
x
(t
)
t
RK4
Condições iniciais: x(0)=0.5, h0=0.1, Erro máx=10
−6
RK4 num. c/ inc. adapt.
RK4 num s/ inc. adapt.
Exato
Figura 2: Soluc¸a˜o atrave´s do me´todo RK4 utilizando incremento adaptativo.
Ao analisarmos as Figuras 1 e 2, nota-se que com um incremento adap-
4
tativo aparecem mais pontos calculados em regio˜es curvas, enquanto que
nas regio˜es em que a func¸a˜o mais se aproxima de uma reta o incremento
aumenta, permitindo passos maiores entre um ponto e outro. Tambe´m e´
vis´ıvel que a quantidade de pontos necessa´rios para andar sobre a curva e´
bastante menor para o me´todo RK4 em relac¸a˜o ao RK2, o que significa que
o RK4 na˜o necessitou adaptar o incremento tantas vezes quanto o RK2,
indicando que o me´todo RK4 tem maior precisa˜o.
Ao fazermos uma interpolac¸a˜o com os pontos calculados utlizando incre-
mentos adaptativos espera-se uma curva muito pro´xima da soluc¸a˜o exata.
Caso deseje-se um nu´mero maior de pontos calculados, e´ poss´ıvel evitar
o aumento do incremento, diminuindo-o e mantendo-o constante ate´ que
seja necessa´rio diminu´ı-lo novamente. No entanto, se a func¸a˜o calculada e´
bastante varia´vel, apresentando curvatura em algumas regio˜es enquanto em
outras na˜o, uma quantidade considera´vel de processamento pode ser evitada
ao aumentarmos o incremento quando poss´ıvel.
Para analisarmos o impacto do erro permitido pelo usua´rio, variou-se
este paraˆmetro e foram plotados os resultados. Abaixo seguem os gra´ficos.
 0
 1
 2
 3
 4
 5
 6
 7
 0 0.5 1 1.5 2 2.5 3
x
(t
)
t
RK2
Condições iniciais: x(0)=0.5, h0=0.1
Erro máx=10−6
Erro máx=10−3
Erro máx=10−1
Exato
Figura 3: RK2 com diferentes erros ma´ximos permitidos.
5
 0
 1
 2
 3
 4
 5
 6
 7
 0 0.5 1 1.5 2 2.5 3
x
(t
)
t
RK4
Condições iniciais: x(0)=0.5, h0=0.1
Erro máx=10−6
Erro máx=10−3
Erro máx=10−1
Exato
Figura 4: RK4 com diferentes erros ma´ximos permitidos.
Veˆ-se nitidamente nas Figuras 3 e 4 que quando o erro ma´ximo permitido
aumenta, os resultados obtidos fogem da soluc¸a˜o exata. Com o erro ma´ximo
da ordem de 10−1 a diferenc¸a entre a soluc¸a˜o exata e a calculada e´ gritante
e o nu´mero de pontos calculados no intervalo e´ muito pequeno, pois o limite
definido permite saltos enormes entre um ponto e outro sem que o incremento
precise ser ajustado.
2.1 Me´todo de Fehlberg
Outro modo de estimar o erro do me´todo nume´rico de Runge-Kutta e´ uti-
lizar um me´todo de ordem m para avaliar um de ordem n < m. Como os
me´todos de maior ordem possuem maior precisa˜o devido ao truncamento em
termos de ordem mais elevada,podemos utilizar o me´todo RK2 para avaliar
a precisa˜o do RK1, ou o RK4 para avaliar a precisa˜o de ambos os anterio-
res, e assim por diante. Enta˜o, para avaliarmos o me´todo Runge-Kutta de
ordem 4, precisaremos de um me´todo de ordem maior. Para os me´todos
de Runge-Kutta em que a ordem e´ m > 4 o valor da func¸a˜o ao final do
intervalo precisa de um nu´mero de pontos maior do que o nu´mero da ordem
do me´todo, levando a uma inconvenieˆncia ao buscarmos me´todos de ordens
maiores. Fehlberg, no entanto, descobriu um me´todo de quinta ordem que
necessita do ca´lculo de 6 posic¸o˜es da func¸a˜o dentro do intervalo h e cuja
combinac¸a˜o destes mesmos 6 pontos pode ser feita de modo a resultar no
me´todo de Runge-Kutta de ordem 4. Assim, em apenas 6 posic¸o˜es dentro
6
do intervalo h obteˆm-se o RK4 e o RK5. Utilizando o me´todo de Fehlberg
podemos aplicar um controle de erro no RK4. O co´digo que implementa
esta te´cnica adjunta do incremento adaptativo discutido anteriormente se-
gue abaixo:
//AUTOR: Fabio Rasera
#include <stdio.h>
#include <math.h>
double f(double t, double x);
int main(void){
int i, k;
double t,e,Ec,x0,x,h,xm;
double m1, m2, m3, m4, k1, k2, k3, k4, k5, k6, a1, a2, a3, a4, a5, a6, b21, b31, b32, b41, b42, b43, b51, b52, b53, b54, b61, b62, b63, b64, b65, b6, c1, c2, c3, c4, c5, c6, cm1, cm2, cm3, cm4, cm5, cm6;
//Coeficientes de Cash e Karp
a2=1.0/5.0; a3=3.0/10.0; a4=3.0/5.0; a5=1.0; a6=7.0/8.0;
b21=1.0/5.0; b31=3.0/40.0; b41=3.0/10.0; b51=-11.0/54.0; b61=1631.0/55296.0;
b32=9.0/40.0; b42=-9.0/10.0; b52=5.0/2.0; b62=175.0/512.0;
b43=6.0/5.0; b53=-70.0/27.0; b63=575.0/13824.0;
b54=35.0/27.0; b64=44275.0/110592.0;
b65=253.0/4096.0;
c1=37.0/378.0; c2=0.0; c3=250.0/621.0; c4=125.0/594.0; c5=0.0; c6=512.0/1771.0;
cm1=2825.0/27648.0; cm2=0.0; cm3=18575.0/48384.0; cm4=13525.0/55296.0; cm5=277.0/14336.0; cm6=1.0/4.0;
FILE *p;
p = fopen("fehlberg.txt","w");
//Condic¸~oes iniciais
x0=0.5;
x=x0;
t=0.0;
h=0.5;
e=pow(10,-8); //Erro ma´ximo permitido
i=1;
//RK4 e RK5 pelo me´todo de Fehlberg
do{
k1=h*f(t,x0);
k2=h*f(t+a2*h, x0+b21*k1);
k3=h*f(t+a3*h, x0+b31*k1 + b32*k2);
k4=h*f(t+a4*h, x0+b41*k1 + b42*k2 + b43*k3);
k5=h*f(t+a5*h, x0+b51*k1 + b52*k2 + b53*k3 + b54*k4);
k6=h*f(t+a6*h, x0+b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5);
x= x0 + cm1*k1 + cm2*k2 + cm3*k3 + cm4*k4 + cm5*k5 + cm6*k6; //RK4
xm= x0 + c1*k1 + c2*k2 + c3*k3 + c4*k4 + c5*k5 + c6*k6; //RK5
Ec=fabs(x-xm); //Erro corrente
if(Ec>=e) //Testa se o erro corrente e´ maior que o permitido
{
h=0.9*h*(fabs(pow(e/Ec,0.2))); //Diminui o incremento
}
else //Caso o erro corrente seja menor que o permitido, segue o programa
{
fprintf(p,"%.12e %.12e\n",t,x0);
x0=x;
t+=h;
h=2*h*(fabs(pow(e/Ec,0.2))); //Aumenta o incremento
}
}while(t<10);
fclose(p);
return 0;
}
double f(double t, double x){
return(x-(t*t)+1);
}
Em seguida um gra´fico ilustra a soluc¸a˜o calculada atrave´s do me´todo de
7
Fehlberg com incrementos adaptativos em comparac¸a˜o com o RK4 de passo
fixo e a soluc¸a˜o exata.
 0
 1
 2
 3
 4
 5
 6
 7
 0 0.5 1 1.5 2 2.5 3
x
(t
)
t
Fehlberg
Condições iniciais: x(0)=0.5, h0=0.5, Erro máx=10
−8
RK4 num s/ inc. adapt.
RK4 num. c/ inc. adapt (Fehlberg).
Exato
Figura 5: Me´todo de Fehlberg com incremento adaptativo.
8
	Introdução ao problema
	Controle de erro e incrementos adaptativos
	Método de Fehlberg