CÁLCULO DA TEMPERATURA NO INTERIOR DE UMA PLACA

COM AUXÍLIO DO MÉTODO DAS DIFERENÇAS FINITAS

Adilandri Mércio Lobeiro – alobeiro@utfpr.edu.br

Viviane Colucci Boromelo – colucci@utfpr.edu.br

Departamento Acadêmico de Matemática (DAMAT)

Cleber Eduardo Fernandes Leal – cleber.efl@hotmail.com

Giovani de Madureira Alves Sobrinho – giovani.eng_civil@hotmail.com

Marlon Vieira Passos – mvpassos25@hotmail.com

Mariana Soares Ribeiro – ribeiromaso@gmail.com

Departamento Acadêmico de Construção Civil (DACOC)

Universidade Tecnológica Federal do Paraná (UTFPR),

Via Rosalina Maria dos Santos, nº 1233

CEP 87301-899 – Campo Mourão – Paraná

Liliana Madalena Gramani - gramani@ufpr.br

Eloy Kaviski - eloydhs@ufpr.br

Departamento de Matemática, Universidade Federal do Paraná - UFPR,

Programa de Pós-Graduação em Métodos Numéricos em Engenharia - PPGMNE.

Centro Politécnico - Caixa Postal 19011

CEP 81531-970 - Curitiba - Paraná

Resumo: Este trabalho tem por objetivo analisar a eficiência do Método das Diferenças Finitas

quando aplicado à um problema comum à física e às engenharias. Trata-se da distribuição de

temperatura ao longo de uma placa quadrada e delgada. O artigo conta com a descrição do

método numérico empregado, incluindo sua dedução e forma de aplicação, evidenciando sua

praticidade e facilidade de aplicação, visando incentivar seu uso e disseminação. Consta

também uma introdução à equação de Poisson, cuja solução traduz o valor da temperatura em

cada ponto do interior da placa sob as condições dadas. Um algoritmo foi implementado e

empregado na resolução de um estudo de caso.

Palavras-chave: Distribuição de Temperatura, Método das Diferenças Finitas, Equação de

Poisson.

1. INTRODUÇÃO

O emprego de softwares como ferramenta de apoio na resolução de problemas nas

diversas áreas da ciência e das engenharias tem se tornado algo cada vez mais frequente. O que

não é de surpreender, visto que seu uso traz inúmeros benefícios, fornecendo praticidade e

precisão aos cálculos e projetos elaborados.

Mesmo contando com essas ferramentas o domínio da teoria é indispensável, pois

resultados confiáveis só podem ser obtidos com o conhecimento de todas as variáveis que

envolvem o problema e com a inserção correta de dados no mesmo.

Este trabalho exemplifica uma das formas de aplicação computacional na área das

ciências exatas, trazendo também a teoria de um Método Matemático nem sempre ministrado

ao longo do curso de Engenharia Civil, o Método das Diferenças Finitas (MDF).

Propõe-se por meio de um estudo de caso referente à análise da temperatura em uma

placa, apresentar as vantagens de se estudar o MDF, procurando incentivar seu ensino e uso na

abordagem de outros problemas pertencentes às demais áreas da física e engenharias.

2. REVISÃO BIBLIOGRÁFICA

Esta seção tem por objetivo apresentar conceitos relativos às Equações Diferenciais

Parciais (EDP’s) e sua classificação, ao Método das Diferenças Finitas e à teoria de propagação

de calor e distribuição da temperatura.

2.1. Equações diferenciais parciais de segunda ordem

Seja Ω ⊂ 𝑅2 um conjunto aberto. Consideremos os seguintes espaços vetoriais:

𝐶0(Ω) = 𝑢: Ω ⟶ 𝑅, 𝑢 é 𝑐𝑜𝑛𝑡í𝑛𝑢𝑎; 𝐶2(Ω) = 𝑢: Ω ⟶ 𝑅, 𝑢 é 𝑑𝑢𝑎𝑠 𝑣𝑒𝑧𝑒𝑠 𝑐𝑜𝑛𝑡𝑖𝑛𝑢𝑢𝑎𝑚𝑒𝑛𝑡𝑒 𝑑𝑖𝑓𝑒𝑟𝑒𝑛𝑐𝑖á𝑣𝑒𝑙. Definimos o operador diferencial da Equação (1)

