📊 Métodos Numéricos para Equações Diferenciais Parciais

www.metodosnumericos.com.br --Prof. Doherty Andrade-- Visualização Interativa Completa: Superfícies 3D Animadas com Fatias

📚 Teoria Geral
🌊 Equação da Onda
🔥 Equação do Calor
🎯 Poisson
⭕ Laplace

📖 Fundamentação Teórica das EDPs

1. Classificação das EDPs de Segunda Ordem

Uma equação diferencial parcial linear de segunda ordem a duas variáveis tem a forma:

$$ Au_{xx} + Bu_{xy} + Cu_{yy} + Du_x + Eu_y + Fu = G $$

O discriminante $\Delta = B^2 - 4AC$ determina o tipo da EDP:

DiscriminanteTipoExemploCaracterísticas
$\Delta > 0$Hiperbólica$u_{tt} - c^2 u_{xx} = 0$Onda, propagação
$\Delta = 0$Parabólica$u_t - \alpha^2 u_{xx} = 0$Calor, difusão
$\Delta < 0$Elíptica$u_{xx} + u_{yy} = 0$Laplace, equilíbrio

2. Condições de Contorno e Iniciais

3. Visualização 3D: Superfícies e Fatias

💡 Conceito Fundamental: A solução de uma EDP evolutiva (calor, onda) não é uma única curva 2D, mas sim uma superfície no espaço tridimensional $(x, t, u)$. As curvas que vemos em livros são apenas fatias (seções transversais) desta superfície em instantes de tempo fixos. A reunião de todas essas fatias gera a superfície completa. Neste texto ilustramos as superficies soluções, as fatias e os métodos numéricos utilizados nas soluções.

4. Métodos Numéricos Utilizados

🌊 Equação da Onda Unidimensional

Fundamentação Teórica

$$ \frac{\partial^2 u}{\partial t^2} = \alpha^2 \frac{\partial^2 u}{\partial x^2}, \quad 0 < x < l, \; 0 < t < T $$

com duas condições iniciais (posição e velocidade):

$$ u(x,0) = f(x), \quad \frac{\partial u}{\partial t}(x,0) = g(x) $$
🌟 Caso Especial $\lambda = 1$ (Solução Exata):

Quando $\lambda = \frac{\alpha k}{h} = 1$, o esquema produz a solução exata da equação da onda. Energia perfeitamente conservada.

⚡ Conservação de Energia: $$ E(t) = \frac{1}{2}\int_0^l \left[\left(\frac{\partial u}{\partial t}\right)^2 + \alpha^2\left(\frac{\partial u}{\partial x}\right)^2\right] dx = \text{constante} $$

🔬 Formulação Numérica (Diferenças Finitas Explícitas)

Malha: uniforme com $h = l/m$ (espaço) e $k = T/N$ (tempo). Denotamos $u_i^j \approx u(x_i, t_j)$ com $x_i = i h$ e $t_j = j k$.

Aplicando diferenças centrais em ambas as derivadas:

$$ \frac{u_i^{j+1} - 2u_i^{j} + u_i^{j-1}}{k^2} = \alpha^2 \frac{u_{i+1}^{j} - 2u_i^{j} + u_{i-1}^{j}}{h^2} $$

Definindo o número de Courant $\lambda = \dfrac{\alpha k}{h}$ e isolando $u_i^{j+1}$:

$$ \boxed{\; u_i^{j+1} = 2(1 - \lambda^2)\, u_i^{j} + \lambda^2\left(u_{i+1}^{j} + u_{i-1}^{j}\right) - u_i^{j-1} \;} $$

Primeiro passo ($j=0 \to 1$): como $u^{-1}$ não existe, usamos $u_t(x,0) = g(x)$ e uma aproximação de diferença central para obter:

$$ u_i^{1} = (1 - \lambda^2)\, f(x_i) + \frac{\lambda^2}{2}\left(f(x_{i+1}) + f(x_{i-1})\right) + k\, g(x_i) $$

Condições de contorno: $u_0^j = u_m^j = 0$ (extremidades fixas).

Estabilidade (CFL): o esquema é estável se $\lambda \le 1$. Quando $\lambda = 1$, o esquema é exato — não há erro de truncamento de fase.

Sistema linear? Nenhum — o esquema é explícito: cada $u_i^{j+1}$ é obtido por substituição direta dos níveis $j$ e $j-1$. Não há sistema a resolver, e a solução avança no tempo por simples substituição.

📘 Exemplo 1 (λ=1)
📗 Exemplo 2 (λ=0.5)
📙 Exemplo 3 (Viajante)

📘 Onda Estacionária com λ = 1 Solução Exata

$f(x) = \sin(\pi x)$, $g(x) = 0$, $u(x,t) = \sin(\pi x)\cos(\pi t)$

📗 Onda Estacionária com λ = 0.5 λ=0.5

