Escreva um código em Python usando diferenças finitas centradas para resolver a equação da difusão:
\frac{\partial^{2} C}{\partial x^{2}}+\frac{\partial^{2} C}{\partial y^{2}}=0
com as seguintes condições de contorno:
\begin{gathered} C(x, 1)=1 \\\ C(x, 0)=0 \\\ \left.\frac{d C}{d x}\right|_{x=0}=\left.\frac{d C}{d x}\right|_{x=1}=0 \end{gathered}

Passo 1

O problema pede resolver a equação
\frac{\partial^{2} C}{\partial x^{2}}+\frac{\partial^{2} C}{\partial y^{2}}=0
no quadrado [0,1]\times[0,1] com as condições de contorno:
\begin{gathered} C(x,1)=1,\\\ C(x,0)=0,\\\ \left.\frac{d C}{d x}\right|_{x=0}=\left.\frac{d C}{d x}\right|_{x=1}=0. \end{gathered}
Lista clara do que temos e do que se pede:
  • Domínio: [0,1]\times[0,1].
  • PDE: laplaciana nula (equação de Laplace em 2D).
  • Condições:
    • topo ( y=1): C=1 (Dirichlet).
    • base ( y=0): C=0 (Dirichlet).
    • laterais ( x=0 e x=1): derivada normal zero, ou seja Neumann com \partial C/\partial x=0.
  • Objetivo prático: construir e implementar um esquema de diferenças finitas centradas para aproximar C(x,y) no domínio.

Passo 2

Discretização do domínio (isso é essencial antes de escrever código).
  • Vamos dividir x em N_x pontos e y em N_y pontos igualmente espaçados.
  • Chamamos \Delta x = \frac{1}{N_x-1} e \Delta y = \frac{1}{N_y-1}.
  • Índices discretos: i=0,\dots,N_x-1 representa x e j=0,\dots,N_y-1 representa y.
  • A solução numérica será armazenada numa matriz C_{i,j} que aproxima C(x_i,y_j).
  • Ponto importante: as condições de contorno serão impostas diretamente nos índices das bordas:
    • linha j=0 (base) recebe C_{i,0}=0.
    • linha j=N_y-1 (topo) recebe C_{i,N_y-1}=1.
    • bordas laterais i=0 e i=N_x-1 têm condição Neumann (ver passo seguinte como implementar).

Passo 3

Aproximações por diferenças finitas centradas (a matemática da discretização).
  • Segunda derivada em x com diferença central:
\frac{\partial^2 C}{\partial x^2}\bigg|_{i,j}\approx \frac{C_{i+1,j}-2C_{i,j}+C_{i-1,j}}{\Delta x^2}.
  • Segunda derivada em y com diferença central:
\frac{\partial^2 C}{\partial y^2}\bigg|_{i,j}\approx \frac{C_{i,j+1}-2C_{i,j}+C_{i,j-1}}{\Delta y^2}.
  • Substituindo na equação \partial^2_x C+\partial^2_y C=0 obtemos a equação discreta:
\frac{C_{i+1,j}-2C_{i,j}+C_{i-1,j}}{\Delta x^2}+ \frac{C_{i,j+1}-2C_{i,j}+C_{i,j-1}}{\Delta y^2}=0.
  • Em malha uniforme com \Delta x=\Delta y (caso comum), isso simplifica para a média dos vizinhos:
C_{i,j}=\frac{1}{4}\bigl(C_{i+1,j}+C_{i-1,j}+C_{i,j+1}+C_{i,j-1}\bigr).
  • 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:
C_{i,j}^{\text{novo}}=(1-\omega)C_{i,j}^{\text{antigo}}+\omega\cdot\frac{C_{i+1,j}+C_{i-1,j}+C_{i,j+1}+C_{i,j-1}}{4}.
  • 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:
\left.\frac{\partial C}{\partial x}\right|_{i=0}\approx\frac{C_{1,j}-C_{-1,j}}{2\Delta x}=0 \Rightarrow C_{-1,j}=C_{1,j}.
  • Substituindo o ponto fantasma na fórmula de Laplace para i=0 obtemos:
C_{0,j}=\frac{1}{4}\bigl(C_{1,j}+C_{-1,j}+C_{0,j+1}+C_{0,j-1}\bigr) =\frac{1}{4}\bigl(2C_{1,j}+C_{0,j+1}+C_{0,j-1}\bigr).
  • Analogamente, na fronteira direita i=N_x-1:
C_{N_x-1,j}=\frac{1}{4}\bigl(2C_{N_x-2,j}+C_{N_x-1,j+1}+C_{N_x-1,j-1}\bigr).
  • 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:
  1. Definir N_x,N_y,\Delta x,\Delta y. Inicializar matriz C com um chute inicial (por exemplo zeros ou uma rampa em y).
  2. Impor cond. de contorno:
    • C_{i,0}=0 para todo i.
    • C_{i,N_y-1}=1 para todo i.
  3. 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.
  4. 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

def solve_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]
    if abs(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 in range(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 = 0
    for it in range(1, maxiter + 1):
        max_diff = 0.0

        # Atualiza nós interiores (i = 1..Nx-2, j = 1..Ny-2)
        for i in range(1, Nx - 1):
            for j in range(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
                if abs(omega * diff) > max_diff:
                    max_diff = abs(omega * diff)

        # Bordas laterais com condição Neumann (derivada em x = 0)
        # i = 0 (esquerda)
        i = 0
        for j in range(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
            if abs(omega * diff) > max_diff:
                max_diff = abs(omega * diff)

        # i = Nx-1 (direita)
        i = Nx - 1
        for j in range(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
            if abs(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ência
        if max_diff < tol:
            if verbose:
                print(f"Convergência alcançada em {it} iterações (max_diff = {max_diff:.2e}).")
            break
        if it % 500 == 0 and verbose:
            print(f"Iter {it}: max_diff = {max_diff:.2e}")

    else:
        # loop terminou sem break
        if 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()

Exercícios relacionados