 |
| 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)