$f(x) = \sin(\pi x)$, $g(x) = 0$, $m = 20$, $N = 20$

📙 Onda Viajante com λ = 1 Onda Viajante

$f(x) = \sin(2\pi x)$, $g(x) = -2\pi\sin(2\pi x)$, $u(x,t) = \sin(2\pi(x-t))$

🔥 Equação do Calor Unidimensional

Fundamentação Teórica

$$ \frac{\partial u}{\partial t} = \alpha^2 \frac{\partial^2 u}{\partial x^2}, \quad 0 < x < l, \; 0 < t < T $$

com condição inicial $u(x,0) = f(x)$ e condições de contorno $u(0,t) = u(l,t) = 0$.

🔥 Método de Crank-Nicolson:

Combina os esquemas explícito e implícito através de uma média. Possui precisão $O(k^2 + h^2)$ e é incondicionalmente estável.

🔬 Formulação Numérica (Crank-Nicolson + Algoritmo de Thomas)

Malha: $h = l/m$, $k = T/N$, com $u_i^j \approx u(x_i, t_j)$. Definimos o parâmetro adimensional $\lambda = \dfrac{\alpha^2 k}{h^2}$.

Ideia do Crank-Nicolson: avaliar a segunda derivada espacial como média entre os níveis $j$ (explícito) e $j+1$ (implícito):

$$ \frac{u_i^{j+1} - u_i^j}{k} = \frac{\alpha^2}{2}\left[ \frac{u_{i+1}^{j+1} - 2u_i^{j+1} + u_{i-1}^{j+1}}{h^2} + \frac{u_{i+1}^{j} - 2u_i^{j} + u_{i-1}^{j}}{h^2} \right] $$

Multiplicando por $k$ e reorganizando para separar incógnitas ($j+1$) e conhecidos ($j$):

$$ -\lambda\, u_{i-1}^{j+1} + 2(1+\lambda)\, u_i^{j+1} - \lambda\, u_{i+1}^{j+1} = \underbrace{\lambda u_{i-1}^{j} + 2(1-\lambda)u_i^{j} + \lambda u_{i+1}^{j}}_{d_i} $$

Sistema tridiagonal para cada passo de tempo (incógnitas $u_1^{j+1},\dots,u_{m-1}^{j+1}$):

$$ \begin{pmatrix} b & c & & \\ a & b & c & \\ & \ddots & \ddots & \ddots \\ & & a & b \end{pmatrix} \begin{pmatrix} u_1^{j+1} \\ u_2^{j+1} \\ \vdots \\ u_{m-1}^{j+1} \end{pmatrix} = \begin{pmatrix} d_1 \\ d_2 \\ \vdots \\ d_{m-1} \end{pmatrix}, \qquad a=c=-\lambda,\;\; b=2(1+\lambda) $$

Solução do sistema: Algoritmo de Thomas — eliminação de Gauss adaptada à estrutura tridiagonal, com custo $O(m)$ por passo (muito mais eficiente que SOR neste caso). É um método direto, não iterativo.

Estabilidade: incondicionalmente estável para qualquer $\lambda > 0$ (ao contrário do esquema explícito, que exigiria $\lambda \le 0.5$).

Precisão: $O(k^2 + h^2)$ — segunda ordem em ambos os passos.

📘 Exemplo 1
📗 Exemplo 2

📘 Crank-Nicolson Incondicionalmente Estável

$\alpha = 1$, $l = 1$, $T = 0.5$, $f(x) = \sin(\pi x)$

📗 Teste de Estabilidade com λ = 5

Crank-Nicolson com $\lambda = 5$ — impossível no método explícito!

🎯 Equação de Poisson 2D

Fundamentação Teórica

$$ -\Delta u = f(x,y), \quad (x,y) \in \Omega = (0,1)^2 $$

com condições de contorno $u = g$ em $\partial\Omega$.

🎯 Esquema de 5 Pontos:

Para $h_x = h_y = h$, a discretização resulta em:

$$ 4u_{i,j} - u_{i+1,j} - u_{i-1,j} - u_{i,j+1} - u_{i,j-1} = h^2 f_{i,j} $$

🔬 Formulação Numérica (Esquema de 5 Pontos + SOR)

Malha: $h = 1/(N+1)$ em ambas as direções. Denotamos $u_{i,j} \approx u(x_i, y_j)$ com $x_i = i h$, $y_j = j h$.

Aplicando diferenças centrais a $-\Delta u = f$ nos pontos interiores:

$$ \frac{-u_{i+1,j} + 2u_{i,j} - u_{i-1,j}}{h^2} + \frac{-u_{i,j+1} + 2u_{i,j} - u_{i,j-1}}{h^2} = f_{i,j} $$

Multiplicando por $h^2$:

$$ 4u_{i,j} - u_{i+1,j} - u_{i-1,j} - u_{i,j+1} - u_{i,j-1} = h^2 f_{i,j} $$

