O perfil de velocidade do escoamento de um fluido em um duto quadrado é dado pela equação de Poisson:
\nabla^{2} v=-1
Sujeito à condição de contorno de não-deslizamento em todas as fronteiras (v=0). A solução analítica deste problema é obtida através da série infinita abaixo:
a) Faça um gráfico da solução analítica em Python usando o domínio 0 \leq x \leq 1 e -1 / 2 \leq y \leq 1 / 2.
b) Resolva a equação de Poisson usando método das diferenças finitas centradas. Você pode construir uma malha de discretização com 100 pontos igualmente espaçados para cada direção.
c) Compare as duas soluções (analítica e numérica) deste problema.
Passo 1
O enunciado pede três coisas:
(a) Fazer um gráfico da solução analítica dada pela série
no domínio 0\le x\le 1 e -\tfrac{1}{2}\le y\le \tfrac{1}{2}.
2. (b) Resolver numericamente a equação de Poisson
\nabla^{2} v = -1
com cond. de contorno de não-deslizamento v=0 em todas as fronteiras, usando diferenças finitas centradas em uma malha com 100 pontos igualmente espaçados em cada direção.
3. (c) Comparar a solução analítica com a solução numérica (gráficos e métricas de erro).
Informações relevantes dadas:
A solução analítica é a série acima (somente termos ímpares n=1,3,5,\dots).
A malha pedida para a FDM tem 100 pontos em x e 100 pontos em y (pontos igualmente espaçados).
Passo 2
Escolha do método e resumo da estratégia (curto e prático):
Para (a) vou truncar a série analítica e calcular a soma numérica em todos os pontos da malha desejada (vetorizando em NumPy).
Para (b) vou usar o esquema de diferenças finitas centradas de segunda ordem (5-point stencil) com passo uniforme em x e em y. Como o domínio tem tamanho 1 em x e 1 em y (pois y vai de -1/2 a 1/2), se usarmos N_x=N_y=100 pontos, o espaçamento será
Ou seja, a matriz do sistema é a clássica matriz do stencil 5-pontos e o vetor direito tem entrada h^{2} para cada nó interior (ajustada se existirem termos de contorno não-nulos — aqui contorno é zero, então sem termos extras).
Resolveremos o sistema linear esparso com um solver direto esparso (por exemplo, scipy.sparse.linalg.spsolve) ou iterativo (conjugate gradient para simétrica definida positiva — aqui a matriz é SPD).
Passo 3
Implementação da solução analítica (truncamento da série e detalhes didáticos).
onde n percorre números ímpares 1,3,5,\dots. Note que:
\operatorname{cosh} é a função cosseno hiperbólica.
O fator \dfrac{4}{n^{3}\pi^{3}} faz com que os termos decaiam como 1/n^{3} — isso garante convergência rápida.
Truncamento prático: escolhemos um número de termos ímpares N_{\text{terms}}. Como a série decai como 1/n^{3}, usar N_{\text{terms}}=101 (ou seja, somar os primeiros 101 termos ímpares: n=1,3,\dots,201) costuma dar ótima aproximação. Se quiser maior precisão, aumente N_{\text{terms}}.
Implementação vetorizada em Python (explicação de cada passo no código):
Cria-se uma malha regular x e y usando np.linspace com N_x e N_y.
Monta-se um array 2D para v_{\text{analytic}} e soma-se cada termo da série.
Código Python (cole e execute no seu ambiente; eu coloquei comentários explicativos):
import numpy as np
import matplotlib.pyplot as plt
# Parâmetros da malha
Nx = 100# pontos em x
Ny = 100# pontos em y
x = np.linspace(0.0, 1.0, Nx)
y = np.linspace(-0.5, 0.5, Ny)
X, Y = np.meshgrid(x, y, indexing='xy') # X.shape = (Ny, Nx)# Soma truncada da série analítica
N_terms = 101# número de termos ímpares (ajuste se quiser mais precisão)
v_analytic = np.zeros_like(X)
# vetor de n ímpares: 1,3,5,...
ns = np.arange(1, 2*N_terms, 2) # comprimento N_termsfor n in ns:
coef = 4.0 / ( (n**3) * (np.pi**3) )
term = coef * (1.0 - np.cosh(n * np.pi * Y) / np.cosh(n * np.pi / 2.0)) * np.sin(n * np.pi * X)
v_analytic += term
Observações:
Usei indexing='xy' para que X varie em x (colunas) e Y em y (linhas), o que facilita o plot.
Se o computador demorar, reduza N_terms ou aumente com paciência; a soma é vetorizada, então é razoavelmente rápida.
Passo 4
Fazer o gráfico da solução analítica (dicas de plot).
Sugestões de visualização:
Um mapa de contorno (contourf) ou imagem (imshow) para ver o campo v(x,y) em cores.
Uma superfície 3D (Axes3D) para visualização 3D.
Código para dois tipos de gráfico:
# 1) Contourf
plt.figure(figsize=(6,5))
cs = plt.contourf(X, Y, v_analytic, levels=50, cmap='viridis')
plt.colorbar(cs, label='v (analítica)')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solução analítica (série truncada)')
plt.tight_layout()
plt.show()
# 2) Superfície 3D (opcional)from mpl_toolkits.mplot3d import Axes3D # necessário para algumas versões do matplotlib
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, v_analytic, cmap='viridis', linewidth=0, antialiased=True)
ax.set_xlabel('x'); ax.set_ylabel('y'); ax.set_zlabel('v (analítica)')
ax.set_title('Solução analítica (superfície)')
plt.tight_layout()
plt.show()
Incentivo: roda esse código agora — dá uma satisfação ver o perfil de velocidade aparecendo.
Passo 5
Derivação do esquema de diferenças finitas (passo-a-passo, didático).
Discretização: Seja a malha com N_x pontos em x e N_y em y. Denotamos v_{i,j} o valor aproximado no ponto (x_i,y_j) com i=0,\dots,N_x-1 e j=0,\dots,N_y-1.
Para pontos internos (não na fronteira) a aproximação por diferenças centrais da segunda derivada é:
Para malha uniforme com h_x=h_y=h (aqui, aproximadamente, já que os comprimentos em x e y são ambos 1 — exatamente h_x=1/(N_x-1) e h_y=1/(N_y-1), com N_x=N_y esses são iguais), multiplicamos por h^{2} e obtemos o esquema 5-pontos:
No lado direito deve-se colocar h^{2} porque o termo fonte é constante -1 e já multiplicamos por h^{2}. Se os contornos fossem não nulos, teríamos termos adicionais no lado direito, mas aqui v=0 nas bordas, então nenhuma correção extra é necessária.
Montamos o sistema linear A u = b, onde u é o vetor com os valores v_{i,j} apenas nos nós interiores (ordem por linhas ou colunas), A é a matriz esparsa do stencil e b é um vetor com todas as entradas iguais a h^{2}.
Passo 6
Montagem e solução do sistema linear em Python (implementação prática com scipy.sparse).
Explicação do mapeamento: suponha que temos N_x=100 e N_y=100. O número de nós interiores em cada direção é n_x = N_x-2 e n_y = N_y-2. O número total de incógnitas é N = n_x \times n_y. Usamos mapeamento p = i + j\cdot n_x (com i=0,\dots,n_x-1 e j=0,\dots,n_y-1) para transformar índices 2D em 1D.
Código (usando matriz laplaciana via kronecker ou diags):
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
# Malha (repetindo Nx, Ny do antes)
Nx = 100
Ny = 100
hx = (1.0 - 0.0) / (Nx - 1)
hy = (0.5 - (-0.5)) / (Ny - 1)
# Aqui Nx == Ny, então hx ~ hy; usaremos h = hx = hy
h = hx
n_x = Nx - 2# nós internos em x
n_y = Ny - 2# nós internos em y
N = n_x * n_y
# Construção da matriz 2D Laplaciana (5-pontos) com h uniforme# Matriz 1D tri-diagonal (para uma linha):
main_1d = 4.0 * np.ones(n_x)
off_1d = -1.0 * np.ones(n_x - 1)
T = sp.diags([off_1d, main_1d, off_1d], offsets=[-1, 0, 1], format='csr')
I = sp.eye(n_x, format='csr')
A = sp.kron(sp.eye(n_y, format='csr'), T, format='csr') + sp.kron(sp.diags([-1.0*np.ones(n_y-1), -1.0*np.ones(n_y-1)], offsets=[-1,1], format='csr'), I, format='csr')
# Vetor b: para cada nó interior, b_i = h^2 (porque -1 * h^2 e trocamos de sinal no esquema)
b = (h**2) * np.ones(N)
# Resolver o sistema esparso
u = spla.spsolve(A, b) # retorna vetor de tamanho N# Reconstruir solução 2D incluindo fronteiras zeros
Vnum = np.zeros((Ny, Nx)) # incluir contorno# Colocar valores interioresfor j inrange(n_y):
for i inrange(n_x):
p = i + j * n_x
Vnum[j+1, i+1] = u[p] # j+1 e i+1 porque 0 e Nx-1 são fronteiras
Observações práticas:
A matriz A foi construída com o padrão de coeficientes -1 nas direções vizinhas e 4 no centro (como no esquema deduzido).
Como as fronteiras são zero, não há termos extras no vetor b.
spsolve é apropriado para problemas desta dimensão (aqui N\approx 98\times98\approx9604, o que é pequeno). Para malhas muito maiores, usar métodos iterativos com precondicionador.
Passo 7
Comparação entre analítica e numérica (gráficos e métricas).
Avaliação da solução analítica nos mesmos pontos da malha usada na solução numérica:
Já construímos v_analytic em X, Y (Passo 3) com dimensões (Ny, Nx).
A solução numérica Vnum também tem dimensão (Ny, Nx) com fronteiras zero.
onde N_{\text{pt}}=N_x\cdot N_y.
3. Código para calcular métricas e plotar a diferença:
# Calcula erro ponto a ponto
diff = v_analytic - Vnum
err_inf = np.max(np.abs(diff))
err_rms = np.sqrt(np.mean(diff.ravel()**2))
print("Erro máximo (L_inf):", err_inf)
print("Erro RMS (L2):", err_rms)
# Plot da diferença
plt.figure(figsize=(6,5))
cs = plt.contourf(X, Y, diff, levels=50, cmap='seismic')
plt.colorbar(cs, label='v_analytic - v_num')
plt.xlabel('x'); plt.ylabel('y')
plt.title('Diferença entre solução analítica e numérica')
plt.tight_layout()
plt.show()
Interpretação:
Se o erro for pequeno (por exemplo, ordem 10^{-3} ou menor dependendo do truncamento e malha), então a solução numérica concorda bem com a analítica.
Se quiser reduzir erro numérico, aumentar N_x,N_y ou usar mais termos N_{\text{terms}} na série analítica ajuda. A ordem de convergência deste esquema FDM é O(h^{2}) (porque usamos diferenças centradas de segunda ordem).
Passo 8
Dicas finais, verificação e melhorias:
Verifique que a fronteira em Vnum esteja realmente zero (condição de contorno). Se ver valores não-nulos na borda, algo na montagem do sistema está errado.
Para checar convergência do truncamento da série: calcule v_analytic com N_terms=51,101,201 e veja a diferença entre essas aproximações; se a diferença for pequena, truncamento é suficiente.
Para checar convergência da FDM: resolva com N_x=N_y=50,100,200 e observe como err_inf decresce — espera-se queda proporcional a h^{2}.
Se preferir resolver o sistema iterativamente (escalável), use spla.cg ou spla.bicgstab com precondicionador spla.spilu (para problemas maiores).
Resposta
a) Código Python para gerar o gráfico da solução analítica foi fornecido nos Passos 3 e 4; execute esse código para obter o gráfico.
b) Código Python para resolver a equação de Poisson por diferenças finitas centradas (malha 100\times100) foi fornecido no Passo 6; execute esse código para obter a solução numérica.
c) A comparação (plots da diferença e métricas de erro \text{err}_{\infty} e \text{err}_{2}) está implementada no Passo 7; execute aquele bloco para imprimir os valores de erro e visualizar o mapa de diferença.