𝐿: 𝐶2(Ω) ⟶ 𝐶0(Ω)

𝑢 ⟼ 𝐿𝑢 = 𝐴(𝑥, 𝑦)𝜕2𝑢

𝜕𝑥2(𝑥, 𝑦) + 𝐵(𝑥, 𝑦)

𝜕2𝑢

𝜕𝑥𝜕𝑦(𝑥, 𝑦) + 𝐶(𝑥, 𝑦)

𝜕2𝑢

𝜕𝑦2(𝑥, 𝑦)

(1)

onde 𝐴, 𝐵, 𝐶: Ω ⟶ 𝑅 são funções reais que dependem das variáveis independentes 𝑥 e 𝑦. Além

disso, para todo (𝑥, 𝑦) ∈ Ω pelo menos um dos coeficientes,𝐴, 𝐵 𝑜𝑢 𝐶 é não nulo, ou seja,

𝐴2(𝑥, 𝑦) + 𝐵2(𝑥, 𝑦) + 𝐶2(𝑥, 𝑦) > 0.

Definimos a função dada pela Equação (2)

𝐹: Ω × 𝑅3 ⟶ 𝑅

(𝑥, 𝑦, 𝜉, 𝜂, 𝜍) ⟼ 𝐹(𝑥, 𝑦, 𝜉, 𝜂, 𝜍) (2)

Definição 1 (EDP de Segunda Ordem Quase Linear) Denomina-se equação diferencial

parcial de segunda ordem, quase linear, na incógnita 𝑢(𝑥, 𝑦), a uma equação com a forma da

Equação (3)

𝐿𝑢 = 𝐹 (𝑥, 𝑦, 𝑢,𝜕𝑢

𝜕𝑥,𝜕𝑢

𝜕𝑦), (3)

sendo 𝐿 e 𝐹 definidas anteriormente, onde os coeficientes 𝐴, 𝐵 𝑒 𝐶 das derivadas de segunda

ordem devem, somente depender das variáveis independentes 𝑥 e 𝑦, e satisfazer a condição

dada pela Equação (4)

𝐴2(𝑥, 𝑦) + 𝐵2(𝑥, 𝑦) + 𝐶2(𝑥, 𝑦) > 0 (4)

(MEDEIROS et al, 1999).

A Equação (5) mostra a forma geral desse tipo de EDP.

𝐴(𝑥, 𝑦)𝜕2𝑢

𝜕𝑥2(𝑥, 𝑦) + 𝐵(𝑥, 𝑦)

𝜕2𝑢

𝜕𝑥𝜕𝑦(𝑥, 𝑦) + 𝐶(𝑥, 𝑦)

𝜕2𝑢

𝜕𝑦2(𝑥, 𝑦) + 𝐷(𝑥, 𝑦)

𝜕𝑢

𝜕𝑥(𝑥, 𝑦)

+ 𝐸(𝑥, 𝑦)𝜕𝑢

𝜕𝑦(𝑥, 𝑦) + 𝐹(𝑥, 𝑦)𝑢(𝑥, 𝑦) + 𝐺(𝑥, 𝑦) = 𝑓(𝑥, 𝑦)

(5)

Esta EDP é homogênea se 𝑓(𝑥, 𝑦) = 0, e é não homogênea se 𝑓(𝑥, 𝑦) ≠ 0. Além disso,

pode-se classificar este tipo de equação em elíptica, parabólica ou hiperbólica, dependendo do

valor do discriminante 𝐵2 − 4𝐴𝐶. Para 𝐵2 − 4𝐴𝐶 > 0, a EDP é dita hiperbólica; para 𝐵2 −4𝐴𝐶 = 0 a EDP é considerada parabólica; para 𝐵2 − 4𝐴𝐶 < 0 a EDP é denominada elíptica

(CHAPRA, 2008).

Para se obter a solução desta equação é necessária a utilização das condições de

contorno, que são valores conhecidos da solução ou de sua derivada sobre os limites do domínio

