Elementos de Análise
Numérica
Prof. Carlos Ruberto Fragoso Júnior
Prof. Marllus Gustavo Ferreira Passos das Neves
Solução de problemas de Engenharia

Sem computador
Formulação

Com computador
Formulação
Solução
Solução
Interpretação
Interpretação
Antes  Problemas com equações conhecidas, mas
sem condições de serem trabalhadas
Tópicos






Aproximação ou Ajuste de curvas
Integração numérica
Derivadas numéricas
Raízes de equações
Sistemas de equações lineares
Sistemas de equações não lineares
Aplicações em Recursos Hídricos

Raizes da equação de manning





Canal prismático
Canal com seção dada em tabela
Equação de remanso
Solução da equação para encontrar Dx ideal para
muskingun cunge (propagação de vazões)
Solução da propagação de reservatório usando
Newton
Aproximação ou ajustes de curvas

Três aplicações



Extrair informações de dados  problemas de previsão
de população, por exemplo
Estudo de Leis ou funções que relacionem duas
variáveis ambientais  largura do rio em função da
área da bacia de aporte; área impermeável em função
da densidade habitacional, volume em função da cota
em um reservatório,...
Achar funções mais simples de se trabalhar do que a
função proposta
Aproximação ou ajustes de curvas

Duas classes de métodos

Interpolação  consideramos os dados precisos  a
curva de ajuste coincidirá com os pontos dados
3,5
Dados
3,0
2,5
Interpolação
2,0
1,5
1,0
Quadrados mínimos
0,5
0,0
0,0

1,0
2,0
3,0
4,0
Método dos quadrados mínimos  leva-se em
consideração erros introduzidos na obtenção dos dados
Interpolação linear
A forma mais simples de interpolação é a interpolação
linear, em que dois pontos são unidos por uma linha reta
volume

Aproximação
f ( x)  f ( x0 ) 
f ( x1 )  f ( x0 )
 ( x  x0 )
x1  x0
x
cota
Interpolação quadrática
Encontra uma parábola que aproxima 3 dados
consecutivos
f ( x)  b0  b1  x  x0   b2  x  x0  ( x  x1 )
volume

Aproximação
x
cota
Interpolação

Aproximação
De forma geral



Temos n+1 pontos (x0,y0), ..., (xn, yn), onde x0 ≠ x1 ≠ ... ≠ xn
Conhecemos y0 = f(x0), ..., yn = f(xn)
Gostaríamos de encontrar o polinômio p(x) tal que p(x0) =
f(x0), ..., p(xn) = f(xn)  polinômio interpolador
3,5
9 pontos
3,0
2,5
A função f é
conhecida em todos
eles
2,0
1,5
1,0
0,5
0,0
0,0
1,0
2,0
3,0
4,0
Interpolação

Aproximação
O que a matemática garante?
3,5
3,0
2,5
2,0
1,5
1,0
0,5
0,0
0,0
1,0
2,0
3,0
4,0
Interpolação


Aproximação
O que a matemática garante?
Seja f(x) uma função conhecida nos n+1 pontos
distintos x0, x1, x2,..., xn. Existe um único polinômio
p(x), de grau menor ou igual a n, tal que
p(xi) = f(xi) para i = 0, 1, 2, ..., n
Parábola (n = 2)  mínimo 3 pontos
Polinômio de grau 8  mínimo de 9 pontos
Splines
Aproximação

Funções polinomiais “por partes”  Splines

As “partes” fazem parte de uma partição devido aos
pontos interpolados

Ao se escolher, por exemplo, Splines cúbicos
(ordem 3), fazemos uma “colagem” de polinômios
de grau 3 em cada subintervalo do intervalo que
caracteriza a partição
Splines
Interpolação numérica
Splines

Interpolação numérica
Alguns softwares de planilha usam splines cúbicos
para suavizar linhas de gráficos
48
46
44
42
40
38
36
34
32
30
0

20
40
60
80
100
120
140
Existem rotinas prontas em praticamente qualquer
linguagem para interpolação com polinômios e splines
 Calculadora, Matlab, Excel, etc…
Splines

Interpolação numérica
Os splines cúbicos podem causar alguns problemas.
Quadrados mínimos



Em alguns casos é necessário gerar funções que aproximam
razoavelmente um conjunto de dados.
Ao contrário da interpolação, no ajuste não é necessário
respeitar todos os pontos.
A idéia é minimizar os erros com uma função simples.
3,5
3,0
2,5
2,0
1,5
1,0
0,5
0,0
0,0
1,0
2,0
3,0
4,0
Quadrados mínimos
Ajuste – exemplo em simulação

Relação entre largura de um rio e área de drenagem obtida a
partir de seções transversais em locais de postos
fluviométricos da ANA
300
250
Utilizada para
calcular os
parâmetros do
modelo Muskingum
Cunge em locais sem
dados
200
Largura do rio (m)

150
100
50
0
Brio  3.2466 A
0.4106
bacia
0
5000
10000
15000
Área da bacia (km2)
20000
25000
30000
Quadrados mínimos
Ajuste – exemplo em simulação

Curva chave de um posto pluviométrico é um ajuste de uma
equação pré-determinada aos dados de medição de vazão.
Integração numérica

Quando utilizar?

quando é necessário obter informações de área
molhada e raio hidráulico de uma seção transversal
de um rio, definida por pares de pontos x e y

Também surgem quando é necessário discretizar
uma função analítica contínua, de forma que sua
área seja mantida
Integração numérica
48
46
44
42
40
38
36
34
32
30
0

20
40
60
80
100
120
140
Idéia básica da integração numérica  aproximação da
função por um polinômio
Matlab

Interpolação 1D  função interp1


Métodos: 'nearest' - vizinho mais próximo, 'linear‘,
'spline' - spline cúbico ....
yi = interp1(x,Y,xi)
Integração numérica

