Calculando o volume numericamente: computação científica com Python

Exemplo de uma figura 3D.

Uma tarefa importante é o cálculo de volumes de figuras definidas por alguma função $z = f(x,y)$. Para alguns casos, é possível calcularmos analiticamente o volume desejado, mas em vários outros esse cálculo analítico não é possível. Nesses casos, podemos tentar fazer uma estimativa numérica do problema. Podemos "fatiar" esse volume em pequenos prismas de base retangulares e altura dada pelo valor da função $f(x,y)$ na posição $(x_n,y_m)$ selecionada. Então, teremos vários prismas com volume $dV_{nm} = f(x_n,y_m) dx dy$ e volume pode ser estimado pela soma ponderada do volume desses prismas: $$V_{estimado} \cong h_x h_y \sum_{n}\sum_{m} p(n,m)f(x_n,y_m) $$ onde $p(n,m)$ é peso na posição $(x_n,y_m)$. Esses pesos podem ser obtidos, por exemplo, pela aplicação da primeira ou segunda regra de Simpson repetida. Por exemplo, para uma integral na área de $X$ de 0 a 1 e $Y$ também de 0 a 1, se o eixo $X$ for dividido em 5 pontos (0, 1/4, 2/4, 3/4, 1) e o eixo $Y$ em 4 pontos (0, 1/3, 2/3, 1), podemos usar os pesos mostrados na tabela abaixo:

Pesos para o cálculo de $dV$:
x/y| 0 & 1/4 & 2/4 & 3/4 & 1
------------------------------------------
0  | 1 & 4 &   2 &   4 &   1
1/3| 3 & 12 &  6 &   12 &  3
2/3| 3 & 12 &  6 &   12 &  3
1  | 1 & 4 &   2 &   4 &   1 

O código abaixo é um exemplo de programa Python que calcula o volume da função $f(x,y) = \sin(3e^{(x^2 + y^2)/4})$ (ver figura inicial desta postagem).

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np

dx = 0.25 

dy = 0.25

X = np.arange(-5, 5+dx, dx)
Y = np.arange(-5, 5+dy, dy)
X, Y = np.meshgrid(X, Y)
R = np.exp((-X**2 - Y**2)/4.0)
Z = np.sin(3*R);
Zp = Z;


fig = plt.figure()
ax = fig.gca(projection='3d')

# Plotando a surperfície.

surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm,linewidth=0, antialiased=False)
plt.title('Função: $\sin(3e^{(x² + y²)/4})$');


# Adicionando cores e uma barra de cores

fig.colorbar(surf, shrink=0.5, aspect=5)

plt.show()

P = np.arange(-5, 5+dx, dx)
fim=len(P)-1print('tam = ',fim)
P[0] = 1;
n=1;
m = 0;

while fim>n:
    P[n] = 4;
    P[n+1]=2;
    n=n+2;

P[fim] = 1;
print(P)
n=0; m = 0; kk=1;
fim=fim+1
while fim>n:
    n=n+1;
    m=1;
    while fim>m:
        kk = P[n]*P[m]
        Zp[n,m] = Z[n,m]*kk
        m=m+1 

Zpsum = dx*dy*np.sum(np.sum(Zp))/9;
print('Volume = ',Zpsum)