Método das diferenças finitas

O método das diferenças finitas (MDF) é um método de resolução de equações diferenciais que se baseia na aproximação de derivadas por diferenças finitas. A fórmula de aproximação obtém-se da série de Taylor da função derivada.^[1] Hoje, os MDFs são a abordagem dominante das soluções numéricas de equações diferenciais parciais.^[2]

O operador de diferenças finitas para derivada pode ser obtido a partir da série de Taylor para as seguintes funções:

f(x+h)=f(x)+f'(x)h+{\frac {f''(x)h^{2}}{2}}+{\frac {f'''(x)h^{3}}{6}}+o(h^{4})\,

f(x-h)=f(x)-f'(x)h+{\frac {f''(x)h^{2}}{2}}-{\frac {f'''(x)h^{3}}{6}}+o(h^{4})

Portanto, a derivada primeira pode ser escrita de três formas distintas como uma diferença-quociente mais um termo de erro, obtido ao desprezar-se termos de ordem superior :

f'(x)={\frac {f(x+h)-f(x)}{h}}+o(h)

, que é conhecida como fórmula das diferenças progressivas, ou

f'(x)={\frac {f(x)-f(x-h)}{h}}+o(h)

, que é conhecida como fórmula das diferenças regressivas, ou ainda

f'(x)={\frac {f(x+h)-f(x-h)}{2h}}+o(h^{2})

, que é conhecida como fórmula das diferenças centradas.

Além disso, é possível obter derivadas de ordem superior. A derivada de segunda ordem é obtida a partir de

f(x+h)+f(x-h)=2f(x)+f''(x)h^{2}+o(h^{4})

e é dada por $f''(x)={\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}+o(h^{2})$

Método das diferenças finitas para problemas lineares

A partir das aproximações por diferença-quociente para derivadas de qualquer ordem, é possível transformar equações diferenciais em problemas lineares. Para isso, é necessário ignorar o termo de erro e tornar $h$ um número muito pequeno, mas grande o suficiente para que não cause instabilidades nas aproximações das derivadas.

Resolução de problemas de contorno

Para uma equação diferencial do tipo $y''(x)=f(x)y'(x)+g(x)y(x)+z(x)$ , onde $x$ varia de $a$ até $b$ , $y(a)=\alpha$ e $y(b)=\beta$ .

A equação é aproximada pelo método das diferenças finitas, com um erro de truncamento igual a $o(h^{2})$ , substituindo-se as derivadas pelas suas representações numéricas, que são dadas por:

$y'(x)={\frac {y(x+h)-y(x-h)}{2h}}$

$y''(x)={\frac {y(x+h)-2y(x)+y(x-h)}{h^{2}}}$

Como é possível perceber, necessita-se definir um valor para $h$ . Este valor pode ser definido pela divisão do intervalo em que se está interessado para a resolução do problema em $N$ intervalos menores. Assim, o valor de $h$ é dado por: $h={\frac {b-a}{N}}$ .

As extremidades destes subintervalos são dadas por $x(i)=a+(i-1)h$ , para $i=1,2,...,N$ .

Para a resolução do problema, o mesmo é escrito na forma $y''(x)-f(x)y'(x)-g(x)y(x)=z(x)$ , que após a substituição das derivadas, torna-se:

${\frac {y(x(i)+h)-2y(x(i))+y(x(i)-h)}{h^{2}}}-f(x(i))\left[{\frac {y(x(i)+h)-y(x(i)-h)}{2h}}\right]-g(x(i))y(x(i))=z(x(i))$ , para $i=2,3,...,N$ .

Como $x(i+1)=x(i)+h$ e $x(i-1)=x(i)-h$ , aquela equação pode ser reescrita como

${\frac {y(x(i+1))-2y(x(i))+y(x(i-1))}{h^{2}}}-f(x(i))\left[{\frac {y(x(i+1))-y(x(i-1))}{2h}}\right]-g(x(i))y(x(i))=z(x(i))$ .

Isolando os termos $y(x(i+1))$ , $y(x(i))$ e $y(x(i-1))$ na fórmula acima, obtêm-se

$y(x(i-1))\left[{\frac {1}{h^{2}}}+{\frac {f(x(i))}{2h}}\right]+y(x(i))\left[{\frac {-2}{h^{2}}}-g(x(i))\right]+y(x(i+1))\left[{\frac {1}{h^{2}}}-{\frac {f(x(i))}{2h}}\right]=z(x(i))$

A partir desta equação é possível resolver o sistema linear a partir de uma matriz $A$ de coeficientes que multiplica os valores de $y$ , sendo que a solução deste sistema é dada por $B$ . Esse sistema linear é representado a seguir.

$A={\begin{bmatrix}1&0&0&\cdots &\cdots &\cdots &0\\{\frac {1}{h^{2}}}+{\frac {f(x(2))}{2h}}&{\frac {-2}{h}}-g(x(2))&{\frac {1}{h^{2}}}-{\frac {f(x(2))}{2h}}&0&\cdots &\cdots &0\\0&\ddots &\ddots &\ddots &\vdots &\vdots &0\\0&\ddots &\ddots &\ddots &{\frac {1}{h^{2}}}+{\frac {f(x(N-1))}{2h}}&{\frac {-2}{h}}-g(x(N))&{\frac {1}{h^{2}}}-{\frac {f(x(N+1))}{2h}}\\0&0&\cdots &\cdots &\cdots &\cdots &1\end{bmatrix}}$