da equação. Tais condições dependem do problema considerado, e serão abordadas na próxima

seção.

2.2. Equação de Poisson

Placas são elementos estruturais planos, ou seja, onde duas dimensões predominam

sobre às demais. O estudo da distribuição de temperatura em placas apresenta grande

relevância, uma vez que diversos elementos estruturais possuem esta forma, como as lajes ou

plataformas na construção civil, as quais muitas das vezes estão submetidas às grandes

variações de temperatura.

A distribuição de temperatura ao longo de uma placa fina pode ser calculada por um

caso particular de EDP elíptica, denominada equação de Poisson, a qual possui a forma geral

da Equação (6)

𝜕2𝑢

𝜕𝑥2 (𝑥, 𝑦) +

𝜕2𝑢

𝜕𝑦2 (𝑥, 𝑦) = 𝑓(𝑥, 𝑦), 𝑝𝑎𝑟𝑎 (𝑥, 𝑦) ∈ 𝑅

𝑢(𝑥, 𝑦) = 𝑔(𝑥, 𝑦), 𝑝𝑎𝑟𝑎 (𝑥, 𝑦) ∈ 𝑆

(6)

Nesta equação, 𝑅 = (𝑥, 𝑦) 𝑎⁄ < 𝑥 < 𝑏, 𝑐 < 𝑥 < 𝑑 e 𝑆 denota a fronteira de 𝑅. Caso

𝑓(𝑥, 𝑦) = 0, temos um caso particular da Equação de Poisson, denominada Equação de

Laplace. Para a obtenção da temperatura, é necessária a determinação das condições de

contorno, que são referentes aos valores de temperatura conhecidos ou funções conhecidas

que descrevam essa temperatura nas bordas ou em outros pontos da placa.

Neste caso propõe-se o cálculo da temperatura nos pontos do interior de uma placa

quadrada que possui a temperatura ao longo de suas extremidades conhecidas, sendo duas delas

constantes e as outras duas variáveis em função das dimensões do elemento. Para encontrar a

temperatura nos demais pontos da placa utiliza-se o MDF cuja veracidade será comprovada por

meio da comparação com a solução analítica dada pela Equação (7) (BURDEN, 2003).

𝑢(𝑥, 𝑦) = 𝑠𝑒𝑛(𝜋𝑥), 𝑝𝑎𝑟𝑎 (𝑥, 𝑦) ∈ 𝑆 (7)

2.3. Método de Diferenças Finitas

O Método das Diferenças Finitas (MDF) é um método numérico para cálculo de

problemas de valor de contorno. Muito utilizado, devido à sua simplicidade e facilidade de

implementação computacional. Utiliza a aproximação das derivadas de primeira e de segunda

ordem da função 𝑢(𝑥, 𝑦), solução da equação diferencial, dadas pelas respectivas equações de

diferenças divididas de primeira e de segunda ordem em 𝑥 e 𝑦 (SHIGUE, 2008).

Para a aplicação do método, primeiramente é necessário que o domínio esteja

discretizado. Discretização, neste contexto consiste na divisão do domínio em várias partes, de

modo a trabalhar não mais com o contínuo mas apenas com certos pontos.

Atentando para o fato deste trabalho modelar uma aplicação governada por uma EDP

cuja solução é uma função de duas variáveis, a forma mais simples de discretizar o domínio é

dividir os intervalos nos eixos 𝑥 e 𝑦 em partes iguais.

Por simplicidade, trabalha-se com um domínio quadrado, onde o lado que está orientado

na direção do eixo 𝑥 esteja limitado pelos valores 𝑎 e 𝑏 (𝑎 < 𝑏) e o outro lado, na direção do

eixo 𝑦, seja limitado por 𝑐 e 𝑑 (𝑐 < 𝑑). Então, para discretizar o domínio, deve-se selecionar

dois números inteiros, 𝑚 e 𝑛, de modo a dividir o intervalo do domínio no eixo 𝑥 em 𝑛

subintervalos e o intervalo no eixo 𝑦 em 𝑚 subintervalos (BURDEN, 2003).