Procura-se desenvolver fórmulas de integração do
tipo:
Pesos da fórmula de integração
b
n
 f ( x)dx   w  f ( x )
a
i 0
i
i
a  x0  x1  ...  xn  b
Pontos de integração
Integração numérica


O uso desta técnica decorre do fato de:

por vezes, f(x) ser uma função muito difícil de
integrar, contrariamente a um polinômio;

a única informação sobre f(x) ser um conjunto de
pares ordenados
Fórmulas de Newton-Cotes.



Regra do Trapézio simples, x0=a e xn=b;
Regra do Trapézio composta, x0=a e xn=b;
Regra de Simpson , x0=a e xn=b.
Integração numérica

Fórmulas de Newton-Cotes
 Usam pontos de integração igualmente
espaçados  (a,b) intervalo de integração
b  a xn  x0
h

n
n

Usa-se um polinômio de grau n, escrito pela
fórmula de Lagrange, que interpola os (n+1)
pontos [xi, f(xi)] i = 0, 1, 2, …, n
Integração numérica
b

Fórmulas de Newton-Cotes
n
 f ( x)dx   w  f ( x )
a
i 0
i
i
( x  x0 )  ... ( x  xi 1 )  ( x  xi 1 )  ... ( x  xn )
wi   li ( x)dx  
dx
( x1  x0 )  ... ( xi  xi 1 )  ( xi  xi 1 )  ... ( xi  xn )
a
a
b
b
n = 1  fórmula w  ( x  x0 ) dx
0
a ( x0  x1 )
dos trapézios
b
x1

x0
h
f ( x)dx   f ( x0 )  f ( x1 )
2
( x  x0 )
w1  
dx
( x1  x0 )
a
b
Regra do trapézio simples
f(x)
x1

x0
h
f ( x)dx   f ( x0 )  f ( x1 )
2
f(x1)
f(x0)
I  b  a  
f  a   f b 
2
x0
x1
x
Aproxima a área sob a curva pela área de um trapézio
Regra do trapézio simples

Intervalo [a, b] relativamente pequeno


aproximação do valor da integral é aceitável
Intervalo [a, b] de grande amplitude

aproximação inadequada  pode-se subdividí-lo
em n sub-intervalos, e em cada um a função é
aproximada por uma função linear
Fórmulas compostas ou fórmulas repetidas
Regra do trapézio composta


Intervalo [a, b] de grande amplitude
Soma da área de n trapézios, cada qual
definido pelo seu sub-intervalo
b  a xn  x0
h

n
n
Subintervalos de igual
comprimento h
Regra do trapézio composta

Fórmula:
xm

x0
f ( x)dx 
h
 f ( x0 )  f ( x1 )  h  f ( x1 )  f ( x2 )
2
2
h
 ...   f ( xN 1 )  f ( xN )
2

xN

x0
Só os termos f(x0) e f(xn) não se repetem,
assim:
h
f ( x)dx   f ( x0 )  2 f ( x1 )  f ( x2 )  ...  f ( xN 1 )  f ( xN )
2
Regra do trapézio composta
Regra do trapézio composta
4
Exemplo: Estimar o valor de
 (1  x )
2 1 / 2
dx
0



Regra do Trapézio Simples: 2 pontos (x0=0,0 e
x1=4,0)
I=h/2*(y0+y1)=2x(1,00000+0,24254) = 2,48508
Regra do Trapézio Composta: 3 pontos (x0=0,0,x1
=2,0,x2 =4,0)
I=h/2(y0+2y1+y2)=1x(1,00000+2x0,44722+ 0,24254) =
2,1369
Regra do Trapézio Composta: 9 pontos
I=(0,5/2)x(y0+2y1+2y2+2y3+2y4+2y5+2y6+2y7+y8)
=2,0936
A aproximação para 9 pontos é melhor, dado
que o valor real é 2,0947
x
y=(1+x²)-1/2
0.0
1,00000
0.5
0,89445
1.0
0,70711
1.5
0,55475
2.0
0,44722
2.5
0,37138
3.0
0,31623
3.5
0,27473
4.0
0,24254
Matlab

Regra do trapézio  função trapz


Z = trapz(Y)  calcula uma aproximação para a
integral de Y (com espaçamento unitário)
Z = trapz(Y,X)  calcula uma aproximação para
a integral de Y, definida pelos pares X, Y
4
 (1  x )
0
2 1 / 2
dx
X=0:0.5:4
Y=sqrt(1+X.^2);
Y=Y.^(-1);
Z = trapz(X,Y);
A regra utilizada é composta?
Regra do trapézio composta

Erro
ERRO!
f(x)
E=I–T