$B={\begin{bmatrix}\alpha \\z(x(2))\\z(x(3))\\\vdots \\\beta \end{bmatrix}}$

Onde a aproximação para $y$ é dada pelos pontos que são solução do sistema $Ay=B$ .

Resolução de problemas de valor inicial e o método de Euler

A partir do método das diferenças finitas também é possível obter o método de Euler, que é usado para obter soluções de problemas de valor inicial bem-posto. Leonhard Euler (1707 - 1783) foi o primeiro matemático de sua época a apresentar o uso do método de diferenças finitas para encontrar aproximações de soluções de equações diferenciais. Entretanto, o método de Euler não é usado na prática, pois possui pouca precisão. Alternativamente a este, são utilizados com maior frequência o método de Euler modificado ou o método de Runge-Kutta para solução de problemas de valor inicial.

Para um dado problema de valor inicial bem posto

$y'(t)=f(t,y)$ , $a$ ≤ $t$ ≤ $b$ , $y(a)=\alpha$ .

Divide-se o intervalo $[a,b]$ em $N$ subintervalos e define-se que $t(i)=a+(i-1)h$ , para $i=1,2,...,N$ . Onde $h$ é o espaçamento da malha.

A partir disto, temos que $h={\frac {(b-a)}{N}}=t(i+1)-t(i)$ , para $i=1,2,...,N$ .

Aproximando a equação diferencial pelo método das diferenças finitas, desprezando-se o termo de erro, temos

${\frac {y(t(i+1))-y(t(i))}{h}}=f(t(i),y(i))$

$y(t(i+1))=y(t(i))+hf(t(i),y(i))$ , que é usada para $i=1,2,3,...,N.$

A equação acima é conhecida como equação de diferença associada ao método de Euler.

O sistema linear é inicializado com $y(1)=\alpha$ e é de fácil solução.

Método das diferenças finitas para problemas não lineares

O método das diferenças finitas é análogo ao utilizado para problemas lineares. Entretanto, é utilizado um processo iterativo para a obtenção da solução do problema, que não é linear.

Resolução de problemas de contorno não-lineares

Para um problema de contorno não-linear geral, dado por $y''=f(x,y,y')$ com $x$ variando de $a$ até $b$ , e sendo as condições de contorno $y(a)=\alpha$ e $y(b)=\beta$ ,

há garantia de solução única se as seguintes condições forem satisfeitas.

$f$ e suas derivadas parciais em relação a $y$ e $y'$ são contínuas em $D=$ { $(x,y,y')|a$ ≤ $x$ ≤ $b,$ - $\infty$ < $y$ < $\infty$ , - $\infty$ < $y'$ < $\infty$ };
$f(x,y,y')$ > $\delta$ , para um certo $\delta$ > $0$ ;
Existem constantes $P$ e $Q$ tais que: $P$ é o valor máximo que o módulo da derivada parcial de $f$ em relação a $y$ atinge em $D$ e $Q$ é o valor máximo que o módulo da derivada parcial de $f$ em relação a $y'$ atinge em $D$ .

Como no caso anterior, a aproximação para a equação é obtida quando os termos de erro são desprezados.

Assim, $y''=f(x,y,y')$ torna-se ${\frac {y(x+h)-2y(x)+y(x-h)}{h^{2}}}=f(x,y,{\frac {y(x+h)-y(x-h)}{2h}})$ , que com a mesma divisão em intervalos anterior, é dada por ${\frac {y(x(i+1))-2y(x(i))+y(x(i-1))}{h^{2}}}=f(x(i),y(i),{\frac {y(x(i+1))-y(x(i-1))}{2h}})$ , para $i=2,3,...,N$ .

As condições de contorno são $y(1)=\alpha$ e $y(N+1)=\beta$ .

A partir da equação $-{\frac {y(x(i+1))-2y(x(i))+y(x(i-1))}{h^{2}}}+f(x(i),y(i),{\frac {y(x(i+1))-y(x(i-1))}{2h}})=0$ , para $i=2,3,...,N$ e das condições de contorno, obtemos um sistema não linear que pode ser resolvido via Método de Newton para sistemas não-lineares. Sendo que o sistema terá solução única se $h$ < $2/Q$ . Se a aproximação inicial utilizada no método de Newton for suficientemente próxima da solução e se a matriz Jacobiana do sistema for não-singular, o sistema converge para a solução exata.

Referências

↑ Richard L. Burden; J. Douglas Faires. Análise Numérica, Editora CENGAGE Learning, 8° edição. [S.l.: s.n.] A referência emprega parâmetros obsoletos |coautor= (ajuda)
↑ Christian Grossmann; Hans-G. Roos; Martin Stynes (2007). Numerical Treatment of Partial Differential Equations. [S.l.]: Springer Science & Business Media. p. 23. ISBN 978-3-540-71584-9