Isso nos dá uma fórmula de iteração: cada ponto interior é atualizado como a média dos quatro vizinhos (esquema clássico de Jacobi/Gauss–Seidel). Na prática, Gauss–Seidel com relaxação (SOR) converge mais rápido: atualizamos direto na matriz e aplicamos um fator de relaxação \omega:
Explicando \omega: se \omega=1 é Gauss–Seidel simples; 1<\omega<2 acelera (over-relaxation) em muitos problemas ellipticos; escolha típica \omega\approx 1{,}5.
Passo 4
Como impor a condição de Neumann lateral (\partial C/\partial x=0) numa malha de diferenças finitas.
Derivada em x na fronteira esquerda i=0: \partial C/\partial x|_{x=0}=0.
Usando aproximação por diferença central com um ponto fantasma i=-1:
Implementação prática: durante a iteração atualize primeiro os nós interiores i=1,\dots,N_x-2, j=1,\dots,N_y-2 com a fórmula da média; depois atualize as bordas laterais com as fórmulas acima (ou trate pontos fantasma explicitamente). Não toque nas linhas de Dirichlet j=0 e j=N_y-1 — elas permanecem fixas.
Passo 5
Algoritmo iterativo completo (passo a passo) que vamos implementar:
Definir N_x,N_y,\Delta x,\Delta y. Inicializar matriz C com um chute inicial (por exemplo zeros ou uma rampa em y).
Impor cond. de contorno:
C_{i,0}=0 para todo i.
C_{i,N_y-1}=1 para todo i.
Iterar:
Para cada iteração:
Atualizar pontos interiores 1\le i\le N_x-2,\\\ 1\le j\le N_y-2 com Gauss–Seidel+relaxação.
Atualizar bordas laterais i=0 e i=N_x-1 usando as fórmulas com ponto fantasma (Neumann).
Reimpor as condições de Dirichlet (por segurança).
Calcular a máxima variação absoluta entre novo e antigo; se menor que a tolerância \text{tol}, pare.
Retornar C.
Critério de convergência: medir o máximo da mudança absoluta em uma iteração (norma L^\infty do resíduo de iteração) e comparar com \text{tol} (por exemplo 10^{-6}). Também limitar o número de iterações para evitar loops infinitos.
Passo 6
Dicas práticas e escolhas:
Malha: N_x=N_y=51 é um bom começo para testes (resolução razoável). Aumente se precisar de mais precisão.
Relaxação: \omega=1{,}5 costuma acelerar; experimente entre 1{,}0 e 1{,}9.
Inicialização: uma rampa linear em y (por exemplo C_{i,j}=y_j) dá um chute melhor que zeros, acelerando convergência.
Visualização: use contour/imshow/pcolormesh para inspecionar resultado.
Complexidade: cada iteração é O(N_x N_y). SOR reduz o número de iterações necessárias.
Passo 7
Agora que você entendeu a derivação, a imposição das condições de contorno e o algoritmo iterativo, aqui está o código Python completo que implementa o método descrito. Cole, rode e experimente os parâmetros (N_x, N_y, \omega, \text{tol}) para ver como muda a convergência.
Resposta
import numpy as np
import matplotlib.pyplot as plt
defsolve_laplace_neumann_sor(Nx=51, Ny=51, omega=1.5, tol=1e-6, maxiter=20000, verbose=True):
"""
Resolve Laplace 2D em [0,1]x[0,1] com:
C(x,1) = 1 (top, Dirichlet)
C(x,0) = 0 (bottom, Dirichlet)
dC/dx |_{x=0} = dC/dx |_{x=1} = 0 (laterais, Neumann)
usando diferenças finitas centradas e Gauss-Seidel com relaxação (SOR).
Retorna a matriz C (shape (Nx,Ny)), vetores x, y e número de iterações.
"""# Malha
x = np.linspace(0.0, 1.0, Nx)
y = np.linspace(0.0, 1.0, Ny)
dx = x[1] - x[0]
dy = y[1] - y[0]
ifabs(dx - dy) > 1e-12:
print("Aviso: dx != dy. O código assume dx ~ dy; resultados simplificados dependem disso.")
# Inicialização: chute inicial linear em y (ajuda na convergência)
C = np.zeros((Nx, Ny), dtype=float)
for j inrange(Ny):
C[:, j] = y[j] # rampa entre 0 e 1# Impor Dirichlet nas bordas superior e inferior
C[:, 0] = 0.0
C[:, -1] = 1.0# Iteração SOR (Gauss-Seidel com relaxação)
it = 0for it inrange(1, maxiter + 1):
max_diff = 0.0# Atualiza nós interiores (i = 1..Nx-2, j = 1..Ny-2)for i inrange(1, Nx - 1):
for j inrange(1, Ny - 1):
# média dos quatro vizinhos
new_val = 0.25 * (C[i+1, j] + C[i-1, j] + C[i, j+1] + C[i, j-1])
# SOR update
diff = new_val - C[i, j]
C[i, j] += omega * diff
ifabs(omega * diff) > max_diff:
max_diff = abs(omega * diff)
# Bordas laterais com condição Neumann (derivada em x = 0)# i = 0 (esquerda)
i = 0for j inrange(1, Ny - 1):
# usando C_{-1,j} = C_{1,j} => C_{0,j} = (2*C_{1,j} + C_{0,j+1} + C_{0,j-1})/4
new_val = 0.25 * (2.0 * C[1, j] + C[0, j+1] + C[0, j-1])
diff = new_val - C[0, j]
C[0, j] += omega * diff
ifabs(omega * diff) > max_diff:
max_diff = abs(omega * diff)
# i = Nx-1 (direita)
i = Nx - 1for j inrange(1, Ny - 1):
# usando C_{N_x,j} = C_{N_x-2,j} => C_{N_x-1,j} = (2*C_{N_x-2,j} + C_{N_x-1,j+1} + C_{N_x-1,j-1})/4
new_val = 0.25 * (2.0 * C[Nx - 2, j] + C[Nx - 1, j+1] + C[Nx - 1, j-1])
diff = new_val - C[Nx - 1, j]
C[Nx - 1, j] += omega * diff
ifabs(omega * diff) > max_diff:
max_diff = abs(omega * diff)
# Reimpor Dirichlet nas bordas superior e inferior por segurança
C[:, 0] = 0.0
C[:, -1] = 1.0# Verifica convergênciaif max_diff < tol:
if verbose:
print(f"Convergência alcançada em {it} iterações (max_diff = {max_diff:.2e}).")
breakif it % 500 == 0and verbose:
print(f"Iter {it}: max_diff = {max_diff:.2e}")
else:
# loop terminou sem breakif verbose:
print(f"Atenção: máximo de iterações ({maxiter}) alcançado. max_diff = {max_diff:.2e}")
return C, x, y, it
if __name__ == "__main__":
# Parâmetros (experimente Nx, Ny, omega)
Nx, Ny = 61, 61
omega = 1.5
tol = 1e-6
C, x, y, iters = solve_laplace_neumann_sor(Nx=Nx, Ny=Ny, omega=omega, tol=tol, verbose=True)
# Plot do resultado
X, Y = np.meshgrid(x, y, indexing='ij') # indexing='ij' para combinar com C[i,j]
plt.figure(figsize=(6,5))
pcm = plt.pcolormesh(X, Y, C, shading='auto', cmap='viridis')
plt.colorbar(pcm, label='C(x,y)')
plt.contour(X, Y, C, levels=10, colors='white', linewidths=0.5, alpha=0.7)
plt.xlabel('x')
plt.ylabel('y')
plt.title(f'Solução aproximada de Laplace (it={iters}, Nx={Nx}, Ny={Ny})')
plt.gca().set_aspect('equal', 'box')
plt.show()