UNIVERSIDADE FEDERAL DO AMAPÁ

PRÓ-REITORIA DE ENSINO E GRADUAÇÃO

CURSO DE LICENCIATURA EM MATEMÁTICA

MÉTODO DE DIFERENÇAS FINITAS: APLICAÇÃO NA EQUAÇÃO

DO CALOR

MACAPÁ-AP
2013
 JOELMA DO SOCORRO MIRANDA DO NASCIMENTO
TAMIRES DA SILVA

MÉTODO DE DIFERENÇAS FINITAS: APLICAÇÃO NA EQUAÇÃO

DO CALOR

Trabalho de conclusão de curso apre-

sentado ao colegiado de Matemática da
Universidade Federal do Amapá, como
parte das exigências para a obtenção do
tı́tulo de Licenciatura em Matemática,
sob orientação do Prof. Dr. Guzmán
Eulálio Isla Chamilco

MACAPÁ-AP
2013
 JOELMA DO SOCORRO MIRANDA DO NASCIMENTO
TAMIRES DA SILVA

MÉTODO DE DIFERENÇAS FINITAS: APLICAÇÃO NA EQUAÇÃO

DO CALOR

Trabalho de Conclusão de Curso apresentado como pré-requisito para obtenção do tı́tulo

de Licenciatura em Matemática da Universidade Federal do Amapá, submetida à aprovação
da banca examinadora composta pelos seguintes membros:

Prof. Dr. Guzmán Eulálio Isla Chamilco

Prof. Dr. José Walter Cárdenas Sotil

Prof. Dr. Gilberlandio Jesus Dias

Avaliado em: 07/05/2013

MACAPÁ-AP
2013
 Agradecimentos

Agradecemos à Deus pela dádiva da vida que nos permite aprender com os fracas-
sos e sentir o prazer das conquistas.
À nossa famı́lia, em especial aos nossos pais, Benedito Ferreira Nascimento e Raimunda
do Socorro Miranda do Nascimento, e Tânia Maria da Silva pelas mensagens e gestos de
incentivos para que completássemos com êxito mais essa etapa de nossas vidas.
Aos colegas da graduação, em especial, à Adriane Bastos, Deniel de Almeida, Elias da
Costa, Fábio Dias e Fagner Carmo pelo companherismo nesses anos de estudo e incentivo
na produção deste trabalho.
Agradecemos em especial ao nosso orientador Prof. Dr. Guzmán Eulálio Isla Cha-
milco pelas contribuições teóricas e de materiais para o desenvolvimento desta monografia.
Deixamos registrados toda nossa gratidão, respeito e admiração.
Aos professores do colegiado deixamos nossos sinceros agradecimentos pelas contri-
buições para nosso progresso desde o inı́cio da graduação.
Enfim, a todos aqueles que cruzaram o nosso caminho nesta trajetória, pois de uma
forma ou de outra, são pessoas que fazem parte de um pedaço de nossas vidas e que jamais
esqueceremos.

A todos muito obrigado!

“As raı́zes do estudo são amargas,
mas seus frutos são doces.”

(Aristóteles)
 Resumo

Neste trabalho de conclusão de curso, abordamos o conceito do Método das Dife-

renças Finitas, juntamente com suas propriedades. Inicia-se o trabalho com uma revisão
histórica, enfatizando alguns matemáticos que foram imprescindı́veis para o desenvolvi-
mento dos métodos de diferenças. No aspecto teórico, trabalhamos com Equações Dife-
renciais Parciais de 2a ordem, com os métodos explı́cito e implı́cito aplicados na Equação
do Calor, e com o método de Crank-Nicholson. Por fim, exemplificamos numericamente
em um problema os conceitos e propriedades estudadas, fazendo uma comparação entre
os métodos.

Palavras-Chave: Método das Diferenças Finitas. Equações Diferenciais Parciais de

2a Ordem. Método Explı́cito. Método Implı́cito. Equação do Calor. Método de Crank-
Nicholson.
 Abstract

In this monograph, we discuss the concept of Finite Difference Method, with their
properties. It begins with a historical review, emphasizing some mathematicians who
were essential to the development of difference methods. In the theoretical aspect, we
deal with Partial Differential Equations of second order, explicit and implicit methods
applied in the Heat Equation, and the Crank-Nicholson method. Finally, we numerically
exemplify those concepts and properties in a problem, making a comparison between the
methods.

Key-words:Finite Difference Method. Partial Differential Equations of Second Order.

Explicit Method. Implicit Method. Heat Equation. Crank-Nicholson Method.

1
Sumário

1 Contexto Histórico 6

2 Definições Básicas de Equações Diferenciais 9

2.1 Equação Diferencial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Equação Diferencial Parcial de Segunda ordem . . . . . . . . . . . . . . . . 10
2.3 Equação do Calor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4 Erro de Truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Métodos Numéricos 14
3.1 Discretização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1 Condições de Contorno . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Método das Diferenças Finitas (MDF) . . . . . . . . . . . . . . . . . . . . 16
3.3 Série de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.4 Consistência, Convergência e Estabilidade . . . . . . . . . . . . . . . . . . 19
3.5 Aplicação do Método Explı́cito na Equação do Calor . . . . . . . . . . . . 21
3.6 Aplicação do Método Implı́cito na Equação do Calor . . . . . . . . . . . . 23
3.7 Método de Crank - Nicholson . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Problema Proposto 31

Considerações Finais 47

Referências Bibliográficas 48

2
Lista de Figuras

3.1 Malha computacional unidimensional . . . . . . . . . . . . . . . . . . . . . 14

3.2 Malha computacional bidimensional . . . . . . . . . . . . . . . . . . . . . . 15
3.3 Diferença Finita Progressiva, Central e Regressiva . . . . . . . . . . . . . . 17
3.4 Método Explı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5 Representação dos nós em subintervalos . . . . . . . . . . . . . . . . . . . . 23
3.6 Método Implı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.7 Método Explı́cito e Implı́cito . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.8 Representação dos nós em subintervalos . . . . . . . . . . . . . . . . . . . . 26
3.9 Representação do Método de Crank-Nicholson . . . . . . . . . . . . . . . . 27
3.10 Representação dos nós em subintervalos . . . . . . . . . . . . . . . . . . . . 29

4.1 Distribuição de temperaturas em uma viga de alumı́nio . . . . . . . . . . . 31

4.2 Variação da temperatura em função do espaço . . . . . . . . . . . . . . . . 35
4.3 Solução Numérica pelo Método Explı́cito para o tempo de 0.1s à 0.4s. . . . 35
4.4 Variação da temperatura em relação do espaço, para diferentes momentos
de tempos, calculado pelo Método Explı́cito. . . . . . . . . . . . . . . . . . 36
4.5 Solução Numérica pelo Método Implı́cito para o tempo de 0,1s à 0,4s. . . . 40
4.6 Variação da temperatura em relação do espaço, para diferentes momentos
de tempos, calculado pelo Método Implı́cito. . . . . . . . . . . . . . . . . . 41
4.7 Solução Numérica pelo Método de Crank-Nicholson para o tempo de 0,1s
à 0,4s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.8 Variação da temperatura em relação do espaço, para diferentes momentos
de tempos, calculado pelo Método de Crank-Nicholson. . . . . . . . . . . . 45

3
Introdução

Muitos fenômenos naturais estudados na Fı́sica, Engenharia e Biologia são traduzidos

para uma linguagem matemática, dando origem a modelos matemáticos. Problemas en-
volvendo movimento de fluidos, fluxo de corrente elétrica, dissipação de calor, propagação
de onda sı́smicas e crescimento populacional podem ser descritas através das equações di-
ferenciais. Por esta razão, saber reconhecer uma equação diferencial é muito importante
para a busca da solução desses problemas.
Uma equação diferencial possui uma famı́lia de soluções e não apenas uma solução
que a satisfaça, entretanto, dada uma condição inicial ou uma condição de contorno e sua
equação diferencial, o que chamamos de problema de valor inicial (PVI) ou problema de
valor de contorno (PVC) e sob certas condições temos a existência e unicidade de solução
garantida. Neste caso, é possı́vel obter aproximações numéricas do PVI ou PVC, sem
obter a solução analı́tica do problema.
A grande necessidade de se obter soluções de questões cada vez mais complexas e
o avanço computacional deram origem aos métodos numéricos. Eles têm como ideia
principal a discretização que reduz o problema contı́nuo, com um número infinito de
variáveis, em um problema discreto, com um número finito de variáveis, podendo ser
resolvido por meio de programas computacionais.
Existem diversos tipos de métodos numéricos utilizados para encontrar a solução de
uma equação diferencial. Um desses métodos é o Método de Diferenças Finitas. A ideia
básica dele é transformar o problema de resolver uma equação diferencial num problema de
resolver um sistema de equações algébricas, usando para isso aproximações das derivadas
que aparecem na equação, por diferenças finitas.
Uma forma de aproximação é aplicar os métodos de diferenças finitas, obtidos através
da discretização do contı́nuo e obter, assim, a solução em apenas pontos discretos do
intervalo considerado.