Sistema linear: $A\vec{u} = \vec{b}$ com $N^2$ incógnitas e matriz $A$ esparsa (5-diagonal). Forma matricial:

$$ \begin{pmatrix} 4 & -1 & & \\ -1 & 4 & -1 & \\ & \ddots & \ddots & \ddots \\ & & -1 & 4 \end{pmatrix} \begin{pmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N,N} \end{pmatrix} = \begin{pmatrix} h^2 f_{1,1} + \text{CC} \\ h^2 f_{2,1} + \text{CC} \\ \vdots \\ h^2 f_{N,N} + \text{CC} \end{pmatrix} $$

onde "CC" representa contribuições das condições de contorno (nulas neste exemplo, pois $u=0$ na fronteira).

Método iterativo SOR (Successive Over-Relaxation):

$$ \boxed{\; u_{i,j}^{(t+1)} = (1-\omega)\, u_{i,j}^{(t)} + \frac{\omega}{4}\left[ u_{i+1,j}^{(t)} + u_{i-1,j}^{(t+1)} + u_{i,j+1}^{(t)} + u_{i,j-1}^{(t+1)} + h^2 f_{i,j} \right] \;} $$

Observe o uso misto de $(t)$ e $(t+1)$: os vizinhos já atualizados ($i-1,j$ e $i,j-1$) entram com o valor novo (Gauss-Seidel), e o fator $\omega$ acelera a convergência (SOR).

Parâmetro de relaxação ótimo: para domínio retangular com malha uniforme $h$,

$$ \omega_{\text{opt}} = \frac{2}{1 + \sin(\pi h)} $$

Critério de parada: $\max_{i,j} |u_{i,j}^{(t+1)} - u_{i,j}^{(t)}| < \varepsilon$.

Convergência: garantida para $0 < \omega < 2$ pois $A$ é simétrica definida positiva.

🎯 Poisson com Termo Fonte Senoidal SOR

$-\Delta u = 2\pi^2 \sin(\pi x)\sin(\pi y)$, $u = 0$ em $\partial\Omega$

Solução exata: $u(x,y) = \sin(\pi x)\sin(\pi y)$

⭕ Equação de Laplace 2D

Fundamentação Teórica

$$ \Delta u = 0, \quad (x,y) \in \Omega = (0,1)^2 $$

com condições de contorno $u = g$ em $\partial\Omega$.

⭕ Princípio do Máximo:

O máximo e o mínimo de uma função harmônica são atingidos na fronteira $\partial\Omega$.

🔬 Formulação Numérica (Esquema de 5 Pontos + SOR)

Caso particular de Poisson com $f = 0$. A discretização é:

$$ 4u_{i,j} - u_{i+1,j} - u_{i-1,j} - u_{i,j+1} - u_{i,j-1} = 0 $$

Sistema linear: $A\vec{u} = \vec{b}$ com $N^2$ incógnitas; agora $\vec{b}$ contém as contribuições das condições de contorno não-homogêneas:

$$ \underbrace{\begin{pmatrix} 4 & -1 & & \\ -1 & 4 & -1 & \\ & \ddots & \ddots & \ddots \\ & & -1 & 4 \end{pmatrix}}_{A} \begin{pmatrix} u_{1,1} \\ u_{2,1} \\ \vdots \\ u_{N,N} \end{pmatrix} = \begin{pmatrix} \text{CC}_1 \\ \text{CC}_2 \\ \vdots \\ \text{CC}_{N^2} \end{pmatrix} $$

onde $\text{CC}_{i,j}$ é não nula apenas para pontos adjacentes à fronteira onde $u = \sin(\pi x) \neq 0$ (ou seja, a base $y=0$).

SOR: mesmas iterações do caso Poisson, sem o termo $h^2 f_{i,j}$:

$$ u_{i,j}^{(t+1)} = (1-\omega)\, u_{i,j}^{(t)} + \frac{\omega}{4}\left[ u_{i+1,j}^{(t)} + u_{i-1,j}^{(t+1)} + u_{i,j+1}^{(t)} + u_{i,j-1}^{(t+1)} \right] $$

ω ótimo: $\omega_{\text{opt}} = \dfrac{2}{1 + \sin(\pi h)}$ (mesma fórmula).

Por que SOR aqui? O sistema é esparso mas grande ($N^2$ equações). Métodos diretos (Gauss) seriam $O(N^6)$ — inviável. SOR converge em $O(N^2 \log N)$ iterações com $\omega$ ótimo, tornando-se competitivo.

⭕ Laplace com Condição de Contorno Não-Homogênea SOR

$\Delta u = 0$, $u(x,0) = \sin(\pi x)$, $u = 0$ nos demais lados

Solução exata: $u(x,y) = \sin(\pi x)\dfrac{\sinh(\pi(1-y))}{\sinh(\pi)}$