ISSN 2317-3297
Solução Numérica de uma Equação Diferencial Parcial
Parabólica pelo Método das Diferenças Finitas
Santos D. Miranda Borjas,
Luiz P. Queiroz,
Depto de Ciências Exatas e Naturais, DCEN, UFERSA,
59625-900, Mossoró, RN
E-mail: [email protected], luiz [email protected],
Palavras-chave: Solução numérica, Equação diferencial parcial parabólica, Diferenças finitas.
Resumo: Neste trabalho apresenta-se uma solução numérica para uma equação diferencial parcial parabólica, mediante o método de diferenças finitas. O método das diferenças finitas é
aplicado a um problema acadêmico de transferência de calor por condução, a qual é modelada
pela equação de difusão do calor unidimensional. Tomamos uma malha retangular com nós,
passo espacial e temporal, para representar o valor da função solução em um ponto. A função
solução é aproximada usando três métodos: Euler Explı́cito, Euler Implı́cito e o método de Crank
Nicholson, usamos a medida do erro para avaliar o desempenho desses métodos.
1
Introdução
Muitos problemas na engenharia são modelados por equações diferenciais parciais (EDP). As
soluções analı́ticas das EDPs são conhecidas apenas para casos muito particulares, em geral
não é fácil encontrar essa solução. Este problema é solucionado usando métodos numéricos
(supondo que exista a solução). Existem diferentes métodos numéricos para obter uma solução
aproximada de uma EDP entre estes temos: método dos volumes finitos, método das diferenças
finitas (MDF), método dos elementos finitos, método dos elementos de contorno, método de
Wavelet, método dos momentos, método das camadas finitas e entre outros [7]. À aplicação do
método das diferenças finitas na solução de diferentes problemas em estudo vem sendo tema de
pesquisa nos últimos anos[1], [2], [3], [4] e [5].
Neste trabalho é apresentada uma solução numérica para uma EDP parabólica usando MDF,
aplicada a um problema acadêmico de transferência de calor por condução, a qual é modelada
pela equação de difusão do calor unidimensional:
ut = cuxx 0 ≤ x ≤ L e
u(x, 0) = f (x)
u(0, t) = u(L, t) = 0
t>0
(1)
0≤x≤L
(2)
0≤t
(3)
onde u=u(x, t) representa a temperatura da barra no ponto x no instante t, c representa a condutividade térmica do material da barra e f (x) é uma função dada que descreve a temperatura
nos vários pontos da barra no instante t = 0. Estamos Supondo que as extremidades na barra
estejam isoladas termicamente. Isso quer dizer que os fluxos de calor através de x = 0 e x = L
são nulos. A equação (2) é chamada de condição inicial e a equação (3) é chamada de condição de
fronteira [6]. Para um problema acadêmico supor f (x)=sen(πx), L = 1, c = 1 e a variável tem2
poral t varia entre 0 e 0.1 assim, a solução analı́tica do problema acima é u(x, t)=−sen(πx)eπ t .
O objetivo neste trabalho é comparar o desempenho de três métodos para aproximar a função
u usando MDF.
395
ISSN 2317-3297
2
Metodologia e procedimentos
Tomamos uma malha retangular com nós (xi , tj ), passo espacial h=1/nx e temporal k = 0.1/nt,
onde xi =0 + ih e tj =0 + jk para i=0, .., nx e j=0, .., nt. A notação ui,j = u(xi , tj ) representa
o valor da função u no ponto (xi , tj ). A função u é aproximada usando método Explı́cito de
Euler, Método Implı́cito de Euler e o Método de Crank Nicholson.
No método de Euler Explı́cito as funções ut e uxx são aproximadas por uma diferença progressiva
e por uma diferença central respectivamente [7], substituindo estas na equação (1) resulta:
ui,j+1 = δ(ui+1,j + ui−1,j ) + (1 + 2δ)ui,j
para i = 0, .., nx; j = 0, ., nt
(4)
onde δ =ck/h2 . A equação (4) nos mostra a solução aproximada para um instante tj+1 a partir
de um tj conhecido. Fixando um instante final T e um certo número de subintervalos temporais
nt então a equação (4) permite aproximar a solução em cada instante até chegar no ponto T .
A condição inicial proporciona os valores ui,0 = f (xi ) para i = 0, .., nx. Este método converge
para δ ≤ 1/2 [7].
A condição de convergência δ ≤ 1/2 obriga a tomar passos temporais pequenos e diminuir
o passo espacial da malha, obtendo melhor aproximação da solução, porém um maior esforço
computacional. Essa limitação pode ser evitada considerando o método Implı́cito de Euler. Este
método usa a diferença regreseiva para aproximar a função ut [7] e a função uxx é aproximada
por uma diferença central, substituindo estas na equação (1) resulta:
ui,j−1 = (1 + 2δ)ui,j − δ(ui−1,j + ui+1,j ) para i = 1, .., nx − 1; j = 0, ., nt
(5)
para cada j , a equação (5) representa um sistema tridiagonal que proporciona a solução no
instante tj a partir da solução no instante tj−1 .
O método de Crank Nicholson é uma combinação do método de Euler Explı́cito para frente em
t, e o método de Euler Implı́cito para trás em t + 1. Este método conciste em discretizar em
dois instantes a equação (1), em seguida faz-se uma média. No primeiro instante tj , aproxima
a função ut usando diferença progressiva, e no segundo instante tj+1 , a função ut é aproximada
por diferença regressiva, dessa forma, resulta [7]:
2(1 + δ)ui,j+1 − δ(ui+1,j+1 + ui−1,j+1 ) = 2(1 − δ)ui,j + δ(ui+1,j + ui−1,j )
(6)
para cada j fixo a equação (6) descreve um sistema tridiagonal, cujos termos independentes se
obtém a partir da solução no instante tj . O sistema nos proporciona a solução no instante tj+1 .
3
Resultados
Nesta secção mostra-se através de tabelas e gráficos a comparação dos três métodos usados para
aproximar a função u . Essa comparação é realizada usando a medida do erro entre a solução
analı́tica u e a solução aproximada us para cada método. No primeiro caso consideramos h=1/20
e k=0.1/100, para esses valores claramente vemos que se satisfaz a condição de convergência
δ ≤ 1/2 e os resultados numéricos em relação ao erro são mostrados na Tabela 1.
Métodos
Erro
Tempo (s)
Euler Explı́cito
0,013
0,4992
Euler Implı́cito
0,013
0,1092
Crank Nicholson
0,013
0,078
Tabela 1: Resultados numéricos do desempenho dos métodos para h=1/20 e k=0.1/100,
Analisando-se os valores da tabela 1, todos os métodos tiveram erro igual a 0,013. Verifica-se
que o tempo de procesamento para computar us é maior para o método Euler Explı́cito.
No segundo caso, com o objetivo de obter um menor valor do erro reduzimos o intervalo espacial
396
ISSN 2317-3297
h = 1/27 e o intervalo temporal k = 0.1/125, seus resultados numéricos são mostrados na Tabela
2.
Métodos
Erro
Tempo (s)
Euler Explı́cito
0,0008694
0,9828
Euler Implı́cito
0,0008699
0,1560
Crank Nicholson
0,0008692
0,0468
Tabela 2: Resultados numéricos do desempenho dos métodos para h = 1/27 e k = 0.1/125
Analisando-se os valores da tabela 2, todos os modelos tiveram um erro aproximado de 0,0008
no entanto a condição de convergência para o método de Euler Explı́cito não é satisfeita pois
0.5 ≤ δ. Verifica-se que o tempo de procesamento para obter a função us é menor para o método
Crank Nicholson, motivo pelo qual opto-se por este método para aproximar a função u, como
mostra a Figura 1.
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
0.01
x
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
t
Figura 1: Método de Crank Nicholson para h = 1/27 e k = 0.1/125
4
Análise e conclusão dos resultados
Os resultados obtidos na Tabela 2 mostra que o critério δ ≤ 0.5 para a convergência do método
de Euler não é uma condição necessária, mas se uma condição suficiente. Para nosso problema
em estudo, verica-se que ao refinar a malha numérica o método Euler Explı́cito apresenta um erro
menor em torno de 0.001, mas isto leva a um maior esforço computacional de aproximadamente
50 %. Para este caso em particular o método de Crank Nicholson apresenta menor esforço
computacional em comparação com os outros métodos.
Referências
[1] D. Alexei , J. Sajeev, Finite difference discretization of semiconductor drift-diffusion equations for nanowire solar cells, Computer Physics Communications, 183 (2012) 2128-2135.
[2] J. Gong, J. Nordström, Interface procedures for finite difference approximations of the
advection diffusion equation, Computational and Applied Mathematics, 236 (2011) 602-620.
[3] A. G. Herwig, Finite difference approximations of first derivatives for two dimensional grid
singularities Journal Computational Physics, 217 (2006) 642-657.
[4] L. Jiang, J. D. Moulton, D. Svyatskiy, Analysis of stochastic mimetic finite difference
methods and their applications in single phase stochastic flows , Computer Methods in
Applied Mechanics and Engineering, 217 (2012) 58-76.
397
ISSN 2317-3297
[5] M. Parvazinia, An elemental scale adjustment method for the finite element and finite
difference solutions of diffusion reaction and convection diffusion equations, Finite Elements
in Analysis and Design, 58 ( 2012) 1-9.
[6] V. M. Iório, ”EDP: Um curso de Graduação”, IMPA, Rio de Janeiro, 2005
[7] W. Yang, W. Cao, T. Chung, J. Morris, ”Applied Numerical Methods Using Matlab”,
Wiley Interscience, New Jersey, 2005.
398