4
 Problemas que envolvem transferência térmica causada por uma diferença de tem-
peratura entre duas regiões em um mesmo meio podem ser representados por equações
diferenciais parciais. Essas equações são denominadas equações do calor, e determinam a
distribuição de temperatura em um corpo.
Nesse trabalho utilizaremos o Método Explı́cito, Implı́cito e Crank-Nicholson, com
o objetivo analisar a transferência de calor no decorrer do tempo em um problema que
envolvem condução de calor unidimensional em uma barra.

5
Capı́tulo 1

Contexto Histórico

O estudo das equações diferenciais nasceu com problemas da Mecânica, com a ex-
plicação do movimento dos planetas, a oscilação do pêndulo e outros que foram estudados
por Leonardo da Vinci (1452-1519), Galileu Galilei (1564-1642), Johann Kepler(1571-
1630).
A mais notável realização matemática do século XVII foi a invenção do Cálculo por
Isaac Newton(1642-1727) e Gottfried Wilhelm Leibinz(1646-1716). Esse feito mudou a
forma de resolver os problemas da época e influenciou toda uma geração futura.
A obra “Methodus Incrementorum Directa et Inversa”, publicada pelo matemático
Brook Taylor(1685-1731), em 1715 foi de grande importância, pois apresentou a série:

0 00 x2 000 x3 xn
f (x + a) = f (a) + f (a) x + f (a) + f (a) + . . . + f n (a) ...
2! 3! n!

que hoje é conhecida como Série de Taylor. Ele desenvolveu o teorema de expansão
de Taylor, sendo muito utilizada nos métodos de diferenças finitas.
As publicações de Leonhard Euler(1707-1783) “Introduction in analysin infinitorum”(1748),
e “Institutiones calculi differentialis”(1755) e “Institutiones calculi integralis” contêm o
estudo mais completo de cálculo e análise até aquela época. Entre as contribuições de
Euler estão o uso de fatores integrantes na resolução de equações diferenciais, método
sistemático para resolver equações diferenciais lineares com coeficientes constantes.
Na obra “Traité de Dynamique”(1743), Jean Le Rond D’Alembert (1717-1783) escre-
veu sobre o princı́pio que as ações e reações internas de um sistema de corpos rı́gidos em
movimento, estão em equilı́brio. Em 1747, ao estudar o problema das cordas vibrantes
desenvolveu a equação:

6
 ∂ 2u ∂ 2u
=
∂t2 ∂x2

e deu a solução u = f (x + t) + g (x − t), onde f e g são funções arbitrárias, f, g ∈ C 2 .

Pierre-Simon Laplace (1749-1827) foi o primeiro a estudar a equação de Laplace (1780),
em seu trabalho sobre potencial do campo gravitacional.
Na obra “Théorie Analytique de La Chaleur ”, Joseph Fourier(1768-1830) introduziu
a Equação do Calor (1810-1822). Utilizando a teoria dos antecessores, em 1807 Fourier
submeteu seu primeiro trabalho a Academia Francesa, onde formalizou e solucionou o
problema da condução de calor. Seu trabalho não foi aceito e um concurso foi feito
para premiar quem solucionasse o problema. Em 1811, Fourier submeteu novamente seu
trabalho, mas a banca julgadora mais uma vez resolveu não publicá-lo, alegando falta de
rigor. A publicação dos seus trabalhos só ocorreu mais tarde, quando Fourier tornou-se
secretário da Academia. Assim, a teoria de Fourier foi reconhecida, porém não finalizado,
pois novos problemas surgiram do seu trabalho. Dessa forma as equações difereciais
parciais de 2a ordem - Hiperbólica, Elı́ptica e Parabólica - já eram conhecidas no começo
do século XIX.
O século XX foi marcado pelo desenvolvimento da Matemática Aplicada com o sur-
gimento dos métodos numéricos implementados no computador. Eles se tornaram muito
importante na resolução das EDP’s.
Em 1928, Richard Courant (1888-1972), Kurt Otto Friedrichs (1901-1982)e Hans Lewy
(1904-1988) publicaram um artigo sobre as equações diferenciais parciais de fı́sica - ma-
temática. O uso de aproximações de diferenças finitas para provar a existência de soluções
de EDP’s foi a principal motivação que os levou a escrever esse artigo.
Os avanços tecnológicos dos métodos numéricos foram criados por Jonh Von Neu-
mam (1903-1957). Ele foi responsável pelo ENIAC (Electronic Numerical Integrator and
Computer) e idealizou os conceitos de armazenamento de programas.
Peter David Lax (1926), nasceu em Budapeste e migrou para EUA em 1941. Participou
do projeto Manhattan. Friedrichs foi seu orientador em Sistemas de Equações Diferenciais
Parciais Hiperbólicas não Lineares em duas Variáveis Independentes. Em 2005, ganhou
o prêmio Abel Prize pelas contribuições em Matemática Aplicada. Desenvolveu métodos
de diferenças Lax-Friedrichs e Lax-Wendroff, respectivamente, com Friedrichs e Burton
Wendroff (1930),e ainda, o Teorema de Lax, no qual estabelece as condições em que uma
solução numérica é uma boa aproximação para solução de uma equação diferencial.

7
 John Crank (1916-2006) e Phyllis Nicholson (1917-1968) foram responsavéis pelo de-
senvolvimento do método numérico, conhecido por Método de Crank - Nicholson, em me-
ados do século XX. Ele consiste em resolver numericamente a equação do calor e equações
diferenciais parciais similares, sendo de segunda ordem no tempo e no espaço, implı́cito
no tempo e é numericamente estável.

8
Capı́tulo 2

Definições Básicas de Equações

Diferenciais

2.1 Equação Diferencial

As equações diferenciais são todas equações cujas incógnitas são funções e que contém
ao menos uma derivada ou diferencial destas funções. Elas são divididas em: equações
diferenciais ordinárias (EDO) e equações diferenciais parciais (EDP).
A equação diferencial ordinária envolve derivada de uma função que contém somente
uma variável independente e tem a forma:
dy d2 y dn y
 
F x, y, , 2 , . . . , n = 0
dx dx dx
Onde x é a variável independente, y é a variável dependente e n designa a maior das
ordens das derivadas presentes na equação.
Já a equação diferencial parcial envolve mais de uma variável independente. Sua forma
geral é:
∂u ∂ 2 u ∂ku
 
∂u
F (x1 , . . . , xn ) , u, ,..., , ,..., k = 0
∂x1 ∂xn ∂x21 ∂xn
Onde (x1 , . . . , xn ) são variáveis independentes e u (x1 , . . . , xn ) é a variável dependente.
As equações diferenciais podem ser classificadas ainda quanto à linearidade, a ordem e o
grau.
A equação pode ser linear ou não-linear. Uma equação diferencial é linear se ela for
linear na função desconhecida e em todas as sua derivadas, com coeficientes dependendo
apenas das variáveis independentes.

9
 A ordem é determinada com a ordem da derivada de maior ordem presente na equação.
Exemplo:
3 4
d3 y d2 y d2 y
 
x− 3 −y 2 = 2− 2 (3a ordem)
dx dx dx

O grau é o maior dos expoentes a que está elevada a derivada de mais alta ordem que
existe na equação.
Exemplo:
2 4
d3 y d2 y d2 y
 
3x − 3 − y 2 = 1 − 2 (3a ordem e 2a grau)
dx dx dx

2.2 Equação Diferencial Parcial de Segunda ordem

Muitos fenômenos que ocorrem na Ótica, Eletricidade, Ondulatória, Mecânica, Flui-
dos e Biologia, podem ser descritos através de uma equação diferencial parcial. Algumas
leis fı́sicas como: Leis de Newton para o resfriamento dos corpos, equações de Maxwell,
Navier-Stokes e da Mecânica Quântica de Schrodinger são descritas por equações diferen-
ciais que relacionam o espaço e suas derivadas com o tempo.
As equações diferenciais parciais lineares de segunda ordem têm a seguinte forma:

∂ 2u ∂ 2u ∂ 2u ∂u ∂u
A 2
+ B + C 2
+D +E + F (u) = G (x, y)
∂x ∂x∂y ∂y ∂x ∂y

Onde A, B, C, D, E e G são funções que dependem de x e y. Quando os coeficientes

de A . . . G são constantes tem-se uma equação diferencial parcial de segunda ordem com
coeficientes constantes.
Um problema com condição inicial ou de contorno para uma EDP de 2a ordem da
forma:

A (x, y) uxx + B (x, y) uxy + C (x, y) uyy + G (x, y, u, ux , uy ) = 0

é aquele que visa obter uma solução u (x, y) 6= 0 para a equação dada sobre um
conjunto M ⊂ R2 de modo que a função u = u (x, y) deva satisfazer a algumas condições
iniciais ou de contorno dada por funções conhecidas.
As EDP’s de 2a ordem são classificadas em: Equações Elı́pticas, Equações Parabólicas
e Equações Hiperbólicas. Segundo Franco [5], as equações parabólicas são adequadas para
modelar problemas de difusão, enquanto que as elı́pticas são para problemas de equilı́brio

10
e as hiperbólicas para problemas de convecção. A classificação das EDP’s lineares de 2a
ordem ocorrem em função do valor do discriminante da EDP. Uma equação diferencial
parcial linear é:

∂ 2u ∂ 2u
• Hiperbólica se ∆ = B 2 − 4AC > 0. Por exemplo, a equação da onda: = ;
∂x2 C 2 ∂t2
∂ 2u ∂ 2u
• Elı́ptica se ∆ = B 2 −4AC < 0. Por exemplo, a equação de Laplace: + = 0;
∂x2 ∂y 2

∂u ∂ 2u
• Parabólica se ∆ = B 2 − 4AC = 0. Por exemplo, a equação do calor: = α 2.
∂t ∂x

2.3 Equação do Calor

A equação do calor ou de difusão de calor é uma equação diferencial parcial parabólica de
segunda ordem. Segundo Incropera [6], essa equação determina o campo de temperatura,
ou seja, representa como a temperatura varia com a posição no meio. Ela envolve três
variáveis , duas variáveis independentes x e t, e uma variável dependente u (t, x). É da
forma:

∂u ∂ 2u
=α 2 (2.1)
∂t ∂x
Onde α ∈ R.
Através da equação do calor é possı́vel expressar um equilı́brio fı́sico fundamental, isto
é, a taxa de calor que entra em qualquer parte da barra é igual à taxa de absorção de
calor naquela parte da barra.
A forma geral, em coordenadas cartesianas da equação do calor é:
     
∂ ∂u ∂ ∂u ∂ ∂u ∂u
k + k + k + q = ρs
∂x ∂x ∂y ∂y ∂z ∂z ∂t

Onde k, ρ e s são constantes e representam, respectivamente, a condutividade térmica

(W/mK), densidade (Kg/m3 ) e calor especı́fico do material da barra (J/KgK). O termo
u expressa a temperatura e q a taxa na qual a energia é gerada por unidade de volume
no meio (W/m3 ). Essa equação postula que em qualquer ponto do meio, a taxa lı́quida
de transferência de energia por condução para o interior de um volume unitário somada
a taxa volumétrica de geração de energia térmica deve ser igual à taxa de variação da
energia acumulada no interior deste volume.

11
 O fluxo de calor é sempre no sentido desde um ponto de maior temperatura a pontos
de menor temperatura.
A equação do calor é de uma importância fundamental em numerosos e diversos cam-
pos da ciência. Na estatı́stica, a equação do calor está vinculada com o estudo do mo-
vimento browniano através da equação de Fokker-Planck. A equação de difusão, é uma
versão mais geral da equação do calor e relaciona-se, principalmente, com o estudo de
processos de difusão quı́mica. A equação do calor é usado em probabilidade e descreve
passeios aleatórios e é aplicada em matemática financeira por esta razão.
Equação do calor :



 ut = αuxx



 u(0, t) = 0
(2.2)


 u(L, t) = 0



 u(x, 0) = f (x)

Com 0 ≤ x ≤ L; e t ≥ 0
A solução do problema é dada pela função:
∞
X nπx n2 π2 2 αt
u (x, t) = A0 + An cos e L
n=1
L

Onde:
Z L Z L
2 2 nπx
A0 = f (x)dx e An = f (x)cos dx
L 0 L 0 L

Essa solução foi obtida usando o método de separação de variável e série de Fourier
(HOLETZ, 2001, p.4).

2.4 Erro de Truncamento

A noção de erro está presente em todos os campos do Cálculo Numérico. De um
lado, os dados, em si, nem sempre são exatos e, de outro lado, as operações sobre valo-
res não exatos propagam esses erros a seus resultados. Finalmente, os próprios métodos
numéricos, frequentemente métodos aproximados, buscam a minimização dos erros, pro-
curando resultados o mais próximo possı́vel do que seriam valores exatos.
Erro é a diferença entre o valor exato e o valor apresentado.

12
 O erro de truncamento é proveniente da utilização, por razões práticas de processos
finitos. Estes processos finitos são muito utilizados na avaliação de funções matemáticas,
tais como, exponenciação, logaritmos, funções trigonométricas e várias outras que uma
máquina pode ter.
Exemplo: Uma máquina poderia calcular a função seno e exponencial utilizando as
seguintes técnicas. Sendo:

x 3 x5 x7
sen(x) = x − + − + ...
3! 5! 7!
x2 x 3
ex = 1 + x + + + ...
2! 3!

Faz-se Truncamento:

x3 x5 x7 xn
sen(x) ∼
=x− + − + . . . + (−1)n
3! 5! 7! n!
x2 x3 xn
ex ∼
=1+x+ + + ... +
2! 3! n!

A ideia é a de interromper os cálculos quando uma determinada precisão é atingida.

De uma maneira geral, pode-se dizer que o erro de truncamento pode ser diminuı́do até
chegar a ficar da ordem do erro desejado; a partir desse ponto, não faz sentido diminuir-se
mais, pois o erro de arredondamento (diferença entre a aproximação calculada de um
número e o seu valor matemático exato) será dominante.

13
Capı́tulo 3

Métodos Numéricos

3.1 Discretização
Para tratar computacionalmente um problema diferencial é necessário expressar de
forma adequada a região (domı́nio) onde o problema será resolvido. Como usualmente não
é possı́vel obter soluções numéricas sobre o domı́nio, uma região contı́nua, devido à uma
infinidade de pontos envolvidos, inicialmente o domı́nio é discretizado, isto é, substituı́do
por um conjunto finito de pontos representativos (malha). Somente nesses pontos é que
as soluções serão obtidas.
A figura 3.1 representa uma malha unidimensional uniforme, onde o espaçamento da
malha é igual para cada intervalo e é representado por h.

Figura 3.1: Malha computacional unidimensional

Fonte: MELO, 2011, p.18.

A figura 3.2 representa uma malha cartesiana ortogonal uniforme para um problema
bidimensional. Qualquer ponto (xi , yi ) fica representado na malha por (i, j) e os vizinhos
a esse ponto vem representados por (i ± 1, j ± 1).

14
 Figura 3.2: Malha computacional bidimensional

Fonte: MELO, 2011, p.18.

Uma vez efetuada a discretização do domı́nio do problema, aplica-se o MDF (Método

das Diferenças Finitas) para a determinação das incógnitas. As derivadas, que aparecem
na equação original, são substituı́das (ou aproximadas) por fórmulas discretas de dife-
renças. A aplicação dessas fórmulas aos pontos do domı́nio discretizado gera um sistema
de equações algébricas, cuja solução fornece os valores das incógnitas do problema nesses
pontos discretos.
Daı́, nos métodos baseados em diferenças finitas completas o domı́nio contı́nuo da
equação é substituı́do por uma malha de pontos, as derivadas são aproximadas por dife-
renças finitas, e assim procura-se soluções numéricas na forma de uma tabela de valores
aproximados em determinados pontos do espaço e do tempo.
A exatidão da solução vai depender das dimensões dos passos escolhidos para o espaço
e para o tempo. Substituindo todas as derivadas parciais por diferenças finitas que resulta
num sistema de equações algébricas, o sistema poderá ser linear ou não dependendo do
tipo de EDP subjacente ao problema.
Nos problemas de valor inicial a solução é obtida partindo de valores iniciais e avançando
no tempo passo a passo, gerando uma sucessão de linhas na tabela da solução.
Os procedimentos baseados em passos de tempo podem ser explı́citos ou implı́citos,
dependendo da fórmula do valor solução usa-se ou não apenas informação relativa aos
pontos do passado.