O tamanho destes subintervalos, denominado por passo daqui em diante, recebe os

nomes ℎ no eixo 𝑥, e 𝑘 no eixo 𝑦. Dados os limites destes intervalos, pode-se definir as

expressões para a obtenção do valor dos passos, sendo ℎ = (𝑏 − 𝑎)/𝑛 e 𝑘 = (𝑑 − 𝑐)/𝑚 para

os eixos 𝑥 e 𝑦, respectivamente.

Agora, são estes pontos que serão utilizados para a construção das equações de

diferenças divididas, utilizando a série de Taylor em torno de cada ponto (𝑥𝑖 , 𝑦𝑗). Como o

problema apresenta apenas termos com derivadas segundas, torna-se necessário conhecer o

termo geral da expansão para funções de duas variáveis, mostrado na Equação (8).

𝑓(𝑥, 𝑦) = 𝑓(𝑥0 + 𝛿𝑥, 𝑦0 + 𝛿𝑦) =

= ∑ [1

𝑝![∑ [

𝑝!

𝑞! ∙ (𝑝 − 𝑞)!∙ (𝛿𝑥)𝑞 ∙ (𝛿𝑦)𝑝−𝑞 ∙

𝜕𝑝𝑓

𝜕𝑥𝑞𝜕𝑦𝑝−𝑞(𝑥0, 𝑦0)]

𝑝

𝑞=0

]]

∞

𝑝=0

(8)

Para a montagem da equação de diferenças divididas no eixo 𝑥, serão feitas duas

expansões, uma para 𝛿𝑥 = ℎ e a outra para 𝛿𝑥 = −ℎ. Sendo ambas com 𝛿𝑦 = 0. O processo

consiste em expandir até a terceira derivada, fazendo surgir um erro devido à interrupção da

série. Desde que os valores 𝑚 e 𝑛 escolhidos sejam suficientemente grandes, os passos ℎ e 𝑘

serão pequenos o bastante para que este erro seja pequeno. Ao fazer as duas expansões e somá-

las, obtém-se a expressão da Equação (9), onde o último termo representa o erro causado pela

interrupção da série.

𝜕2𝑢

𝜕𝑥2(𝑥𝑖 , 𝑦𝑗) =

𝑢 (𝑥𝑖+1, 𝑦𝑗) − 2𝑢 (𝑥𝑖, 𝑦𝑗) + 𝑢 (𝑥𝑖−1, 𝑦𝑗)

ℎ2−

ℎ2

12

𝜕4𝑢

𝜕𝑥4(𝜉𝑖, 𝑦𝑗) (9)

Este termo, a princípio desconhecido, pode ser desconsiderado se substituirmos

𝑢 (𝑥𝑖 , 𝑦𝑗) por 𝑤𝑖,𝑗, onde 𝑢 (𝑥𝑖 , 𝑦𝑗) ≈ 𝑤𝑖,𝑗, na soma das expansões. Quando realizada esta

substituição, a Equação (10) é obtida. Esta é a chamada fórmula das diferenças centrais de

segunda ordem, cuja função é aproximar o valor da derivada segunda em 𝑥.

𝜕2𝑢

𝜕𝑥2(𝑥𝑖 , 𝑦𝑗) ≈

𝑤𝑖+1,𝑗 − 2𝑤𝑖,𝑗 + 𝑤𝑖−1,𝑗

ℎ2 (10)

O processo ocorre de forma análoga para o eixo 𝑦, expandindo com os valores para

𝛿𝑦 = 𝑘 e 𝛿𝑦 = −𝑘, ambos com 𝛿𝑥 = 0. Após realizar a soma das duas expansões e a

substituição de 𝑢 (𝑥𝑖 , 𝑦𝑗) pela aproximação 𝑤𝑖,𝑗, tem-se a expressão dada na Equação (11), que

é a fórmula das diferenças centrais de segunda ordem, cuja função é aproximar o valor da

derivada segunda em 𝑦.

𝜕2𝑢

𝜕𝑦2(𝑥𝑖 , 𝑦𝑗) ≈

𝑤𝑖,𝑗+1 − 2𝑤𝑖,𝑗 + 𝑤𝑖,𝑗−1

