quinta-feira, 23 de fevereiro de 2017

Python: integração numérica com a regra de Simpson composta


Na biblioteca Numpy existe a função "trapz" para o cálculo numérico de uma integral definida usando o método dos trapézios. Esse método é bom, mas a regra dos 3/8 de Simpson (ver aqui) composta resulta em um erro menor. O código abaixo implementa essa regra e compara o resultado para algumas funções.
Regra dos trapézios composta:
$$I_T = \frac{h}{2} (y_0 + 2y_1 + 2y_2 + ... + 2y_{m-1} + y_m) $$
Regra dos 3/8 de Simpson composta:
$$ I_T = \frac{3h}{8} (y_0 + 3y_1 + 3y_2 + 2y_3 + ... + y_m) $$ com $m$ múltiplo de 3.

Funções de teste escolhidas: $e^{-x}$, $xe^{-x}$, $x^{1/2}$, no intervalo $\{0, 6\}$, com passo $h = 0,5$.


import numpy as np

dt = 0.5;
x = np.arange(0, 6.01, dt);

y1 = x*np.exp(-x);
y2 = np.sqrt(x);
y3 = np.exp(-x);

# Regra dos trapézios: 
y1t = np.trapz(y1, x);
y2t = np.trapz(y2, x);
y3t = np.trapz(y3, x);

# Valores teóricos: 
ytt1 = 1 - (6 + 1)*np.exp(-6);
ytt2 = 2*(6**(3/2))/3;
ytt3 = 1 - np.exp(-6);

erro1 = y1t - ytt1;
erro2 = y2t - ytt2;
erro3 = y3t - ytt3;

# Regra dos 3/8 de Simpson: 
soma1 = y1[0];
soma2 = y2[0];
soma3 = y3[0];
n = 1;
tam = len(x);
while n
    soma1 = soma1 + 3*y1[n] + 3*y1[n+1] + 2*y1[n+2];
    soma2 = soma2 + 3*y2[n] + 3*y2[n+1] + 2*y2[n+2];
    soma3 = soma3 + 3*y3[n] + 3*y3[n+1] + 2*y3[n+2];
    n = n + 3
 
# foi considerado a mais no laço acima: 
# integral pela regra de Simpson composta:
soma1 = soma1 - y1[tam-1];  
soma1 = 3*soma1 * dt/8; 
soma2 = soma2 - y2[tam-1];  
soma2 = 3*soma2 * dt/8; 
soma3 = soma3 - y3[tam-1]; 
soma3 = 3*soma3 * dt/8

erros1 = ytt1 - soma1;
erros2 = ytt2 - soma2;
erros3 = ytt3 - soma3;

# print(x); print(y);print('Rg. Trap., Rg. Simpson, Erro Trap, Erro Simpson:')
print("%1.6f" % y1t, "  %1.6f" % soma1, "   %1.6f" % erro1, "  %1.6f" % erros1);
print("%1.6f" % y2t, "  %1.6f" % soma2, "   %1.6f" % erro2, "  %1.6f" % erros2);
print("%1.6f" % y3t, "  %1.6f" % soma3, "   %1.6f" % erro3, "  %1.6f" % erros3); 


Resultados:

Rg. Trap., Rg. Simpson, Erro Trap, Erro Simpson:
0.961816   0.980515    -0.020833   0.002134
9.728712   9.763015    -0.069247   0.034944
1.018217   0.998257    0.020696   -0.000735

Pode ser percebido que o erro absoluto obtido com a Regra de Simpson é menor que usando a Regra dos Trapézios.

Nenhum comentário:

Postar um comentário