T - valor da integral
numérica.
I - valor da integral
obtida pela integração
de f(x)
f(x1)
f(x0)
x0
x1
x
Regra do trapézio composta
Erro da Regra do Trapézio Simples
(b  a)3
h3
E( f )  
f ´´( )   f ´´( ),
12
12
paraum certo   ]a, b[
Erro da Regra do Trapézio Composta
h3
Nh3 f ´´(i )
EN ( f )   f ´´(i )  
12
i 1 12
N
O erro final de uma fórmula repetida é obtido pela soma
dos erros parciais
Regra do trapézio composta
1

Exemplo: Seja
I   e dx ,
x
0
calcule uma aproximação para I usando a Regra dos Trapézios
Simples. Estime o erro cometido.
h  b  a  1 0  1
1
x1

x0
h
f ( x)dx   f ( x0 )  f ( x1 )
2
I   e x dx 
0
1

1 0
e e
2

I   e x dx  1,859141
0
12:02
Regra do trapézio composta
Estimativa do erro cometido:
ETR
(1)3 

e ,   (0,1)
12
Portanto :
ETR
1
 máx e x  0,226523
12 x[ 0 ,1]
e1  máx e x
x[ 0 ,1]
12:03
Integração numérica
n = 2  fórmula de Simpson
( x  x0 )  ... ( x  xi 1 )  ( x  xi 1 )  ... ( x  xn )
wi   li ( x)dx  
dx
( x1  x0 )  ... ( xi  xi 1 )  ( xi  xi 1 )  ... ( xi  xn )
a
a
b
b
( x  x1 )  ( x  x2 )
h
dx 
( x0  x1 )  ( x0  x2 )
3
a
b
w0  
( x  x0 )  ( x  x2 )
4h
dx 
( x1  x0 )  ( x1  x2 )
3
a
b
w1  
( x  x0 )  ( x  x1 )
h
w2  
dx 
( x2  x0 )  ( x2  x1 )
3
a
b
Regra de Simpson

Fórmula
f(x)
x2

x0
h
f ( x )dx   f ( x0 )  4 f ( x1 )  f ( x2 )
3
f(x1)
f(x2)
f(x0)
x0
x1
x2
Aproxima pela área de um polinômio de grau 2
x
Regra de Simpson composta

xn

x0
Considerando n sub-intervalos (n deve ser
um número par):
h
f ( x)dx   f ( x0 )  4 f ( x1 )  2 f ( x2 )   2 f ( xn  2 )  4 f ( xn 1 )  f ( xn )
3
Regra de Simpson composta
Regra de Simpson
Exemplo: Estimar o valor de



1
dx
0 1  x
Dividindo [0,1] em seis subintervalos, temos:
h=1/6
Regra de simpson
S =1/18.[1+4(6/7+2/3+6/11)+2. (3/4+3/5)+1/2] = 0,69317
Valor da integral
I = ln(2) = 0,69315
x
y=(1+x)-1
0.0
1,00000
1/6
6/7
2/6
3/4
3/6
2/3
4/6
3/5
5/6
6/11
1
1/2
Regra de Simpson- Erro
Erro da Regra de Simpson
5
h
IV
E ( f )    f ( ),
90
para um certo  ]a, b[
Diferenciação numérica

Idéia básica da diferenciação numérica
Aproximar a derivada real em um ponto utilizando
diferenciais pequenos.
df
Δf Δf f(x  Δx)  f(x)
 lim


dx Δx0 Δx Δx
Δx

Utilizando principalmente na discretização de
equações diferenciais
Diferenciação numérica
f
df
dx x  x1
Df f1  f0

Dx x1  x0
x0
x1
x
Diferenciação numérica  Erros de
truncamento

As derivadas numéricas são apenas uma
aproximação razoável das derivadas analíticas
df Δf

dt Δt

É possível avaliar o erro cometido nesta
aproximação utilizando as séries de Taylor
Séries de Taylor

A série de Taylor permite estimar o valor de uma
função num ponto a partir do valor da função e das
suas derivadas em um ponto próximo.
f(xi ) 2 f(xi ) 3
f(xi1 )  f(xi )  f(xi )  h 
h 
 h  ...  Rn
2!
3!



Onde h é a diferença entre xi+1 e xi.
A série de Taylor é infinita.
A aproximação da derivada numérica é finita
Séries de Taylor

O resto
f(xi ) 2 f(xi ) 3
f(xi1 )  f(xi )  f(xi )  h 
h 
 h  ...  Rn
2!
3!
O resto é dado por
f n1 ( ) n1
Rn 
h
(n  1)!
onde fn+1 é a derivada de ordem n+1 e
xi+1 e xi

é um valor entre
Séries de Taylor e derivadas
f ( xi ) 2 f ( xi ) 3
f ( xi 1 )  f ( xi )  f ( xi )  h 
h 
 h  ...  Rn
2!
3!
f ( xi 1 )  f ( xi )  f ( xi )  h  R1
f ( xi 1 )  f ( xi ) R1
f ( xi ) 

h
h
A derivada numérica tem erro de
truncamento dado por Rn/h
O valor do erro R1/h é da ordem
de h  O(h)  pode-se
expressar
f ( xi ) 
f ( xi 1 )  f ( xi )
 O ( h)
h
Erro da ordem de h  quanto menor o passo (incremento),
menor o erro da aproximação
Erros de arredondamento x truncamento


Erro de arredondamento  soma das incertezas
associadas à representação do sistema de
numeração na máquina  o computador utiliza
uma representação binária com um número finito
de bytes para representar os números reais
Erro de truncamento  aquele associado ao
truncamento de um processo infinito  como o
processo infinito não se conclui, somos forçados a
adotar uma aproximação obtida após a execução
de alguns passos
Tipos de derivadas numéricas
f(xi ) 
f(xi1 )  f(xi )
 O(h)
h
f(xi ) 
f(xi )  f(xi1 )
 O(h)
h
f(xi ) 
f(xi1 )  f(xi1 )
 O(h2 )
2 h
Progressiva forward
Regressiva backward
Centrada Centered
Considerando que h é pequeno, o erro de truncamento da
derivada numérica centrada é menor do que os outros.
Tipos de derivadas numéricas
Derivada segunda:
f ( xi ) 2 f ( xi ) 3
h 
 h  ...  Rn
2!
3!
f ( xi ) 2 f ( xi ) 3

f ( xi 1 )  f ( xi )  f ( xi )  h 
h 
 h  ... Rn
2!
3!
f ( xi 1 )  f ( xi )  f ( xi )  h 
f ( xi 1 )  2 f ( xi )  f ( xi 1 )
2
f ' ( xi ) 

O
(
h
)
2
h
12:16
Tipos de derivadas numéricas
progressiva
f
analítica
regressiva
x0
centrada
12:17
x1
x2
x
Exemplo derivada numérica


A celeridade cinemática de
propagação de perturbações no
escoamento é calculada por
onde c é a celeridade, Q é a
vazão e A é a área da seção
transversal
dQ
c
dA
Exemplo derivada numérica

Considerando uma seção prismática regular
h
A R  S
Q
n
2
3
1
2
dQ
c
dA
Qh Dh   Qh 
Dh
c
Ah Dh   Ah 
Dh
dQ
c
dh
dA
dh
Exemplo derivada numérica

Considerando uma seção qualquer
dQ
48
dQ
c
dA
46
44
42
40
h
38
c
dh
dA
dh
36
34
32
30
0
20
40
A R  S
Q
n
2
3
1
2
60
80
100
Tabelas de
A; R e Q em
função de h
120
140
interpolação
Qh   Qh  Dh 
Dh
c
Ah   Ah  Dh 
Dh
Raízes de equações


Recursos hídricos  surgem muitas equações
de difícil solução analítica, com termos
implícitos e não lineares
Os métodos aqui apresentados são iterativos
 estabelecemos uma expressão (função de
iteração) que, aplicada repetidas vezes, a
partir de uma aproximação inicial conhecida,
produz uma sequencia de aproximações que
convergem para a solução do problema.
Raízes de equações

Determinação da aproximação inicial para o
caso de uma única função  do cálculo
diferencial e integral:


Se y = f(x) é uma função contínua e muda de sinal no
intervalo [a,b], isto é, se f(a) . f(b) < 0  existe pelo
menos um ponto c E [a,b] tal que f(c) = 0. Se, além
disso, f’(x) não muda de sinal em [a,b]  c é a única
de f(x) neste intervalo
O ponto médio do intervalo pode ser uma
aproximação inicial  método da bisseção
Métodos numéricos para encontrar raízes
de equações




Bissecção
Falsa posição
Newton-Raphson
Secantes
f(x)
raiz
x
Raízes de equações
Método de bissecção


No método de bissecção é necessário fornecer
duas estimativas iniciais (limites do intervalo) de
valor de x que “cercam” a raiz
Dadas as duas estimativas iniciais xu e xl, uma
primeira estimativa para a raiz é dada por:
F(x)
x
Raízes de equações
xu  xl
xr 
2
Método de bissecção
Raízes de equações
Supõe-se que a raiz esteja exatamente
entre xu e xl
F(x)
xu  xl
xr 
2
Se f(xr).f(xl) negativo, então
Busca entre xr e xl
Se não, busca entre xr e xu
Busca entre xr e xu
x
Busca termina de acordo
Com critério de parada
Método de bissecção

Raízes de equações
Critérios de parada


Incremento de x menor que um dado limite
Diferença entre f(x) no ponto testado e zero é
menor do que um dado limite
Método de falsa posição
F(x)
Raízes de equações
Supõe-se que a raiz esteja onde estaria a raiz de
uma linha reta unindo os dois pontos
x
Método de falsa posição
F(x)
Raízes de equações
Supõe-se que a raiz esteja onde estaria a raiz de
uma linha reta unindo os dois pontos
x
Método de falsa posição
F(x)
Raízes de equações
Supõe-se que a raiz esteja onde estaria a raiz de
uma linha reta unindo os dois pontos
x
Método de falsa posição
F(x)
Raízes de equações
Supõe-se que a raiz esteja onde estaria a raiz de
uma linha reta unindo os dois pontos
x
Problemas dos métodos
anteriores
Raízes de equações

Bissecção e falsa posição sempre encontram
a raiz, mas podem ser demorados

Além disso, exigem que sejam dadas duas
tentativas iniciais com sinais contrários da
função
Raízes de equações
Método de Newton-Raphson

Raízes de equações
Combina duas ideias básicas muito comuns em
aproximações numéricas:
 Linearização  substituir (numa certa
vizinhança) um problema complicado por sua
aproximação linear que, por via de regra, é mais
facilmente resolvida
 Iteração  um processo iterativo, ou
aproximações sucessivas  repetição
sistemática de um certo procedimento até que
seja atingido um grau de aproximação desejado
Método de Newton-Raphson

Raízes de equações
Linearização de uma função  valor de f em x3
Quanto mais próximo eu tomo um ponto de x3, mais a
reta se aproxima da curva
Coef. angular da reta L(x) L(x 3 )  f(x 2 )
que passa em x2 :
Δx
f(x)
Este coef. angular também
é dado por f’(x2)
f(x2)
f(x3)
L(x3)
x1
x2
x3
x
Método de Newton-Raphson
Raízes de equações
L(x3 )  f(x2 )  f'(x2 )  Δx
f(x2)
L(x3 )  f(x2 )  f'(x2 )  Δx
f(x3)
f(x3 )  f(x2 )  f'(x2 )  Δx
Dx
L(x3)
x2
x3
Para que x3 seja a raiz
0  f(x2 )  f'(x2 )  Δx
L(x 3 )  f(x 2 )
f' (x 2 ) 
Δx
7:36
f(x2 )
x3  x2 f'(x2 )
Método de Newton-Raphson

Raízes de equações
Pela série de Taylor
f ( xi ) 2 f ( xi ) 3
f ( xi 1 )  f ( xi )  f ( xi )  h 
h 
 h  ...  Rn
2!
3!
se
h  xi 1  xi

Supondo que
f ( xi 1 )  f ( xi )  f ( xi )  xi 1  xi 
f ( xi 1 )  0
(xi+1 é a raiz)
f ( xi )
xi 1  xi 
f ( xi )
Método de Newton-Raphson
Raízes de equações
Supõe-se que a raiz pode ser encontrada seguindo uma
linha reta dada pela derivada da função no ponto inicial
F(x)
x
Tentativa inicial
Método de Newton-Raphson
Raízes de equações
Supõe-se que a raiz pode ser encontrada seguindo uma
linha reta dada pela derivada da função no ponto inicial
F(x)
derivada
Tentativa inicial
x
Método de Newton-Raphson
Raízes de equações
Supõe-se que a raiz pode ser encontrada seguindo uma
linha reta dada pela derivada da função no ponto inicial
F(x)
derivada
Tentativa inicial
x
Método de Newton-Raphson
Raízes de equações
Supõe-se que a raiz pode ser encontrada seguindo uma
linha reta dada pela derivada da função no ponto inicial
F(x)
x
derivada
Método de Newton-Raphson
Raízes de equações
Supõe-se que a raiz pode ser encontrada seguindo uma
linha reta dada pela derivada da função no ponto inicial
F(x)
x
Método de Newton-Raphson

Raízes de equações
Critérios de parada  quando f(xi) for
suficientemente próximo de zero ou quando a
diferença de dois iterados torna-se muito pequena
Problemas do método de
Newton-Raphson

Raízes de equações
É melhor que a primeira estimativa não esteja longe
demais da raiz
x
Problemas do método de
Newton-Raphson

Raízes de equações
É melhor que a primeira estimativa não esteja longe
demais da raiz
x
Método das Secantes
Raízes de equações

Um possível problema do método de NewtonRaphson, especialmente em recursos hídricos, é
que pode ser difícil estimar a derivada da função

Neste caso é possível utilizar uma aproximação
numérica para a derivada, gerando o método
das secantes
f ( xi 1 )  f ( xi )
f ( xi ) 
xi1  xi 
Método das Secantes
Raízes de equações
Semelhança dos triângulos abaixo
f ( xi ) 
xi 1  xi 
f ( xi 1 )  f ( xi )
xi1  xi 
f(x)
f ( xi )  xi 1  xi 
f ( xi 1 )  f ( xi )
x
Tentativa inicial
secante
Método das Secantes
f ( xi ) 
xi 1  xi 
f ( xi 1 )  f ( xi )
xi1  xi 
Raízes de equações
f(x)
f ( xi )  xi 1  xi 
f ( xi 1 )  f ( xi )
x
Tentativa inicial
secante
Método das Secantes
Raízes de equações
f(x)
f ( xi )  xi 1  xi 
xi 1  xi 
f ( xi 1 )  f ( xi )
x
Tentativa inicial
secante
Método das Secantes
Raízes de equações
Comparação de métodos
Raízes de equações

Newton-Raphson é mais rápido, seguido do
método das secantes, da falsa posição e
finalmente bissecção

Newton-Raphson e Secantes podem divergir

Secantes pode ser aplicado para funções em
que é difícil obter derivadas (comuns em
simulação hidrológica)

Mas podemos usar derivadas numéricas
Comparação de métodos
Raízes de equações
Comparação de métodos
Raízes de equações
Comparação de métodos
Raízes de equações
Comparação de métodos

Raízes de equações
Mesmo exemplo no excel  Newton
Comparação de métodos

Raízes de equações
Mesmo exemplo no excel  Newton com
derivadas numéricas
Comparação de métodos

Raízes de equações
Mesmo exemplo no excel  Secante
Matlab
Exemplo

Raízes de equações
Calcule o nível da água h se:
A R  S
Q
n
2
h
3
1
2
B
Q=15 m3/s
S=0,001 m/m
n=0,02
B=8 m
Raízes de equações
A R  S
G ( h)  Q 
n
2
3
1
2
0
Exemplo
Raízes de equações
Calcule o nível da água h se:

A R  S
Q
n
2
1
3
1
2
h
m
B
Q=15 m3/s
S=0,001 m/m
n=0,02
B=8 m
m=1,5
Raízes de equações
A R  S
G ( h)  Q 
n
2
3
1
2
0
Exemplo

Raízes de equações
Calcule a vazão de um vertedor
h


Q2

Q  C  L   h 
2
2
2h  L  g 

g=9,81 m/s2
h=20 cm
L=10 m
C=2
3
2
Exemplo

Raízes de equações
Calcule o nível h para uma dada vazão Q
Q = 15 m3/s
S = 0,001 m/m
n = 0,02
48
46
44
42
40
h
38
36
34
A R  S
G ( h)  Q 
n
2
32
30
0
20
40
60
A R  S
Q
n
2
3
80
1
2
100
120
140
3
1
2
0
Tabelas de A; R e Q em função de h
Simples busca e interpolação da tabela
Outro exemplo: balanço hídrico
de reservatório com vertedor
dV
 I O
dt
I  f (t )
O  f ( h)
V  g ( h)
Raízes de equações
Q  C  L  h  hs 
3/ 2
Equação de vertedor
Supondo um reservatório
V  200 h  30  h 0,3856

 200 h  30 h 
Dt
 200 h  30 h   I  C  L  h
dV
200 h t 1  30  h t 1

dt
200 h
t 1
 30  h t 1
0 , 3856
Dt
0 , 3856
t 0 , 3856
t
t
Raízes de equações
t 0 , 3856
t 1
 hs

3/ 2
Como tornar o termo de h no tempo t+1 explícito?
Raízes de equações

 C  L  h t  hs
2

3/ 2
Como encontrar raízes de
equações implícitas
200 h
t 1
0,3856
 30 ht 1
 200 h  30 h   I  C  L  h
t
Dt
t 0,3856
Raízes de equações
t 1
 hs

3/ 2

 C  L  ht  hs
2
Método de bissecção
Método de Newton-Raphson
Método das secantes
E se houver operação de comportas durante uma cheia?

3/ 2
Exemplo

Raízes de equações
Na aplicação do método de Muskingum-Cunge para a simulação da
propagação de vazão em rios, utilizam-se sub-trechos, cujos
comprimentos ideais podem ser encontrados resolvendo a equação
abaixo:
Q0
0 ,8
0, 2
Dx 
 0,8  c0  Dt   Dx 
B  S0  c0






Aplique considerando:
Q0=100 m3/s
c0=1,0 m/s
B = 30 m
S0=0,001 m/m
Dt = 1 hora (3600 s)
Use a equação abaixo para
a estimativa inicial
2,5  Q0
Dx 
B  S0  c0
Solver do Excel



Raízes de equações
O solver pode ser utilizado para encontrar raízes
de equações
Não está claro que método que Solver utiliza
Chute inicial deve estar relativamente próximo da
raiz
Raízes de equações
Sistemas de equações - Introdução



Problema comum em engenharia;
A utilização do método está liga a dois
condicionantes: (a) matriz de coeficientes, (b)
eficiência da solução;
Classificação:



Quanto ao tipo: (a) linear, (b) não linear;
Quanto ao tipo de solução: (a) direta (ex. Gauss),
(b) iterativa (ex. Gauss-Seidel);
Quanto à solução: (a) compatível e determinada;
(b) compatível e indeterminada; (c) incompatível.
Sistemas de equações lineares

Pode ser definido como:
a11 x1  a12 x2  a13 x3    a1n xn  b1
a21 x1  a22 x2  a23 x3    a2 n xn  b2
a31 x1  a32 x2  a33 x3    a3 n xn  b3





an1 x1  an 2 x2  an 3 x3    ann xn  bn
Sistemas de equações lineares

Em forma matricial:
A X   B
 a11
a
 21
 a31

 
an1
a12
a13
a22
a32
a23
a33


an 2
an 3
 a1n 
 a2 n 
 a3 n 


 
 ann 
Matriz do coeficientes
 x1 
x 
 2 
 x3 

 
 xn 
b1 
b 
 2 
b3 

 
bn 
Vetor das incógnitas
ou vetor solução
Vetor das
constantes
Sistemas de equações lineares

Classificação quanto à solução:

Possível e determinado → Possui uma única
solução.



Possível e indeterminado → Possui infinitas
soluções


Solução trivial → Det(A) ≠ 0 e B = 0;
Solução não trivial → Det(A) ≠ 0 e B ≠ 0
Det(A) = 0 e B = 0 ou B é múltiplo de uma coluna de A
Impossível → Não possui soluções

Det(A) = 0 e B ≠ 0 e B não é múltiplo de nenhuma
coluna de A
Soluções de sistemas de equações lineares

Método de Gauss (direto)


Método direto  fornecem a solução do sistema
após a realização de um n° finito de passos. Os
erros são basicamente de arredondamento da
máquina
Método de Gauss-Seidel (iterativo)

Métodos iterativos  baseiam-se na construção
de sequências de aproximações; em cada passo
valores calculados anteriormente são usados
para melhorar a aproximação
Método de Gauss
a11x1  a12 x2  a13 x3 ....... a1n xn  b1
~
~
~
~
a22 x2  a23 x3 ....... a2 n xn  b2
~
~
~
a33 x3 ....... a3n xn  b3
...........
~
~
ann xn  bn

Consiste em transformar a matriz A em uma
matriz triangular equivalente através das
seguintes operações:


Subtração de uma linha por outra multiplicada por
uma constante;
Formação de uma matriz diagonal superior.
Método de Gauss
Considere,
onde:
e,
2 x1  4 x2  2 x3  5

4 x1  9 x2  2 x3  7
 2 x  2 x  9 x  3
1
2
3

4  2
2
 x1 
5
 
 
A   4
9  2 , X   x2  e B  7 
x 
3
 2  2 9 
 3
 
4  2 5
2
A   4
9  2 7 
 2  2 9 3
Método de Gauss

1o passo: Definir um multiplicador para cada
linha baseado na primeira


m2 = a21/a11; m3 = a31/a11
2o passo: Subtrair o produto do multiplicador
da 2a e 3a linha pela 1a linha

a’i,j=ai,j- mi . ai-1,j , onde i = 2,3 e j = 1,2,3
Método de Gauss
O multiplicadores são: m2 = a21/a11 = 4/2 = 2 e m3 = a31/a11 = -2/2 = -1
-1)
2 x1  4 x2  2 x3  5 (x 2)