𝑘2 (11)

Tendo em mãos as Equações (10) e (11), que são as aproximações para as duas derivadas

presentes na equação diferencial a ser solucionada, podemos substituir estas duas equações na

equação diferencial a ser resolvida. Atribuindo valores de i e j na expressão resultante, obtemos

um sistema de (𝑛 − 1)(𝑚 − 1) equações que, se resolvido, fornece os valores de 𝑤𝑖,𝑗 em todos

os pontos da quadrícula.

3. METODOLOGIA

Nesta seção, serão apresentados o problema a ser analisado e resolvido, seu método de

solução e os fatores referentes à implementação do algoritmo criado para solução de equações

como a que foi vista na Equação (6) via MDF.

3.1. Estudo de caso Uma placa metálica, quadrada, delgada (espessura desprezível) e com dimensões 0,5 m

por 0,5 m foi aquecida até que a temperatura em todos os seus pontos não variasse com o passar

do tempo. Sabe-se que duas de suas bordas são mantidas a 0º C e em seus outros dois limites o

valor da temperatura aumenta linearmente de 0º C, em um canto, para 100º C no lugar onde

ambos os lados se encontram. O objetivo é determinar a distribuição de temperatura deste

estado estável em cada ponto ao longo da extensão da placa.

Se retratarmos as bordas da placa com temperaturas conhecidas como as condições de

contorno ao longo dos eixos 𝑥 e 𝑦, obteremos um Problema de Valor de Contorno (PVC) dado

pela Equação (12)

𝜕2𝑢

𝜕𝑥2(𝑥, 𝑦) +

𝜕2𝑢

𝜕𝑦2(𝑥, 𝑦) = 0,

𝑢(0, 𝑦) = 0, 𝑢(𝑥, 0) = 0, 𝑢(𝑥, 0,5) = 200𝑥, 𝑢(0,5, 𝑦) = 200𝑦.

(12)

para (𝑥, 𝑦) no conjunto 𝑅 = (𝑥, 𝑦)|0 < 𝑥 < 0,5, 0 < 𝑦 < 0,5.

3.2. Solução da equação

Substituindo as aproximações das Equações (10) e (11) na EDP considerada na Equação

(12), obtém-se a Equação (13), que é o ponto de partida para a aplicar o MDF.

𝑤𝑖+1,𝑗 − 2𝑤𝑖,𝑗 + 𝑤𝑖−1,𝑗

ℎ2+

𝑤𝑖,𝑗+1 − 2𝑤𝑖,𝑗 + 𝑤𝑖,𝑗+1

𝑘2= 0 (13)

Multiplicando toda a Equação (9) por −ℎ2 e agrupando os termos semelhantes,

encontra-se a Equação (14), que é a chamada equação de diferenças para o problema descrito

2 [1 + (ℎ

𝑘)

2

] 𝑤𝑖,𝑗 − (𝑤𝑖+1,𝑗 + 𝑤𝑖−1,𝑗) − (ℎ

𝑘)

2

(𝑤𝑖,𝑗+1 + 𝑤𝑖,𝑗+1) = 0 (14)

Atribuindo valores para 𝑖 e 𝑗 monta-se o sistema de equações que fornece os valores de

𝑤𝑖,𝑗 necessários. Como, para este problema, ℎ = 𝑘, a Equação (14) assume uma forma mais

simples, dada na Equação (15).

4𝑤𝑖,𝑗 − 𝑤𝑖+1,𝑗 − 𝑤𝑖−1,𝑗 − 𝑤𝑖,𝑗+1 − 𝑤𝑖,𝑗+1 = 0 (15)

Após resolver o sistema, os valores de 𝑤𝑖,𝑗 são os valores aproximados da temperatura

nos pontos considerados. Devido ao pequeno valor de ℎ e 𝑘, os cálculos manuais não são

recomendados, dado o grande número de equações obtido. Por isso, foi elaborado um algoritmo

para a solução deste problema, sendo que seus detalhes são descritos na próxima seção.

4. IMPLEMENTAÇÃO DO ALGORITMO

