www.metodosnumericos.com.br --Prof. Doherty Andrade-- Visualização Interativa Completa: Superfícies 3D Animadas com Fatias
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:
| Discriminante | Tipo | Exemplo | Caracterí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 |
com duas condições iniciais (posição e velocidade):
$$ u(x,0) = f(x), \quad \frac{\partial u}{\partial t}(x,0) = g(x) $$Quando $\lambda = \frac{\alpha k}{h} = 1$, o esquema produz a solução exata da equação da onda. Energia perfeitamente conservada.
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.
$f(x) = \sin(\pi x)$, $g(x) = 0$, $u(x,t) = \sin(\pi x)\cos(\pi t)$
$f(x) = \sin(\pi x)$, $g(x) = 0$, $m = 20$, $N = 20$
$f(x) = \sin(2\pi x)$, $g(x) = -2\pi\sin(2\pi x)$, $u(x,t) = \sin(2\pi(x-t))$
com condição inicial $u(x,0) = f(x)$ e condições de contorno $u(0,t) = u(l,t) = 0$.
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.
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.
$\alpha = 1$, $l = 1$, $T = 0.5$, $f(x) = \sin(\pi x)$
Crank-Nicolson com $\lambda = 5$ — impossível no método explícito!
com condições de contorno $u = g$ em $\partial\Omega$.
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} $$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.
$-\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)$
com condições de contorno $u = g$ em $\partial\Omega$.
O máximo e o mínimo de uma função harmônica são atingidos na fronteira $\partial\Omega$.
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.
$\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)}$
A solução decai exponencialmente à medida que $y$ aumenta, ilustrando o princípio do máximo. A solução exibida é a numérica obtida por SOR.