Método de diferenças finitas

Considere o problema de valor de contorno, não linear e de segunda ordem \begin{align} \label{pvc} y'' &= f(x,y,y'), \quad a \lt x \lt b,\\[2mm] y(a) &= \alpha, \nonumber\\ y(b) &= \beta. \nonumber \end{align}

Queremos aproximações para os valores da função $y$ sobre uma malha regular para o intervalo $[a,b].$ Mais claramente, dado $n\in\mathbb{N},$ seja $h=(b-a)/n$ e defina $x_j = a + jh.$ Queremos encontrar $y_j,$ tal que \begin{equation}\label{yj} y_j \approx y(x_j), \quad j=0,1,\ldots,n. \end{equation} Conhecidos os valores $y_j,$ se realmente quisermos uma função contínua que aproxime a função $y,$ podemos resolver um problema de interpolação por partes. Se quisermos uma função diferenciável, podemos realizar a interpolação por splines cúbicas.

Como temos $(n+1)$ incógnitas a determinar, os valor $y_j,$ precisamos montar um sistema com $(n+1)$ equações. As condições de contorno de Dirichlet fornecem duas equações que imediatamente determinam os valores de $y$ nos extremos do intervalo. Com efeito, \begin{align*} y_0 & \equiv \alpha = y(x_0),\\ y_n & \equiv \beta = y(x_n). \end{align*} Nos resta conseguir outras $(n-1)$ equações envolvendo $y_j,$ $j=1,\ldots,(n-1).$ Como a equação diferencial \eqref{pvc} vale para todo $x\in(a,b),$ em particular vale nos pontos $x_j$ da malha, ou seja, \begin{equation}\label{disc1} y''(x_j) = f(x_j,y(x_j),y'(x_j)), \quad j=1,\ldots,(n-1). \end{equation}

Com as duas equações obtidas diretamente das condições de contorno e as $(n-1)$ equações acima, temos as $(n+1)$ equações que precisamos. Entretanto, nestas equações as incógnitas são $y(x_j)$, $y'(x_j)$ e $y''(x_j)$, ao passo que queremos ficar apenas com incógnitas sobre os valores da função $y.$ Para superar esta dificuldade, a estratégia é impor que as equações acima sejam satisfeitas não com $y(x_j),$ $y'(x_j)$ e $y''(x_j),$ mas sim com aproximações para estas quantidades, todas escritas apenas em termos de $\{y_0,y_1,\ldots,y_n\}.$

No lugar de $y(x_j)$ usamos diretamente a definição de $y_j$, dada em \eqref{yj}. No lugar de $y'(x_j)$ e $y''(x_j)$ usamos as aproximações de segunda ordem para derivadas primeira e segunda, já estudada na aula sobre diferenças finitas, a saber \begin{align} y'(x_j) &\approx {y_{j+1} - y_{j-1}\over 2h}, \\[6pt] y''(x_j) &\approx { y_{j-1}-2y_j+y_{j+1}\over h^2}. \end{align}

Desta forma, ficamos com o sistema não linear em $\{y_0,y_1,\ldots,y_n\},$ dado por \begin{align*} y_0 & = \alpha, \\ y_0 - 2y_1 + y_2 & = h^2 f\left(x_1,y_1, {(y_2-y_0)\over 2h} \right), \\ y_1 - 2y_2 + y_3 & = h^2 f\left(x_2,y_2, {(y_3-y_1)\over 2h} \right), \\ & \vdots \\ y_{n_2} - 2y_{n-1} + y_n & = h^2 f\left(x_{n-1},y_{n-1}, {(y_n-y_{n-2})\over 2h}\right), \\ y_n & = \beta. \end{align*}

Esta estratégia de construir um sistema de equações pelo emprego de aproximações de diferenças finitas para as derivadas é conhecida como o método de diferenças finitas. Este é o método numérico mais importante para a aproximação da solução de um PVC.

Exemplo: pêndulo com haste rígida

Considere como exemplo o PVC que descreve o movimento de um pêndulo com haste rígida de comprimento $L$, sujeito apenas à força da gravidade, \begin{equation}\label{exemplo} y'' = -{g\over L} \sin(y), \qquad 0 \lt t \lt 1, \qquad y(0) = -{\pi\over 3}, \quad y(1)={\pi\over 3}, \end{equation} onde $g$ é a aceleração da gravidade e $y(t)$ representa o ângulo que a haste do pêndulo faz com a vertical, no tempo $t.$

Discretizando o intervalo $[0,1]$ nos 5 pontos regularmente espaçados $\{0,{1\over 4},{1\over 2},{3\over 4}, 1\}$ $(h = 1/4),$ ficamos com o sistema \begin{align*} y_0 & = -{\pi\over 3}, \\ y_0 - 2y_1 + y_2 & = h^2 \left[ -{g\over L} \sin(y_1) \right],\\ y_1 - 2y_2 + y_3 & = h^2 \left[ -{g\over L} \sin(y_2) \right],\\ y_2 - 2y_3 + y_4 & = h^2 \left[ -{g\over L} \sin(y_3) \right],\\ y_4 & = {\pi\over 3}. \end{align*} A solução deste sistema pelo método de Newton para sistemas não lineares foi exibida no vídeo. Na figura abaixo, está a solução do mesmo PVC, mas em uma malha com 21 pontos. Esta solução foi obtida tomando como aproximação inicial para o método de Newton a reta que une os pontos dados na condição de contorno. Apenas 3 iterações do método de Newton foram necessárias para satisfazer o critério de parada.

Exemplo: Sistema massa–mola com vibrações forçadas

Considere o PVC abaixo, que surge do modelamento de um sistema massa–mola com amortecimento, sujeito a uma vibração forçada periódica \begin{equation}\label{exemplo2} my''+cy'+ky= A\cos(\omega x), \quad 0\lt x\lt 2, \qquad y(0) = y(2) =0, \end{equation} onde $m$ é a massa do objeto preso à mola, $c$ é coeficiente de amortecimento, $k$ é a constante elástica da mola, $A$ é amplitude de oscilação forçada e $\omega_0$ é frequência de oscilação.

Vamos discretizar o intervalo $[0,2]$ em $(n+1)$ pontos igualmente espaçados por $h = (2-0)/n.$ Ficamos com o seguinte sistema \begin{align*} y_0 & = 0, \\ y_0 - 2y_1 + y_2 & = h^2 \left[ -{c \over m} \left({y_2 - y_0 \over 2h} \right) - {k\over m} y_1 + {A\over m}\cos(\omega x_1 )\right],\\ y_1 - 2y_2 + y_3 & = h^2 \left[ -{c \over m} \left({y_3 - y_1 \over 2h} \right) - {k\over m} y_2 + {A\over m}\cos(\omega x_2 )\right],\\ & \;\; \vdots \\ y_{n-2} - 2y_{n-1} + y_n & = h^2 \left[ -{c \over m} \left({y_n- y_{n-2} \over 2h} \right) - {k\over m} y_{n-1} + {A\over m}\cos(\omega x_{n-1} )\right],\\ y_n & = 0, \end{align*} onde $x_j = 0 + h\cdot j,$ para $j=0,1,2\ldots,n.$

O sistema acima, apesar de parecer complicado, é na verdade linear nas incógnitas $\{y_0,y_1,\ldots,y_n\}.$

Na próxima aula vamos nos concentrar no caso particular de um PVC linear, como o deste último exemplo.

1. Considere o exemplo \eqref{exemplo}, da seção aprofunde e seja $y=(y_0,\ldots,y_4)^T,$ tomando $L=1$m.

• Escreva o sistema não linear do exemplo, na forma $F(y) = 0.$
• Em sendo um método iterativo, o método de Newton precisa de uma aproximação inicial para dar início. Como já sabemos que $y_0 = -{\pi\over 3} = - y_4$, defina os demais valores de $y_j$ de maneira que variem linearmente entre estes extremos.
• Aplique o método de Newton para sistemas não lineares e estime a solução daquele PVC.

2. Consisdere o exemplo \eqref{exemplo2}, ao final da seção aprofunde e seja $y=(y_0,\ldots,y_{10})^T$, $h = 2/10$ s. Tome $m = 1$ kg, $c = 0.5$ kg/s, $k=1$ kg/s², $A=2$ kg/s² e $\omega = 5$ s⁻¹.

• Escreva como fica a equação diferencial discretizada, em um ponto $x_j$ interno ao intervalo.
• Escreva o sistema linear do exemplo, na forma $My = b.$
• Usando o Octave, resolva o sistema linear e exiba o gráfico para a aproximação obtida para $y$ no intervalo $[0,2].$