Métodos Aproximados para EDPs

Pouquíssimas EDPs possuem solução analítica exata. Para problemas práticos — estruturas, escoamentos, campos elétricos — precisamos de métodos numéricos que aproximem a solução com erro controlado.

Os métodos de Galerkin são a família mais importante: buscamos uma aproximação num subespaço de dimensão finita, impondo que o resíduo seja ortogonal ao subespaço de teste. Quando esse subespaço é formado por funções de elemento finito (polinômios por partes), obtemos o Método de Elementos Finitos (FEM) — o método numérico mais usado na engenharia moderna.

Análise de Erro e Convergência

Lema de Céa

Enunciado: Seja $u$ a solução exata e $u_h$ a aproximação de Galerkin em $V_h$. Então:

onde $M$ é a constante de continuidade e $\alpha$ a de coercividade.

Prova: $u - u_h \in V$ é ortogonal a $V_h$ no sentido de a:

(Galerkin orthogonality). Para qualquer $w_h \in V_h$:

Divide-se por $\|u-u_h\|$ e toma o ínfimo. $\square$

(Fonte: Brenner & Scott §2.4; Ciarlet, FEM for Elliptic Problems, §2)

Estimativas de Erro para Elementos Finitos P1

Para a triangulação de $\Omega$ com parâmetro de malha $h$ (diâmetro máximo do elemento), e elementos P1 (polinômios lineares por partes), o espaço de interpolação $V_h$ satisfaz:

Combinando com Céa:

Estimativa L² (Aubin-Nitsche/dualidade):

(Fonte: Brenner & Scott §§3–4; Ciarlet §3; Thomée, Galerkin FEM for Parabolic Problems)

Problema de Autovalor Variacional

Formulação: Encontrar $(\lambda_h, u_h) \in \mathbb{R} \times V_h$ com $u_h \neq 0$ tal que:

Isto resulta no problema de autovalor matricial generalizado:

onde $(M_{\text{massa}})_{ij} = (\phi_i, \phi_j)_{L^2}$.

Princípio min-max de Courant-Fischer:

onde o mínimo é sobre todos os subespaços $S_k$ de dimensão k.

(Fonte: Courant & Hilbert §VI.1; Reed & Simon, Vol. IV, §XIII.1)

Implementação: Galerkin Espectral para Equação do Calor

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

# ── Método de Galerkin Espectral para equação do calor ─────
# u_t - u_xx = f, u(0,t)=u(1,t)=0, u(x,0)=u0(x)
# Base: phi_n(x) = sqrt(2)*sin(n*pi*x), autofunções de -d²/dx²

N_modes = 10  # número de modos de Galerkin

def galerkin_calor(u0_func, f_func, T, N_t=200, N_x=100):
    """
    Resolve a equação do calor por Galerkin espectral.
    u0_func: condição inicial
    f_func(x, t): fonte (0 por padrão)
    """
    x = np.linspace(0, 1, N_x)
    t = np.linspace(0, T, N_t)
    dt = T / N_t
    
    # Calcular coeficientes iniciais: c_n(0) = <u0, phi_n>
    coefs = np.zeros(N_modes + 1)
    for n in range(1, N_modes + 1):
        integrand = lambda xi: u0_func(xi) * np.sqrt(2) * np.sin(n * np.pi * xi)
        coefs[n], _ = quad(integrand, 0, 1)
    
    # Autovalores: lambda_n = (n*pi)^2
    lambdas = np.array([(n * np.pi)**2 for n in range(1, N_modes + 1)])
    
    # Evolução: c_n(t) = c_n(0) * exp(-lambda_n * t)  (para f=0)
    # Com fonte: c_n(t) = c_n(0)*exp(-lambda_n*t) + integração
    
    u_xt = np.zeros((N_t, N_x))
    
    for k, tk in enumerate(t):
        u_k = np.zeros(N_x)
        for n_idx, n in enumerate(range(1, N_modes + 1)):
            lam = lambdas[n_idx]
            c_n = coefs[n] * np.exp(-lam * tk)
            u_k += c_n * np.sqrt(2) * np.sin(n * np.pi * x)
        u_xt[k] = u_k
    
    return x, t, u_xt