15
3.1.1 Condições de Contorno

A seguir são apresentadas as condições de contorno que são comumente encontradas

na solução de equações diferenciais parciais:

• Condição de contorno de Dirichlet: especifica o valor da função no contorno;

• Condição de contorno de Neumann: especifica uma derivada normal da função no

contorno;

• Condição de contorno de Cauchy: especifica uma combinação dos dois tipos ante-
riores.

3.2 Método das Diferenças Finitas (MDF)

Utiliza-se o Método das Diferenças Finitas (MDF) como uma forma alternativa para
obtenção da solução aproximada de uma equação diferencial parcial. Tem como principal
objetivo transformar em um sistema de equações algébricas o resultado de uma equação
diferencial, substituindo as derivadas por diferenças.
No Método das Diferenças Finitas e de modo geral nos métodos numéricos, permitem-
se observações importantes em termos computacionais para matrizes de coeficientes que
serão produzidas quando o MDF for empregado. Para a programação de um software os
métodos numéricos são importantes, pois como já dito, não se restringem a nenhum caso
particular, podendo assim ser empregados de forma segura e precisa, dependendo do grau
de refinamento do cálculo.
Na engenharia, dificilmente se conhece a solução matemática analı́tica dos fenômenos
fı́sicos. As equações diferenciais que regem esses fenômenos são muitas vezes complicadas
e em geral não lineares. O Método das Diferenças Finitas consiste em resolver a equação
diferencial em pontos discretos. Estes pontos são igualmente espaçados, ou seja, a malha é
regular. O uso da técnica de Diferenças Finitas procura escrever os operadores diferenciais
em sua forma discreta, ou seja, em função de valores pontuais da solução.
Por definição, a derivada de um função u (x) em um ponto xi é dada por:

du  u (xi + h) − u (xi )
= lim (3.1)
dx x=xi h→ 0 h


16
 Onde h = ∆x . De forma aproximada, utilizando-se um incremento h , porém finito,
pode-se escrever:

• Diferença Finita Progressiva:

0 ui+1 − ui
ui ≈ (3.2)
h

• Diferença Finita Central:

00 ui−1 − 2ui + ui+1
ui ≈ (3.3)
h2

• Diferença Finita Regressiva:

0 ui − ui−1
ui ≈ (3.4)
h

Onde:
xn − x1
h= , xi = x1 + (i − 1) h
n−1
ui+1 = u (xi + h), ui = u (xi ), ui−1 = u (xi − h)

Figura 3.3: Diferença Finita Progressiva, Central e Regressiva

Fonte: ALVES, 2007, p.33.

O conhecimento da solução, mesmo que de forma aproximada, em alguns pontos dá

uma boa ideia da solução contı́nua, à medida que essa nuvem de pontos é adensada o
valor da resposta numérica se aproxima do valor real.
A formulação básica da série de Taylor dá inı́cio ao Método das Diferenças Finitas,
através da aproximação das derivadas.

17
3.3 Série de Taylor
Este método utiliza como técnica de solução de equações diferenciais, a substituição
das derivadas por formas de diferenças finitas que são obtidas pela expansão em série de
Taylor e truncamento ao nı́vel da ordem do erro desejada.

00
0 f (a) (x − a)2 f n−1 (a)n−1
f (x) = f (a) + f (a) f (x − a) + + ... + + Rn
2! (n − 1)!

O principal objetivo é determinar, de acordo com a função, uma boa aproximação

através de um polinômio a ela. Não será utilizado pontos de sua imagem, mas, sim, de
sua derivada de ordem n. Com os resultados, pode-se aproximar tão bem quanto se quer
a função por meio do polinômio.
Teorema (Expansão em Série de Taylor) Seja uma função f : D → Rinfinitamente
diferenciável f (x) qualquer. Então, tem-se:
00
0 f (a) 2 f n (a)
f (x) = f (a) + f (a)(x − a) + (x − a) + ... + (x − a)n + E(x)
2 n!

Onde E(x) é a função que representa o erro de truncamento e f (n) (x) é a n-ésima
derivada de f. Esta representação de f é dita centrada em a, pois o polinômo que aparece
antes do erro de truncamento assume valor igual a f (a), ou seja, o erro de a é igual a
zero. Uma maneira alternativa de representar a função, iterando nosso processo o quanto
quisermos, é:
∞
X f (n) (a) (x − a)n
f (x) =
n=0
n!

Demonstração: A intenção é achar um polinômio P (x)tal que este molde suficien-

temente bem a função, daı́:

p(a) = f (a)

0 0
p (a) = f (a)

00 00
p (a) = f (a)

000 000
p (a) = f (a)

..
.

18
 pn (a) = f n (a)

então tem-se que:

A1 (x − a)n + A2 (x − a)n−1 + ... + An−1 (x − a) + An ≈ f (x)

Deve-se determinar o coeficiente independente de P de modo que a função se aproxime

no ponto a. É natural tomar este coeficiente como f (a). Os outros é possı́vel conseguir
através de sucessivas derivações: Derivando até tornar o termo independente, obtém-se:
f (n) (a)
A1 =
n!
Integrando esta última equação, com o devido coeficiente lı́der em seu lugar, e fazendo
x = a, encontra-se A2 :
f (n−1) (a)
(n − 1)!A2 = f (n−1) (a) ⇔ A2 =
(n − 1)!
Dessa forma:
f (n−k) (a)
Ak =
(n − k)!
Logo, o método de aproximação polinomial que melhor molda a função em uma vizi-
nhança de um ponto a é:
00
0 f (a) f n (a)
f (x) = f (a) + f (a)(x − a) + (x − a)2 + ... + (x − a)n + E(x)
2 n!

3.4 Consistência, Convergência e Estabilidade

Para se garantir que a solução numérica fornecida por um esquema numérico repre-
sente uma aproximação razoável da solução exata do problema matemático, é necessário
que o esquema utilizado apresente propriedades de consistência, convergência e estabili-
dade. Essas propriedades estão inter-relacionadas na solução numérica e são funções dos
erros envolvidos.

• Consistência:

Um esquema de diferenças finitas é dito consistente quando, ao refinarem-se as apro-

ximações (por diferenças finitas), no limite, as equações aproximadas se tornam matema-
ticamente equivalentes às equações diferenciais originais. Assim, se ∆x → 0, ∆t → 0,
∆y → 0,  → 0 ( Onde  é o erro de truncamento da equação aproximada).

19
 • Convergência:

A convergência é uma qualidade importante e difı́cil de ser demonstrada diretamente,

em geral, utiliza-se uma técnica indireta para demonstrar a convergência de um método
numérico, baseada no Teorema de Equivalência de Lax: para a solução de um problema
linear de valor inicial bem posto, por uma discretização consistente, a estabilidade do
método numérico é condição necessária e suficiente para a convergência. Sintetizando:

Consistência + Estabilidade −→ Convergência

Desta forma, como foi visto, todas as discretizações são consistentes com as EDP’s e,
para os métodos estáveis e para os que têm critérios de estabilidade, se obedecidos, todos
os métodos numéricos serão convergentes.

• Estabilidade:

A estabilidade de métodos numéricos está intimamente associada ao erro numérico.

Um método de diferenças finitas é estável se os erros produzidos em um passo de tempo
do cálculo não provocam um aumento dos erros à medida que os cálculos avançam. Nos
esquemas que tratam de problemas evolutivos, a estabilidade condicional proporciona uma
solução numérica limitada. Os erros mais comuns introduzidos na solução são provenientes
de condições de fronteira ou iniciais aproximadas de forma incorreta, e também, devido
ao acúmulo de erros de arredondamento cometidos pelo computador durante os cálculos.
Em relação à estabilidade dos métodos numéricos, eles podem ser classificados em:

1. Condicionalmente estável : depende de certos parâmetros para que os erros per-

maneçam limitados e não instabilizem. Se os erros diminuı́rem ou chegarem a de-
saparecer, o método numérico é dito como sendo estável. Métodos explı́citos, em
geral, são desse tipo;

2. Incondicionalmente estável : não necessitam satisfazer quaisquer critérios de estabili-

dade para produzirem soluções estáveis. Em geral, métodos implı́citos (incluindo os
métodos semi-implı́citos) e alguns métodos explı́citos muito particulares pertencem
a essa categoria;

20
 3. Instável : se os erros aumentarem com o tempo, a solução numérica irá divergir