Visando minimizar o encargo dos cálculos manuais, os quais seriam demorados e

repetitivos, a automatização dos mesmos foi realizada por meio de um algoritmo implementado

no software MATLAB®. Passou-se a ser necessária, para a obtenção da solução do problema,

apenas a inserção dos dados no programa.

Os dados fornecidos ao algoritmo foram advindos do estudo de caso e os passos de

discretização adotados foram ℎ = 𝑘 = 0,01 𝑚. O algoritmo, que aplica o Método Iterativo de

Gauss-Seidel para resolver o sistema linear produzido pelo uso do método das Diferenças

Finitas, encontra-se descrito abaixo, e os resultados obtidos serão mostrados e analisados na

seção seguinte.

function [] = Eliptica_Poisson()

syms x y

func = inputdlg('f(x, y)= ','Dados'); %edp

limites = inputdlg('a: ', 'b: ', 'c: ', 'd: ', 'm: ', 'n: ','Limites'); %a<x<b / c<y<d

contorno = inputdlg('u(a,y) = ', 'u(b,y) = ', 'u(x,c) = ', 'u(x,d) = ','PVC');

parada = inputdlg('Número de iterações: ', 'Tolerância: ');

a = str2num(limites1); %x0

b = str2num(limites2); %xn

c = str2num(limites3); %y0f

d = str2num(limites4); %ym

m = str2num(limites5);

n = str2num(limites6);

h = (b-a) / n; %passo em x

k = (d-c) / m; %passo em y

uF = matlabFunction(sym(func),'vars',[x y],'file','fxy');

uA = matlabFunction(sym(contorno1),'vars',[x y],'file','Uay');

uB = matlabFunction(sym(contorno2),'vars',[x y],'file','Uby');

uC = matlabFunction(sym(contorno3),'vars',[x y],'file','Uxc');

uD = matlabFunction(sym(contorno4),'vars',[x y],'file','Uxd');

w = zeros(n+1,m+1); %cria a matriz nula NxM

xVector = zeros(1, n+1); %discretiza x

i = 1;

for j = a : h : b

w(i,1) = uC(j,c);

w(i,m+1) = uD(j,d);

xVector(i) = j;

i = i+1;

end

i = 1;

yVector = zeros(1, m+1); %discretiza y

for j = c : k : d

w(1,i) = uA(a,j);

w(n+1,i) = uB(b,j);

yVector(i) = j;

i = i+1;

end

lam = (h^2) / (k^2);

u = 2*(1 + lam);

l = 1;

while l <= str2num(parada1) %iterações gauss-seidel

z = (((-h^2)*(uF(xVector(2),yVector(m))))+uA(a,yVector(m))+(lam*uD(xVector(2),d))+(lam*w(2,m-

1))+w(3,m))/u;

norm = abs(z - w(2,m));

w(2,m) = z;

for i = 3:n-1

z = (((-h^2)*(uF(xVector(i),yVector(m))))+(lam*uD(xVector(i),d))+w(i-1,m)+w(i+1,m)+(lam*w(i,m-

1)))/u;

if abs(w(i,m)-z) > norm

norm = abs(w(i,m)-z);

end

w(i,m) = z;

end

z = (((-h^2)*(uF(xVector(n),yVector(m))))+uB(b,yVector(m))+(lam*uD(xVector(n),d))+w(n-

1,m)+(lam*w(n,m-1)))/u;

if abs(w(n,m)-z) > norm

norm = abs(w(n,m)-z);

end

w(n,m) = z;

for j = m-1: -1 : 3

z = (((-h^2)*(uF(xVector(2),yVector(j))))+uA(a,yVector(j))+(lam*w(2,j+1))+(lam*w(2,j-1))+w(3,j))/u;

if abs(w(2,j)-z) > norm

norm = abs(w(2,j)-z);

end

w(2,j) = z;

for i = 3:n-1

z = (((-h^2)*(uF(xVector(i),yVector(j))))+w(i-1,j)+(lam*w(i,j+1))+w(i+1,j)+(lam*w(i,j-1)))/u;

if abs(w(i,j)-z) > norm