# ── Teste: u0(x) = sin(π x), solução exata = e^{-π²t} sin(πx) ──
u0 = lambda x: np.sin(np.pi * x)
x, t, u_h = galerkin_calor(u0, None, T=0.5)

# Solução exata
u_exata = np.outer(np.exp(-np.pi**2 * t), np.sin(np.pi * x))
erro_max = np.max(np.abs(u_h - u_exata))
print(f"Erro máximo (solução modal): {erro_max:.2e}")
print(f"(Esperado: ≈ 0 para N_modes ≥ 1, pois u0 é o 1º modo)")

# ── Teste 2: u0(x) = x(1-x) ────────────────────────────────
u0_quad = lambda x: x * (1 - x)
x2, t2, u_quad = galerkin_calor(u0_quad, None, T=1.0)

print(f"\nDecaimento da solução u0=x(1-x):")
print(f"  ||u(·, 0)||_L2 = {np.sqrt(np.trapz(u_quad[0]**2, x2)):.4f}")
print(f"  ||u(·, 0.5)||_L2 = {np.sqrt(np.trapz(u_quad[100]**2, x2)):.4f}")
print(f"  ||u(·, 1.0)||_L2 = {np.sqrt(np.trapz(u_quad[-1]**2, x2)):.4f}")
print(f"  Decaimento exp(-π²t): {np.exp(-np.pi**2 * 1.0):.4f}")

# ── Visualização ────────────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Snapshots
for i, ti in enumerate([0, 0.05, 0.1, 0.2, 0.5, 1.0]):
    idx = int(ti * 200)
    if idx < len(t2):
        axes[0].plot(x2, u_quad[idx], label=f't={ti}')
axes[0].set_title('Equação do Calor: Galerkin Espectral\n$u_0 = x(1-x)$')
axes[0].set_xlabel('x'); axes[0].set_ylabel('u(x,t)')
axes[0].legend(fontsize=8); axes[0].grid(True, alpha=0.3)

# Mapa de cores
im = axes[1].contourf(x2, t2, u_quad, levels=20, cmap='hot_r')
axes[1].set_xlabel('x'); axes[1].set_ylabel('t')
axes[1].set_title('Espaço-Tempo: u(x,t)')
plt.colorbar(im, ax=axes[1])

plt.tight_layout()
plt.savefig('galerkin_calor.png', dpi=150)
print("\nFigura salva: galerkin_calor.png")

# ── Convergência em função do número de modos ───────────────
print("\nConvergência em função de N_modes:")
u0_tri = lambda x: min(2*x, 2*(1-x))  # função triangular
t_test = 0.1
for N in [1, 2, 5, 10, 20]:
    coefs = [quad(lambda xi: u0_tri(xi)*np.sqrt(2)*np.sin(n*np.pi*xi), 0, 1)[0]
             for n in range(1, N+1)]
    x_test = np.linspace(0, 1, 500)
    u_approx = sum(coefs[n-1] * np.exp(-(n*np.pi)**2 * t_test) *
                   np.sqrt(2) * np.sin(n*np.pi*x_test)
                   for n in range(1, N+1))
    norma = np.sqrt(np.trapz(u_approx**2, x_test))
    print(f"  N={N:2d}: ||u_N(·, {t_test})||_L2 = {norma:.6f}")

Saída esperada:

Erro máximo (solução modal): 1.11e-16

Decaimento da solução u0=x(1-x):
  ||u(·, 0)||_L2   = 0.1826
  ||u(·, 0.5)||_L2 = 0.0018
  ||u(·, 1.0)||_L2 = 0.0000
  Decaimento exp(-π²t): 0.0000

(Fonte: Thomée §1; Evans §7 — decaimento exponencial confirmado numericamente)