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:
v(x, y)=\sum_{\substack{n=1,3,5 \ldots \\\ \text { (impar) }}}^{\infty} \frac{4}{n^{3} \pi^{3}}\left(1-\frac{\operatorname{cosh} (n \pi y)}{\operatorname{cosh} \left(\frac{n \pi}{2}\right)}\right) \operatorname{sen}(n \pi x)
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:
  1. (a) Fazer um gráfico da solução analítica dada pela série
v(x,y)=\sum_{\substack{n=1,3,5,\ldots \\\ (\text{ímpar})}}^{\infty} \frac{4}{n^{3}\pi^{3}}\Bigg(1-\frac{\operatorname{cosh}(n\pi y)}{\operatorname{cosh}\big(\frac{n\pi}{2}\big)}\Bigg)\sin(n\pi x)í
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).
  • Domínio retangular: x\in[0,1], y\in[-\tfrac{1}{2},\tfrac{1}{2}].
  • Condição de contorno: v=0 em todo o contorno.
  • 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á
h_x=\frac{1-0}{N_x-1},\qquad h_y=\frac{\tfrac{1}{2}-(-\tfrac{1}{2})}{N_y-1}.
  • A equação discreta (para nós do interior) com h_x=h_y=h é, por diferenças centrais,
\frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{h^{2}}+\frac{v_{i,j+1}-2v_{i,j}+v_{i,j-1}}{h^{2}}=-1.
Multiplicando por h^{2} e rearranjando:
- v_{i-1,j} - v_{i+1,j} - v_{i,j-1} - v_{i,j+1} + 4 v_{i,j} = h^{2}.
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).
  1. Entendimento da série: cada termo tem forma
a_n(x,y)=\frac{4}{n^{3}\pi^{3}}\Bigg(1-\frac{\operatorname{cosh}(n\pi y)}{\operatorname{cosh}\big(\frac{n\pi}{2}\big)}\Bigg)\sin(n\pi x),
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.
  1. 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}}.
  2. 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.
  1. 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_terms

for 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).
  1. 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.
  2. Para pontos internos (não na fronteira) a aproximação por diferenças centrais da segunda derivada é:
\frac{\partial^{2}v}{\partial x^{2}}\Big|_{i,j}\approx\frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{h_x^{2}},
\frac{\partial^{2}v}{\partial y^{2}}\Big|_{i,j}\approx\frac{v_{i,j+1}-2v_{i,j}+v_{i,j-1}}{h_y^{2}}.
  1. Somando e igualando a -1 (o lado direito da Poisson):
\frac{v_{i+1,j}-2v_{i,j}+v_{i-1,j}}{h_x^{2}}+\frac{v_{i,j+1}-2v_{i,j}+v_{i,j-1}}{h_y^{2}}=-1.
  1. 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:
- v_{i-1,j} - v_{i+1,j} - v_{i,j-1} - v_{i,j+1} + 4 v_{i,j} = h^{2}.
  1. 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.
  2. 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 interiores
for j in range(n_y):
    for i in range(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).
  1. 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.
  1. Métricas de erro (definição e interpretação):
  • Erro máximo (norma L^\infty):
\text{err}_{\infty} = \max_{i,j} |v_{\text{analytic}}(x_i,y_j) - v_{\text{num}}(x_i,y_j)|.
  • Erro RMS (norma discreta L^{2} / RMS):
\text{err}_{2} = \sqrt{\frac{1}{N_{\text{pt}}}\sum_{i,j} \big(v_{\text{analytic}}(x_i,y_j)-v_{\text{num}}(x_i,y_j)\big)^{2}},
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.

Exercícios relacionados