Text Material Preview
INTRODUÇÃO
AOS MÉTODOS
NUMÉRICOS
Celme Torres F. da Costa
UNIVERSIDADE FEDERAL DO CEARÁ
Campus Cariri
2010
Introdução aos Métodos Numéricos
2
Sumário
1. Introdução aos Métodos Numéricos ................................................................... 4
1.1. Soluções Analíticas .................................................................................... 5
1.2. Método das diferenças Finitas ................................................................... 6
1.3. Método dos elementos finitos .................................................................... 6
1.4. Método dos momentos (MOM) .................................................................. 7
1.5. Método dos elementos analíticos ............................................................... 7
1.6. Método das camadas finitas ...................................................................... 7
1.7. Método de elementos de contorno (MEC) ................................................. 8
1.8. Método dos volumes finitos........................................................................ 8
2. Condições de Contorno .................................................................................... 10
2.1. Condições de contorno naturais .............................................................. 11
2.2. Condições de contorno geométricas ........................................................ 11
3. Sistemas de Equações Diferenciais ................................................................. 12
3.1. Sistemas de Equações Diferenciais ......................................................... 12
3.1.1. Sistemas homogêneos ...................................................................... 12
3.1.2. Sistemas não homogêneos ............................................................... 14
3.2. Conceito de autovalor e autovetor ........................................................... 15
3.3. Diagonalização de matrizes ..................................................................... 18
3.4. Solução geral de sistemas homogêneos ................................................. 19
3.5. Solução geral de sistemas não homogêneos .......................................... 20
3.5.1. Exercícios e Aplicações de sistemas de equações diferenciais ........ 21
3.6. Solução de sistemas lineares .................................................................. 22
3.6.1. Exercícios .......................................................................................... 27
4. Métodos dos Resíduos Ponderados ................................................................ 28
4.1. Método de colocação pontual .................................................................. 29
4.2. Método dos subdomínios ......................................................................... 32
4.3. Método dos mínimos quadrados .............................................................. 34
4.4. Método de Galerkin .................................................................................. 36
4.5. Exercícios ................................................................................................ 38
5. Introdução aos Métodos Variacionais............................................................... 39
5.1. Cálculo variacional ................................................................................... 39
5.2. Máximo e mínimo de funções .................................................................. 39
5.3. Funcional ................................................................................................. 43
5.3.1. Funcional de uma função .................................................................. 44
5.3.2. Determinação do Funcional - Equação de Euler ............................... 47
5.3.3. Exercícios .......................................................................................... 51
6. Métodos Variacionais ....................................................................................... 53
6.1. Método de Ritz ......................................................................................... 53
6.1.1. Exercícios .......................................................................................... 58
6.2. Método de Rayleigh-Ritz .......................................................................... 59
7. Métodos das Diferenças Finitas ....................................................................... 66
7.1. Derivadas de ordem superior ................................................................... 68
7.2. Exercícios ................................................................................................ 73
8. Métodos dos Elementos Finitos ....................................................................... 75
Profª Celme Torres
3
8.1. Sistemas discretos e funções de base ..................................................... 78
8.1.1. Funções de base ................................................................................ 78
8.1.2. Exemplos de sistemas discretos ........................................................ 81
8.1.2.1. Mola elástica ................................................................................ 81
8.1.2.2. Torção em barra circular .............................................................. 82
8.1.2.3. Fluxo de um fluido em um tubo .................................................... 83
9. Métodos dos Elementos Finitos Bidimensional ................................................. 94
9.1. Modelos de elementos finitos ................................................................... 95
9.1.1. Funções de forma .............................................................................. 95
9.1.2. Elemento triangular linear................................................................... 96
9.1.3. Teorema de Green-Gaus ................................................................... 98
9.1.4. Imposição das condições de contorno ............................................. 100
9.1.5. Determinação da matriz de rigidez ................................................... 101
9.1.6. Determinação do vetor de carga nodal............................................. 102
Introdução aos Métodos Numéricos
4
1. Introdução aos Métodos Numéricos
Os modelos são ferramentas projetadas para representar uma versão simplificada de um
problema real. É uma tentativa para compreensão dos processos físicos, químicos e bioló-
gicos traduzidos em termos matemáticos. O objetivo da modelagem é o de prever ou predi-
zer cenários onde estão envolvidos variáveis desconhecidas, como, por exemplo, a varia-
ção da carga hidráulica ou distribuição de concentrações de espécies químicas em um sis-
tema aqüífero no tempo e no espaço (BEDIENT et al., 1994), variação de temperaturas,
análise de estruturas, entre outras.
No desenvolvimento de um modelo, o primeiro passo é o desenvolvimento de um modelo
teórico ou conceitual, que consiste na descrição dos processos envolvidos e no conheci-
mento do comportamento do sistema analisado. O próximo passo é traduzir esse modelo
teórico em termos matemáticos, o chamado modelo matemático, que é o conjunto de equa-
ções diferenciais associada ao conjunto de equações de condições iniciais e de condições
de contorno. Finalmente, a solução das equações sujeitas às condições de contorno impos-
tas é obtida utilizando métodos analíticos ou numéricos.
Quando os modelos consideram as variáveis de entrada e/ou os parâmetros dependentes de
funções aleatórias (probabilidade) esses modelos são definidos como modelos ESTO-
CÁSTICOS ou PROBABILÍSTICOS, do contrário, são DETERMINÍSTICOS. Os mo-
delos baseados em resultados experimentais, são ditos modelos EMPÍRICOS (COSTA,
2000).
Uma vez que o modelo teórico é formulado e os parâmetros apropriados tenham sido de-
terminados,o próximo passo é achar uma solução para as equações que governam o pro-
blema, considerando as condições iniciais e de contorno. As técnicas mais utilizadas na
solução destas equações são: solução analítica, método dos momentos, método das cama-
das finitas, elementos de contorno, diferenças finitas e elementos finitos, descritas a seguir.
Profª Celme Torres
5
1.1.1. Soluções Analíticas
Qualquer função f definida em algum intervalo I, que, quando substituída na equação dife-
rencial reduz a equação a uma identidade, é chamada de solução analítica da equação di-
ferencial no intervalo. Porém, nem toda equação diferencial possui necessariamente uma
solução analítica.
Quando resolvemos uma equação diferencial, normalmente obtemos uma família de curvas
ou funções contendo um parâmetro arbitrário tal que, cada membro da família é uma solu-
ção da equação diferencial.
As soluções das equações diferenciais são divididas em Explícitas ou Implícitas. Uma so-
lução explícita é qualquer função y(x) que verifique a equação num intervalo a x b .
Ou seja, uma solução na qual a variável dependente é expressa somente em termos da vari-
ável independente e das constantes. Para nossos propósitos, vamos pensar que uma solução
explicita seja da forma y f (x) , a qual pode ser manipulada, calculada e diferenciável por
meio das regras padrão. Isto é, a variável dependente pode ser isolada e igualada a uma ex-
pressão, a qual é função apenas da variável independente, não ambígua.
Uma solução implícita é uma relação g(x, y) 0 que verifique a equação. Dizemos que
uma relação g(x, y) 0 é uma solução implícita de uma EDO em um intervalo I, se ela
define uma ou mais soluções explícitas em I. Isto é, a variável dependente (função) não
pode ser isolada e igualada a uma expressão que dependa apenas da variável independente,
ou quando isto for possível então a expressão será ambígua.
As soluções analíticas são ideais para um rápido cálculo preliminar, são válidas para estu-
dos de sensibilidade e úteis para checar resultados de análises mais complexas. Toda base
teórica para determinação das soluções analíticas de equações diferenciais foram descritas
e assimiladas em disciplinas anteriores envolvendo equações diferencias ou matemática
aplicada.
Introdução aos Métodos Numéricos
6
1.1.2. Método das diferenças Finitas
O método das diferenças finitas substitui a equação diferencial por um conjunto de equa-
ções lineares algébricas que são resolvidas de forma aproximada (McDONALD; HAR-
BAUGH, 1988). As idéias fundamentais foram estabelecidas e usadas por matemáticos do
século XVIII, tais como Taylor e Lagrange (GOMES, 2000). No método das Diferenças
Finitas, uma equação diferencial (de natureza contínua) é substituída por uma série de e-
quações algébricas, chamadas de diferenças finitas, em pontos discretos. Esse método é o
mais antigo e o mais divulgado e foi o primeiro a ser usado para a solução sistemática de
problemas de água subterrânea.
Figura 1. Malha de diferenças finitas.
1.1.3. Método dos elementos finitos
O método dos elementos finitos substitui a equação diferencial por uma formulação vari-
acional do problema, que leva a resolução aproximada de um sistema de equações lineares
algébricas (ZIENKIEWICZ, 1977). Baseia-se na divisão do domínio do problema em sub-
domínios, o que permite a representação de domínios complexos como um conjunto de
subdomínios mais simples – chamados Elementos Finitos. Cada elemento finito é conecta-
do aos elementos vizinhos através de nós. O conjunto dos elementos e dos nós é chamado
de malha. A Figura 1.0 mostra uma malha de elementos finitos.
i, ji-1, j i+1, j
i, j+1
i, j-1i-1, j-1 i+1, j-1
i-1, j+1 i+1, j+1 i+2, j+2
i+2, j-2
x
y
x
y
Profª Celme Torres
7
Figura 2 – Malha de elementos finitos.
1.1.4. Método dos momentos (MOM)
O método dos momentos (MOM) tem a sua origem – como muitas outras idéias em esta-
tística - no trabalho precursor de Lord Karl PEASON (1857). Pearson havia construído o
seu sistema de famílias de distribuições de probabilidade e necessitava um procedimento
automático para, uma vez escolhida a família apropriada, estimar os parâmetros da particu-
lar curva (ou função de densidade) que melhor se ajustaria aos dados. Como tais parâme-
tros eram, geralmente, relacionados aos momentos da distribuição através de funções algé-
bricas não muito complicadas, ele propôs que se igualasse os valores dos momentos amos-
trais aos momentos teóricos, conforme descritos por tais funções. Desse modo, indo-se até
o momento de ordem igual ao número de parâmetros a estimar, se obtinha um sistema de
tantas equações quanto incógnitas – os parâmetros – a resolver.
1.1.5. Método dos elementos analíticos
O método dos elementos analíticos é baseado no princípio da superposição de várias fun-
ções analíticas. Este princípio consiste da adição de soluções individuais para cada elemen-
to. A solução do problema é obtida através da adição de todas as influências dos elementos
analíticos individuais, representados pelas funções analíticas, que correspondem às caracte-
rísticas do sistema em estudo.
1.1.6. Método das camadas finitas
O método das camadas finitas, descrito por ROWE e BOOKER (1985 a, b; 1987), se a-
plica a situações onde a hidroestratigrafia pode ser idealizada como composta de camadas
Introdução aos Métodos Numéricos
8
horizontais, com as propriedades do solo sendo as mesmas em qualquer local na camada. O
método das camadas finitas é muito utilizado em problemas de dimensionamento de pavi-
mentos.
1.1.7. Método de elementos de contorno (MEC)
O Método dos elementos de contorno ou Boundary Element Method – BEM, é um méto-
do computacional para a solução de sistemas de equações diferenciais, formuladas em
forma integral. É aplicado em diversas áreas da engenharia, como em mecânica dos flui-
dos, acústica, eletromagnetismo e estudo de fraturas.
O método possui melhor desempenho que o método dos elementos finitos em certas cir-
cunstâncias, como por exemplo, quando o domínio de estudo for infinito ou semi-infinito.
No método de elementos de contorno (MEC), primeiramente ocorre a transformação da
equação diferencial parcial que descreve o comportamento da incógnita no interior e no
contorno do domínio, em uma equação integral que envolve somente incógnitas em seu
contorno. A seguir, faz-se a discretização do contorno em elementos de superfície e, por
fim, encontra-se a solução do sistema de equações algébricas resultantes. As principais
vantagens do MEC estão no fato de que somente o contorno (ou contornos) do domínio
deve ser dividido em sub-regiões (discretizado). Nos outros métodos (Elementos Finitos e
Diferenças Finitas), todo o domínio deve ser discretizado. Desta forma, a dimensão do
problema é efetivamente reduzida em uma dimensão (BREBBIA e SKERGET, 1984).
1.1.8. Método dos volumes finitos
O método dos volumes finitos é um método de resolução de equações às derivadas parci-
ais baseado na resolução de balanços de massa, energia e quantidade de movimento a um
determinado volume de meio contínuo.
Este método evoluiu das diferenças finitas e não apresenta problemas de instabilidade ou
convergência, por garantir que, em cada volume discretizado, a propriedade em questão
(por exemplo, a massa) obedece à lei da conservação de massa. É um método largamente
Profª Celme Torres
9
utilizado na resolução de problemas envolvendo transferência de calor ou massa e em me-
cânica dos fluidos.
Figura 1.1 - Diagrama de métodos de solução de um problema físico.
Introdução aos Métodos Numéricos
10
2. Condições de Contorno
O objetivo da maioria das análises numéricas é determinar funções desconhecidas,chama-
das variáveis dependentes, que satisfaçam um dado conjunto de equações diferenciais em
um domínio ou região e algumas condições de contorno na borda do domínio. Um domínio
é uma coleção de pontos no espaço com propriedades que se P é um ponto no domínio, en-
tão todos os pontos próximos a P pertencem ao mesmo domínio. Esta definição implica
que um domínio consiste somente de pontos internos. Os símbolos Ω são utilizados para
denotar um domínio arbitrário e Γ para determinar seu contorno.
Γ (contorno)
Ω (domínio)
Figura 1.2 – Esquema de apresentação do domínio e seu contorno.
Em matemática, no ramo de equações diferenciais, um problema de valor de contorno ou
de fronteira é uma equação diferencial provida de um conjunto de restrições adicionais,
as chamadas condições de contorno. Uma solução para um problema de valor de contorno
é aquela que seja solução da equação diferencial e satisfaça as condições de contorno.
Problemas de valor de contorno surgem em diversos ramos, assim como toda equação dife-
rencial também o terá. Problemas envolvendo a equação de onda, bem como a determina-
ção dos modos normais, são frequentementes classificados como problemas de valor de
contorno. Um problema de valor de contorno deve ser bem determinado. Isto é, dado certas
condições para o problema, haverá então solução única, que depende continuamente das
condições citadas.
A solução de qualquer equação diferencial dependente do tempo (transiente) requer a espe-
cificação de condições na fronteira do sistema estudado (condições de contorno) e as con-
Profª Celme Torres
11
dições no início do processo físico (condição inicial). A definição precisa das condições de
contorno e inicial é a parte mais importante para os processos de modelagem.
Quando as variáveis dependentes são funções de uma variável independente (x), o domínio
é um segmento de reta, dito, nesse caso, unidimensional, e os pontos limites são chamados
pontos de contorno. Quando as variáveis dependentes são funções de duas variáveis inde-
pendentes (x, y), o domínio é bidimensional, sendo dado por uma superfície ou área, ge-
ralmente um plano.
2.1.1. Condições de contorno naturais
As condições de contorno naturais se referem à derivada da função que está sendo analisa-
da.
No caso de escoamento corresponde aos valores de fluxo no contorno do domínio, assim
como os pontos de perda ou ganho (poços, vazões, carga hidráulica). No caso de análise de
estruturas as condições de contorno naturais representam as forças concentradas e distribu-
ídas existentes no sistema físico analisado.
2.1.2. Condições de contorno geométricas
As condições de contorno geométricas são os valores da função ao longo do contorno ou
domínio. No caso de análise de estruturas correspondem aos apoios e as restrições do sis-
tema.
Introdução aos Métodos Numéricos
12
3. Sistemas de Equações Diferenciais
3.1.1. Sistemas de Equações Diferenciais
As equações diferenciais ordinárias simultâneas envolvem duas ou mais equações que con-
tém derivadas de duas ou mais funções incógnitas de uma única variável independente. Se
x, y e z são funções da variável t, então:
2
2
2
2
d x
4 5x y
dt
d y
2 3x y
dt
e,
2
x ' 3x y ' z ' 5
x ' y ' 2z ' t
x y ' 6z ' t 1
São exemplos de sistemas de equações diferenciais.
Soluções de um Sistema
Uma solução de um sistema de equações diferenciais é um conjunto de funções diferenciá-
veis x f (t) , y g(t) , z h(t) , etc., que satisfaz cada equação do sistema em algum in-
tervalo I.
Os sistemas de equações diferenciais podem ser resolvidos por:
(i) Método dos operadores, através da eliminação algébrica sistemática;
(ii) Método da transformada de Laplace
(iii) Método de autovetores e autovalores de matrizes
Sistemas homogêneos
Um sistema de equações diferenciais é um conjunto de n ED, com uma variável indepen-
dente e n variáveis dependentes, que podem ser escritos da seguinte forma:
' ' '1
1 1 2 n 1 2 n
' ' '2
2 1 2 n 1 2 n
' ' 'n
n 1 2 n 1 2 n
dy
F y , y ,..., y , y , y ,..., y , x
dx
dy
F y , y ,..., y , y , y ,..., y , x
dx
dy
F y , y ,..., y , y , y ,..., y , x
dx
Profª Celme Torres
13
onde F1, F2, ..., Fn são quaisquer funções de (2n+1) variáveis reais, que definem o sistema.
O sistema de EDO lineares pode ser escrito numa forma simplificada, usando a notação
vetorial,
onde y é um vetor com n componentes, cada uma delas função de x, e F é um vetor com n
componentes, funções de y, x e y’.
Um caso especial na teoria de ED é quando o vetor de funções F, no sistema de equações
tem a forma:
onde A é uma matriz quadrada n x n de funções que dependem unicamente de x, f é um
vetor
com n componentes dependentes de x.
Exemplo 3.1
1. Escreva em forma matricial os sistemas
(a)
1 1 2
2 1 2
y ' x x
y ' 4x x
Na forma vetorial, o sistema pode ser escrito como:
O vetor f é nulo e a matriz A é igual a
1 1
A
4 1
dy dy
F y, x,
dx dt
F Ax f
1 1
2 2
y ' x1 1
y ' 4 1 x
Princípio da Superposição
Seja 1 2 nx , x ,..., x um conjunto de vetores solução de um sistema de EDO homogêneo
em um intervalo I. Então a combinação linear
1 1 2 2 n nx c x c x ... c x
Onde, c1, c2 e cn são constantes arbitrárias, é também solução do sistema no intervalo.
Introdução aos Métodos Numéricos
14
(b)
t
t
dx
3x 4y e sin 2t
dt
dy
5x 9y 4e cos 2t
dt
A forma matricial temos:
t
t
x ' 3 4 x e sin 2t
y ' 5 9 y 4e cos 2t
t
t
3 4 e sin 2tdX
X
5 9dt 4e cos 2t
Onde,
x
X
y
2. Verifique que
2t
1 2t
e
x
e
e
6t
2 6t
3e
x
5e
são soluções do sistema de EDO’s
1 3
X ' X
5 3
, no intervalo ( , ) , onde
x
X
y
.
Temos que,
2t
1 2t
2e
x '
2e
2t 2t 2t 2t
1 12t 2t 2t 2t
1 3 2e e 3e 2e
Ax x '
5 3 2e 5e 3e 2e
6t
2 6t
18e
x '
30e
6t 6t 6t 6t
2 26t 6t 6t 6t
1 3 3e 3e 15e 18e
Ax x '
5 3 5e 15e 15e 30e
Logo, x1 e x2 são soluções da EDO.
Sistemas não homogêneos
Para sistemas não homogêneos, uma solução particular xp em um intervalo I é qualquer
vetor, sem parâmetros arbitrários, cujos elementos são funções que satisfazem o sistema,
dX
A(t)X F(t)
dt
Profª Celme Torres
15
Exemplo 3.2
Verifique que o vetor
p
3t 4
X
5t 6
é uma solução particular do sistema não-homogêneo
1 3 12t 11
X ' X
5 3 3
, no internalo ( , ) .
Solução
Temos que,
p
3
X '
5
p
1 3 12t 11 1 3 3t 4 12t 11
X
5 3 3 5 3 5t 6 3
12t 14 12t 11 3
2 3 5
Logo, o vetor
p
3t 4
X
5t 6
é solução particular do sistema homogêneo.
3.1.2. Conceito de autovalor e autovetor
Dado o sistema:
8 2 1 14
2 11 3 35
Teorema
Seja 1 2 nx , x ,..., x um conjunto de vetores solução de um sistema de EDO homogêneo
em um intervalo I e seja Xp um vetor arbitráriosolução do sistema não homogêneo no
mesmo intervalo. Então a combinação linear
1 1 2 2 n n px c x c x ... c x X
É também uma solução do sistema não-homogêneo no intervalo, para quaisquer valores
das constantes c1, c2 e cn .
vetor de entrada
vetor de saída
Introdução aos Métodos Numéricos
16
Tendo,
O vetor de saída está na mesma direção e sentido que o vetor de entrada. Neste caso:
Sendo um escalar e v um vetor.
Logo, v é chamado de AUTOVETOR de A e um AUTOVALOR de A associado ao ve-
tor v.
Considere a matriz quadrada . Se for um número complexo e v um vetor com-
plexo não nulo de dimensões n, satisfazendo a identidade:
8 2 1 12 1
12
2 11 2 24 2
Av v
n nA R
Av v
vetor de saída vetor de entrada
vetor de saída
vetor de entrada
Teorema de Cayley-Hamilton
Qualquer matriz quadrada satisfaz sua própria equação característica. Isto é, se e somen-
te se,
n n 1 2
n n 1 2 1 odet(A I) b b ... b b b
Então, n n 1 2n n 1 2 1 ob A b A ... b A b A b I 0
Profª Celme Torres
17
Então é um autovalor de A e v é um autovetor de A associado ao autovalor de .
Se A é uma matriz n x n sobre R e I é a matriz identidade de mesma ordem de A, defini-
mos o polinômio característico de A como:
f ( ) det( I A)
Os autovalores de A podem ser obtidos a partir da existência de escalares e vetores não-
nulos v para os quais Av v .
Este sistema pode ser escrito como:
Av Iv
(A I)v 0
Sendo I é a matriz identidade. Para que a solução não seja nula, o determinante deve ser
igual a zero (equação característica), logo:
det( I A) 0
OBS: O número de autovalores é sempre igual a ordem da matriz.
Exemplo 3.3
(1) Calcule o autovalor da EDO: 2y'' 6y ' 4y 0
Pelo teorema de Cayley-Hamilton
2
1
2
2 6 4 0
logo,
2
1
(2) Calcule os autovalores e autovetores da matriz
3 4
A
1 7
Pelo teorema de Cayley-Hamilton
2 2
1 2
3 4
det(A I) ( 5) 0
1 7
5
Para encontrar os autovetores correspondentes ao autovalor, fazemos:
Introdução aos Métodos Numéricos
18
1
2
1 2
1 2
(A I)k 0
(A 5I)k 0
k3 4 1 0
5 0
1 7 0 1 k
2k 4k 0
k 2k 0
Por esse sistema temos que 1 2k 2k . Escolhendo a solução trivial 2k 1 , temos o úni-
co autovetor do sistema, dada por:
1
2
K
1
3.1.3. Diagonalização de matrizes
Seja C uma matriz qualquer
a b
C
c d
Se multiplicarmos C por
1
0
a b 1 a
C
c d 0 c
1ª coluna de C
Multiplicando C por
0
1
a b 0 b
C
c d 1 d
2ª coluna de C
Suponha uma matriz A2 x 2 cujos autovetores e autovalores são:
1 1
2 2
v
v
Se uma matriz E 2 x 2 é formada pelos autovetores de A e se,
1
a
v
c
e 2
b
v
d
Profª Celme Torres
19
a b
E
c d
Considere a matriz resultante de 1E AE , multiplicada por
1
0
1
1
a
v
c
1
E A E
0
.
Temos, 1 1E Av
Onde,
1 1 1Av v
1 11 1 1 1E v E v
Se,
1
1
E v
0
1
1
1
E v
0
Assim,
11
1 1E v
0
1ª coluna de 1 1E Av
Do mesmo modo se multiplicarmos
1E AE por
0
1
, obtemos a 2ª coluna de 1 1E Av
Portanto,
11
2
0
E AE
0
, diagonaliza a matriz A.
3.1.4. Solução geral de sistemas homogêneos
Seja 1 2 n, ,..., n autovalores reais distintos de uma matriz de coeficientes A do sistema
X' AX , onde
x
X
y
, e sejam 1 2 nk ,k ,...,k os autovetores correspondentes. Então a
solução geral do sistema no intervalo ( , ) é dada por:
1 2 nt t t
1 1 2 2 n nX c k e c k e ... c k e
Exemplo 3.4
Resolva
dx
2x 3y
dt
dy
2x y
dt
A equação característica é
Introdução aos Métodos Numéricos
20
22 3det(A I) 3 4 0
2 1
Ao autovalores do sistema são: 1 1 e 2 4
Para 1 1 , temos para (A I)K 0
1 2
1 2
3k 3k 0
2k 2k 0
Resolvendo o sistema: 1 2k k
Quando 2k 1 , o autovetor correspondente é 1
1
k
1
Para 2 4 , temos para (A I)K 0
1 2
1 2
2k 3k 0
2k 3k 0
Resolvendo o sistema: 21
3k
k
2
Quando 2k 2 , o autovetor correspondente é 2
3
k
2
Como a matriz A de coeficientes é uma matriz 2 x 2, achamos duas soluções linearmente
independentes na forma:
t1
1
X e
1
e 4t2
3
X e
2
A solução geral do sistema de EDO homogêneo é dada por:
1 2t t t 4t1 1 2 2 1 2
1 3
X c k e c k e c e c e
1 2
t 4t
1 2
t 4t
1 2
x(t) c e 3c e
y(t) c e 2c e
3.1.5. Solução geral de sistemas não homogêneos
Seja Xp uma solução do sistema não-homogêneo X' A(t)X F(t) em um intervalo I, e
denotemos por
p 1 1 2 2 n nX c X c X ... c X
A solução geral no mesmo intervalo, do sistema homogêneo X' A(t)X correspondente.
Define-se solução geral do sistema não-homogêneo no intervalo como c pX X X .
Profª Celme Torres
21
Exemplo 3.5
Verificamos que o vetor
p
3t 4
X
5t 6
é uma solução particular do sistema não- homo-
gêneo
1 3 12t 11
X ' X
5 3 3
, no internalo ( , ) .
Para a parte homogênea da EDO temos
dx
1x 3y
dt
dy
5x 3y
dt
Resolvendo aplicando autovetor e autovalores de matriz, temos:
2t
1
1
X e
1
e 6t
2
3
X e
5
A solução não homogênea Xn é dada por
2t 6t
n
1 3
X e e
1 5
A solução geral do sistema de EDO’s pela definição é n pX X X , logo:
2t 6t
1 3 3t 4
X e e
1 5 5t 6
2t 6t
2t 6t
x(t) e 3e 3t 4
y(t) e 5e 5t 6
Exercícios e Aplicações de sistemas de equações diferenciais
(1) Encontre a solução geral dos sistemas de equações diferenciais
a)
dx
x 2y
dt
dy
4x 3y
dt
b)
0,1 0,075 1,5
Q' Q
0,1 0,2 3
c)
dx
2x y
dt
dy
5x 4y
dt
Introdução aos Métodos Numéricos
22
d)
1 5 sin t
X ' X
1 1 2cos t
(2) Sejam T1 = T1(t) e T2 = T2(t) as temperaturas, no instante t, nos ambientes 1 e 2, res-
pectivamente. Admita que as temperaturas estejam variando nas seguintes taxas
dT1
T1 4T2
dt
dT2
2T1 3T2
dt
Suponha que T1(0) = 20 e T2(0) = 10.
Determine as temperaturas, no instante t, nos ambientes 1 e 2.
(3) O tanque A contém 50 litros de água que foram dissolvidos 25 gramas de sal. Um se-
gundo tanque B, contém 50 litros de água pura. Bombeia-se o líquido para dentro e para
fora dos tanques as taxas indicadas na figura. O sistema de equações diferenciais que des-
creve o problema de tanques de mistura é
1
1 2
2
1 2
dx 2 1
x x
dt 25 50
dx 2 2
x x
dt 25 25
onde as condições iniciais são dadas por 1x (0) 25 e 2x (0) 0 . Aplicando autovetor e
autovalor de matriz encontre a solução do sistema de equações diferenciais
3.1.6.Solução de sistemas lineares
Vamos considerar a solução do sistema linear: usando o fato de que qualquer
matriz quadrada pode ser expressa como um produto de uma matriz triangular infe-
rior [L] e uma matriz triangular superior [U], temos:
(3.1)
Considerando n = 3,
(3.2)
[A]{x} {B}
n n[A]
[A] [L][U]
11 12 13 11 12 13
21 22 23 21 22 23
31 32 33 31 32 33
a a a l 0 0 1 u u
a a a l l 0 0 1 u
a a a l l l 0 0 1
Profª Celme Torres
23
Resolvendo o sistema para as variáveis e teremos:
(3.3)
(3.4)
(3.5)
(3.6)
Uma vez que decomposição da matriz A é feita, a resolução do sistema é feita da seguinte
forma:
(3.7)
(3.8)
Definindo um vetor [Z], como:
(3.9)
Portanto,
(3.10)
O sistema pode ser resolvido por substituição:
onde, i = 1,..., n (3.11)
Encontrados os valores de {Z} podemos calcular {X} usando:
(3.12)
Ou,
(3.13)
A solução é dada então por:
ijl iju
11 11l a 21 21l a 31 31l a
12
12
11
a
u
l
22 22 21 12l a l u
32 32 31 12l a l u
j 1
ij ij ik kj
k 1
i 1
ij ij ik kj
k 1ii
l a l u se i j
1
u a l u se i j
l
[A]{X} {B}
[L][U]{X} {B}
[U]{X} {Z}
[L]{Z} {B}
i 1
i ik k
k 1
i
ii
b l z
z
l
[U]{X} {Z}
1 12 2 13 3 1n n 1
2 23 3 2n n 2
3 3n n 3
n n
x u x u x ... u x z
x u x ... u x z
x ... u x z
x z
Introdução aos Métodos Numéricos
24
(3.14)
A solução de sistemas lineares pode ser resolvida também pelo Método de Gauss.
Exemplo 3.6
Seja o sistema:
Para eliminar x1 da segunda linha, basta multiplicar a primeira linha por 5/10 , isto é, por
a21/a11 , o que torna o coeficiente de x1 na primeira linha igual ao de x1 na segunda linha.
Em seguida basta subtrair a primeira linha modificada da segunda, o que zerará o coefi-
ciente de x1, eliminando-o da segunda equação.
Repetir essa operação para a terceira linha, multiplicando a primeira por a31/a11 , antes
de subtraí-la da terceira. O mesmo para a quarta linha, quinta linha etc...
Na enésima linha multiplica-se a primeira por an1/a11 e subtrai-se da enésima, eliminando-
se, dessa linha, a variável x1.
O mesmo raciocínio é repetido para a variável x2 da terceira equação em diante, para a
variável x3 da quarta em diante, até que a penúltima variável, xn-1 seja eliminada da última
equação.
Nesse momento o sistema estará triangularizado.
A partir desse ponto começa a segunda fase, “backward”, quando se calcula, na última
equação, a última variável, xn ; leva-se xn à penúltima equação e se calcula a penúltima
variável, xn-1 ; leva-se xn e xn-1 à antepenúltima equação e se calcula a antepenúltima vari-
ável, xn-2 , etc... até ser calculada a primeira variável na primeira equação.
Sendo Li a linha i , a eliminação da variável xj dessa linha se dará pela operação:
Nova Li = Li – aij/ajj . Lj , variando-se j de 1 a n-1 e i de j+1 até n.
No sistema apresentado:
n 1
i i i,k 1 k 1
k 1
x z u x
1 2 3
1 2 3
1 2 3
10x 5x 3x 17
5x 8x 2x 19
2x 3x 8x 0
1 2 3
1 2 3
1 2 3
10x 5x 3x 17
5x 8x 2x 19
2x 3x 8x 0
Profª Celme Torres
25
Observe que após eliminarmos x1 da segunda à última equação, esquecendo-se de L1, fica-
se com um sistema de 2 equações a 2 incógnitas, isto é, diminuiu-se de 1 a ordem do sis-
tema anterior; no caso passamos a ter um sistema de duas equações a duas incógnitas.
Para se eliminar x2 da terceira equação, trabalha-se da mesma maneira:
Tem-se agora um sistema de 1 equação a 1 incógnita.
O sistema completo fica sendo:
O sistema está triangularizado e começa, nesse momento, a segunda parte, a solução
propriamente dita do sistema.
2 2 1L L 5/10L
2 1 2 3L 0x (8 2,5)x (2 1,5)x 19 8,5
2 1 2 3L 0x 5,5x 0,5x 10,5
3 1 2 3L 0x 2,0x 7,4x 3,4
2 3
2 3
5,5x 0,5x 10,5
2,0x 7,4x 3,4
3 3 2L L 2,0 / 5,5L
3 2 3L 0x (7,4 0,182)x 3,4 3,818
37,218x 7,218
1 2 3
1 2 3
1 2 3
10x 5x 3x 17
0x 5,5x 0,5x 10,5
0x 0x 7,218x 7,218
Introdução aos Métodos Numéricos
26
Solução:
Algumas Considerações:
Numa dada matriz, quando se subtrai de uma linha outra linha multiplicada por uma
constante, o determinante da matriz não se altera. Como no processo de eliminação,
o que se fez foi sempre subtrair de uma linha outra linha multiplicada por uma cons-
tante, conclui-se que ao ser triangularizada, a matriz original manteve constante o
valor do determinante.
O determinante de matriz triangular é dado pelo produto dos elementos da diagonal
principal. Dessa maneira, o determinante da matriz original, antes de ser triangulari-
zada pelo método de Eliminação de Gauss, será dado pelo produto dos elementos da
diagonal principal da matriz triangular resultante.
No processo de eliminação da variável xj da linha i, faz-se a transformação
ij
i j
jj
a
L L
a
. No caso de ajj ser zero, deve-se renumerar as linhas para se conseguir ajj
diferente de zero.
Na verdade, o ideal seria conseguir uma linha j tal que, além de ter ajj diferente de
zero, tivesse esse elemento com o maior módulo possível e os demais elementos da
linha j o menor possível para que alterasse o mínimo a linha i, da expressão
ij
i j
jj
a
L L
a
. Dessa maneira diminuem os erros propagados de uma linha para outra no
processo de eliminação. Levar em conta a escolha adequada da linha que será a linha
j, para minimizar a propagação dos erros aumentando a precisão do método, é a
chamada Condensação Pivotal.
3
2
1
x 1,0
x (10,5 0,5( 1,0)) / 5,5 2,0
x (17 5.2,0 3.( 1,0)) /10 1,0
1x 1,0 2x 2,0 3x 1,0
Profª Celme Torres
27
Exercícios
Resolva os sistemas abaixo aplicando o método de Gauss
a) b)
1 2 3
1 2 3
1 2 3
2x 6x x 7
x 2x x 1
5x 7x 4x 9
x 3y 2z 7
4x y 3z 5
2x 5y 7z 19
Introdução aos Métodos Numéricos
28
4. Métodos dos Resíduos Ponderados
São métodos de aproximação utilizados para resolver equações diferenciais. Estes métodos
têm diversos procedimentos, dentre os quais se destacam o Método de Galerkin, o Método
da Colocação, o Método do Subdomínio e o Método dos Mínimos Quadrados.
A solução aproximada, para o caso de problemas unidimensionais, é do tipo:
(4.1)
Onde são as funções que satisfazem as condições de contorno do problema e são
coeficientes que serão determinados a partir da imposição de que o erro resultante do uso
desta solução aproximada seja mínimo. A seguir essa solução é substituída diretamente na
equação diferencial.
Como a solução aproximada não satisfaz a equação, haverá como resultado da
substituição um erro ou resíduo .
Esse erro é então multiplicado por uma função ponderadora e a integral desse pro-
duto é igualada a zero:
(4.2)
Onde D é o domínio do problema e o número de funções ponderadoras é igual ao
número de coeficientes desconhecidos .
Os métodos dos resíduos ponderados mais conhecidos são o método de colocação pontual,
método dos subdomínios, método dos mínimos quadrados e método de Galerkin. Veremos
a seguir cada um deles.
n
n i i
i 1
h (x) a N (x)
iN (x) ia
nh (x)
R(x)
iW (x)
i
D
W (x)R(x)dx 0
iW (x)
ia
Profª Celme Torres
29
4.1.1. Método de colocação pontual
Para este método,as funções ponderadoras são funções ponderadoras do tipo de
Delta-Dirac .
Uma das propriedades da função Delta-Dirac é:
(4.3)
Sendo assim, se a função aproximada tem n coeficientes , então:
(4.4)
Portanto,
(4.5)
Ou ainda,
(4.6)
Expandindo a expressão acima teremos:
(4.7)
iW (x)
ix x
i if (x) x x dx f (x )
ia
1 2 nR R x,a ,a ,...,a
i 1 2 n
D
W (x)R x,a ,a ,...,a dx 0
i 1 2 n
D
x x R x,a ,a ,...,a dx 0
1 1 2 n
2 1 2 n
3 1 2 n
n 1 2 n
R x ,a ,a ,..., a 0
R x ,a ,a ,..., a 0
R x ,a ,a ,..., a 0
R x ,a ,a ,..., a 0
Introdução aos Métodos Numéricos
30
Em outras palavras, este método impõe que o resíduo seja zero em determinados pontos do
domínio . O número de pontos n é o mesmo número de coeficientes da solu-
ção aproximada.
Exemplo 4.1
Vamos considerar nossa já conhecida equação diferencial:
Sujeita as condições de contorno
e
Uma solução aproximada para esta equação pode ser:
Observe que satisfaz as condições de contorno no problema.
Usando o método de colocação pontual, nos vamos resolver o problema de forma que o
resíduo R(x) seja zero para , assim:
Logo,
Impondo a condição de R(x) = 0
Deste modo,
1 2 nx , x ,..., x
2
2
2
d h
10x 0
dx
0 x 1
h(0) 0 h(1) 0
1 1h (x) a x 1 x
1h (x)
1x 0,5
2
2
2
2
1 1 1
1 1
2
12
d h
10x R(x)
dx
h a x(1 x) a x a x
dh
a 2a x
dx
d h
2a
dx
2
1 1R(x,a ) 2a 10x
1
2
1
1
R(x,a ) 0
2a 10x 0
a 1,25
1h (x) 1,25x 1 x
Profª Celme Torres
31
Exemplo 4.2
Considerando a mesma equação diferencial, com o objetivo de verificar o acréscimo de
precisão da solução ao se aumentar o número de termos da função de aproximação, va-
mos agora considerar uma solução aproximada do tipo:
Assim,
Logo,
Simplificando,
Como agora existem temos dois coeficientes a determinar, vamos impor que
em dois pontos do domínio do problema. Podemos escolher esses pontos
aleatoriamente.
Escolhendo e , temos:
Resolvendo o sistema acima para a1 e a2,
Assim ,
22 1 2h (x) a x 1 x a x 1 x
2
2
2
2 2 2 3
1 2 1 1 2 2
2
1 1 2 2
2
1 2 22
d h
10x R(x)
dx
h a x 1 x a x 1 x a x a x a x a x
dh
a 2a x 2a x 3a x
dx
d h
2a 2a 6a x
dx
2
1 2 1 2 2R(x,a ,a ) 2a 2a 6a x 10x
2
1 2 1 2 2R(x,a ,a ) a a 3a x 5x
1 2R(x,a ,a ) 0
1x 0,3 2x 0,6
2
1 2 2
2
1 2 2
a a 3a (0,3) 5(0,3) 0
a a 3a (0,6) 5(0,6) 0
1
2
a 0,6
a 1,5
22h (x) 0,6x 1 x 1,5x 1 x
Introdução aos Métodos Numéricos
32
Comparação gráfica entre e e a solução exata.
Figura 5.2 – Comparação gráfica entre as soluções encontradas
utilizando o Método de Colocação Pontual e a solução exata.
4.1.2. Método dos subdomínios
Neste caso, todas as funções ponderadoras são a unidade . Isto equivale a reque-
rer que a integral do resíduo desapareça em determinados intervalos do domínio.
O número desses intervalos é igual ao número de coeficientes indeterminados da solução
aproximada.
Vamos aplicar esse método à equação diferencial dos exemplos anteriores.
Exemplo 4.3
Vamos considerar nossa já conhecida equação diferencial:
Sujeita as condições de contorno
1h (x) 2h (x)
0,00
0,05
0,10
0,15
0,20
0,25
0,30
0,35
0,40
0,45
0 0,2 0,4 0,6 0,8 1 1,2
Solução h1
Solução h2
Solução Exata
iW (x) 1
2
2
2
d h
10x 0
dx
0 x 1
Profª Celme Torres
33
e
Considerando a solução aproximada,
Como há apenas um coeficiente a ser determinado, o intervalo de integração será o pró-
prio domínio do problema .
Vimos que,
Assim,
Portanto,
No caso de escolhermos uma função do tipo:
Vimos que,
Como agora existem dois coeficientes a serem determinados, usaremos dois intervalos
e .
h(0) 0 h(1) 0
1 1h (x) a x 1 x
0 x 1
2
1 1R(x,a ) 2a 10x
1
1
0
R(x,a )dx 0
1
3
1
0
10
2a x x 0
3
1
5
a
3
1
5
h (x) x 1 x
3
22 1 2h (x) a x 1 x a x 1 x
2
1 2 1 2 2R(x,a ,a ) 2a 2a 6a x 10x
0 x 0,5 0,5 x 1,0
0,5
1 2
0
1,0
1 2
0,5
R(x,a ,a )dx 0
R(x,a ,a )dx 0
Introdução aos Métodos Numéricos
34
O que dá origem ao seguinte sistema de equações:
Resolvendo o sistema, temos:
e
Assim,
Comparação gráfica entre e e a solução exata utilizando o método dos
subdomínios.
Figura 5.3 – Comparação gráfica entre as soluções encontradas
utilizando o Método dos Subdomínios e a solução exata.
4.1.3. Método dos mínimos quadrados
Este método requer que a integral do quadrado do resíduo seja minimizado. A integral do
quadrado do resíduo é dada por:
(4.8)
1 2
1 2
a 0,25a 0,4167
a 1,25a 2,917
1
5
a
6
2
5
a
3
22
5 5
h (x) x 1 x x 1 x
6 3
1h (x) 2h (x)
0,00
0,05
0,10
0,15
0,20
0,25
0,30
0,35
0,40
0,45
0,50
0 0,2 0,4 0,6 0,8 1 1,2
Solução h1
Solução h2
Solução Exata
2
1 2 n
D
Q R(x,a ,a ,...,a ) dx
Profª Celme Torres
35
Minimizando essa integral, temos:
(4.9)
Ou,
(4.10)
Ou seja, por esse método, as funções ponderadoras são:
(4.11)
2
1 2 n
D
1 1
R(x,a ,a ,...,a ) dx
Q
0
a a
2
1 2 n
D
2 2
2
1 2 n
D
n n
R(x,a ,a ,..., a ) dx
Q
0
a a
R(x,a ,a ,..., a ) dx
Q
0
a a
2
1D
R dx 0
a
1D
R
Rdx 0
a
2
2D
R dx 0
a
2D
R
Rdx 0
a
2
nD
R dx 0
a
nD
R
Rdx 0
a
1
1
2
2
n
n
R
W (x)
a
R
W (x)
a
R
W (x)
a
Introdução aos Métodos Numéricos
36
Exemplo 4.4
Aplicando o método dos mínimos quadrados a nossa já conhecida equação diferencial:
Sujeita as condições de contorno
e
Usado a função aproximada,
Vimos que,
Logo,
Portanto,
Que é o mesmo resultado obtido usando o método dos subdomínios. Desta forma, se usar-
mos uma função de aproximação com dois coeficientes, obteremos o mesmo resultado do
método dos subdomínios.
É importante ressaltar que essa igualdade dos resultados entre os métodos acontece para a
equação diferencial usada nos problemas e não se estende genericamente para qualquer
tipo de equação diferencial.
4.1.4. Método de Galerkin
Neste método, as funções ponderadoras são iguais as funções de aproximação
. Assim,
2
2
2
d h
10x 0
dx
0 x 1
h(0) 0 h(1) 0
1 1h (x) a x 1 x
2
1 1R(x,a ) 2a 10x
1
1
R
W (x) 2
a
1 1
2
1
10 0
R
Rdx a 5x dx 0
a
1
5
a
3
iW (x)
iN (x)
Profª Celme Torres
37
(4.12)
Sendo,
o conjunto , com . Podendo
o mesmo ser escrito ma forma:
(4.13)
Exemplo 4.5
Vamos aplicar o método de Galerkin para resolver a equação diferencial:
Sujeita as condições de contorno
e
Usado como função de aproximação,
Vimos que,
e,
Portanto,
Logo,
Usando agora uma função de aproximação com dois coeficientes, temos:
E repetindo o mesmo procedimento, teremos: e,
i i
D D
W (x)R(x)dx N (x)R(x)dx 0
1 2 n[N] N , N ,..., N i
D
N (x)R(x)dx 0 i 1,2,...,n
T
D
N R(x)dx 0
2
2
2
d h
10x 0
dx
0 x 1
h(0) 0 h(1) 0
1 1h (x) a x 1 x
1N (x) x(1 x)
2
1R(x) 2a 10x
1
2
1
0
x(1 x)( 2a 10x )dx 0 1
3
a
2
22 1 2h (x) a x 1 x a x 1 x
1
3
a
2
2
5
a
3
Introdução aos Métodos Numéricos
38
4.1.5. Exercícios
(1) Aplique o método de colocação pontual e o método dos mínimos quadrados para
encontrar a solução aproximada da equação diferencial no intervalo de 0 x 20 .
Utilize um programa gráfico (Maple ou Excel) para representar graficamente as
soluções indicando qual mais se aproxima da solução exata da EDO.
y'' 0,25x 0 y(0) 2 e y'(0) 2
(2) Aplique o método de Galerkin e o método dos subdomínios utilizando dois
coeficientes para encontrar a solução aproximada das equações diferenciais no
intervalo de 0 x 1 . Utilize um programa gráfico (Maple ou Excel) para representar
graficamente as soluções indicando qual mais se aproxima da solução exata da EDO.
2
2
d y dy
2x
dx dx
y(0) 1 e y'(0) 3
(3) Encontre a solução aproximada da equação diferencial no intervalo indicando
utilizando os métodos dos resíduos ponderados, colocação pontual, mínimos
quadrados, subdomínios e o método de Galerkin. Represente graficamente as soluções
aproximadas comparando com a solução exata da EDO.
(a) 216y'' 8y ' 0,150x 0 x 10 y(0) 2 e y'(0) 1
(b) ty'' 3y ' 4y 2e 0 x 1 y(0) 1 e y'(0) 0
(c)
2
2
d y
4y 3sin 2t
dt
0 x 2 y(0) 2 e y'(0) 1
(4) Considere a equação diferencial que rege a distribuição da temperatura T(x), em uma
barra trocando calor com o ambiente, cuja temperatura é de 25ºC.
Sujeita as condições de contorno:
T(0) 150
dT
(2) 0
dx
Proponha uma solução com um termo e outra com dois termos e ache os valores dos coefi-
cientes dessas soluções usando o método de colocação pontual, o método dos subdomínios,
dos mínimos quadrados e o método de Galerkin.
Compare a graficamente a solução com dois termos de cada método com a solução analíti-
ca e determine qual método mais preciso, neste caso.
2
2
d T
3(T 25) 0
dx
0 x 2
Profª Celme Torres
39
5. Introdução aos Métodos Variacionais
A análise de problemas de engenharia envolve o isolamento de um elemento diferencial da
região ou domínio onde o fenômeno está ocorrendo.
MODELO MATEMÁTICO RESULTANTE EQUAÇÃO DIFERENCIAL
Tipos de solução de uma Equação Diferencial
1. Solução analítica ou Exata
2. Solução aproximada ou solução numérica.
A solução analítica é a ideal, pois ela se constitui de uma expressão (equação) que satisfaz
a equação diferencial e suas condições de contorno, portanto para obtermos o valor exato
da variável dependente da equação diferencial em qualquer ponto do domínio do problema,
basta substituirmos as coordenadas desse ponto na equação solução.
5.1.1. Cálculo variacional
O cálculo variacional é uma parte da análise matemática que trata de problemas de máxi-
mos e mínimos de tipos especiais de funções, chamadas FUNCIONAIS.
Vários são os motivos que tornam convenientes o estudo do cálculo variacional. O princi-
pal deles é que este se constitui em uma alternativa que permite definir com mais clareza
os conceitos básicos que correspondem a problemas de engenharia governados por siste-
mas de equações diferenciais.
O cálculo variacional é especialmente aplicável a problemas de mecânica dos meios contí-
nuos. Através do uso de princípios de energia, é possível formular os problemas na forma
variacional.
5.1.2. Máximo e mínimo de funções
Suponha a seguinte função apresentada no gráfico.
Introdução aos Métodos Numéricos
40
O domínio de F(x) é a bx x x . A relação funcional implica que, para cada valor da va-
riável x, um único valor será obtido para a função F(x), dentro do domínio.
O menor valor de F(x) dentro do domínio é o mínimo absoluto da função F(xa) e o máxi-
mo absoluto é dado por F(xb). Além disso, podem existir valores mínimos e máximos rela-
tivos.
Por exemplo; a função F(x) terá mínimo relativo ou um máximo relativo, para um certo
valor da variável independente x, quando todos os valores de F(x) na vizinhança infinite-
simal de x: x x , com x tão pequeno quanto se queira, seja maior ou menor, respecti-
vamente, que F(x). Para estudar as condições que determinam os valores extremos relat i-
vos, é conveniente expandir F(x) em uma Série de Taylor em torno de x.
2 3
2 3
2 3
x
dF 1 d F 1 d F
F(x x) F(x) x x x ...
dx 2! 3!dx dx
(4.1)
O incremento F será dado por:
2 3
2 3
2 3
x
dF 1 d F 1 d F
F(x x) F(x) x x x ...
dx 2! 3!dx dx
(4.2)
Observe que haverá um mínimo relativo quando:
F(x) é sempre menor que F(x x) , não interessando o valor de x , para F 0
F positivo;
Para F 0 F negativo máximo relativo.
F(x)
x
xb x2 x1 xa
Profª Celme Torres
41
Como x é infinitamente pequeno, x é sempre maior que
2x e que 3x e assim su-
cessivamente.
Portanto,
Se
x
dF
0
dx
, o sinal de F será indefinido, pois depende do sinal de x ;
Se
x
dF
0
dx
, não teremos ponto de mínimo e de máximo.
Então, para se garantir que a função terá pontos de mínimo e de máximo relativo a primei-
ra derivada obrigatoriamente será zero
x
dF
0
dx
.
Os pontos que satisfazem essa condição são chamados de pontos estacionários. Neste ca-
so, como
2x é sempre positivo, a característica de um ponto estacionário dependerá do
sinal da 2ª derivada de F(x):
2
2
x
d F
dx
.
Se
2
2
x
d F
0
dx
mínimo relativo
Se
2
2
x
d F
0
dx
máximo relativo
Se a derivada segunda for nula, haverá um ponto neutro. Neste caso, devem ser analisados
os sinais das derivadas de ordem superior, para se determinar a natureza do ponto estacio-
nário. Se todas as derivadas até a ordem n 1 forem nulas, a natureza do ponto estacioná-
rio dependerá do sinal da enésima derivada.
Se n for par, será aplicado para o último termo da Série de Taylor o mesmo critério que se
usa para a 2ª derivada.
Se n for impar, o termo
n
n
n
x
1 d F
x
n! dx
, vaiará com o sinal de x e não se terá um máximo
ou mínimo relativo, e sim um ponto indefinido.
Introdução aos Métodos Numéricos
42
No caso da 2ª derivada ser nula, tais pontos são denominados de pontos de sela ou pontos
de inflexão.
Uma análise prévia pode ser estendida a funções de mais de uma variável independente.
Por exemplo, funções do tipo: F(x, y) , com duas variáveis independentes. Ao se expandir
esta função em uma Série de Taylor, temos:
2 2 2
2 2
2 2
x,y x,y x,y x,y x,y
F F(x x, y y) F(x, y)
F F 1 F F F
F x y x 2 x y y ...
x x 2! x yx y
(4.4)
Como visto anteriormente, haverá um ponto estacionário quando:
x,y x,y
F F
0
x y
E a característica do ponto estacionário será determinada, na forma:
2 2
2
2 2 2
x,y x,y2 2
2 2 2 2
2x,y x,y x,y
2
x,y x,y
F F
x yx xF F F
x 2 x y y x y
yx yx y F F
y
x y y
(4.5)
A fórmula acima é quadrática e o sinal depende da natureza da matriz das derivadas de 2ª
ordem, assim:
Quando a matriz é positiva definida, com autovalores positivos, tem-se um ponto
de mínimo.
Quando a matriz é negativa definida, com autovalores negativos, tem-se um ponto
de máximo.
Exemplo 5.1
Considere a função
3 2
3
f (x) x x 6x
2
Os pontos estacionários dessa função são dados por:
Profª Celme Torres
43
2
df
3x 3x 6 0
dx
A qual será definida para
1
2
x 1
x 2
Ao se verificar a segunda derivada :
2
2
d f
6x 3 0
dx
Pode verificar que f(1) é um ponto de mínimo relativo e f(-2) é um ponto de máximorelati-
vo.
Exemplo 5.2
Considere a função
2 2f (x, y) x 2xy 4y 2x 8y
Haverá um ponto estacionário quando:
f
2x 2y 2 0
x
f
2x 8y 8 0
y
Reescrevendo temos:
2x 2 2
2x 8y 8
A matriz da 2ª derivada é dada por:
2 2
2
2 2
2
F F
x y 2 2x
2 0F F
x y y
Cujos autovalores são positivos ponto estacionário é um ponto de mínimo.
5.1.3. Funcional
A seção anterior tratava de determinar pontos estacionários de funções. O problema do cál-
culo variacional também se refere à determinação de pontos estacionários, mas definidos
para funções especiais, chamadas funcionais.
Introdução aos Métodos Numéricos
44
Por exemplo, seja a seguinte função: f (x, y,z,...)
Ao se dar valores às variáveis independentes, um valor numérico é obtido para a função f.
Determina-se assim, um grupo de valores de x, y, z, ..., para o qual define-se um ponto es-
tacionário para f.
O FUNCIONAL é um tipo especial de função, no qual as variáveis independentes não são
simples números, mas funções desconhecidas. Ou seja, um funcional é um grupo (ou uma
função) que depende de outras funções. Funcionais são chamados de F ou .
1 2
2
0
1
x xx
0
3df 8d f
F (2f )dx
dx dx
F G(g,g ,g , x)dx
Onde, x
df
f
dx
e
2
xx 2
d f
g
dx
Tanto F quanto G são considerados FUNCIONAIS. Quando se adota diferentes valores
para as funções desconhecidas f, se obtém diferentes valores para os funcionais.
Em geral um funcional pode depender de várias funções, as quais podem depender de uma
ou mais variáveis independentes.
5.1.4. Funcional de uma função
Considere o funcional no domínio
a bx x x
b
a
x
x
x
F G(g,g , x)dx
Condições de contorno para g:
a a
b b
g(x ) g
g(x ) g
Profª Celme Torres
45
Considere agora que g é uma função que faz o funcional F ser estacionário.
Definindo um grupo de funções:
g( ,x) g(x) (x)
onde, é um escalar infinitesimal e (x) é uma função conhecida com as seguintes condi-
ções:
a
b
(x ) 0
(x ) 0
Podendo escrever então:
f g g
onde, g (x)
Sendo o símbolo variacional e significa a 1ª variação de g.
No cálculo variacional α corresponde a x da Série de Taylor. Dessa forma, o funcional é
dado por:
b
a
x
x
xx
G G
F g g dx
g g
(4.6)
5.1.5. Propriedades de um Funcional:
d( f ) (df )
(f g) f g
xa xb
𝑛(𝑥)
𝑓(𝑥)
𝑔 𝑥 = 𝑓 𝑥 + 𝑛(𝑥)
𝛼𝑛(𝑥)
Introdução aos Métodos Numéricos
46
(f.g) f g fg
2
f f g fg
g g
n n 1(f ) nf f
Pretende-se agora determinar as condições para que o funcional seja estacionário através
do estudo da variação de f f . Isto é semelhante ao caso de funções onde a existência de
um ponto estacionário é determinada para um certo valor de x através do estudo variacio-
nal de f(x) no ambiente x x .
Se considerarmos,
b
a
x
x x
x
(g) F(f , f , x)dx (4.8)
Onde f e são conhecidos, portanto o funcional só depende de α.
Expandindo em torno de f através de uma Série de Taylor, temos:
2 3
2 3
2 3
0 0 0
d 1 d 1 d
(g) (f ) ...
d 2! 3!d d
(4.9)
2 3
2 3
2 3
0 0 0
d 1 d 1 d
(g) (f ) ...
d 2! 3!d d
(4.10)
Dessa forma é possível definir a variação do funcional , de acordo com a ordem da equa-
ção diferencial:
1ª ordem
0
d
d
2ª ordem
2
2 2
2
0
d
d
Profª Celme Torres
47
3ª ordem
3
3 3
3
0
d
d
A equação para variacional total pode ser escrita na forma:
2 3
1 1
...
2! 3!
(4.11)
De maneira semelhante aos caso de funções, a natureza de ponto estacionário, isto é, se é
um ponto de máximo, de mínimo ou de inflexão, dependerá do sinal da derivada de 2ª
ordem.
5.1.6. Determinação do Funcional - Equação de Euler
Considere uma integral da forma:
(4.12)
Com,
(4.13)
Onde o integrando é uma dada função de x, h e hx. Tanto I(h) como F são
chamados de funcional. A palavra funcional significa função de funções.
Considere o funcional . Uma mudança de em h, onde é uma constante e
w é uma função, é chamada de variação de h e é denominada:
(4.14)
Onde é chamado de símbolo variacional.
A primeira variação de F em h é definida por:
(4.15)
b
x
a
I(h) F(x,h,h )dx
h h(x) x
dh
h
dx
xF(x,h,h )
xF(x,h,h ) w
h w
x
x
F F
F h h
h h
Introdução aos Métodos Numéricos
48
O símbolo age como um operador diferencial com respeito às variáveis dependentes,
com a seguinte propriedade:
(4.16)
Considere agora o problema de se determinar o mínimo do funcional:
(4.17)
A condição necessária para minimizar esse funcional é que a primeira variação de I deve
ser zero, ou:
(4.18)
Usando a condição para o funcional ser estacionário, temos:
(4.19)
Podemos ainda escrever:
(4.20)
Assim,
(4.21)
Usando integral por partes ( ) com, e
(4.22)
b b
a a
dh d( h)
dx dx
h(x)dx h(x)dx
b
x
a
I(h) F(x,h,h )dx
b
x
a
I F(x,h,h )dx 0
b
x
x
a
F F
I h h dx 0
h h
x
dh d( h)
h
dx dx
b
x
a
F F d( h)
I h dx 0
h h dx
udv uv vdu v h
x
F
u
h
b b
x x a
a
F d F F
I h h dx h
h dx h h
Profª Celme Torres
49
(4.23)
O que só pode ser verdade para qualquer variação de se:
(4.24)
com, Essa expressão é gerada pela condição de contorno do problema.
Este resultado é uma equação diferencial com h e hx como pseudovariáveis independentes.
A equação 4.24 é chamada de Equação de Euler ou Euler-Lagrange.
É importante salientar que, a função h(x) que minimiza o funcional I é a mesma função que
satisfaz a equação diferencial de Euler.
Exemplo 5.3
Determine o funcional da equação diferencial:
Sujeita as condições de contorno
e
Como h é uma constante nos extremos do intervalo , a variação de é zero
nestes pontos. Este tipo de condição de contorno é denominada essencial ou geométrica.
Neste caso,
b b
x x a
a
F d F F
I h dx h 0
h dx h h
h
x
F d F
0
h dx h
b
x a
F
0
h
2
2
2
d h
10x 0
dx
h(0) 0 h(1) 0
0 x 1 h
1
x 0
F
h 0
h
Introdução aos Métodos Numéricos
50
Determinação do funcional F usando a equação de Euler:
Assim, usando a equação (10), temos:
Integrando a 2ª integral por partes e considerando e , temos:
(E.1)
Analisando novamente a 2ª integral da expressão acima, usando d( h) (dh) , a seguinte
identidade é considerada:
Se, na equação (E.1) assumirmos que,
, teremos:
Como a variação de x
2
é zero, e usando , temos:
Portanto, o funcional F associado à equação diferencial do problema é:
2
2 2
2
F d F d h d dh
10x 10x
h dx x dx dxdx
1
2
0
1 1
2
0 0
d dh
I 10x hdx 0
dx dx
d dh
I 10x hdx hdx 0
dx dx
u h
dh
v
dx
1 11
2
0
0 0
d hdh dh
I 10x hdx h dx 0
dx dx dx
1 1 1 2
0 0 0
d hdh dh dh 1dh
I dx dx dx
dx dx dx dx 2 dx
1
0
dh
h 0
dx
1 2
2
0
1 dh
I 10x h dx
2 dx
b b
a a
Fdx Fdx
1 2
2
0
1 dh
I 10x h dx
2 dx
Profª Celme Torres
51
É bom enfatizar que a função h(x) que minimiza o funcional acima é também a solução da
ED que rege o problema. O funcional acima é chamado de forma variacional ou forma
fraca da ED. Lembrando que I não é uma solução, apenas o funcional da ED.
5.1.7. Exercícios
(1) Determine o funcional associado às equações diferenciais sujeita as condições de con-
torno apresentadas.
a) sendo,
x x
h(0) 0 h(5) 100
h (0) 1 h (5) 0
b)
2
2 2
2
d f
9f 5x
dx
onde, 0 x 5 sendo,
f (0) 9 e f (5) 0
c)
3
3
d f
5f 0
dx
onde, 0 x 4 sendo,
df (0)
0
dx
e
f (4) 0
(2) Determine a equação de Euler associada ao funcional
Sabendo que:
1 2
2
0
1 dh
I 10x h dx
2 dx
4 2
4 2
d h d h
8 4h 10
dx dx
0 x 5
b
x xx
a
I F(x,h,h ,h )dx 0
b
x xx a
F d F
h 0
h dx h
Introdução aos Métodos Numéricos
52
b
x
x a
F
h 0
h
Profª Celme Torres
53
6. Métodos Variacionais
6.1.1. Método de Ritz
O primeiro passo do método de Ritz é assumir uma solução aproximada para a equação
diferencial na forma de uma série finita:
(6.1)
Onde as constantes ai são chamadas coeficientes de Ritz e as funções Ni(x) são conhecidas
como funções de aproximação.
Como hn é uma solução da equação diferencial, hn deve satisfazer as condições de contorno
essenciais. Qualquer condição de contorno natural já está incluída na formulação variacio-
nal.
O passo seguinte no Método de Ritz usa o princípio de que a função que minimiza o
funcional também é a solução da equação diferencial associada ao funcional, assim
como:
(6.2)
Após a substituição de em , este funcional se torna função dos parâmetros
ou e a condição necessária para minimizar I será:
(6.3)
n
n i i
i 1
h a N (x)
h(x)
I(x)
n
n i i
i 1
h a N (x)
nh I(h)
1 2 na ,a ,...,a 1 2 nI(a ,a ,...,a )
1
2
n
I
0
a
I
0
a
I
0
a
Introdução aos Métodos Numéricos
54
O que resultará em um sistema de n equações com n coeficientes a determinar
.
Uma vez conhecidos os valores desses coeficientes (conhecidos como coeficientes de
Ritz), então a solução aproximada será conhecida.
Exemplo 6.1
Considere a equação diferencial:
Sujeita as condições de contorno
e
Como vimos o funcional associado a essa equação diferencial é:
Usando o método de Ritz, vamos assumir a seguinte solução aproximada:
Observe que h1 satisfaz as condições essenciais e de contorno, portanto:
i(a ,i 1,2,...,n)
nh
2
2
2
d h
10x 0
dx
0 x 1
h(0) 0 h(1) 0
1 2
2
0
1 dh
I 10x h dx
2 dx
1 1h a x(1 x)
1 1
2
1 1
1
1
h a x(1 x)
h a (1x x )
dh
a (1 2x)
dx
Profª Celme Torres
55
Logo,
Dessa forma:
Minimizando então I:
Portanto,
Vamos agora usar a seguinte a solução aproximada:
Usando os dois resultados acima, obtemos:
Minimizando agora I com respeito a e obteremos o seguinte sistema de equações:
2 2
21 1
1
dh a1
a (1 2x)
2 dx 2
1
2
23 4 1
1
0
a
I 10a x x 1 2x dx
2
1 14 5
2 2 3
1 1
0 0
2
1 1
x x 4
I 10a a x 2x x
4 5 3
a a
I
2 6
1
I
0
a
1a1 0
2 3
1
3
a
2
1
3
h x 1 x
2
2
2 1 2
22
1 2
h a x(1 x) a x (1 x)
dh
a (1 2x) a (2x 3x )
dx
2 2
1 1 2 1 2 2a a a a a aI
6 2 3 6 15
1a 2a
Introdução aos Métodos Numéricos
56
Resolvendo o sistema acima para e , obtemos:
e,
Dessa forma,
A seguir faremos uma comparação gráfica, usando o Maple, entre e e a so-
lução exata:
Figura 6.1 – Comparação gráfica entre as soluções encontradas
utilizando o Método de Ritz e a solução exata.
Verificamos que a função aproximada encontrada ao se usar o Método de Ritz através do
método variacional coincidiu com a solução obtida aplicando o Método de Galerkin. Tal
fato não foi uma coincidência. Considerando a equação resultante da minimização do fun-
cional I,
1 2
1
1 2
2
a aI 1
0
a 2 3 6
a 2aI 1
0
a 3 6 15
1a 2a
1
3
a
2
2
5
a
3
22
2 5
h x 1 x x 1 x
3 3
1h (x) 2h (x)
410h(x) x x
12
0,00
0,05
0,10
0,15
0,20
0,25
0,30
0,35
0,40
0,45
0 0,2 0,4 0,6 0,8 1 1,2
Solução h1
Solução h2
Solução Exata
Profª Celme Torres
57
Usando a solução aproximada,
Onde,
e,
Neste caso, como as funções de aproximação são invariáveis,
Substituindo a expressão acima na equação do funcional, temos:
Como, não depende de x,
Considerando que, é a equação diferencial associada ao funcional I, e que
a solução aproximada foi usada,
Teremos,
Para a expressão acima ser verdadeira para qualquer valor de ,
b b
x x a
a
F d F F
I h dx h 0
h dx h h
n
T
m i i
i 1
h h (x) a N (x) [a][N]
1 2 n[a] [a ,a ,...,a ]
1 2 n[N] [N , N ,..., N ]
T[N]
Th [a][N]
b b
T T
x x a
a
F d F F
I [a][N] dx [a][N] 0
h dx h h
[a]
b b
T T
x x a
a
F d F F
I [a] [N] dx [a][N]
h dx h h
x
F d F
h dx h
mh (x)
m x
F d F
dx R(x)
h dx h
b b
T T
x a
a
F
I [a] [N] R(x)dx [N]
h
[a]
Introdução aos Métodos Numéricos
58
O primeiro termo da expressão acima é o princípio básico do Método de Galerkin. O se-
gundo termo refere-se às condições de contorno naturais do problema.
Podemos então concluir que havendo um princípio variacional para o problema, a solução
através do uso do método de Galerkin coincide com a solução do método variacional. Essa
é a razão pelo qual o Método de Galerkin é o escolhido como fundamento para o Método
dos Elementos Finitos (MEF).
6.1.2. Exercícios
(1) Considere a seguinte equação diferencial
Sujeitas as condições de contorno:
a) Ache o funcional associado a essa equação diferencial
b) Pode a seguinte equação 1
x
T (x) sin
10
ser usada na solução aproximada? Se puder,
usando o Método de Ritz, ache a solução aproximada usando a função acima e compare
graficamente com a solução exata.
c) Usando a função aproximada , ache os coeficientes e
.
d) Faça um gráfico comparando as soluções (b) e (c) com a solução exata.
(2) Dada a seguinte equação diferencial
Sujeitas as condições de contorno:
b b
T T
x a
a
F
[N] R(x)dx [N] 0
h
2
2
d T
100 0
dx
0 x 10
T(0) 0
T(10) 0
2 1 2
x 3 x
T (x) a sin a sin
10 10
1a
2a
2
2
d y
6y 10x
dx
0 x 2
Profª Celme Torres
59
y(0) 1
y(2) 0
Encontre o funcional associado à equação diferencial e encontre uma solução aproxi-
mada usando o Método de Ritz.
6.1.3. Método de Rayleigh-Ritz
No método de Ritz apresentado anteriormente, as funções de Ni(x) eram válidas em todo o
domínio do problema. No método de Rayleigh-Ritz, essas funções de aproximação são vá-
lidas apenas em trechos do problema. Portanto, o domínio do problema é dividido em ele-
mentos e asfunções de aproximação são específicas para esses elementos.
Vamos analisar agora como encontrar essas funções de aproximação Ni(x), também co-
nhecidas como funções de forma.
Considere um problema cujo domínio é o intervalo que será dividido em M in-
tervalos ou elementos finitos. Nos limites de cada elemento, há dois nós – denominados
genericamente de i e j. Como, nesse caso, cada elemento e tem dois nós, uma interpolação
polinomial linear deve ser usada como função de interpolação no elemento, portanto:
De acordo com a Figura 5.4, podemos constatar que:
(6.4)
Portanto, (6.5)
a x b
eh (x) cx d
e
i i
e
j j
x x h h
x x h h
i i
j j
h cx d
h cx d
Introdução aos Métodos Numéricos
60
Figura 6.2 – Divisão de um domínio qualquer em elementos lineares.
Resolvendo o sistema para achar c e d, temos:
e (6.6)
Assim,
(6.7)
Onde e são as funções de forma.
Logo,
(6.8)
h (x)
e
hi
hj
xjxi
h(x)
x
x = a x = b
1 2 3 i j ne
1 e
2 e
m
e
j i
j i
h h
c
x x
i j j i
j i
h x h x
d
x x
je i
i j i i j j
j i j i
x x x x
h (x) h h N (x)h N (x)h
x x x x
iN (x) jN (x)
j i
j i j i
x x x x
[N]
x x x x
Profª Celme Torres
61
Figura 6.3 – Interpolação linear de um elemento.
Exemplo 6.2
Considere a já conhecida equação diferencial:
Sujeita as condições de contorno
e
Como vimos o funcional associado a esse equação diferencial é:
Se dividirmos o domínio em vários trechos (elementos), temos:
Sendo,
e,
Minimizando o funcional do elemento, temos:
h =N (x)h +
e
i i N (x)h j j
hi
hj
x
i je
1,0 1,0N (x)i N (x)j
2
2
2
d h
10x 0
dx
0 x 1
h(0) 0 h(1) 0
1 2
2
0
1 dh
I 10x h dx
2 dx
j
j
i
i
x 2x e
e e 2 e
x
x
dh 1 dh
I h 10x h dx
dx 2 dx
n
e
e 1
I I
eeh N a
ie
j
a
a
a
Introdução aos Métodos Numéricos
62
Mas como, podemos reescrever o termo da seguinte forma:
Sendo a matriz transposta de .
Da mesma forma podemos simplificar o segundo termo da integral, onde:
Portanto,
Ou ainda,
O que pode ser escrito de forma simplificada:
Sendo, [K]
e
a matriz de rigidez do elemento; {f}
e
é conhecido como vetor de carregamento
nodal do elemento e {a}
e
é o vetor dos valores aproximados a função no nó do elemento.
Depois que as matrizes [K]
e
e os vetores {f}
e
e {a}
e
de todos os elementos serem conheci-
dos, encontra-se o sistema global do problema, na forma:
j
j
i
i
x 2xe e e e
2
e e e e
x
x
dI dh dh dh 1 d dh
10x dx 0
dx 2 dxd a d a d a d a
eeh N a
e
T
e
dh
N
d a
T
N N
2
e e e e
e e e
e e
e
e e
e
dh 1 d dh 1 dh d dh
(2)
2 dx 2 dx dxd a d a d a
d dh dh
dx dxd a
d N d Nd
a a
dx dxd a
T
ed N d N
a
dx dx
j j
j
i
i i
x xx Te
T T e2
e
x
x x
d N d NdI dh
N 10x N dx a dx 0
dx dx dxd a
j j
j
i
i i
x xxT
e T T2
x
x x
d N d N dh
dx a N 10x N dx
dx dx dx
e e e
K a f
K a f
Profª Celme Torres
63
Agora as condições de contorno podem ser impostas. Voltando ao nosso exemplo. Se divi-
dirmos o domínio em dois elementos, temos dois elementos e três nós, conforme a Figura
5.6.
Figura 6.4 – Elementos lineares
Para cada elemento, temos:
Para os dois elementos do problema:
Analisando o termo,
Analisando agora a integral,
j
i
x T
e
j i
x
1 1d N d N 1
K dx
1 1dx dx x x
1 2 2 2K K
2 2
j
j
i
i
xx
T Te 2
x
x
dh
f N 10x N dx
dx
j
i
x e e
i jT i ij i
j ij jx
e e
j i
e
i
e
j
N (x ) N (x )dh (x ) dh (x )dh
N
N (x )N (x )dx dx dx
dh (x )0 1 dh (x )
1 0dx dx
dh (x )
dx
dh (x )
dx
e
1 e
2
1 2 3
0,5 0,5
Introdução aos Métodos Numéricos
64
Assim, para o elemento 1, temos:
Para o elemento 2,
Fazendo a montagem do sistema global:
Levando em consideração a equação da continuidade, podemos impor, que:
Dessa forma o sistema global se reduz a forma:
j
i
x 4 3 4
j j i iT2
4 4
j i j i j ix
x 4x x 3x10
10x N dx
10 x x 3x 4x x x
1
1
1
1
2
dh (x )
0,104
dx
f
dh (x )
0,313
dx
2
2
2
2
3
dh (x )
1,146
dx
f
dh (x )
1,771
dx
1
1
1 1 2
2 2
2
3 2
3
dh (x )
0,104
dx
h2 2 0
dh (x ) dh (x )
2 2 2 2 h 1,459
dx dx
0 2 2 h
dh (x )
1,771
dx
1 2
1 2dh (x ) dh (x )
dx dx
Profª Celme Torres
65
Usando a 2ª linha do sistema encontramos o valor para h2
A solução gráfica apresenta-se da seguinte forma:
Figura 6.5 – Comparação gráfica entre as soluções encontradas
utilizando o Método de Rayleigh-Ritz e a solução exata.
1
1
1
2
2
3 3
dh (x )
0,104
h2 2 0 dx
2 4 2 h 1,459
0 2 2 h dh (x )
1,771
dx
2h 0,365
0,00
0,05
0,10
0,15
0,20
0,25
0,30
0,35
0,40
0,45
0 0,25 0,5 0,75 1 1,25
Solução Exata
Rayleigh-Ritz
Introdução aos Métodos Numéricos
66
7. Métodos das Diferenças Finitas
Como mencionado anteriormente, a idéia básica do método das diferenças finitas (MDF)
consiste na aproximação das derivadas parciais de uma equação diferencial através de fór-
mulas de diferenças, e satisfação dessa equação aproximada de diferenças em determina-
dos pontos do domínio (região). Isto conduz a um sistema de equações algébricas cuja so-
lução fornece o valor da função desejada nos pontos pré-determinados.
Por definição, a derivada de uma função u(x) num ponto é dada por:
(7.1)
De forma aproximada, pode-se tomar:
(7.2)
Onde h é um incremento pequeno, porém finito. Esta aproximação é chamada de diferença
progressiva, já que no seu cálculo aparecem os valores de u nos pontos e
. De maneira análoga, pode-se definir uma aproximação de diferenças regressi-
vas na forma:
(7.3)
E um esquema de diferença central:
(7.4)
As aproximações acima podem ser visualizadas na Figura 6.1.
1x x
i
i i
h 0
x x
u(x h) u(x )du
lim
dx h
i
i i
x x
u(x h) u(x )du
dx h
1x x
1x x h
i
i i
x x
u(x ) u(x h)du
dx h
i
i i
x x
u(x h) u(x h)du
dx 2h
Profª Celme Torres
67
Uma forma alternativa de se obter fórmulas aproximadas de diferenças finitas é através de
séries de Taylor. Desta maneira, é possível estimar o erro cometido em cada tipo de apro-
ximação.
A expansão em série de Taylor do valor de u em , em torno do valor de u em
, é dada por:
(7.5)
Figura 7.1 – Interpretação gráfica das aproximações por diferenças finitas
Que pode ser reescrita na forma:
(7.6)
Tornando-se:
(7.7)
ix x h
ix x
i i i i
2 2 3 3 n n
i i 2 3 n
x x x x x x x x
du h d u h d u h d u
x(x h) u(x ) h...
dx 2 6 2ndx dx dx
x +hixi
y
x
x - hi
h h
y= u(x)
tg
i i i i
2 2 3 n n
i i
2 3 n
x x x x x x x x
u(x h) u(x )du h d u h d u h d u
...
dx h 2 6 2ndx dx dx
i
i i
x x
u(x h) u(x )du
dx h
Introdução aos Métodos Numéricos
68
Que é a aproximação de diferença progressiva. Pode-se notar que foram desprezados os
termos relativos à derivada de segunda ordem e de ordem superior. Como h é um valor pe-
queno, maior termo desprezado é igual a uma constante multiplicada por h, ou seja, é da
ordem de h.
Analogamente,
(7.8)
Que dá a aproximação:
(7.9)
Para se obter a aproximação de diferença central, subtrai-se (5.8) de (5.5), chegando-se a:
(7.10)
Ou,
(7.11)
O erro cometido pelo esquema de diferença central é menor que o de diferença progressiva
ou regressiva.
7.1.1. Derivadas de ordem superior
Para se obter aproximações por diferenças finitas para derivadas de ordem superior, pode-
se empregar séries de Taylor ou uma maneira mais intuitiva, através da aplicação repetida
da aproximação (6.2) a (6.4). Tomando-se, por exemplo, a derivada de segunda ordem, po-
de-se somar as expressões (6.5) e (6.8), obtendo-se:
(7.12)
i i i i
2 2 3 3 n n
i i 2 3 n
x x x x x x x x
du h d u h d u h d u
x(x h) u(x ) h ...
dx 2 6 2ndx dx dx
i
i i
x x
u(x ) u(x h)du
dx h
i i i
3 3 n n
i i 3 n
x x x x x x
du h d u h d u
x(x h) u(x h) 2h ...
dx 3 ndx dx
i
i i
x x
u(x h) u(x h)du
dx 2h
i i i
2 4 4 n n
2
i i i 2 4 n
x x x x x x
d u h d u h d u
u(x h) u(x h) 2u x h ...
12 3ndx dx dx
Profª Celme Torres
69
Reescrevendo a expressão (7.12), desprezando-se as derivadas de ordem superior, temos:
(7.13)
Que é a aproximação de diferença central para a derivada de segunda ordem.
Outra opção é aproximar a derivada segunda como uma diferença entre derivadas de pri-
meira ordem:
(7.14)
(7.15)
Cujo resultado é igual a equação (7.13). Os dois processos acima podem ser estendidos pa-
ra a determinação de qualquer ordem de derivada.
Exemplo 7.1
Seja o problema descrito pela equação e condições de contorno abaixo:
Resolver a equação diferencial pelo método das diferenças finitas, com aproximação de
diferenças centrais.
A equação aproximada de diferenças centrais é dada por:
i
2
i i i
2 2
x x
u(x h) 2u(x ) u(x h)d u
dx h
i i
i
2
x x h / 2 x x h / 2
2
x x
du du
dx dxd u
hdx
i
i i i i
2
2
x x
u x h u(x ) u x u(x h)
d u h h
hdx
2
2
d u
u x 0 (1 x 0)
dx
0 em x 0
u
0 em x 1
Introdução aos Métodos Numéricos
70
Reescrevendo os termos, temos:
(i)
Adotando-se , o problema irá apresentar 4 incógnitas, correspondentes ao valor
de u nos pontos:
É necessário então, gerar 4 equações, que são obtidas aplicando a equação (i) nesses 4
pontos.
Introduzindo as condições de contorno e arrumando em forma matricial:
Dessa forma, obtemos os valores de u nos pontos determinados:
A solução exata do problema é dada por:
Comparando graficamente a solução dada por diferenças finitas com aproximação central
com a solução exata da ED, temos:
i i i
i i2
u x h 2u(x ) u(x h)
u(x ) x 0
h
2 2i i i iu x h (h 2)u(x ) u(x h) h x
h 0,2
x 0, 2
x 0, 4
x 0,6
x 0,8
i
i
i
i
x 0,2 0,4u 1,96u(0,2) 0,0u 0,008
x 0,4 0,6u 1,96u(0,4) 0,2u 0,016
x 0,6 0,8u 1,96u(0,6) 0,4u 0,024
x 0,8 1,0u 1,96u(0,8) 0,6u 0,032
1,96 1 0 0 0,2u 0,008
1 1,96 1 0 0,4u 0,016
0 1 1,96 1 0,6u 0,024
0 0 1 1,96 0,8u 0,032
u(0,2) 0,0362
u(0,4) 0,0630
u(0,6) 0,0713
u(0,8) 0,0527
sinx
u
sin1 x
Profª Celme Torres
71
Figura 7.2 – Comparação gráfica da solução exata e da solução por diferenças finitas
com aproximação central.
Exemplo 7.2
Resolver a Equação Diferencial do exemplo anterior, com a condição de contorno natural,
dada por:
O problema agora possui 5 incógnitas, já que também não se conhece o valor de u no pon-
to . É necessário, então, obter mais uma equação pela aproximação de diferença re-
gressiva da condição de contorno natural.
Substituindo esta aproximação e a condição de contorno em no sistema de equa-
ções, temos:
du
1 em x 1
dx
x 1
x 1
du 1,0u 0,8u
1 u(1,0) 0,2 u(0,8)
dx 0,2
x 0
Introdução aos Métodos Numéricos
72
Uma alternativa à solução acima é a de se aplicar a equação (i) também no ponto ,
envolvendo o valor de u no ponto fictício e utilizar uma aproximação de diferen-
ças centrais para a condição de contorno natural. Nesta caso, o sistema final terá 5 equa-
ções e 5 incógnitas, na forma:
A solução da ED considerando o primeiro sistema através da condição de contorno natu-
ral com diferença finita regressiva é:
Considerando a condição natural com diferenças finitas aproximadas, temos:
A solução exata do problema é dada por:
Comparando graficamente a solução dada por diferenças finitas com aproximação regres-
siva com a solução dada por diferenças centrais e a solução exata da ED:
1,96 1 0 0 0,2u 0,008
1 1,96 1 0 0,4u 0,016
0 1 1,96 1 0,6u 0,024
0 0 1 0,96 0,8u 0,232
x 1
x 1,2
0,2u 0,0081,96 1 0 0 0
0,4u 0,0161 1,96 1 0 0
0,6u 0,0240 1 1,96 1 0
0,8u 0,0320 0 1 1,96 1
1,0u 0,4400 0 0 2 1,96
u(0, 2) 0,4415
u(0, 4) 0,8573
u(0,6) 1, 2229
u(0,8) 1,5155
u(1,0) 1,7155
u(0, 2) 0,5423
u(0, 4) 1,0550
u(0,6) 1,5094
u(0,8) 1,8794
u(1,0) 2,1422
2sinx
u
cos1 x
Profª Celme Torres
73
Figura 7.3 – Comparação entre a solução exata e a solução por diferenças finitas com
aproximação regressiva e aproximação central.
7.1.2. Exercícios
(1) Seja o problema descrito pela EDO sujeito as condições de contorno apresentadas:
2
2
2
d H
10x 0
dx
0 x 1
H(0) 0
H(1) 0
Resolva a EDO pelo método das diferenças finitas com aproximação de diferença central.
(2) Considere a seguinte equação diferencial de transferência de calor:
Sujeitas as condições de contorno:
2
2
d T
100 0
dx
0 x 10
T(0) 0
T(10) 0
Introdução aos Métodos Numéricos
74
Encontre a solução numérica pelo método das diferenças finitas com aproximação central,
aproximação regressiva e aproximação progressiva. Compare graficamente as soluções a-
proximadas com a solução exata da EDO.
(3) Aplique o método das diferenças finitas com aproximação progressiva para encontrar
a solução aproximada da equação diferencial no intervalo. Mostre graficamente as so-
luções exata e aproximada utilizando um programa gráfico (maple ou Excel).
y' 4 2x y(1) 5 0 x 1 considere h = 0,25
(4) Aplique o método das diferenças finitas com aproximação regressiva para encontrar a
solução aproximada da equação diferencial no intervalo. Mostre graficamente as solu-
ções exata e aproximada utilizandoum programa gráfico (maple ou Excel).
y' cos(x) 1 y(0) 1 0 x 120 considere h =
(5) Aplique o método das diferenças finitas com aproximação central para encontrar a so-
lução aproximada da equação diferencial no intervalo. Mostre graficamente as solu-
ções exata e aproximada utilizando um programa gráfico (maple ou Excel).
2y' yx y y(0) 1 0 x 2 considere h =
(6) Aplique o método de diferenças finitas com aproximação central para encontrar a so-
lução aproximada das equações diferenciais nos intervalos definidos.
a) y'' 2y ' y x 0 x 1 y(0) 2 e y'(1) 0 (use h =
0,25)
b) y'' 2y ' x 0 x 6 y(0) 4 e y'(6) 2 (use h =
1,50)
Profª Celme Torres
75
8. Métodos dos Elementos Finitos
O Método dos Elementos Finitos (MEF) é uma técnica de análise numérica para obtenção
de soluções aproximadas em problemas de valor de contorno, sendo aplicados uma grande
variedade de problemas de engenharia. Uma das principais publicações sobre MEF ilustra
sua aplicação na análise de tensões em estruturas aeronáuticas, mas sua utilização esten-
deu-se rapidamente para outros campos da engenharia (Huebner, 1995).
O MEF consiste em dividir (discretizar) o domínio de soluções do contínuo em uma quan-
tidade finita de subdomínios simples, denominados elementos finitos. Esses elementos
são representados matematicamente quase sempre por expansões polinomiais, denomina-
das funções de forma do elemento. O conjunto desses elementos, que podem assumir as
mais variadas formas, é denominado malha. A solução global é obtida pelo somatório das
soluções locais de cada elemento.
Um elemento finito pode ser unidimensional, bidimensional ou tridimensional. A Figura
8.1 ilustra alguns tipos de elementos finitos.
Elemento unidimensional linear de dois nós.
Elemento triangular de três nós.
Elemento quadrilátero de quatro nós.
Elementos quadrilátero de oito nós.
Introdução aos Métodos Numéricos
76
Elemento tridimensional de oito nós
Figura 8.1 – Tipos de elementos finitos
A fundamentação do método dos elementos finitos conta com três etapas distintas: sua
concepção matemática, aplicações e análise física e sua utilização na engenharia.
A formulação do MEF requer a existência de uma equação integral, de modo que seja pos-
sível substituir a integral sobre um domínio complexo (de volume V) por um somatório de
integrais de subdomínios de geometria simples (de volume Vi).
(8.1)
Onde se pressupõe que:
(8.2)
Se for possível calcular todos as integrais dos subdomínios Vi, basta efetuar o somatório
correspondente ao segundo membro de (1) para se obter a integral estendida a todo o do-
mínio. Cada subdomínio (Vi) corresponde a um elemento finito de geometria simples (e.g.,
segmento de reta, triângulo, quadrilátero, tetraedro, paralelepípedo).
Neste módulo apresentaremos problemas de engenharia em uma dimensão, nos quais o
método dos elementos finitos pode ser aplicado. Como nesses casos só existe uma variável
independente, as equações que regem estes problemas são equações diferenciais ordinárias.
Uma vez conhecida a equação diferencial que rege o problema e discretizado o domínio
(ou geometria) em elementos finitos, nós devemos encontrar funções que se aproximem da
solução verdadeira (geralmente desconhecida) dentro de cada elemento finito. Essas fun-
ções de aproximação são também conhecidas como funções de forma.
i
n
i 1v v
fdV fdV
n
i
i 1
V V
Profª Celme Torres
77
Para exemplificar, vamos considerar um problema unidimensional no qual h(x) é a variável
dependente. Então as funções de forma são usadas para interpolar o valor da variá-
vel h entre os nós dos elementos, da seguinte maneira:
(8.3)
Onde o superscrito e significa que a função só é válida dentro de um elemento e não
em todo o domínio. são os valores que a variável h assume nos nós do elemento e n é o
número de nós de cada elemento.
Escrevendo a equação acima em forma matricial:
(8.4)
Onde é a matriz linha que contém as funções de aproxi-
mação e é o vetor que contém o valor da variável h nos nós do elemento.
A maior dificuldade que os métodos aproximados estudados apresentam, como é o caso do
método dos elementos finitos, é a escolha da função de aproximação a ser empregada em
um problema complexo. Essas funções, em geral, têm que satisfazer exatamente as condi-
ções de contorno essenciais do problema e, além disso, representar adequadamente a geo-
metria e propriedades físicas do meio. Em termos computacionais, para poder-se usar esses
conceitos de forma sistemática, empregam-se funções de aproximação locais ao invés das
globais.
Desta maneira, discretiza-se (divide-se) inicialmente o domínio contínuo em um certo nú-
mero de elementos e aplica-se uma aproximação local da variável do problema em cada
um desses elementos. Fazendo-se, posteriormente, a junção de todos os elementos, obtém-
se um sistema de equações relativo ao modelo discreto, com um número finito de incógni-
tas.
iN (x)
n
e e
i i
i 1
h (x) N (x)h
eh (x)
e
ih
e eh (x) [N]{h}
1 2 n[N] [N (x) N (x) ...N (x)]
1
2e
3
h
h
{h}
h
Introdução aos Métodos Numéricos
78
8.1.1. Sistemas discretos e funções de base
8.1.2. Funções de base
As funções de base representam, em princípio, funções que são definidas e possuem valo-
res diferentes de zero em todo o domínio. A formulação clássica não fornece um método
sistemático para definição das funções de base, existindo para cada caso uma infinidade de
funções possíveis. Levando-se em conta que a precisão e qualidade dos resultados depende
fortemente da escolha das funções de base, isto se constitui uma séria desvantagem. A situ-
ação é ainda mais complicada quando se considera domínios em 2 e 3 dimensões com ge-
ometrias complexas. Uma escolha inadequada das funções de base pode levar a um sistema
de equações mal-condicionado e de difícil solução ou até mesmo impossível de resolver
numericamente. Por estes motivos, o Método de Galerkin possui uso restrito na solução de
casos práticos.
O MEF pode ser considerado uma modificação do Método de Galerkin onde as funções de
base são definidas de forma sistemática e levando a sistemas de equações numericamente
estáveis e fáceis de resolver, evitando desta forma as dificuldades mencionadas. No MEF
são utilizadas funções de base que só tem valores diferente de zero em uma pequena parte
do domínio. Conforme será visto, a escolha deste tipo de funções simplifica consideravel-
mente a obtenção da solução aproximada do problema. Devido à sua forma particular, elas
são também conhecidas como funções piramidais. O uso deste tipo particular de funções
constitui-se, desta forma, numa das características principais do MEF, sendo que a sua efi-
cácia está ligada à esta forma particular das funções de base. Exemplos de funções de base
piramidais em uma dimensão são mostrados na Figura 7.2.
Samiria
Highlight
Samiria
Highlight
Samiria
Highlight
Samiria
Highlight
Profª Celme Torres
79
Figura 8.2 –Subdivisão do domínio e funções de base: (a) subdivisão do domínio, (b)
funções de base, (c) derivada das funções de base e (d) forma da função de aproxima-
ção.
Com o objetivo de definir conjunto de funções piramidais, como as mostradas na Figura
8.2, é necessário primeiramente dividir o intervalo [0,1] em um conjunto de n subinterva-
los, os quais por simplicidade serão todos do mesmo tamanho. Este conjunto de intervalos
será denominado genericamente de malha. Considerando, portanto, intervalos de igual ta-
manho resulta para a dimensão de cada um dos intervalos a relação, Figura 8.2(a):(8.5)
Adotando-se a nomenclatura comum no MEF, os subintervalos são denominados de ele-
mentos e os pontos comuns de cada subintervalo são denominados de nós. Assim, cada e-
lemento é limitado pelos seguintes nós:
1 0 1
h
n n
Introdução aos Métodos Numéricos
80
(8.6)
De acordo com a Figura 8.2(b), as funções de base são definidas como segue:
(8.7)
As funções lineares definidas por (7.7) são chamadas de funções de base definidas por tre-
chos. A sua característica principal é o fato de que elas são diferentes de zero apenas em
um pequena parcela do domínio, no caso considerado o domínio é o intervalo [0,1]. Por
este motivo elas são por vezes denominadas de funções de base local. Esta característica
particular simplifica enormemente a solução do problema por métodos numéricos. As deri-
vadas das funções de base são mostradas na figura 1c e dadas matematicamente por:
(8.8)
Com as definições anteriores, a solução aproximada será composta de segmentos de
retas, ou seja ela será to tipo linear por trechos, conforme mostra a figura 7.2(d). Observa-
se ainda que função de aproximação é contínua, não havendo saltos nos pontos de transi-
ção entre um elemento e outro. Esta característica é em geral desejada na solução. Haverá
entretanto uma descontinuidade na derivada da função, conforme mostra a figura 7.2(e).
Isto não causa, todavia, problemas. Caso seja necessário continuidade, também na derivada
de primeira ordem,da função de aproximação, poderão ser escolhidas funções de base de
ordem maior, por exemplo um polinômio de segunda ordem. O procedimento delineado
kx k h k 0,1,2,3,...n
k 1
k 1
k 1 k
j
k 1
k k 1
k 1
0 0 x x
x x
x x x
h
(x)
x x
x x x
h
0 x x 1
k 1
k 1 k
j
k k 1
k 1
0 0 x x
1
x x x
d (x) h
1dx
x x x
h
0 x x 1
nu (x)
Profª Celme Torres
81
não se altera ao ser alterada a função de base. As matrizes do sistema no entanto possuirão
um número maior de elementos diferentes de zero.
8.1.3. Sistemas discretos
8.1.3.1. Mola elástica
Uma mola elástica é um elemento linear discreto onde a força-deslocamento é expresso
pela relação:
(8.9)
(8.10)
(8.11)
onde F é a força (N), é o deslocamento (m) e k é a constante elástica da mola (N/m).
A relação entre as forças ) e os deslocamentos de um elemento elástico tí-
pico (Figura 7.3) pode ser considerado um sistema discreto.
Figura 7.3 –Sistema linear elástico discreto
A força aplicada no nó 1é igual a constante de deslocamento multiplicada pelo deslo-
camento relativo do nó 1 com respeito ao nó 2, :
(8.12)
Igualmente, a força aplicada no nó 2 é igual a:
F ma
2 2
2 2
dv d x d x
a m k
dt dt dt
F k
e e
1 2(F ,F
e e
1 2( , )
e
1F
e e
1 2
e e e e e1 e 1 2 e 1 e 2F k k k
Introdução aos Métodos Numéricos
82
(8.13)
Observe que a força de equilíbrio, , é automaticamente satisfeita. A equação
acima pode ser escrita em forma matricial, como:
(8.14)
A equação 7.12 é aplicada para qualquer elemento elástico onde a relação força-
deslocamento é linear.
8.1.3.2. Torção em barra circular
Outro problema que pode ser considerado como um sistema discreto é a torção de uma bar-
ra circular. Os vetores momento, T e T' designam-se momentos de torção ou momentos
torsores (torques).
Figura 7.4 – Torção em uma barra cilíndrica
Equação de equilíbrio de momentos em torno do eixo, é dada por :
(8.15)
A distribuição de tensões tangenciais é indeterminada, é por isso necessário analisar as de-
formações. A distribuição de tensões tangenciais não pode considerar-se uniforme no eixo,
ao contrário do que acontecia na barra com carregamento axial.
e e e e e2 e 2 1 e 1 e 2F k k k
e e
1 2F F 0
e e
1 1
e e e
2 2
F1 1
k
1 1 F
T dF dA
e e1 1T ,
e e
2 2T ,
he
Profª Celme Torres
83
No caso de mecânica dos sólidos, o ângulo de torção da seção circular de um corpo ci-
líndrico é descrita pelo torque T em torno do eixo, sendo expresso na forma:
(8.16)
Onde J é o momento, L é o comprimento e G é o momento de torção. A equação 6.14 pode
ser usada para descrever as relações entre os torques e os ângulos de um
elemento cilíndrico de comprimento . Escrevendo na forma matricial, temos:
(8.17)
8.1.3.3. Fluxo de um fluido em um tubo
Outro exemplo é dado para o caso de um fluido incompressível através de tubo circular. A
velocidade sendo laminar o fluxo do fluido viscoso através de um tubo circular é dada por:
(8.18)
Onde é o gradiente de pressão, d é o diâmetro do tubo e é a viscosidade do flui-
do.
GJ
T
L
e e1 2T ,T e e1 2,
eh
e ee e
1 1
e e
e 2 2
T1 1G J
1 1h T
2
x
1 dP 2r
v 1
4 dx d
dP / dx
ee 4
e
128 h
R
d
Introdução aos Métodos Numéricos
84
Figura 8.5 – Fluxo viscoso através de um tubo cilíndrico
O volume da taxa de fluxo, Q, á obtida através da integração de em relação a área toda
seção transversal do tubo. Assim a relação entre Q e o gradiente de pressão é dada
pela equação:
(8.19)
O sinal negativo indica que o fluxo ocorre na direção do gradiente de pressão negativo.
A equação 7.19 pode ser aplicada para relação entre o valor nodal da taxa de volume de
fluxo, , a pressão, , em um tubo de comprimento e diâmetro . A
taxa de volume que entra no nó 1 é dada por:
(8.20)
Igualmente a taxa de volume que entra no 2 é dada por:
(8.21)
Logo, na forma matricial, temos:
A constante é chamada de resistência do tubo.
he1 2
xv
dP / dx
4d dP
Q
128 dx
e e1 2Q ,Q e e1 2P ,P eh ed
4
e e ee
1 2 1
e
d
Q P P
128 h
4
e e ee
2 1 2
e
d
Q P P
128 h
e e4
1 1
e e
e 2 2
P Q1 1d
1 1128 h P Q
e
e 4
e
128 h
R
d
e1P
e
2P
e1Q
e
2Q
Profª Celme Torres
85
Exemplo 8.1
Seja o problema descrito pela equação e condições de contorno abaixo:
Condições de contorno:
Resolver a equação diferencial pelo método dos elementos finitos.
O domínio será discretizado em 5 elementos de mesma dimensão ( ), conectados
por meio de nós onde o valor da incógnita u será calculada.
Onde no elemento 1, considerou-se implicitamente que , condição de contorno do
problema.
A equação básica pelo Método de Galerkin para a equação diferencial, pode ser escrita
na seguinte forma discreta:
(ii)
O passo seguinte é a substituição das variações de u e , em cada elemento, no lado es-
querdo da equação (ii).
Tomando o elemento 2, temos:
2
2
d u
u x 0 (1 x 0)
dx
u 0 em x 0
du
1 em x 1
dx
h 0,2
u 0
i
i
x 15
x 1
i 1 x
du d u
(u x) u dx u
dx dx
u
Introdução aos Métodos Numéricos
86
Efetuando as integrações acima, chega-se a:
Repetindo-se as operações acima para todos os elementos e agrupando termo a termo,
chega-se à equação:
onde:
Como as variações são arbitrárias, a identidade acima só é válida se os termos que
multiplicam cada , dos dois lados da equação, forem idênticos simultaneamente. In-
troduzindo o valorde , temos o sistema global, dado por:
A solução do sistema é dada por:
Comparando graficamente as soluções:
2h
3 2 3 2 3 2 3 22
h
1 x x x x
u u u u 1 u 2 u x 1 u 2 u dx
h h h hh
2
3 2 3 2 3 3 2 2 2 3 3 2 3 2
1 h h h
u u u u u u u u u u u u 5 u 4 u
h 3 6 6
2 3 2
2 3 4 3
3 4 5 4
4 5 6 5
5 6 6 6
au bu c u
bu au bu 2c u
bu au bu 3c u
bu au bu 4c u
a 7
u u c u u
2 3
2
2
3
a 4h 12
b h 6
c 6h
iu
iu
h 0,2
2
3
4
5
6
u 0,04811,84 6,04 0 0 0
u 0,0966,04 11,84 6,04 0 0
u 0,1440 6,04 11,84 6,04 0
0,1920 0 6,04 11,84 6,04 u
1,3120 0 0 6,04 11,84 u
u(0,2) 0,5335
u(0,4) 1,0379
u(0,6) 1,4851
u(0,8) 1,8494
u(1,0) 2,1085
Profª Celme Torres
87
Introdução aos Métodos Numéricos
88
Note que a matriz do sistema é simétrica e em banda, o que permite a utilização, em com-
putador, de algoritmos eficientes de armazenamento e solução.
Em termos computacionais, para uma melhor sistematização do programa, deve-se consi-
derar um elemento genérico (ao invés de se definirem funções de interpolação específicas
para cada elemento) e introduzir as condições de contorno do problema e fazer a expansão
das matrizes dos elementos para obtenção da matriz global. Conforme apresentado no Mé-
todo de Rayleigh-Ritz.
Exemplo 8.2
A parede apresentada na figura consiste de três materiais. De um lado da parede a tempe-
ratura é de 200 ºC e do outro lado a temperatura é de 50 ºC, com um coeficiente de con-
vecção de ( ). A equação que governa o problema é dada por:
Determine a temperatura na parede.
Condições de contorno:
100w 2m K
2
2
d T
kA 0
dx
0 x L
h1 h2 h3
1
h1
2 43
Material 1 Material 2 Material 3
T = 200 C0
o T = 50 Cf
o
1
2
3
1
2
3
2
k 70W /(m K)
k 40W /(m K)
k 20W /(m K)
h 2 cm
h 2,5 cm
h 4 cm
10W /(m K)
Profª Celme Torres
89
Onde A é a área da seção transversal (que pode ser considerada igual a 1) e k é a condu-
tividade.
Solução
Para uma malha com 3 elementos lineares, considerando a matriz de rigidez dos elemen-
tos, temos:
A matriz global do sistema é dada por:
As condições de contorno e as condições de equilíbrio, são dadas por:
Omitindo a primeira linha e a primeira coluna, temos:
o
f
x L
T(0) T
dT
kA A T T 0
dx
1 1
1
2 2
2
3 3
3
1 1 1 1 3500 3500k A 70.1
[K]
1 1 1 1 3500 3500h 0,02
1 1 1 1 1600 1600k A 40.1
[K]
1 1 1 1 1600 1600h 0,025
1 1 1 1 500 500k A 20.1
[K]
1 1 1 1 500 500h 0,04
[K] U Q
1
11
1 2
2 2 1
2 3
3 2 1
3
4
2
QU3500 3500 0 0
U Q Q3500 3500 1600 1600 0
U0 1600 1600 500 500 Q Q
0 0 500 500 U Q
1
2 1
1 2
2 3
2 1
3
2 4 f
U 200
Q Q 0
Q Q 0
Q A U T
Introdução aos Métodos Numéricos
90
A solução é dada por:
ºC
ºC
ºC
ºC
A quantidade de calor por unidade de área é dada por:
W/m
2
na parte externa (lado esquerdo)
W/m
2
na parte interna (lado direito)
Exemplo 8.3
A estrutura representada na figura é unidimensional, tem quatro nós (1 a 4) e quatro bar-
ras (A a D). Cada barra tem as suas características, nomeadamente, o módulo de Young
(E), a área da secção transversal (A) e o comprimento (L). Em cada nó existe um único
grau de liberdade. Em correspondência com os quatro graus de liberdade existem quatro
deslocamentos nodais (a) e quatro forças nodais equivalentes à ação exterior (F). Cada
barra tem dois graus de liberdade (um em cada extremidade).
2
3
44
U5100 1600 0 3500 200
1600 2100 500 U 0
0 500 510 10U 500U
5
2
3
4
U5100 1600 0 70 10
1600 2100 500 U 0
0 500 510 500U
1U 200
2U 199,58
3U 198,67
4U 195,76
1
1 1 2Q 3500U 3500U 1457,6
32 4Q 10 1 U 50 1457,6
2,0 m 2,5 m 2,5 m
Dados:
E = 2500 Pa;
A = 20cm x 40cm
Forças: F1: 4 N
F2: 3 N
F3: 5 N
F4: 2 N
Profª Celme Torres
91
Usando o MEF encontre os deslocamentos nodais (a1, a2, a3 e a4) das barras da estrutura
sabendo que a Equação diferencial que rege o problema é dada por:
d du
EA 0
dx dx
Solução
# Matriz de rigidez das barras:
Barra A: A
1 1 1 1 100 100EA 2500 0,08
[K]
1 1 1 1 100 100L 2,0
Barra B: B
1 1 1 1 80 80EA 2500 0,08
[K]
1 1 1 1 80 80L 2,5
Barra C: C
1 1 1 1 80 80EA 2500 0,08
[K]
1 1 1 1 80 80L 2,5
Barra D: D
1 1 1 1 40 40EA 2500 0,08
[K]
1 1 1 1 40 40L 5,0
# Montagem da matriz de rigidez global:
RELEMBRANDO: Mecânica dos Sólidos
O módulo de Young ou módulo de elasticidade é um parâmetro mecânico que proporciona uma
medida da rigidez de um material sólido. Obtém-se da razão entre a tensão (ou pressão) exerci-
da e a deformacão unitária sofrida pelo material. Isto é,
Tensão F/ A FL
E
Deformação a / L Ax
onde (em unidades do SI):
E é o módulo de Young, medido em pascal (1 Pa ≡ 1 N/m²);
F é a força medida em Newton;
A é a secção através da qual é exercida a tensão, e mede-se em metros quadrados;
a é a deformação, medido em metros;
L é o comprimento natural medido em metros.
Introdução aos Métodos Numéricos
92
A A
11 12
A A B D B D
21 22 11 11 12 12
B B C C
21 22 11 12
D C C D
21 21 22 22
K K 0 0
K K K K K K
[K]
0 K K K K
0 K K K K
100 100 0 0 100 100 0 0
100 100 80 40 80 40 100 220 80 40
[K]
0 80 80 80 80 0 80 160 80
0 40 80 80 40 0 40 80 120
A formulação do MEF, uma vez que a relação envolvendo todos os graus de liberdade da
estrutura, é dada por:
[K]{a} {F}
# Montagem do vetor de força {F}
O vetor de força é a soma algébrica de todas as forças que atuam no nós da estrutu-
ra,sendo portanto:
1
2
3
4
F 4,0
F 3,0
{F}
F 5,0
F 2,0
# Montagem do vetor dos deslocamentos em todos os graus de liberdade da estrutura {a}
1
2
3
4
a
a
{a}
a
a
# Solução do sistema
1
2
3
4
a100 100 0 0 4,0
a100 220 80 40 3,0
a0 80 160 80 5,0
a0 40 80 120 2,0
1
2
3
4
a 0,66
a 0,70
a 0,64
a 0,64
Profª Celme Torres
93
Aplicando o Maple temos:
A matriz de rigidez global do sistema é dada por:
O vetor de carga nodal é:
Os deslocamentos são dados pelo vetor a, onde:
Introdução aos Métodos Numéricos
94
9. Métodos dos Elementos Finitos Bidimensional
A análise bidimensional de elementos finitos envolve alguns conceitos básicos, assim co-
mo foi descritopara problemas unidimensionais. A análise do problema é algumas vezes
complicada pelo fato dos problemas bidimensionais serem descritos por equações diferen-
ciais parciais sobre uma região geometricamente complexa. O contorno () de um domínio
() bidimensional é em geral uma curva. Assim, em um problema bidimensional nos não
só procuramos uma solução aproximada de um dado problema em um domínio, como tam-
bém aproximamos o domínio de uma malha de elementos finitos adequada. Conseqüente-
mente, iremos ter erros de aproximação devido a aproximação da solução, como também,
erros de discretização devido a aproximação do domínio por elementos finitos.
A malha de elementos finitos (discretização) consiste de simples elementos bidimensio-
nais, tal como, triângulos, retângulos ou quadriláteros, que permite uma única derivação
das funções de interpolação.
No desenvolvimento de um modelo de elementos finitos bidimensional, é necessário seguir
alguns passos:
Discretização do domínio em um conjunto de elementos finitos;
Formulação ponderada da equação diferencial;
Encontrar as funções de interpolações dos elementos finitos;
Desenvolvimento do modelo usado às funções de ponderação;
Agrupamento dos elementos finitos para obtenção de um sistema global de equa-
ções algébricas;
Imposição das condições de contorno;
Solução das equações;
Pós-processamento da solução graficamente.
Profª Celme Torres
95
9.1.1. Modelos de elementos finitos
A formulação que descreve o modelo de elementos finitos bidimensional é dada por:
n
e e e e
ij j i i
j 1
K u f Q
(i 1,2,3,...,n) (8.1)
onde,
e
e e e ee e
j j j je e ei i
ij 11 12 21 22 00 i j
N N N NN N
K a a a a a N N dxdy
x x y y x y
(8.2)
e
e e
i if fN dxdy
(8.3)
e
e e
i n iQ q N ds
(8.4)
Em notação matricial, a equação (8.4) fica na forma:
e e e eK u f Q (8.5)
Note que
e e
ij jiK K (i.e., [K
e
] é uma matriz simétrica de ordem n n , somente quando
12 21a a
9.1.2. Funções de forma
Se u(x, y) é a variável dependente de um determinado problema genérico bidimensional,
então as funções de forma iN (x, y) são usadas para interpolar o valor da variável u entre
os nós dos elementos da seguinte forma:
n
e e
i i
i 1
u (x, y) N (x, y)u
(8.6)
Introdução aos Métodos Numéricos
96
Onde o supercrito e significa que a função eu (x, y) só é válida dentro do elemento e não
em todo domínio; eiu são os valores que a variável u assume nos nós do elemento e n é o
número de nós de cada elemento.
A equação (8.6) pode ser escrita em forma matricial, como segue:
e eu (x, y) N u (8.7)
onde, 1 2 n[N] N (x, y) N (x, y) ... N (x, y) é a matriz linha que contém as funções de
forma, ou funções de aproximação; e e{u } é o vetor que contém os valores da variável u
nos nós do elemento, sendo dado por:
1
2e
n
u
u
{u }
u
(8.8)
9.1.3. Elemento triangular linear
Considerando o elemento triangular linear apresentado na figura 8.1.
Figura 8.1 – Elemento triangular linear
h
y
x
ui
uk
uj
i
(xi, yi)
k
(xk, yk)
j
(xj, yj)
Profª Celme Torres
97
O polinômio de interpolação é dado por:
1 2 3u x y (8.9)
O polinômio de interpolação deve obedecer as seguintes condições nodais:
i i i
j j i
k k k
u u x x ; y y
u u x x ; y y
u u x x ; y y
(8.10)
O qual pode ser resolvido para 1 2 3, e :
j k k j i k i i k j i j j i k
1
x y x y u x y x y u x y x y u
2A
(8.12)
j k i k i j i j k
2
y y u y y u y y u
2A
(8.13)
k j i i k j j i k
1
x x u x x u x x u
2A
(8.14)
Com,
i i
j j
k k
1 x y
2A 1 x y
1 x y
(8.15)
Onde A é a área do elemento triangular.
Substituindo 1 2 3, e no polinômio e rearranjando os termos:
i i j j k kh N u N u N u (8.16)
Onde, Ni, Nj e Nk são as funções de forma, dadas por:
i i i
i
j j j
j
k k k
k
a b x c y
N
2A
a b x c y
N
2A
a b x c y
N
2A
(8.17)
Os coeficientes a, b e c foram definidos nas equações (8.12), (8.13) e (8.14), sendo:
Introdução aos Métodos Numéricos
98
i j k k j
i j k
i k j
a x y x y
b y y
c x x
j k i i k
j k i
j i k
a x y x y
b y y
c x x
k i j j i
i i j
i j i
a x y x y
b y y
c x x
(8.18)
Para derivação e exemplificação do método dos elementos finitos aplicado a problemas
bidimensionais no estado permanente, será usada a seguinte equação diferencial:
x y
u u
g g Q 0
x x y y
(8.19)
onde, u é a variável desconhecida, gx e gy são coeficientes que regem o problema descrito
pela ED(constantes) nas direções x e y e Q é um termo que representa forças ou cargas a-
plicadas.
De acordo com o Método dos Resíduos Ponderados de Galerkin, para problemas bidimen-
sionais, temos:
e e
T T T
x y
A A A
u u
[N] R(x, y)dA [N] g g dA [N] QdA 0
x x y y
(8.20)
9.1.4. Teorema de Green-Gaus
Segundo o Teorema de Green-Gaus em 2D, dado um vetor q que é o produto de um escalar
por um vetor p, e uma área A, cujo contorno delimitador é C, o qual tem um vetor unitá-
rio normal n, como mostrado na Figura 8.2.
Figura 8.2 – Teorema de Green-Gaus.
Área A
Contorno C
j
i
x
y
dA
dC
q
n
Profª Celme Torres
99
A C
q dA q n dC (8.21)
onde,
i j
x y
e x yn n i n j (8.22)
sendo, i e j os vetores unitários nas direções x, y.
Mas, 𝑞 = 𝛽𝑝, portanto,
A C
p dA p ndC (8.23)
Sabendo que:
( p) p p (8.24)
Logo,
A A C
pdA pdA p ndC (7.25)
ou
A C A
pdA p ndC pdA (8.26)
A equação (8.26) define o Teorema de Green- Gauss.
Voltando a integral (8.20), o primeiro termo a direita, dado por:
e e
T
x y
A
u u
[N] g g dA
x x y y
(8.27)
Pode ser escrito na forma:
e
x
T
eA
y
u
g
x
[N] dA
u
g
y
(8.28)
Aplicando o Teorema de Green-Gaus
Introdução aos Métodos Numéricos
100
e
x
T
e
y
u
g
x
[N] p
u
g
y
(8.29)
Substituindo no teorema de Green-Gaus, temos:
e e
x xe e
T T T
x x y ye e
A C A
y y
(a)
(b)
u u
g g
u ux x
[N] dA [N] g n g n dC [N] dA
x yu u
g g
y y
(8.30)
Analisando o primeiro termo da equação (8.30), identifica-se que o termo (a) está relacio-
nado com a forças ou cargas aplicadas no contorno do domínio, assim o termo (a) da equa-
ção (8.30) fica:
T
C
[N] qdC (8.31)
Sendo, q o módulo ou valor das forças ou cargas aplicadas no contorno do elemento. Por
convenção, é adotado o valor positivo para as forças ou cargas aplicadas entrando no do-
mínio.
A integral (8.31) é válida para qualquer elemento bidimensional e pode ser determinada
desde que as funções de forma [N] estejam determinadas. Esta integral representa a con-
tribuição das condições de contorno naturais do problema.
9.1.5. Imposição das condições decontorno
Existem duas condições de contorno, as chamadas condições de contorno essenciais ou
geométricas e as condições de contorno naturais.
As condições de contorno naturais se referem à derivada da função que está sendo anali-
sada e as condições de contorno geométricas ou essenciais, são valores especificados da
função ao longo do contorno ou domínio. As condições de contorno geométricas são in-
corporadas diretamente no vetor das incógnitas {u}, como será mostrado a seguir.
Profª Celme Torres
101
9.1.6. Determinação da matriz de rigidez
Voltando a equação (8.31), temos:
eT
qC
[N] qdC f (8.32)
que é chamado de vetor de carga do elemento.
Analisando o termo (b) da equação (8.30)
e
x
T
eA
y
u
g
x
[N] dA
u
g
y
, temos:
e e
x
x
e e
y
y
u u
g
g 0x x
0 gu u
g
y y
(8.33)
Como, eh (x, y) [N(x, y)]
e e1 2 n{h } N N ... N {h } (8.34)
e
1 2 n
e e
e
1 2 n
N N Nu
...
x x x x
{h } [N]{h }
N N Nu ...
y y yy
(8.35)
Chamando,
x
y
g 0
[D]
0 g
, temos:
e
x
T T e
eA A
y
u
g
x
[N] dA [N] [D] [N]dA {f }
u
g
y
(8.36)
A integral T e
A
[N] [D] [N]dA {f }
é chamada Matriz de Rigidez, portanto:
Introdução aos Métodos Numéricos
102
e
x
T e e
n n n 1eA
y
u
g
x
[N] dA [K ] {f }
u
g
y
9.1.7. Determinação do vetor de carga nodal
Analisando o segundo termo do lado direito da equação (8.20), temos:
T
A
[N] QdA (8.38)
Essa integral representa a contribuição das forças e cargas no Vetor de Carga Nodal, lo-
go:
T eQ
A
[N] QdA {f } (8.39)
Portanto, a equação matricial (8.5) para solução de problemas bidimensionais através do
método de elementos finitos é composta por uma matriz de rigidez [K] um vetor de incóg-
nitas e o vetor de carga nodal que incorpora, em sua derivação, as condições de contorno
do problema.
e e e eq QK u f f 0 (8.40)
9.1.8. Formulação de elementos finitos
Fazendo,
e e eq Qf f f (8.41)
Temos,
e e eK u f (8.42)
O sistema apresentado na equação (8.42) é um sistema de equações lineares no qual {u
e
} é
o vetor das incógnitas, resultante da formulação do problemas transcrito por uma equação
diferencial parcial bidimensional.
Profª Celme Torres
103
Exemplo 9.1
Elementos finitos aplicados à análise bidimensional de problemas de fluxo hídrico
subterrâneo.
A equação diferencial que rege o problemabidimensional de fluxo hídricos subterrâneo em
estado permanente é dada por:
x y
h h
k k 0
x x y y
Considere um aqüífero no qual o potencial hidráulico é conhecido no contorno do domí-
nio, sendo dado por h 10 m no lado esquerdo, h 5 m no lado direito e h 0,0 m no
contorno impermeável. Sabendo que a condutividade hidráulica do aqüífero é
5
x yk k 25 10
m/s, determine o potencial hidráulico no interior do aqüífero.
Solução
O domínio foi dividido em 4 elementos e 5 nós.
Área dos elementos: A 1,0 m2
# Elemento 1 i = 1; j = 2; k = 5
1
1
1
1
x 0
y 0
1
2
1
2
x 0
y 0
1
5
1
5
x 1
y 1
1
1
b 0 1 1
c 1 2 1
2
2
b 1 0 1
c 0 1 1
5
5
b 0 0 0
c 2 0 2
Matriz de rigidez do elemento 1:
//////////////////////////////////////////
/////////////////////////////////////////
2,0 m h = 10 m h = 5 m
h = 0 m
h = 0 m
impermeável
2,0 m
y
x
5
3
2
4
1
1
2
4
3
Introdução aos Métodos Numéricos
104
1 5
12,5 0 12,5
[K ] 0 12,5 12,5 10
12,5 12,5 25
# Elemento 2 i = 1; j = 5; k = 4
1
1
b 1 2 1
c 0 1 1
5
5
b 2 0 2
c 0 0 0
4
4
b 0 1 1
c 1 0 1
Matriz de rigidez do elemento 2:
2 5
12,5 12,5 0
[K ] 12,5 25 12,5 10
0 12,5 12,5
# Elemento 3 i = 2; j = 5; k = 3
3 2 5
12,5 12,5 0
[K ] [K ] 12,5 25 12,5 10
0 12,5 12,5
# Elemento 4 i = 4; j = 5; k = 3
4 1 5
12,5 0 12,5
[K ] [K ] 0 12,5 12,5 10
12,5 12,5 25
# Montagem da matriz global do sistema
5
25 0 0 0 25
0 25 0 0 25
[K] 100 0 25 0 25
0 0 0 25 25
25 25 25 25 100
# Vetor de carga nodal e vetor das incógnitas
Profª Celme Torres
105
O fluxo especificado no contorno do domínio é q = 0, devido ao contorno impermeável. O
nó 5 é um nó interno e deve atender a equação da continuidade, sendo portanto q5 = 0,
logo o fluxo que entra, necessariamente, é igual ao fluxo que sai do sistema.
Os valores das incógnitas (h), potencial hidráulico nos nós 1, 2, 3 e 4 são conhecidos, logo
o sistema global é descrito por:
1
2
5
3
4
5 5
q1025 0 0 0 25
q00 25 0 0 25
[K] 10 0 q0 0 25 0 25
100 0 0 25 25 q
h25 25 25 25 100 q
Tomando a última linha do sistema encontramos o valor da incógnita h5 5h 5,0 m
Graficamente, temos a variação do potencial hidráulico, ou a superfície piezométrica, do
exemplo em estudo:
O mapa dos vetores de fluxo pode ser observado na figura abaixo, indicando que o movi-
mento da água ocorre do maior para o menor potencial hidráulico. No contorno superior
e inferior, como não existe fluxo (h = 0) os vetores são adjacentes ao contorno.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
4.4
4.8
5.2
5.6
6
6.4
6.8
7.2
7.6
8
8.4
8.8
9.2
9.6
10
Carga hidráulica h (m)
Introdução aos Métodos Numéricos
106
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Profª Celme Torres
107
Referências Bibliográficas
AKIM, J.E. Finite Element for Analysis and Design. Editora Academic Press. 1994.
REDDY, J.N. An Introduction to the Finite Element Analysis. Editora Mc Graw Hill.
1984.
COOK, R.D. Concepts and Aplications os the Finite Element Analysis. Editora Jonh Wiley
& Song. 1989.
ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The Finite Element Method: Its Basis
And Fundamentals.
ASSAN, A. E. Método dos Elementos Finitos: Primeiros Passos. 2ª Edição. Editora Uni-
camp. 2003.More content from this subject
3FDED8C3-7E4D-48C7-9633-113171AE4761- AE1 - Atividade de Estudo 1
- AE2 - Atividade de Estudo 2
- AE3 - Atividade de Estudo 3
- AE4 - Atividade de Estudo 4
- Planejamento Urbano: Áreas dos Lotes
- AVF - Avaliação
- AE4 - Atividade de Estudo 4_
- Introdução Algebra (28)
- SAP Certified Application Associate - SAP Enable Now C_SEN_2305 Dumps
- SAP Certified Application Associate - SAP S4HANA Sales 2022 C_TS462_2022 Dumps