em relação à realidade (solução exata) e então o método numérico é dito como
sendo instável. Não existem valores de ∆t que permitem a eles fornecerem soluções
estáveis. Em princı́pio, não devem ser utilizados.

A estabilidade de métodos numéricos pode ser averiguada pela análise de estabilidade

de von Neumann. Para problemas dependentes do tempo, a estabilidade garante que o
método numérico produza uma solução limitada sempre que a solução da equação diferen-
cial for limitada. Estabilidade em geral, pode ser dificilmente averiguada, especialmente
se a equação em questão for não-linear.

3.5 Aplicação do Método Explı́cito na Equação do

Calor
O princı́pio do método de diferenças finitas, para obtenção da solução de uma dada
equação diferencial, é aproximar as derivadas desta equação por um conjunto de valores
da função em um determinado número de pontos.
Os MDF podem ser classificados em explı́citos e implı́citos. Os esquemas explı́citos de
diferenças finitas são simples e de fácil implementação. Nestes esquemas, a determinação
do passo de tempo da discretização temporal é dependente da condição de estabilidade
(condicionalmente estável) ka fim de obter uma solução estável.
A equação de condução de calor requer aproximações para a segunda derivada no
espaço e a primeira derivada no tempo. Utiliza-se a diferença central para representar a
aproximação da segunda derivada no espaço, da seguinte maneira:

∂ 2u uki+1 − 2uki + uki−1

= (3.5)
∂x2 (∆x)2
A primeira derivada no tempo pode ser representado utilizando a Diferença Progres-
siva, expressa da seguinte forma:

∂u uk+1 − uki
= i (3.6)
∂t ∆t
Substituindo essas equações em 2.1, obtém-se:

21
 uk+1 − uki uk − 2uki + uki−1
i
= α i+1 (3.7)
∆t (∆x)2

∆t
Considerando β = α , daı́:
(∆x)2

uk+1
i = βuki+1 + (1 − 2β) uki + βuki−1 (3.8)

Assim, a equação 3.8 fornece uma forma explı́cita para calcular os valores de cada nó,
baseado nos valores atuais e em seus vizinhos.

Figura 3.4: Método Explı́cito

Fonte: SOUZA, 2003, p.3.

Exemplo: Considere a seguinte equação diferencial parcial:

∂u ∂ 2u
=
∂t ∂x2

Com as condições de fronteiras:

u (t, 0) = u (t, 1) = 0, para t > 0

e Condição inicial:

u (0, x) = x (1 − x)2 , para 0 ≤ x ≤ 1

Usando intervalos na variável de espaço de 0.25 e na variável de tempo 0.01. Calcule

pelo Método Explı́cito o valor da solução numérica nos 5 nós para t = 0.01.
Solução:
Como x ∈ [0, 1], e tomando o valor de ∆x, tem-se:
Da condição inicial, obtém-se:

22
 Figura 3.5: Representação dos nós em subintervalos

Fonte: Elaborado pelas autoras

u00 = u04 = 0

u01 = u (0, 0.25) = 0.25 (0.75)2 = 0.140625

u02 = u (0, 0.5) = 0.5 (0.5)2 = 0.125

u03 = u (0, 0.75) = 0.75 (0.25)2 = 0.046875

Do enunciado temos ∆x = 0.25,∆t = 0.01 e α = 1, obtemos β = 0.16. Substituindo

na equação (3.8), para k = 0 e i = 1, 2, 3, temos:

i = 1; ⇒ u11 = 0.16u02 + 0.68u01 ⇒ u11 = 0.115625

i = 2; ⇒ u12 = 0.16u03 + 0.68u02 + 0.16u01 ⇒ u12 = 0.115

i = 3; ⇒ u13 = 0.68u03 + 0.16u02 ⇒ u13 = 0.051875

Portanto:



 u10 = 0


u1 = 0.115625


 1


u12 = 0.115 (3.9)


u13 = 0.051875






u14 = 0



3.6 Aplicação do Método Implı́cito na Equação do

Calor
O Método Implı́cito é incondicionalmente estável, ou seja, o passo de tempo é inde-
pendente da resolução espacial da grade adotada. Deve-se salientar que a estabilidade

23
garantida pelos métodos implı́citos não garante a precisão dos resultados, pois esta de-
pende da discretização espacial, do passo de tempo e das condições de contorno adotadas.
Além disso, eles geram um amortecimento nos resultados.

Figura 3.6: Método Implı́cito

Fonte: SOUZA, 2003, p.3.

No método das diferenças finitas implı́cito, a primeira derivada de tempo é aproximada

pela equação da diferença progressiva, enquanto que a aproximação da 2a derivada no
espaço é representada da seguinte maneira:

∂ 2u uk+1 k+1
i+1 − 2ui + uk+1
i−1
2
= 2 (3.10)
∂x (∆x)
Quando é substituı́do a expressão 3.6 e 3.10 na equação 2.1, obtém-se a equação a
seguir:

uk+1 − uki uk+1 − 2uk+1 + uk+1

i
= α i+1 i i−1
(3.11)
∆t (∆x)2
∆t
Considerando β = α , então:
(∆x)2

uki = −βuk+1 k+1

i+1 + (1 + 2β) ui − βuk+1
i−1 (3.12)

Essa equação possui várias incógnitas, por isso não pode ser resolvida por uma reor-
ganização algébrica simples, como a do método explı́cito.
As fórmulas implı́citas, juntamente com as condições de contorno resultam em um
conjunto de equações lineares algébricas com o mesmo número de incógnitas. Logo, o
método se reduz a solução de um conjunto de equações simultâneas em cada ponto de
tempo.

24
 A fundamental diferença entre os métodos explı́cito e implı́cito está representada na
figura a seguir:

Figura 3.7: Método Explı́cito e Implı́cito

Fonte: SOUZA, 2003, p.4.

Para a forma explı́cita, a derivada espacial é aproximada no nı́vel k de tempo. Já para
a forma implı́cita, a derivada espacial aproximada em k + 1, nı́vel de tempo avançado.
Exemplo: Considere a seguinte equação diferencial parcial:

∂u ∂ 2u
=
∂t ∂x2

Com as condições de fronteiras:

u (t, 0) = u (t, 1) = 0, para t > 0

e Condição inicial:

u (0, x) = x (1 − x)2 , para 0 ≤ x ≤ 1

Usando intervalos na variável de espaço de 0.25 e na variável de tempo 0.01. Calcule

pelo Método Implı́cito o valor da solução numérica nos 5 nós para t = 0.01.
Solução:
Como x ∈ [0, 1], e tomando o valor de ∆x, tem-se:
Da condição inicial, obtém-se:

u00 = u04 = 0

u01 = u (0, 0.25) = 0.25 (0.75)2 = 0.140625

u02 = u (0, 0.5) = 0.5 (0.5)2 = 0.125

25
 Figura 3.8: Representação dos nós em subintervalos

Fonte: Elaborado pelas autoras

u03 = u (0, 0.75) = 0.75 (0.25)2 = 0.046875

Do enunciado temos ∆x = 0.25,∆t = 0.01 e α = 1, obtemos β = 0.16. Substituindo

na equação (3.11), para k = 0 e i = 1, 2, 3, temos:

i = 1; ⇒ 1.32u11 − 0.16u12 = 0.140625

i = 2; ⇒ −0.16u11 + 1.32u12 − 0.16u13 = 0.125

i = 3; ⇒ −0.16u12 + 1.32u13 = 0.046875

Na forma Matricial tem-se:

     
1.32 −0.16 0 u11 0.140625
     
 −0.16 1.32 −0.16  .  u12  =  0.125 (3.13)
     

     
1
0 −0.16 1.32 u3 0.046875

Da resolução da matriz tridiagonal e das condições de contorno, tem-se:




 u10 = 0


u1 = 0.1205


 1


u12 = 0.1153 (3.14)


u13 = 0.495






u14 = 0



26
3.7 Método de Crank - Nicholson
Na análise numérica, o método de Crank Nicholson é um método de diferenças finitas
usado para resolver numericamente a equação do calor e equações diferenciais parciais
similares. É um método de segunda ordem no tempo e no espaço, implı́cito no tempo e é
numericamente estável. O método foi desenvolvido por John Crank e Phyllis Nicholson na
metade do século XX. Nele é proposto um esquema implı́cito alternativo para chegar mais
próximo do valor real da derivada de segunda ordem no espaço e primeira ordem no tempo
da equação de difusão. Para alcançar esse objetivo, aproximações para as diferenças são
implementadas no ponto médio do incremento de tempo.

Figura 3.9: Representação do Método de Crank-Nicholson