norm = abs(w(i,j)-z);

end

w(i,j) = z;

end

z = (((-h^2)*(uF(xVector(n),yVector(j))))+uB(b,yVector(j))+w(n-1,j)+(lam*w(n,j+1))+(lam*w(n,j-1)))/u;

if abs(w(n,j)-z) > norm

norm = abs(w(n,j)-z);

end

w(n,j) = z;

end

z = (((-

h^2)*(uF(xVector(2),yVector(2))))+uA(a,yVector(2))+(lam*uC(xVector(2),c))+(lam*w(2,3))+w(3,2))/u;

if abs(w(2,2)-z) > norm

norm = abs(w(2,2)-z);

end

w(2,2) = z;

for i = 3:n-1

z = (((-h^2)*(uF(xVector(i),yVector(2))))+(lam*uC(xVector(i),c))+w(i-1,2)+(lam*w(i,3))+w(i+1,2))/u;

if abs(w(i,2)-z) > norm

norm = abs(w(i,2)-z);

end

w(i,2) = z;

end

z = (((-h^2)*(uF(xVector(n),yVector(2))))+uB(b,yVector(2))+(lam*uC(xVector(n),c))+w(n-

1,2)+(lam*w(n,3)))/u;

if abs(w(n,2)-z) > norm

norm = abs(w(n,2)-z);

end

w(n,2) = z;

format long;

if (norm <= str2num(parada2))

for i = 1:n-1

for j = 1:m-1

h = plot3(xVector(i),yVector(j),w(i,j),'r.');

hold on;

end

end

xlabel('x');

ylabel('y');

zlabel('w');

title('Equação Elíptica');

h = surf(xVector,yVector,w);

saveas(h,'Grafico','fig');

break;

end

l = l +1;

end

if l > str2num(parada1)

warndlg('Numero de iterações excedido. Abortar.');

else

file = fopen('resultados.txt','wt');

solAnalitica = inputdlg('u(x, y) = ','Solução Analítica', 1, '0');

funcAnalitica = sym(solAnalitica);

if ~(strcmp(char(funcAnalitica),'0'))

funAnalitica = matlabFunction(funcAnalitica,'vars',[x y],'file','fun_analitica');

fprintf(file,'i \t j \t xi \t\t yj \t\t wij \t\t\t u(i,j) \t\t Erro absoluto \t\t Erro percentual\n');

for i = 1 : n+1

for j = 1 : m+1

u(i,j)= funAnalitica(xVector(i),yVector(j));

fprintf(file,'%d \t %d \t %.4f \t %.4f \t %.10f \t %.10f \t %.10e \t %.10e\n',i-1, j-1, xVector(i),

yVector(j), w(i,j), u(i,j), abs(u(i,j) - w(i,j)), 100*abs(u(i,j) - w(i,j))/(u(i,j)));

end

end

else

fprintf(file,'i \t j \t xi \t\t yj \t\t wij\n');

for i = 1 : n+1

for j = 1 : m+1

fprintf(file,'%d \t %d \t %.4f \t %.4f \t %.10f\n', i-1, j-1, xVector(i), yVector(j), w(i,j));

end

end

end

end

fclose(file);

5. RESULTADOS E DISCUSSÃO

Após a implementação do estudo de caso no algoritmo, obtêm-se os valores de 𝑤𝑖,𝑗 nos

pontos discretizados. Como são muitos pontos, apenas alguns dos valores obtidos são

mostrados na Tabela 1, a fim de permitir a interpretação dos valores aproximados bem como

da ordem de grandeza dos erros obtidos.

Tabela 1 – Valores Analíticos e Numéricos para alguns dos pontos discretizados.

𝑖 𝑗 𝑥𝑖 𝑦𝑗 𝑤𝑖,𝑗 𝑢(𝑥𝑖 , 𝑦𝑗) Erro absoluto Erro percentual

25 0 0.2500 0.0000 0.0000000000 0.0000000000 0.0000000000e+00 0.000000000e+00

25 5 0.2500 0.0500 4.9999999925 5.0000000000 7.4975421427e-09 1.4995084285e-07