4 x1  9 x2  2 x3  7
 2 x  2 x  9 x  3
1
2
3

2 x1  4 x2  2 x3  5
 x2  2 x3  3
 2 x2  7 x3  8
((-)
Método de Gauss
Os multiplicadores são: m2 = a21/a11 = 4/2 = 2 e m3 = a31/a11 = -2/2 = -1
a21
'
a
a

a

m
a

a

a11  0
2 linha:
21
21
2 11
21
a11
a'22  a22  m2 a12  9  2  4  1
a'23  a23  m2 a13  2  2   2   2
b2'  b2  m2b1  7  2  5  3
3a linha:
a31
a  a31  m3 a11  a31 
a11  0
a11
'
31
'
a32
 a32  m3 a12  2   1  4  2
'
a33
 a33  m3 a13  9   1   2   7
b3'  b3  m3b1  3   1  5  8
Método de Gauss
Após estes passos, a matriz aumentada fica da seguinte forma:
a11
A   0
 0
a12
a13
a' 22
a'32
a' 23
a'33
b1  2 4  2 5 
b' 2   0 1 2  3
b'3  0 2 7
8 
Repentindo os passos de 1 a 3, só que agora tomando como
base a linha 2:
Método de Gauss
Calculando os novos multiplicadores: m’3 = a’32/a22=2/1=2
2 x1  4 x2  2 x3  5

 x2  2 x3  3 (x 2)


 2 x2  7 x3  8

2 x1  4 x2  2 x3  5
 x2  2 x3  3
 3 x3  14
(-)
Método de Gauss
Calculando os novos multiplicadores: m’3 = a’32/a22=2/1=2
3a linha:
'
a
''
'
'
a32
 a32
 m3' a'22  a32
 '32 a'22  0
a22
''
'
a33
 a33
 m3' a'23  7  2  2  3
b3''  b3'  m3' b2'  8  2   3  14
Após estes passos, a matriz aumentada agora tem a seguinte forma:
a11
A   0
 0
a12
a13
a'21
0
a'23
''
a33
b1  2 4  2 5 
b2'   0 1 2  3
b3''  0 0 3 14 
Método de Gauss
Equivalente a:
2 x1  4 x2  2 x3  5

 x2  2 x3  3


 3 x3  14

Resolvendo o novo sistema, obtem-se:
 x3  4 ,67

 x2  12,33
 x  31,83
 1
Método de Gauss


Método diretos como o de Gauss tem a
vantagem de fornecer a solução após um n°
finito de passos e não dependem de
condições de convergência
Podem ser inviáveis quando o sistema é
muito grande ou mal condicionado
Método iterativo de Gauss-Seidel



É um dos métodos mais comum e simples de
ser programado;
O método converge somente sob certas
condições e normalmente conduz a um
número maior de operações quando
comparado com métodos diretos
Como qualquer método iterativo 
convenientes para sistemas grandes e
esparsos que aparecem após discretização de
EDPs
Método iterativo de Gauss-Seidel
A equação utilizada para iterações é a seguinte:
k 1
i
x
1

ai ,i
i 1
n

k 1
 bi   ai , j  x j   ai , j  x kj

j 1
j i 1





Pode-se utilizar um coeficiente para acelerar o processo
de convergência:
k 1
i
x

 
i 1
bi   ai , j  x

ai ,i 
j 1
k 1
j

  ai , j  x 
j i 1

n
k
j
Método iterativo de Gauss-Seidel
Seja o sistema de equações:
a11 x1  a12 x2  a13 x3    a1n xn  b1
a21 x1  a22 x2  a23 x3    a2 n xn  b2
a31 x1  a32 x2  a33 x3    a3 n xn  b3





an1 x1  an 2 x2  an 3 x3    ann xn  bn
Método iterativo de Gauss-Seidel
Obtemos o valor de x1 a partir da primeira equação, o valor de x2 a
partir da segunda equação e assim sucessivamente:
 

x1k 1 
1
b1  a12 x2k  a13 x3k  a14 x4k    a1n xnk
a11
x2k 1 
1
b2  a21 x1k 1  a23 x3k  a24 x4k    a2 n xnk
a22
x3k 1 
1
b3  a31 x1k 1  a32 x2k 1  a34 x4k    a3 n xnk
a33

xnk 1
 
 






1

bn  an1 x1k 1  an 2 x2k 1  an 3 x3k 1    an n 1 xnk11
ann
 

Método iterativo de Gauss-Seidel

Ponto de partida


Conjunto de valores iniciais  na falta de melhores
informações, podemos usar x1 = x2 = ... = xn = 0
Critério de parada


Número de iterações excedeu um determinado
valor m;
A seguinte condição atenta uma precisão adotada:
n
n
 b  a
j 1
i
i 1
i, j
xj  
Método iterativo de Gauss-Seidel

Convergência do método:


Existe um critério de convergência, através de um
teorema, que envolve autovalores de matrizes, o
que nem sempre é trivial
Este teorema, no entanto, permite estabelecer
outras condições de convergência de verificação
mais simples


O método converge se a matriz A é diagonalmente
dominante
O método converge se a matriz A é uma matriz positiva
definida
Método iterativo de Gauss-Seidel

Convergência do método:

matriz de coeficientes seja positiva definida

Inspeção da diagonal principal (necessária):
ai , j  0 , para i  j

Domínio da diagonal (suficiente):
ai ,i   ai , j e ai ,i   a j ,i
j i

j i
Método dos menores principais (necessária e suficiente):
Ri  0
Método iterativo de Gauss-Seidel
Considere
2 x1  4 x2  2 x3  5

4 x1  9 x2  2 x3  7
 2 x  2 x  9 x  3
1
2
3

k 1
1
x
Aplicando o método, tem-se
x
k 1
2
x3k 1


1
 5  4 x2k  2 x3k
2
1
k 1
k
 7  4 x1  2 x3
9
1
 3  2 x1k 1  2 x2k 1
9




Método iterativo de Gauss-Seidel
Considerando o ponto de partida com Xk=(x1, x2, x3)=(0, 0, 0),
a primeira iteração fica: k 1 1
x1
x2k 1
x3k 1
Adotando ɛ = 0.0001, após
244 iterações a solução
converge para:

5  4  0  2  0  2 ,5
2
1
 7  4  2 ,5  2  0   0 ,333
9
1
 3  2  2 ,5  2   0 ,333  0 ,815
9
 x1  31,83

 x2  12,33
 x  4 ,67
 3
Método iterativo de Gauss-Seidel
Exercício para casa:
- Desenvolver um algoritmo para resolução de sistemas lineares pelo
método iterativo de Gauss-Seidel.
Sistemas de equações não lineares

Pode ser definido como:
f 1  x1 , x2 , x3 ,...,xn   0
f 2  x1 , x2 , x3 ,...,xn   0
f 3  x1 , x2 , x3 ,...,xn   0

f n  x1 , x2 , x3 ,...,xn   0
onde f é uma função não linear em função de x1,x2,…,xn.
Sistemas de equações não lineares

Método iterativo de Newton


Se baseia no método Newton-Rapson para
solução de equações não lineares
Transforma o sistema não linear em um sistema
linear (linearização), este resolvido a cada uma
das várias iterações de modo que a solução do
linear se aproxime daquela esperada (não linear)
Método iterativo de Newton

Um sistema de equações não lineares:
f 1  x1 , x2 , x3 ,...,xn   0
f 2  x1 , x2 , x3 ,...,xn   0
f 3  x1 , x2 , x3 ,...,xn   0

f n  x1 , x2 , x3 ,...,xn   0
pode ser expandido para série de Taylor de primeira ordem:
Método iterativo de Newton

Resultando em um sistema de equações lineares:
f 1 x1 , x2 , x3 ,...,xn 
k 1
f 2  x1 , x2 , x3 ,...,xn 
k 1
f 3 x1 , x2 , x3 ,...,xn 
k 1
 f 1 x1 , x2 , x3 ,...,xn  
k
 f 2 x1 , x2 , x3 ,...,xn  
k
 f 3 x1 , x2 , x3 ,...,xn 
k
k
k
k
k
k
k
k
k
k
k
k
k
f 1
f
f
Dx1  1 Dx2    1 Dxn  0
x1
x2
xn
f 2
f
f
Dx1  2 Dx2    2 Dxn  0
x1
x2
xn
f
f
f
 3 Dx1  3 Dx2    3 Dxn  0
x1
x2
xn

f n  x1 , x2 , x3 ,...,xn 
k 1
 f n x1 , x2 , x3 ,...,xn  
onde Δxi = xik+1- xik
k
f n
f
f
Dx1  n Dx2    n Dxn  0
x1
x2
xn
Método iterativo de Newton

Em forma matricial:
J   X 
k
 f 1
 x
 1
 f 2
 x1
 f
 3
 x1
 
 f n
 x
 1
f 1
x2
f 2
x2
f 3
x2

f 2
x2
f 1
x3
f 2
x3
f 3
x3

f 3
x2
f 1 
xn 

f 2 

xn 
f 3 


xn 

 
f n 

xn 

Jacobiano (k)
k 1
 B
 x1 
x 
 2 
 x3 

 
 xn 
Vetor das incógnitas
ou vetor solução (k+1)
k












f 1
f
f
x1  1 x2    1
x1
x2
xn
f
f
f
f 2  2 x1  2 x2    2
x1
x2
xn
f
f
f
f 3  3 x1  3 x2    3
x1
x2
xn

f
f
f
f n  n x1  n x2    n
x1
x2
xn
f1 
Vetor das
Constantes (k)

xn 

xn 



xn 


xn 

Método iterativo de Newton

Ponto de partida


Conjunto de valores iniciais
Critério de parada


Número de iterações excedeu um determinado
valor m;
Verifique se a seguinte condição atenda uma
precisão adotada:
n

i 1
f i k 1  f i k  
Método iterativo de Newton

Convergência do método:

É necessário que a matriz de coeficientes seja
positiva definida

Inspeção da diagonal principal (necessária):
ai , j  0 , para i  j

Domínio da diagonal (suficiente):
ai ,i   ai , j e ai ,i   a j ,i
j i

j i
Método dos menores principais (necessária e suficiente):
Ri  0
Método iterativo de Newton
Exercício para casa:
- Desenvolver um algoritmo para resolução de sistemas
não lineares pelo método iterativo de Newton.
Trabalho




7:36
Desenvolver uma iteração, manualmente, do
sistema não linear resultante do chamado
problema dos três reservatórios a seguir
Verifique as condições de convergência
(matriz A diagonalmente dominante e
positiva definida)
Utilize o método de Gauss-Sidel após a
linearização do sistema não-linear
Prazo: 1 semana após esta aula
Trabalho

o sistema abaixo é composto por 3 reservatórios.
Não se sabe quais os valores de vazão nos trechos
nem a cota piezométrica (CP) no ponto de
convergência dos trechos. Para Determiná-los,
desprezando as perdas de carga localizadas e as
cargas cinéticas
trecho
AD
DB
DC
100m
90m
A
B
D
7:36
80m
C
L(m) D(mm)
f
300
400 0,03
300
400 0,03
900
500 0,02
Trabalho



Para resolver este problema, faz-se a
hipótese de que a CPD = 90, o que equivale
a dizer que QDB = 0
Depois testa-se a hipótese
Do resultado do teste, ou o problema
acaba ou se monta um sistema de
equações não-lineares com 4 incógnitas.
A seguir o resumo do processo
7:36
Trabalho
Completando a tabela  b=0,0826f
trecho L(m) D(m) f
b
AD
300 400 0,03 0,00248
DB
300 400 0,03 0,00248
100m
DC
90m
A
B
D
80m
C
900 500 0,02 0,00165
Trabalho
Hipótese  CPD= 90m  QDB=0 m3/s
Calcular QAD e QDC. Por exemplo,
100m
QAD
90m
A
B
D
80m
C
 ΔHADDAD
 
β
L
AD
AD

m




1
n
Trabalho
Hipótese  CPD= 90m  QDB=0 m3/s
trecho
AD
DB
DC
100m
90m
A
B
D
DH Q (m3/s)
10
0,37
0
0,00
10
0,46
80m
C
QAD< QDC
QDB≠ 0
Trabalho
100 - CPD  72,62Q
Sistema de
equações
2
AD
2
DB
2
DC
90 - CPD  72,62Q
CPD  80  47,59Q
QAD  QDB  QDC
Resultado  CPD=89,63m,
QAD=0,38m3/s, QDB=0,07m3/s e
QDC=0,45m3/s
Trabalho
Resultado
100m
A
0,38 m3/s
D
90m
89,63m
B
80m
0,07 m3/s
0,45m3/s
C
Trabalho

Para facilitar, chame:




CPD de x1
QAD de x2
QDB de x3
QDC de x4
100 - x1  72,62x
2
2
2
3
2
4
90 - x1  72,62x
x1  80  47,59x
x2  x3  x4
7:36