Fonte: SOUZA, 2003, p.5.

As soluções aproximadas, encontradas por meio desse método, podem ainda conter
oscilações significativas caso a razão entre o passo de tempo e o quadrado do passo de
espaço for grande (usualmente maior que 1/2). Por essa razão, sempre que grandes passos
de tempo forem tomados, o método menos preciso de euler implı́cito é frequentemente
utilizado, o qual é estável e imune à oscilações.
O método de Crank-Nicholson é baseado em diferenças centradas no espaço, e na regra
trapezoidal no tempo, é de segunda no tempo e no espaço. Por exemplo, para um caso
unidimensional, se a equação diferencial parcial for:

∂u ∂ 2 u
 
∂u
=F u, x, t, , 2 (3.15)
∂t ∂x ∂x
Em seguida, fazendo: u (i∆x, k∆t) = uki . A equação para o método de Crank -
Nicholson é a combinação do método de Euler explı́cito em k e do método de Euler
implı́cito em k + 1 (deve-se notar, contudo, que o método por si só não é simplesmente a

27
média desses dois métodos, já que a equação tem uma dependência implı́cita na solução:

• Euler Explı́cito:

uk+1 − uki ∂u ∂ 2 u
 
i k
= Fi u, x, t, , 2 (3.16)
∆t ∂x ∂x

• Euler Implı́cito:

uk+1 − uki ∂u ∂ 2 u
 
i k+1
= Fi u, x, t, , 2 (3.17)
∆t ∂x ∂x
De 3.16 e 3.17, tem-se o Método de Crank - Nicholson:

uk+1 − uki ∂u ∂ 2 u ∂u ∂ 2 u
    
i 1 k+1 k
= (Fi u, x, t, , 2 + Fi u, x, t, , 2 (3.18)
∆t 2 ∂x ∂x ∂x ∂x

A função F deve ser discretizada no espaço por diferenças centradas. Note que este
é um método implı́cito. Para conseguir o valor posterior de u no tempo, um sistema de
equações algébricas deve ser resolvido. Se a equação diferencial parcial for não-linear, a
discretização também deverá ser não-linear para que o avanço no tempo envolva a solução
do sistema algébrico não-linear de equações, embora que linearizações sejam possı́veis.
O método de Crank-Nicholson é frequentemente aplicado a problemas da equação do
calor. Como por exemplo:

∂u ∂ 2u
=α 2 (3.19)
∂t ∂x
A discretização de Crank-Nicholson é dada por:

uk+1
i − uki α  k k k

k+1 k+1 k+1


= ui+1 − 2u i + u i−1 + ui+1 − 2u i + u i−1 (3.20)
∆t 2 (∆x)2

α∆t
Fazendo β = , obtém-se:
2 (∆x)2

uk+1 = uki + β uki+1 − 2uki + uki−1 + uk+1 k+1

+ uk+1



i i+1 − 2ui i−1 (3.21)

Então:

28
 −βuk+1 k+1
i+1 + (1 + 2β) ui − βuk+1 k k k
i−1 = βui+1 + (1 − 2β) ui + βui−1 (3.22)

Exemplo:
Considere a seguinte equação diferencial parcial:

∂u ∂ 2u
=2 2
∂t ∂x

Com as condições de fronteiras:

u (t, 0) = u (t, 1) = 0, para t > 0

e Condição inicial:

u (0, x) = 4 (1 − x) x, para 0 ≤ x ≤ 1

Usando intervalos na variável de espaço de 0.25 e na variável de tempo 0.0625. Calcule

pelo método de Crank-Nicholson o valor da solução numérica nos 5 nós para t = 0.0625.
Solução:
Como x ∈ [0, 1], e tomando o valor de ∆x, tem-se:

Figura 3.10: Representação dos nós em subintervalos

Fonte: Elaborado pelas autoras

Da condição inicial, obtém-se:

u00 = u04 = 0

u01 = u (0, 0.25) = 4 (0.75) 0.25 = 0.75

u02 = u (0, 0.5) = 4 (0.5) 0.5 = 1

u03 = u (0, 0.75) = 4 (0.25) 0.75 = 0.75

29
 Do enunciado temos ∆x = 0.25,∆t = 0.0625 e α = 2, obtemos β = 1. Substituindo
na equação (3.19), para k = 0 e i = 1, 2, 3, temos:

i = 1; ⇒ −u12 + 3u11 = 0, 25

i = 2; ⇒ −u13 + 3u12 − u11 = 0, 5

i = 3; ⇒ 3u13 − u12 = 0, 25

Na forma Matricial tem-se:

     
3 −1 0 u11 0, 25
     
 −1 3 −1  .  u12  =  0, 5 (3.23)
     

     
1
0 −1 3 u3 0, 25

Da resolução da matriz tridiagonal e das condições de contorno, tem-se:




 u10 = 0


u1 = 0.182


 1


u12 = 0.287 (3.24)


u13 = 0.179






 u1 = 0

4

30
Capı́tulo 4

Problema Proposto

0
Considere uma viga, com comprimento de 10 cm e os seguintes valores: k = 0, 49cal/
(s.cm.o C), ∆x = 2cm e ∆t = 0, 1s. A t = 0, a temperatura da viga é igual a zero e as
condições de contorno são fixas para todos os tempos a T (0) = 100o C e T (10) = 50o C.
Observe que a viga é de alumı́nio com C = 0, 2174cal/(g o C) e ρ = 2, 7g/cm3 . Portanto,
α = 0, 835cm2 /s e β = 0, 020875. Resolva numericamente a distribuição da temperatura
nessa viga utilizando os métodos: Explı́cito, Implı́cito e de Crank-Nicholson. Comente os
resultados.

Figura 4.1: Distribuição de temperaturas em uma viga de alumı́nio

Fonte: Elaborado pelas autoras

Observação:
A temperatura u = u(x, t) da viga de comprimento L = 10cm, satisfaz a equação do
calor:
∂u ∂ 2u
=α 2
∂t ∂x
Com as condições de fronteira:

u (0, t) = 100 e u (10, t) = 50 para t > 0

E com condições iniciais:

u(x, 0) = 0 para 0 < x < 10

31
 A difusividade térmica dada por é α = 0, 835cm2 /s e a solução analı́tica deste problema
e expressa por:

  ∞
T2 − T1 −n2 π 2 αt
 nπx 
X  
u (x, t) = x + T1 + Cn e L2 sen (4.1)
L n=1
L

Onde:
Z L
2 h x i  nπx 
Cn = f (x) − (T2 − T1 ) − T1 sen dx
L 0 L L
Sendo L é o comprimento da barra. Como os resultados precisos podem ser obtidos,
geralmente, utilizando-se apenas alguns termos da série, é possı́vel calcular a solução
analı́tica do problema usando a série acima (HOLETZ, 2001, p.98).
Para o cálculo da solução do problema pelo método das diferenças finitas, utiliza-se a
premissa de que o número de Fourier β deve ser igual ou menor que 0, 5 para o método
convergir.
A seguir será apresentado a solução do problema por meio dos métodos: Explı́cito,
Implı́cito e Crank-Nicholson.

32
 • Método Explı́cito

Solução:
Aplicando a equação 3.8 para ∆x = 2cm, ∆t = 0.1s, i = 1, 2, 3, 4 e β constante.

uk+1
i = βuki+1 + (1 − 2β) uki + βuki−1

Para t = 0.1s e k = 0, o valor do nó em x = 2cm é dado por:

u11 = 0, 020875u02 + (1 − 2.0, 020875) u01 + 0, 020875u00

u11 = 0, 020875.100 = 2, 0875

Os valores incógnitos da função são aqueles correspondentes aos nós em x = 4, 6 e

8cm, daı́:

u12 = 0

u13 = 0

u14 = 1, 04375

Para t = 0.2s e k = 1, os valores para os quatro nós interiores são:

u21 = 4, 0878

u22 = 0, 0436

u23 = 0, 02179

u24 = 2, 0439

Para t = 0.3s e k = 2, os valores para os nós no interior são:

u31 = 6, 0055

u32 = 0, 1276

u33 = 0, 06445

u34 = 3, 0028

33
 Para t = 0.4s e k = 3, tem-se:

u41 = 7, 8449

u42 = 0, 2489

u43 = 0, 1271

u44 = 3, 9225

Após utilizar a equação do método explı́cito, obtém-se o comportamento da tempera-

tura para os diferentes momentos de tempo e comprimento da barra como mostrado na
tabela:

Tempo(s) Comprimento (cm)

0 2 4 6 8 10
0 100 0 0 0 0 50
0,1 100 2,0875 0 0 1,04375 50
0,2 100 4,0878 0,0436 0,02179 2,0439 50
0,3 100 6,0055 0,1276 0,06445 3,0028 50
0,4 100 7,8449 0,2489 0.1271 3,9225 50
Fonte: Elaborado pelas autoras

Através da tabela observou-se que o problema não possui simetria em relação ao eixo
médio da viga. Após a obtenção das soluções foram desenvolvidos gráficos que facilitam
o estudo da distribuição de temperatura no decorrer do tempo e comprimento da viga.
Para representar a condição inicial do problema, onde as extremidades são mantidas a
100o C e 50o C, respectivamente, e o comprimento da viga a 0o C, foi desenvolvido o gráfico
a seguir, que relaciona a distribuição de temperatura no decorrer do comprimento da viga.

34
 Figura 4.2: Variação da temperatura em função do espaço

Fonte: Elaborado pelas autoras

Os gráficos abaixo mostram a variação de temperatura no comprimento da viga cal-

culado pelo método explı́cito para os tempos de: 0.1s, 0.2s, 0.3s e 0.4s, respectivamente.

Figura 4.3: Solução Numérica pelo Método Explı́cito para o tempo de 0.1s à 0.4s.

Fonte: Elaborado pelas autoras

Em suma, é possivél representar o comportamento da temperatura no decorrer do

tempo na viga, em apenas um gráfico como mostrado a seguir:

35
Figura 4.4: Variação da temperatura em relação do espaço, para diferentes momentos de
tempos, calculado pelo Método Explı́cito.

Fonte: Elaborado pelas autoras

• Método Implı́cito

Solução:
Utilizando a equação 3.12 para ∆x = 2cm, ∆t = 0.1s e β = 0.020875, i = 1, 2, 3, 4
tem-se:

uki = −βuk+1 k+1

i+1 + (1 + 2β) ui − βuk+1
i−1

Para t = 0.1s e k = 0, o valor do nó em x = 2cm é dado por:

1, 04175u11 − 0, 020875u12 = 2, 0875

−0, 020875u11 + 1, 04175u12 − 0, 020875u13 = 0

−0, 020875u12 + 1, 04175u13 − 0, 020875u14 = 0

−0, 020875u13 + 1, 04175u14 = 1, 04375

Na forma matricial:

     
1, 04175 −0, 020875 0 0 u112, 0875
     
  1  
 −0, 020875 1, 04175 −0, 020875
 
0   u2   0 
 . = 
  1  
−0, 020875 1, 04175 −0, 020875
 
 0   u3   0 
     
0 0 −0, 020875 1, 04175 u14 1, 04375
(4.2)

36
Logo:

u11 = 2, 0046

u12 = 0, 0406

u13 = 0, 0209

u14 = 1, 0023

Para t = 0.2s e k = 1 tem-se:

1, 04175u21 − 0, 020875u22 = 4, 0922

−0, 020875u21 + 1, 04175u22 − 0, 020875u23 = 0, 0406

−0, 020875u22 + 1, 04175u23 − 0, 020875u24 = 0, 0209

−0, 020875u23 + 1, 04175u24 = 2, 04605

Na forma matricial:

     
1, 04175 −0, 020875 0 0 u21 4, 0922
     
  2  
 −0, 020875 1, 04175 −0, 020875
 
0   u2   0, 0406 
.
  2 =
   
−0, 020875 1, 04175 −0, 020875
 
 0   u3   0, 0209 
     
0 0 −0, 020875 1, 04175 u24 2, 04605
(4.3)
Logo:

u21 = 3, 9306

u22 = 0, 1189

u23 = 0, 0618

u24 = 1, 9653

37
Para t = 0.3s e k = 2 tem-se:

1, 04175u31 − 0, 020875u32 = 6, 0181

−0, 020875u31 + 1, 04175u32 − 0, 020875u33 = 0, 1189

−0, 020875u32 + 1, 04175u33 − 0, 020875u34 = 0, 0618

−0, 020875u33 + 1, 04175u34 = 3, 00905

Na forma matricial:

     
1, 04175 −0, 020875 0 0 u316, 0181
     
  3  
 −0, 020875 1, 04175 −0, 020875
 
0   u2   0, 1189 
.
  3 =
   
−0, 020875 1, 04175 −0, 020875
 
 0   u3   0, 0618 
     
0 0 −0, 020875 1, 04175 u34 3, 00905
(4.4)
Logo:

u31 = 5, 7816

u32 = 0, 2324

u33 = 0, 1219

u34 = 2, 8909

Para t = 0.4s e k = 3:

1, 04175u41 − 0, 020875u42 = 7, 8691

−0, 020875u41 + 1, 04175u42 − 0, 020875u43 = 0, 2324

−0, 020875u42 + 1, 04175u43 − 0, 020875u44 = 0, 1219

1, 04175u44 − 0, 020875u43 = 3, 93465

38
 Na forma matricial:

     
1, 04175 −0, 020875 0 0 u41 7, 8691
     
  4  
 −0, 020875 1, 04175 −0, 020875
 
0   u2   0, 2324 
.
  4 =
   
−0, 020875 1, 04175 −0, 020875
 
 0   u3   0, 1219 
     
4
0 0 −0, 020875 1, 04175 u4 3, 93465
(4.5)
Logo:

u41 = 7, 561386

u42 = 0, 3786179

u43 = 0, 2003662

u44 = 3, 7809769

A seguir é dado o comportamento da temperatura para diferentes momentos de tempo

e comprimento da viga.

Tempo(s) Comprimento (cm)

0 2 4 6 8 10
0 100 0 0 0 0 50
0,1 100 2,0047 0,0406 0,0209 1,0023 50
0,2 100 3,9306 0,1189 0,0618 1,9653 50
0,3 100 5,7816 0,2324 0,1219 2,8909 50
0,4 100 7,5614 0,3787 0,200366 3,78098 50
Fonte: Elaborado pelas autoras

39
 Analisando o gráfico para t = 0, 1s, pode-se observar que há uma pequena variação da
solução númerica nos pontos entre 2 e 3cm de comprimento. Essa variação será atenuada
com a resolução do problema, onde o resultado numérico apresenta uma mı́nima oscilação,
conforme mostrado nos gráficos a seguir:

Figura 4.5: Solução Numérica pelo Método Implı́cito para o tempo de 0,1s à 0,4s.

Fonte: Elaborado pelas autoras

O gráfico (4.5),a seguir, descreve como ocorre o aumento da temperatura no compri-

mento da barra com o decorrer do tempo. Nota-se que no instante t = 0, 2s, a região
central da viga está com uma temperatura de, aproximadamente, 10o C. Com o passar
do tempo essa região terá um aumento de temperatura. Em t = 0, 4s observa-se que a
temperatura na região central aproxima-se de 20o C.

40
Figura 4.6: Variação da temperatura em relação do espaço, para diferentes momentos de
tempos, calculado pelo Método Implı́cito.

Fonte: Elaborado pelas autoras

• Método de Crank-Nicholson

Solução:
Aplicando a equação 3.22 para ∆x = 2cm, ∆t = 0.1s, β = 0.0104375 e i = 1, 2, 3, 4

−βuk+1 k+1
i+1 + (1 + 2β) ui − βuk+1 k k k
i−1 = βui+1 + (1 − 2β) ui + βui−1

Para t = 0.1s e k = 0, o valor do nó em x = 2cm é dado por:

1, 020875u11 − 0, 0104375u12 = 2, 0875

−0, 0104375u11 + 1, 020875u12 − 0, 0104375u13 = 0

−0, 0104375u12 + 1, 020875u13 − 0, 0104375u14 = 0

−0, 0104375u13 + 1, 020875u14 = 1, 04375

Na forma matricial:

     
1, 020875 −0, 0104375 0 0 u11 2, 0875
     
  1  
 −0, 0104375 1, 020875 −0, 0104375
 
0   u2   0 
 . = 
  1  
−0, 0104375 1, 020875 −0, 0104375
 
 0   u3   0 
     
1
0 0 −0, 0104375 1, 020875 u4 1, 04375
(4.6)

41
 Logo:

u11 = 2, 0450294

u12 = 0, 0210176

u13 = 0, 0106692

u14 = 1, 0225163

Para t = 0.2s e k = 1, tem-se:

1, 020875u21 − 0, 0104375u22 = 4, 090058782475

−0, 0104375u21 + 1, 020875u22 − 0, 0104375u23 = 0, 0420352117375

−0, 0104375u22 + 1, 020875u23 − 0, 0104375u24 = 0, 02133836553125

1, 020875u24 − 0, 0104375u23 = 2, 0450326320125

Na forma matricial:

     
1, 020875 −0, 0104375 0 0 u21 4, 090058782475
     
  2  
 −0, 0104375 1, 020875 −0, 0104375
 
0   u2   0, 0420352117375 
.
  2 =
   
−0, 0104375 1, 020875 −0, 0104375
 
 0   u3   0, 02133836553125 
     
0 0 −0, 0104375 1, 020875 u24 2, 0450326320125
(4.7)
Logo:

u21 = 4, 007269

u22 = 0, 0825781

u23 = 0, 0422318

u24 = 2, 0036473

42
 Para t = 0.3s e k = 2, tem-se:

1, 020875u31 − 0, 0104375u32 = 6, 01197916854375

−0, 0104375u31 + 1, 020875u32 − 0, 0104375u33 = 0, 1231209467625

−0, 0104375u32 + 1, 020875u33 − 0, 0104375u34 = 0, 0631251887875

1, 020875u34 − 0, 0104375u33 = 3, 006011957025

Na forma matricial:

     
1, 020875 −0, 0104375 0 0 u31 6, 01197916854375
     
  3  
 −0, 0104375 1, 020875 −0, 0104375
 
0   u2   0, 1231209467625 
.
  3 =
   
−0, 0104375 1, 020875 −0, 0104375
 
 0   u3   0, 0631251887875 
     
0 0 −0, 0104375 1, 020875 u34 3, 006011957025
(4.8)
Logo:

u31 = 5, 890904

u32 = 0, 1817915

u33 = 0, 0938081

u34 = 2, 9455037

Para t = 0.4s e k = 3, tem-se:

1, 020875u41 − 0, 0104375u42 = 7, 85732882778125

−0, 0104375u41 + 1, 020875u42 − 0, 0104375u43 = 0, 24046203498125

−0, 0104375u42 + 1, 020875u43 − 0, 0104375u44 = 0, 1244909995625

1, 020875u44 − 0, 0104375u43 = 3, 92874543230625

43
 Na forma matricial:

     
1, 020875 −0, 0104375 0 0 u41 7, 85732882778125
     
  4  
 −0, 0104375 1, 020875 −0, 0104375
 
0   u2   0, 24046203498125 
.
  4 =
   
−0, 0104375 1, 020875 −0, 0104375
 
 0   u3   0, 1244909995625 
     
4
0 0 −0, 0104375 1, 020875 u4 3, 92874543230625
(4.9)
Logo:

u41 = 7, 6998913

u42 = 0, 3159515

u43 = 0, 1645393

u44 = 3, 8500921

Obteve-se o comportamento da temperatura para os diferentes momentos de tempo e

comprimento da viga, como mostrado a seguir:

Tempo(s) Comprimento (cm)

0 2 4 6 8 10
0 100 0 0 0 0 50
0,1 100 2,0450294 0,0210176 0,0106692 1,0225163 50
0,2 100 4,00726 0,0825781 0,0422318 2,0036473 50
0,3 100 5,890904 0,1817915 0,0938081 2,9455037 50
0,4 100 7,6998913 0,3159515 0,1645393 3,8500921 50
Fonte: Elaborado pelas autoras

44
 Os gráficos da figura (4.6) mostram a temperatura variando no decorrer da viga nos
tempos: 0.1s, 0.2s, 0.3s, e 0.4s. Pode-se representar o comportamento da temperatura em
função do comprimento da viga para as variações do tempo, em apenas um gráfico, como
mostrado na figura (4.7).

Figura 4.7: Solução Numérica pelo Método de Crank-Nicholson para o tempo de 0,1s à
0,4s.

Fonte: Elaborado pelas autoras

Figura 4.8: Variação da temperatura em relação do espaço, para diferentes momentos de

tempos, calculado pelo Método de Crank-Nicholson.

Fonte: Elaborado pelas autoras

45
 Representação de diferentes soluções para u em x = 2cm e t = 10s.

∆t β Explı́cito Implı́cito Crank-Nicholson

10 2,0875 208,75 53,01 79,77
5 1,04375 -9,13 58,48 64,79
2 0,4175 67,12 62,22 64,87
1 0,20875 65,91 63,49 64,77
0,5 0,104375 65,33 64,12 64,74
0,2 0,04175 64,97 64,49 64,73
Fonte: Elaborado pelas autoras

A tabela apresenta resultados numéricos para a temperatura em x = 2 e t = 10.

Note que uma gama de intervalos são utilizados, os resultados indicam uma série das
propriedades dos métodos numéricos.
Em primeiro lugar, percebe-se que o método explı́cito é instável para os valores
elevados de β. Essa instabilidade não é manisfestada por qualquer abordagem implı́cita.
Em segundo lugar, o método de Crank-Nicholson converge mais rapidamente, pois β é
reduzido e proporciona resultados moderadamente precisos.
Finalmente, se ∆x e ∆t forem diminuı́dos como β foi diminuido (isto é, forem utilizados
segmentos mais espaciais), a solução numérica se aproxima do resultado analı́tico.

46
Considerações Finais

Este trabalho, teve como objetivo principal analisar a transferência de calor no decor-
rer do tempo, em um problema que envolve a condução de calor unidimensional em uma
viga, utilizando os métodos: Explı́cito, Implı́cito e de Crank-Nicholson.
Na questão apresentada, observou-se que é possı́vel obter resultados coerentes e pre-
cisos utilizando o método das diferenças finitas, seguindo o critério de convergência dele.
E fazendo a discretização, verificou-se que os termos das formulações geram matrizes
tridiagonais, para os quais a resolução dos sistemas é muito fácil de ser efetuada compu-
tacionalmente, utilizando programas como o SCILAB.
Outro programa utilizado foi Microsoft Excel. Com ele determinou-se as simulações
dos gráficos, e a partir destes notou-se que houve uma mı́nima variação na temperatura
com relação ao comprimento, no decorrer dos tempos empregados.
Verificou-se que o comportamento dos gráficos, tanto dos métodos de Euler quanto o de
Crank-Nicholson, se assemelham devido os valores obtidos serem bastante aproximados.
Ao término deste trabalho, foi possı́vel observar que as equações diferenciais parciais
representam um modelo fı́sico, por isso é necessário adicionar condições auxiliares de modo
a caracterizar melhor a situação modelada. Essas condições são denominadas, dependendo
do problema, de condições iniciais e de contorno, e a sua adequada implementação é
fundamental para se obter uma solução numérica de boa qualidade.

47
Referências Bibliográficas

[1] ALVES,L.M. Introdução aos Métodos Aproximados em Engenharia: Álgebra Linear,

Geometria Analı́tica, Cálculo e Equações Diferenciais. 2007, 287p. Programa de Pós-
Graduação em Métodos Numéricos em Engenharia- Universidade Federal do Paraná,
Curitiba.

[2] BILLO, E.J. Excel for Scientists and Engineers. Manhattan: Wiley, 2007.

[3] BOYER, C.B. Tópicos de História da Matemática para Uso em Sala de Aula. São
Paulo: Atual, 1992.

[4] CHAPRA, S. C.; CANALE, R. P. Numerical Methods for Engineers. 6 edição. New
York: McGraw-Hill,2010.

[5] FRANCO, N. B. Cálculo numérico. São Paulo: Person Prentice Hall, 2006.

[6] HOLETZ, M. S. Método de Fourier para a resolução de problemas de valores ini-

cial e fronteira para a Equação do Calor. 2001, 105p. Centro de Ciências Fı́sicas e
Matemática - Universidade Federal de Santa Catarina.

[7] INCROPERA, F. P. Z. Fundamentos de Transferência de calor e de massa. Tradução

e revisão técnica Eduardo Queiroz, Fernando Luiz Pellegrni Pessoa. Rio de Janeiro :
LT C, 2008.

[8] JUSTO, D.A.R. Métodos Numéricos Para Equações Diferenciais. Rio Grande do Sul:
2008.

[9] LEVEQUE, R. J. Finite Diference Methods for Diferential Equations. Washington:

University of Washington, 1998.230p.

[10] SOARES, A.A.B. O Método das Diferenças Finitas Aplicado à Teoria das Vigas.
2010, 142p. Universidade da Amazônia, Belém.

48
[11] SOUSA, N. Método de Diferenças Finitas: Conceitos e Interpretações. 2009, 89p.
Mestrado Profissional em Matemática - Instituto de Matemática, Estatı́stica e Com-
putação Cientı́fica, Campinas.

49