25 10 0.2500 0.1000 9.9999999856 10.0000000000 1.4402715465e-08 1.4402715465e-07

25 15 0.2500 0.1500 14.9999999800 15.0000000000 2.0020392455e-08 1.3346928303e-07

25 20 0.2500 0.2000 19.9999999762 20.0000000000 2.3768969015e-08 1.1884484508e-07

25 25 0.2500 0.2500 24.9999999748 25.0000000000 2.5240218804e-08 1.0096087522e-07

25 30 0.2500 0.3000 29.9999999758 30.0000000000 2.4243128394e-08 8.0810427979e-08

25 35 0.2500 0.3500 34.9999999792 35.0000000000 2.0827116032e-08 5.9506045805e-08

25 40 0.2500 0.4000 39.9999999847 40.0000000000 1.5281969468e-08 3.8204923669e-08

25 45 0.2500 0.4500 44.9999999919 45.0000000000 8.1139432950e-09 1.8030985100e-08

25 50 0.2500 0.5000 50.0000000000 50.0000000000 0.0000000000e+00 0.0000000000e+00

Além dos valores de 𝑤𝑖,𝑗, o algoritmo fornece um gráfico, exibido na Figura 1.

Figura 1 – Plotagem da Solução Numérica

Nota-se que os valores da solução numérica apresentam grande semelhança com a

solução analítica. Além disso os erros obtidos foram pequenos (na ordem de 10−8) e, portanto,

podem ser considerados satisfatórios.

6. CONCLUSÕES

Mediante os resultados apresentados, é possível concluir que o Método das Diferenças

Finitas pode se tornar uma preciosa ferramenta na busca de soluções de problemas que

envolvam EDP’s. A facilidade de assimilação, de aplicação do método bem como a eficiência

de sua aproximação comprovam sua importância na resolução não só desta, mas de inúmeras

aplicações a serem exploradas nas áreas da física e engenharias, como o comportamento dos

fluidos, distribuição de um campo elétrico em um circuito, entre outras.

Outra ferramenta utilizada no trabalho e de grande importância foi o software

matemático MATLAB®. A facilidade em manipulá-lo e o algoritmo implementado permitiram

uma automação eficiente dos cálculos o que representou uma grande contribuição para a

solução do problema apresentado.

Agradecimentos

Agradecemos à UTFPR pela concessão de bolsas de apoio à produção científica e também

pela disponibilização do espaço físico, pelo incentivo aos alunos-autores e pela manutenção,

durante a graduação, do projeto de pesquisa que originou este trabalho.

7. REFERÊNCIAS

Livros: BURDEN, Richard L.; FAIRES, John Douglas. Análise Numérica. São Paulo: Pioneira

Thomson Learning, 2003.

CHAPRA, Steven C.; CANALE, Raymond P. Métodos Numéricos para Engenharia. 5. ed. São

Paulo: McGraw-Hill, 2008.

MEDEIROS, Luis Adauto; FERREL, Juan Limaco; BIAZUTTI, Angela Cássia. Métodos

Clássicos em Equações Diferenciais Parciais. 2. ed. Rio de Janeiro: UFRJ-IM, 1999. 163p.

Internet: SHIGUE, Carlos Y. Equações Diferenciais Parciais. 2008. Disponível em:

<http://www.demar.eel.usp.br/metodos/mat_didatico/Equacoes_diferenciais_parciais.pdf>.

Acesso em: 14 jun. 2014.

CALCULATION OF TEMPERATURE WITHIN A PLATE WITH THE

AID OF THE METHOD OF FINITE DIFFERENCES

Abstract: This paper addresses the analysis of the efficiency of the Finite Difference Method

when applied to a common problem of physics and engineering: the temperature distribution

along a thin, square plate. The paper contains a description of the numerical method, including

its deduction and way of application, highlighting its praticality and ease of use, aiming to

incentivate its use and dissemination. Futhermore, there is an introduction to the Poisson’s

equation, whose solution represents the value of the temperature at each point in the interior

of the plate under certain conditions. An algorithm was implemented and applied in the

resolution of a case study.

Keywords: Distribution of Temperature, Finite Difference Method, Poisson equation.