sábado, 28 de junho de 2014

Calculando derivadas numericamente


Podemos calcular a derivada de uma função por uma das três fórmulas abaixo:


Em alguns casos é necessário calcular a derivada de alguma função de forma numérica. Em outras situações não temos a função, mas somente uma tabela descrevendo a função. Nesses dois casos podemos calcular as derivadas usando uma das fórmulas acima. Existe alguma fórmula mais indicada? Sim. É possível mostrar que última fórmula (c) apresenta menor erro numérico. E quanto ao valor de "h"? Em princípio, quanto menor o valor para "h", melhor a precisão do resultado numérico. Porém, quando usamos o computador (ou uma calculadora) para efetuar esses cálculos existem erros de arredondamento que aumentam quando o valor de "h" é muito pequeno. Então, vai existir um valor ótimo para o valor de "h". Esse valor ótimo vai depender da precisão numérica do software/hardware usado e também da função. O exemplo numérico a seguir mostra o cálculo desse valor ótimo para três funções e para quatro valores distintos de "x". Dos gráficos abaixo, concluímos que h = 0,00001 é um valor ótimo, isto é, o erro entre a derivada numérica e a derivada analítica é (próximo) do mínimo.

Funções:

Valor ótimo de h para alguns valores de x:

O valor ótimo varia com a função e com x, mas para h ~ 10^-5 ou 10^-6 o erro é mínimo.

Eixos em escala log.

Código Scilab:

function rr=fx1(x)
    rr = cos(2*x).*exp(x);
endfunction

function rrd=fxd(x)
    rrd = (cos(2*x) - 2*sin(2*x)).*exp(x);
endfunction

function rr=fx2(x)
    rr = x.*x.*log(x);
endfunction

function rrd=fxd2(x)
    rrd = x.*(2*log(x) + 1);
endfunction

function rr=fx3(x)
    rr = x.*x.*x - (2.0)./x;
endfunction

function rrd=fxd3(x)
    rrd = 3*x.*x + 2/(x.*x);
endfunction

clc;
close; close; close;

xa = 0.5:0.01:2.5;
f1 = fx1(xa);
f2 = fx2(xa);
f3 = fx3(xa);
plot(xa,f1,'b',xa,f2,'r',xa,f3,'k');
title('Funções de teste');
xlabel('x'); ylabel('fi(x)');
legend('cos(2*x).*exp(x)','x.*x.*log(x)','x.*x.*x - (2.0)./x');

figure;
for px=2:9;
h=0.01; xx=px/5;
erro1 = zeros(1,22);
hk = erro1;
erro2 = erro1;
erro3 = erro1;
for k=1:22
    da = fxd(xx);
    dn = (fx1(xx+h) - fx1(xx-h))/(2*h);
    erro1(k) = da - dn;
    da = fxd2(xx);
    dn = (fx2(xx+h) - fx2(xx-h))/(2*h);
    erro2(k) = da - dn;
    da = fxd3(xx);
    dn = (fx3(xx+h) - fx3(xx-h))/(2*h);
    erro3(k) = da - dn;
    hk(k) = h;
    h = h / 2;
end

if px==6 then
    figure;
end
if px<6 p="" then="">    subplot(2,2,px-1);
    plot(log10(hk),log10(abs(erro1)),'b-o',...
    log10(hk),log10(abs(erro2)),'r-x',...
    log10(hk),log10(abs(erro3)),'k-s');
    legend('Erro função 1','Erro função 2','Erro função 3');
    title('x = ' + string(xx));
else
    subplot(2,2,px-5);
    plot(log10(hk),log10(abs(erro1)),'b-o',...
    log10(hk),log10(abs(erro2)),'r-x',...
    log10(hk),log10(abs(erro3)),'k-s');
    legend('Erro função 1','Erro função 2','Erro função 3');
    title('x = ' + string(xx));
end;
end;



Nenhum comentário:

Postar um comentário