matematica computacional

Matemática Computacional, MEMec, LEAN, MEAer

Representação de números - Conversão de base b para base 10

Números em base 10

0, 1, 2, … , 8, 9

( )1 1 10 02 1... , ...n nN d d d d d d− −−=

unidadesdezenas

centenas10 dígitos que

constituem a base

Valor depende da posição dos dígitos

( ) 1 010

3 24579 4000 500 70 9 4 10 5 10 7 10 9 10N = = + + + = × + × + × + ×Exemplo1:

Generalizando

( ) 2 110

0 1 2684,75 6 10 8 10 4 10 7 10 5 10N − −= = × + × + × + × + ×Exemplo2:

1 21

1 1 01 1 0 210 10 ...10 10 ... 10 10n n

n nd d d dd d − −−

−−−= × + × + + × + + × + × +×

0 9 10 1id =≤ −≤

, 10, 11, 12, … , 19 , 100, 101, … , 472, … , 999, 1000, 1001, …, 20, 21, … , 99


Representação de números - Conversão de base b para base 10

Números em base b ≠10, por exemplo base 3 (ou seja b=3)

0, 1, 2

( )1 311 0 2.. , .. ..n nN d d d d d d− −−=

3 dígitos da base

Generalizando

1 1 01 1

1 21 203 3 3 3... ..3 .3n n

n n d dd d d d−−

− −− −= × + × + + × + + × + × +× 20 3 1id≤ ≤ = −

3

4

5

6

7

8

9

10

11

Exemplo:4 2 0

3 12

1

3( 2 0 1 0 2 , 2 1 )N↑ ↑ ↑ ↑

−↑−

↑ ↑

= 3 1 104 2 20 3 0 32 3 1 3 2 2 33 1 3− −= × + × + + ×× ++ ×× + ×

1 12 81 1 9 2 2 13 9

173,777.

2 0 3

.

0

.

7= × + × + + × + ×

=

+ × + ×

, 10, 11, 12 , 100, 101, 102, 110, 111, 112, 120, 121, 122, 200, …, 20, 21, 22


Representação de números - Conversão de base 10 para base b

Números inteiros

Exemplo:3 02 1

023 1

4( 1 3 3 2 )d dd d

N↑↑

↑ ↑

=

2 10

3 011 4 3 4 3 4 2 4 126 (126)== × + × + × + × =

126 4

2 31

o resto da divisão inteira de 126 por 4 é o dígito da posição d0 do número 126 escrito em base 4

Dividindo 126 por 4 resulta,

0

3 2 1 0126 1 4 3 4 3 4 2 44 4

d

× + × + × + ×=

Ou seja,

rest

2 1

o

3

0

1

21 4 3 4 3 44

= × + × + × +

0resto

2314

d=

= +



31 4

3 7

o resto da divisão inteira de 31 por 4 é o dígito da posição d1 do número 126 escrito em base 4

1

2 1 031 1 4 3 4 3 44 4

d

× + × + ×=

Ou seja,

Dividindo o resultado da anterior divisão (= 31) por 4 resulta,

Ou seja, efectuando divisões sucessivas por 4, os restos das divisões vão ser os dígitos do número escrito em base 4

126 431 4

7 41 4

23

31 0 4

0 0

4 4126 ( 1 3 30 ... 0 2 ) (1 3 3 2 )= =

1resto

374

d=

= +7

r

0

e

1

sto

31 4 3 44

= × + × +



Números fraccionários puros

Exemplo:1 2

24

1

( , 1 3 )d d

x− −

↑−↑

−

=

1 2101 4 3 4 0,4 (0375 ,4375)− −= × + × = =

Multiplicando 0,4375 por 4 resulta, 1 2

1 20,4375 4 ( 1 4 3 4 ) 4d d− −

− −× = × + × ×

Retirando a parte inteira ao resultado anterior e multiplicando novamente por 4 resulta,

22

10,75 4 ( 3 4 ) 4 3

dd −−

−× = × × =

Ou seja, efectuando multiplicações sucessivas por 4, as partes inteiras vão ser os dígitos do número escrito em base 4

0,4375 x 4 = 1,750,75 x 4 = 3

13 4 40,4375 ( ,1 3 ) ( 0, )0 1 3= =

1 1

0,75

11 3 4 1 ,75

d d− −

−= + × =


Sistema de ponto flutuante

Exemplo: FP(10,4,2,A) e FP(10,4,2,T)

↑ ↑↑↑

=

2 dígit4 dígitos na mantissa

base 10 A=ArredondamentoT=Truncat

os no exp uraoente

( , , ,_) ( 10 , ,4 _2 , )pP FqF b P

Formato normalizado – com excepção da representação do número zero, d‒1≠0

−±

− − −±= ± × → = ± ×

2 dígitos(base 10)

1 01 2 3 4

4 dígitos (da base

( )

10)

(0, ) 10 t ttx d dx dbm d

−≤ ≤ → ≤ < → ≤ <1pelo que 0,1000 0,9999 0,1 1 1m m b m



Exemplo: Representar x = 805,174 em FP(10,4,2,T) e em FP(10,4,2,A)

3805,1 0,805174 74 10x = = × 03(10,4,2, ) : ( ) (0,8051) 10FP fl xT x += = + ×

03(10,4,2, ): ( ) (0,8052) 10FP fl xA x += = + ×

805,174

805,0 805,1 805,2 805,3

805,15 805,25



Exemplo: Representar x = 805,174 em FP(10,4,2,T) e em FP(10,4,2,A)

3805,1 0,805174 74 10x = = × 03(10,4,2, ) : ( ) (0,8051) 10FP fl xT x += = + ×

03(10,4,2, ): ( ) (0,8052) 10FP fl xA x += = + ×

Erro absoluto: E x x= −

3 3 3(10,4,2, ) : 0,80 74 751 10 0,805 41 10 0,0000 10 0,074FP E x xT = − = × − × = × =

3 3 374(10,4,2, ): 0,8051 10 26 20,8052 10 0,0000 10 0, 60FP E x xA = − = × − × = − × = −

Erro relativo: x x Eex x−= =

50, 70(10,4,2, ) : 9,2 108

47405,1

EFPx

T e −= = = ×

50, 20(10,4,2, ): 3,2 108

67405,1

EF AP ex

−−= = = − ×

39,2 10 %−×

33,2 10 %−− ×

x100

x100


Unidade de arredondamento, u

Com truncatura, FP(b,p,q,T)

Unidade de arredondamento, u – majorante do erro relativo (possível de ser cometido) narepresentação dum número

− × − × −= = =× t t

t

x x m b m b m mex m b m

1 1min( ): (0,10...0) 1 1 min( )b b m bm m m− −≤ < → ≤ < → =

− − −→ = = ≤ < max( )max ( ) max max

min( )x x m m m me ux m m

111

max( )min( )

p ppb b b

bbm m

u um −−

− −−− ≤

→ = ==

− max( ):m m

(0,10... 0)bp

↑−

(0,1...0 1)bp

↑−

u

Com truncatura, a unidade de arredondamento é a maior distância relativa entre dois números consecutivos (representados nesse sistema)

( 1) ( 2)...(0,10... 0 )

(0,10... 0)

bp

b

p

p

pdm

m

d↑

−

↑

− + − +

−

= =

( 1) ( 2)...(0,0... 0 ) (0,0...0 1)p pp

bpp

bd dm bm − + − +−

↑↑−−

− = < =


Unidade de arredondamento, uUnidade de arredondamento, u – majorante do erro relativo (possível de ser cometido) narepresentação dum número

− × − × −= = =× t t

t

x x m b m b m mex m b m

− − −→ = = ≤ < max( )max ( ) max max

min( )x x m m m me ux m m

Com arredondamento, FP(b,p,q,A) – a unidade de arredondamento é metade do valor da unidade de arredondamento com truncatura

112

pu b −=

(0,10... 0)bp

↑−

(0,1...0 1)bp

↑−

uCom arredondamento, a unidade de arredondamento é metade da maior distância relativa entre dois números consecutivos (representados nesse sistema)


Norma IEEE754 – formatos simples e duplo (base 2)

Formato simples 32 bits = 4 bytes S Expoente Mantissabits -> 1 8 23

Formato duplo 64 bits = 8 bytes S Expoente Mantissabits -> 1 11 52

Número normalizado ou desnormalizado

• Número normalizado – no caso do expoente não ser nem todo “zeros” nem todo “uns”

• Número desnormalizado• Se o expoente for todo “zeros” – representação do número zero

ou representação de underflow

• Se o expoente for todo “uns” – representação de overflow (infinito ou NaN)


Norma IEEE754 – formato simples


Número normalizado24 bits

22127

1 3( 1) (1, ) 2S ed dx −−

−= − × ×

Exemplo: Valor correspondente ao conjunto de bits (em formato simples IEEE754)1 01011001 0010100…0

O número é normalizado (porque o expoente armazenado não é nem todo “zeros” nem todo “uns”) 24 bits

22127

1 3( 1) (1, ) 2S ed dx −−

−= − × ×

6 4 32

3 52

1

(01011001) 2 2 2 1 89

( ,001010...0) 2 2 1,11 51 625

s

e

m − −

== = + + + == = + + =

127 31 8 8 129( 1) 2 1,15625 2 4,20641... 11,15625 0x − − −= − × × = − × = − ×




Números desnormalizados

Exemplo: Valor correspondente ao conjunto de bits (em formato simples IEEE754)a) 0 00000000 0010100…0 b) 1 00000000 00…000

Expoente todo “zeros”24

1

bit

12623 2

s

( 1) ( , ) 20S d dx −− −= − × ×

3 52

0

( ,001010...0) 2 2 0,156250

s

m − −

== = + =

120 396( 1) 20 1,8367... 10,15625x − −= − × × = ×

(i) Expoente todo “zeros” – representação do número zero ou de underflow gradual 24

1

bit

12623 2

s

( 1) ( , ) 20S d dx −− −= − × ×

2

1( ,00...0) 00

sm

== =

261 1( 1) 2 00x −= − × × = −

a)

b)




Números desnormalizados

Exemplo: Valor correspondente ao conjunto de bits (em formato simples IEEE754)a) 0 11111111 0010100…0 b) 1 11111111 00…000

Expoente todo “uns” – representação de overflow

(ii) Expoente todo “uns” – representação de overflow

a)

b)

0 NaNm ⎯⎯→≠

0 I fm n⎯⎯→ ±= Como 1, o resultado é s Inf= −




Formato normalizado −− −= − × ×

1||

2

1270

4 bits

1 23 2( 1) ( , ) 2S ex d d d

≤ ≤ ⇔ ≤≤ −≤ ⇔ ≤ −2 2Expoente: (00000001) (1111111 1260) 1 2 1 754 127 2ee e

− −× = × = − × × 254 127 127 23 127 128 382 2Limite de overflow: (1,11 11) 2 (1,11 11) 2 (2 2 ) 2 2 3,4 10

− − − −× = × = × 1 127 126 126 382 2Limite de underflow: (1,00 00) 2 (1,00 00) 2 2 1,2 10

− − − − −

↑−

× = × = × 126 23 126 149 452

23

Limite de underflow gradual: (0,00 0 1 ) 2 2 2 2 1,4 10

− − − −= = = ×1 1 24 23 7Unidade de arredondamento c/ truncatura: 2 2 1,2 10pu b

− − − −= × = × = × ×1 1 24 23 71 1 1Unidade de arredondamento c/ arredondamento: 2 2 0,6 10

2 2 2pu b


Norma IEEE754 – formato duplo

64 bits = 8 bytes S Expoente Mantissabits -> 1 11 52

Formato normalizado

≤ ≤ ⇔ −⇔ ≤ ≤ − ≤≤2 2Expoente: (00000000001) (1111111 10221110) 1 20 1046 23 1023ee e

+ − + + +× = − × × 1023 52 1023 1024 3082Limite de overflow: (1,11 11) 2 (2 2 ) 2 2 1,8 10

− − − −× = × = × 1 1023 1022 1022 3082 2Limite de underflow: (1,00 00) 2 (1,00 00) 2 2 2,2 10

− − − − −

↑−

× = × = × 1022 52 1022 1074 3242

52

Limite de underflow gradual: (0,00 0 1 ) 2 2 2 2 4,9 10

− − − −= = = ×1 1 53 52 16Unidade de arredondamento c/ truncatura: 2 2 2,2 10pu b

− − − −= × = × = × ×1 1 53 52 161 1 1Unidade de arredondamento c/ arredondamento: 2 2 1,1 10

2 2 2pu b

−− −= − × × 1023

0 1 52

1||

3 bits

2

5

( 1) ( , ) 2S ex d d d

Formato duplo


Exemplo 1)

Operações elementares em ponto flutuante (FP)Passos a seguir:1) Decomposição dos operandos nas mantissas e expoentes2) No caso de soma e subtracção, alinhamento das mantissas3) Operações com mantissas e com expoentes4) Normalização da mantissa5) Arredondamento da mantissa

= + → × + ×3 1123,4 4,321 0,1234 10 0,4321 10y

Exemplos em FP(10,4,2,T)

×+ ×

×

3

3

3

0,12340,0043210,127 217

101010

= = × 3( ) 0,1277 10fl y y


Exemplo 2)


= − → × − ×3 1427,3 2,183 0,4273 10 0,2183 10y


×− ×

×

3

3

3

0,42730,0021830,425 171

101010

= = × 3( ) 0,4251 10fl y y

Nota: se não existirem dígitos de guarda

×−

30,42730,00

1021 83 ×

×3

30,42521010

= = × 3( ) 0,4252 10fl y y


Exemplo 3)


− −×= → = × = ××

11 3 2

3

8,475 0,8475 10 0,8475 10 5,478345 10154,7 0,1547 10 0,1547

y


−= = × 1( ) 0,5478 10fl y y



Notas:1) Existindo dígitos de guarda, a simulação duma operação elementar em FP corresponde a escrever o resultado obtido no formato em FP, arredondando o resultado para o número de dígitos existentes na mantissa.

2) As operações com os expoentes são operações com números inteiros pelo que não introduzem aproximações (operações exactas).

3) As operações em FP, em geral, não respeitam as propriedades comutativas, distributiva e associativa da aritmética exacta.


Erros nas operações elementares em ponto flutuante (FP)Nota: as operações com os expoentes são exactas, os erros provêm das operações com as mantissas

= + → = + ⋅ = ⋅ += = ⋅

( )( ) (1 )

fl x x Efl x x x e x eEe E x e

x

[ ]= + = ⋅ + + ⋅ + ⋅ + 1 2 1 1 2 2 3( ) (1 ) (1 ) (1 )y fl x x x e x e e

Soma:

= + + + + + +1 2 1 1 2 2 3 1 2( )y

x x e x e x e x x ϑ

arredondamentodo argumento

arredondamentodo resultado

termos de ordem superior

= +1 2y x x

= − = + + + + 1 1 2 2 3 1 2( )E y y e x e x e x x ϑ

= ≤ +2Ee e uy

ϑ

(x1 e x2 têm o mesmo sinal)

→ ≤ ⋅ + ⋅ + ⋅ + + 21 2 1 2 ( )E u x u x u x x uϑ → ≤ ⋅ + +

21 22 ( )

y

E u x x uϑ


Erros nas operações elementares em ponto flutuante (FP)Nota: as operações com os expoentes são exactas, os erros provêm das operações com as mantissas

= + → = + ⋅ = ⋅ += = ⋅

( )( ) (1 )

fl x x Efl x x x e x eEe E x e

x

[ ]= ⋅ = ⋅ + ⋅ ⋅ + ⋅ + 1 2 1 1 2 2 3( ) (1 ) (1 ) (1 )y fl x x x e x e e

Multiplicação:

= = ⋅ + + + +1 2 1 1 2 2 1 2 3 1 2...y

x x e x x e x x e x x ϑ




= ⋅1 2y x x

= − = + + + 1 1 2 2 1 2 3 1 2E y y e x x e x x e x x ϑ

Analogamente se conclui para a divisão: +≤ 3e u ϑ

→ ≤ ⋅ ⋅ +1 23y

E u x x ϑ

= ≤ +3Ee e uy

ϑ


Erros nas operações elementares em ponto flutuante (FP)

[ ]= − = ⋅ + − ⋅ + ⋅ + 1 2 1 1 2 2 3( ) (1 ) (1 ) (1 )y fl x x x e x e e

Subtracção:

= − + − + − +1 2 1 1 2 2 3 1 2( )y

x x e x e x e x x ϑ




= −1 2y x x

= − = − + − + 1 1 2 2 3 1 2( )E y y e x e x e x x ϑ

Se |x1‒x2| for “muito pequeno”, o erro relativo pode ser muito grande

-> cancelamento subtractivo

= + + + − + + +1 1 3 2 2 3(1 ) (1 )x e e x e eϑ ϑ

erro absoluto “pequeno”

(em relação à grandeza dos argumentos)

(x1 e x2 têm o mesmo sinal)

→ ≤ ⋅ + ⋅ + ⋅ − +1 1 2 2 3 1 2| | | | | ( )|E e x e x e x x ϑ

( )→ ≤ ⋅ + + 21 22 ( )E u x x uϑ( )→ ≤ ⋅ + ⋅ + ⋅ + + 2

1 2 1 2 ( )E u x u x u x x uϑ

+= = ≤

− −+1 2

1 2 1 2

2x xE Ee e u

y x x x xϑ


Exemplo de cancelamento subtractivo

2 212,46485 12,45012 0,1246 10 0,1245 10y = − → × − ×Cálculo em FP(10,4,2,A):

2

2

2

0,12460,12450,000

1010101

×−×

×

2( ) 0,0001 10fl y y= = ×


Valor exacto:

12,46485 12,45012y = −

12,46485 12,45012 0,01473y = − =

Erro absoluto: 20,01473 0,0001 10 0,01473 0,01 0,00473E y y E= − = − × = − =

Erro relativo: 100

0,00473 0,32 32%0,01473

y y Ee ey y ×

−= = = = ⎯⎯→



Soma: = +1 2y x x ϑ≤ + 22 ( )e uu(x1 e x2 têm o mesmo sinal)

Multiplicação e divisão: = ⋅ =1 2 1 2, /y x x y x x ϑ≤ + 23 ( )e uu

Subtracção: = −1 2y x x ϑ≤−

++1 2

1

2

2

| | | |2| |

( )x xex

uux


Processos que podem originar acumulação de erros

Somatório:=

=1

n

i

i

y x

ϑ+→ ≤ + 2( 1) ( )n u ue

(xi – números positivos e negativos)

Algoritmo: Inicialização: s0=0para i=1 até n fazer

si=si–1 +xifim do ciclo iy=sn

No caso de os xi possuírem o mesmo sinal é possível estimar um majorante do erro relativo

Notar que a ordem pelo qual o cálculo é efectuado não é indiferente

Para minimizar o erro, a variável “auxiliar” si pode ser declarada com precisão acrescida. Se não ocorrer cancelamento subtractivo, o erro raramente ultrapassa uma unidade de arredondamento (independentemente do valor de n)


Processos que podem originar acumulação de erros

Produto interno (de vectores):=

= ⋅ = ⋅

1

n

i i

i

s x yx y

ϑ+→ ≤ + 2( 2) ( )n u ue

Algoritmo: Inicialização: s0=0para i=1 até n fazer

si=si–1 +xi . yifim do ciclo iy=sn

No caso dos termos (xi yi) possuírem o mesmo sinal é possível encontrar um majorante do erro relativo

Tal como no caso do somatório, para minimizar o erro, a variável “auxiliar” si pode ser declarada com precisão acrescida


Número de condição

Avaliar a propagação de erros: análise directa vs. análise indirecta

Dedução de número de condição

→

− −= →− −

( ) ( ) ( ) ( )

'( ) lim '( )x x

f x f x f x f xf x f xx x x x

Análise indirecta – número de condição

xx ( )f x

( )f xsituação 1

( )f xsituação 2

perturbação de x

situação mal condicionada

situação bem condicionada

( )− − ( ) ( ) '( )f x f x f x x x

− ⋅ − ⋅

cond ( )

( ) ( ) '( )( ) ( )

xfefe x

f x f x x f x x xf x f x x


Número de condição

Ou seja,⋅= × '( )

cond ( ) , cond ( )( )f x

x f xe f x e f xf x

cond f(x) representa o factor de ampliação entre o erro relativo do argumento x e o erro dovalor da função f(x)

Se cond f(x) for grande, então uma perturbação no valor do argumento x é muito ampliada

Se cond f(x) ≈ 1 (valor pequeno) – função é bem condicionada

Se cond f(x) ≈ 106 (valor “grande” (?)) – função é mal condicionada

Nota 1: Se uma função for bem condicionada (num ponto), então deverá existir algoritmoque permita calcular (nesse ponto) o valor da função com precisão. Contudo, podem existiralgoritmos que originem imprecisões no cálculo da função.


Número de condiçãoNota 2: cond f(x) é “grande” ou “pequeno” dependendo do nosso objectivo e da incertezados argumentos

Considerar, por hipótese,− −= → ≈ × 3 3cond ( ) 10 10y y x xf xy x

a) se os erros dos argumentos forem da ordem da representaçãodos números em computador (por exemplo em formato simples) −− ≤ ≈

710x x ux

− −− ≈ × = 3 7 410 10 10y yy

erro inferior a 0,01% erro pequeno (?)(depende da aplicação)

b) se os erros dos argumentos forem erros de leitura numa escala (temperatura,distância, velocidade, etc), por exemplo se os erros forem inferiores a 10–4

− −− ≈ × = 3 4 110 10 10y yy

erro inferior a 10% erro grande (?)(depende da aplicação)


Problema geral de interpolação== = ) ( )( , 0,1,2, , , 0,1, ,( ) i

j ji ijx ny i mfEncontrar p(x) que verifique as condições:

Exemplo: encontrar p(x) que verifique:

= == == =

1 , (1) 43 , (3 , '(3) 0

3) 2x f

x fx f

x

p(x)

1 2 3

123

4

– se todos os mi = 0(não há informação de derivadas)

− ≠ se 0, todos iguaisi im m

Interpolação de Lagrange

Interpolação de Hermite

derivada de ordem j

valores nodaisnós

– no caso contrário Interpolação de Birkhoff


Definições

∞ ∈=

[ , ]( ) max ( )

x a bf x f x

Polinómio: combinação linear de monómios xi

x

f(x)

∞( )f x

a b

xf(x)

∞( )f x

ab

x

f(x)

∞( )f x

a b

|f(x)|

= + + + +21 2( ) n

o np x a a x a x a x

– Cálculo com número finito de operações elementares (+, –, x, /)

– grau do polinómio, { }= ≠deg ( ) max , 0ip x i a

Norma do máximo (ou norma do infinito):


Algoritmo de Horner com centros

Algoritmo:

( ) ( )( ) ( )( )( )

( ) ( ){ }( )

= + − + − − + − − −

= − + − + − +

0 1 1 2 1 2 3 1 2 3

3 3 2 2 1 1 0

( )

( )

p x a a x c a x c x c a x c x c x c

p x a x c a x c a x c a

( )+

== −= × − +

=

1

Para até fa 1 0

(

zer

Fim do ciclo)

n

i i

y ai n

y y x c a

p x y

y

y

y

y


Teoremas

∞

Ω >∀ ∈ Ω ∃ − <

intervalo finito, 0

( ) ( ), : ( ) ( )f x C p f x p x

εε

Teorema de Weierstrass:

ou seja, para qualquer função contínua existe um polinómio tão próximo quanto se queira

Teorema: Se z1, z2, … , zk forem zeros distintos do polinómio p(x) então

= − ⋅ − − ⋅1 2( ) ( ) ( ) ( ) ( )kp x x z x z x z r x

Se o grau de p(x) for n então o grau de r(x) será n – k

Exemplo:

↑↑

= + − − = = − ⋅ = − ⋅ + +

grau 2grau

3 2 2

3

( ) 3 4 8 4 tem raiz 2 pelo que, ( ) ( 2) ( 2)( ) (3 5 2)p x x x x z p x xxr x xx

Teorema da unicidade: p(x) e q(x) (de grau ≤ n) interpolam os pontos(x0,y0), (x1,y1), …, (xn,yn), então p(x) – q(x) = 0, ou seja p(x) = q(x)

ou seja, o polinómio de grau ≤ n que interpola os pontos (x0,y0), (x1,y1), …, (xn,yn) é único


Teoremas

Teorema: Se z1, z2, … , zk forem respectivamente zeros de multiplicidade m1, m2, … , mkentão

1 21 2( ) ( ) ( ) ( ) ( )km m m

kp x x z x z x z r x= − ⋅ − − ⋅

Se o grau de p(x) for n então o grau de r(x) será n – (m1 + m2 + … + mk )

Teorema da unicidade: p(x) e q(x) (de grau ≤ n) verificam as seguintes condições(n+1 condições para os k+1 pontos):

ou seja, o polinómio de grau ≤ n que interpola as n+1 condições indicadas é único

( ) ( )( ) , ( ) , ..., ( ) 0,1,..., 1, 0,1,..., ,I I j ji i i i i i i i if x y f x y f x y j m i k= = = = − =

0

1k

i

i

m n=

= +então p(x) – q(x) = 0, ou seja p(x) = q(x)

Definição: Diz-se que z é um zero de multiplicidade m se( 1) ( )( ) 0, '( ) 0, ''( ) 0, , ( ) 0 , ( ) 0m mp z p z p z p z p z−= = = = ≠

Se m=1 o zero diz-se simples, se m=2 o zero diz-se duplo, etc.

( n+1 condições)


Formula de Lagrange

x

p(x)

x1

1

y1

x0

y0

L1(x)

y1L1(x)

y0L0(x)

L0(x)

−=−

10

0 1

( )( )

( )x xL xx x

−=−

01

1 0

( )( )

( )x xL xx x

= ≠= → = =

( ) 0 se ( )

( ) 1 se k i

k i kik i

L x i kL x

L x i kδ

Propriedade:

Polinómios de Lagrange:

Polinómio interpolador:

= ⋅ + ⋅0 0 1 1( ) ( )y L x y L x

= +− ⋅−−−⋅ 1

00 1

01

1 0

((

( ))) )(

)(

p x x x xxxx

yx

yx

== 0 10 1,( ) ( )p xx y ypEncontrar p(x) que verifique as condições:


Formula de Lagrange= = =10 2 210 ( ) ,( ),) (pxy xyp x ypEncontrar p(x) que verifique as condições:

x

p(x)

x1

1

y1

x0

y0

L1(x)

y1L1(x)

y0L0(x)

L0(x)

− −=− −

1 20

0 1 0 2

( )( )( )

( )( )x x x xL x

x x x x− −=− −

0 21

1 0 1 2

( )( )( )

( )( )x x x xL x

x x x x


= ≠= → = =

( ) 0 se ( )

( ) 1 se k i

k i kik i

L x i kL x

L x i kδ


= ⋅ + ⋅ + ⋅0 0 1 1 2 2( ) ( ) ( )y L x y L x y L x

− −=− −

0 12

2 0 2 1

( )( )( )

( )( )x x x xL x

x x x x

x2

L2(x)

y2L2(x)

y2

Propriedade:

= + − −⋅

+ − −⋅

−

− −

− −⋅− − −

1 20

0 1 0 2

0 21

1 0 1 2

0 12

2 0 2 1

( )( )

( )( )( )

( )( )( )(

( )

)

)

(

(

))

(x x x x

x x x xyx x

x x x xyx

x

x xy

xx x xp

x

xx


Formula de Lagrange= = =0 0 1 1( ) , ( ) , ... , ( )n np x y p x y p x yEncontrar p(x) que verifique as condições:

− +

− +

− − − −=− − − −

1 1

1 1

0

0

( ) ... ( )( ) ... ( )( )

( ) ... ( )( ) ... ( )k k

kk k k k

n

nkk

x x x x x x x xL xx x x x x x x x


= ≠= → = =

( ) 0 se ( )

( ) 1 se k i

k i kik i

L x i kL x

L x i kδ


= ⋅ + ⋅ + + ⋅0 0 1 1( ) ( ) ( ) ... ( )n np x y L x y L x y L x=

= ⋅0

( ) ( )n

k k

k

p x y L x

Propriedade:

=≠

−=−∏

1

( )( )

( )

n

ik

k iki k

x xL xx x


Método de Vandermonde2 3

0 1 2 3( ) ... nnp x a a x a x a x a x= + + + + +

0 0 02

0 1 22

0 1 22

0 1 2

2

0 0

2 2 2 2 2

1 1 1 1 1

0 1 2

( ) ( ) ... ( )

( ) ( ) ... ( )

( ) ( ) ... ( )

( ) ( ) ... ( )n n

nn

nn

nn

n n nn

n

p a a a a

p a a a a

p a a a a

x x x x

x x x x y

p a a a ax x

y

x

x y

x

x

x x y

= + + + + =

= + + + + == + + + + =

= + + + + =...

O cálculo dos coeficientes ai resulta na resolução dum sistema de equações

que se pode escrever na forma matricial

Sistema com matriz “cheia” – resolução com número de operações da ordem de n3

0 0

11 1

2

300

311

322

3

02

02

12

2

2

22

( )( )( )

(

( )( )

( )( )(

111 )

))

( )

()1 ( n

n

n

n

nnn n nn

yxyxy

x xax x

x x

y

xx

ax

ax

a

x

x

xx

=

10 210 2( , ) , ( , ) , ( , ) , ... , ( , )n nx y x xyx y y

Nota: Com o aumento do valor de n,a matriz de Vandermonde torna-secada vez mais mal condicionada


Forma de Newton

0 1 2

1

0 00 1 2

( ) ( ) ( )

1 10 2

0 21

)

3

(

1

( ) ( ) ( )( ) ( )( )( ) ...

... ( )( )( )...( )n

W x W x W x

n n

W x

p x a a x a x x a x x x

a x x x x

x xx

x

x

x

x x

x x−

−

= + − + − − + − − − +

+ − − − −

10 210 2( , ) , ( , ) , ( , ) , ... , ( , )n nx y x xyx y y

0 1 2 1( ) ( )( )( )...( )k kW x x x x x x x x x −= − − − −

Wk(x) designa-se por polinómio nodal

Nota: Wk(x) tem grau k+1

x

x0 x2 x3x1

W4(x)

x4


Forma de Newton

0 0 00 1( ) (p a a xx x= + − 20

0 0) ( xxa=

+ − 10

0 0 0)( ) ... (nx xx a x=

− + + − 1 2 1

0 1 0

0

2 0 1

0 0 0 0

1 1 1 1

)( )( )...( )

( ) ( ) ( )(

nx x x

p a a xx x x

x x

a x xx

x y−=

− − − =

= + − + − − 10

10 1) ... ( )(n x xa x x=

+ + − − 2 1

0 1 0 2 0 1 0 12 2 2 2 22 2 2

1 10

1)( )...( )

( ) ( ) ( )( ) ... ( )( )(

n

n

x x

p a a x ax x x x x x x

x x

x a x x

y

x x=

−− − =

= + − + − − + + − − − 1

0 1 0 2 0 1 0 1 2

2

1

20

)...( )

( ) ( ) ( )( ) ... ( )( )( ) ... ( )n n n n nn nn n n

n

nx x

x

p a a x a x x a x xx

x y

x x x xx x yx

=−

−

− =

= + − + − − + + − − − − =...

O cálculo dos coeficientes ai resulta na resolução dum sistema de equações

0 1 2

1

0 00 1 2

( ) ( ) ( )

1 10 2

0 21

)

3

(

1

( ) ( ) ( )( ) ( )( )( ) ...

... ( )( )( )...( )n

W x W x W x

n n

W x

p x a a x a x x a x x x

a x x x x

x xx

x

x

x

x x

x x−

−

= + − + − − + − − − +

+ − − − −

10 210 2( , ) , ( , ) , ( , ) , ... , ( , )n nx y x xyx y y


Forma de Newton

0 0

2 2 2 2 2

0

0 1 0

0 1 0 2 0 1

0 1 0 2 0 1 0 1 2

1 1 1

1

( )( ) ( )( ) ( ) ( )( )

( ) ( ) ( )( ) ... ( )( )( ) ... ( )n n n n n n n nn nn

x x y

x x x x x

p ap a a xp a a x a x x

p a a x a x x a x x x xx x

x y

x x

x y

x x y

−

= == + − == + − + − − =

= + − + − − + + − − − − =...

Ou seja, resulta na resolução do sistema de equações

que se pode escrever na forma matricial

−

− −

= − −

−−

− − − −

1 0 1 1

2 0

0 0 1 1

2 0 2 1 2 2

0

0 0

1

0 0( ) 0( ) 0

( ) ( )( )

111

00

(

1

)( )

..( ) . )

000

(( )n nn n n nn n n

x x x x a y

x

x x a yx x

x x x x x

a y

x ax x x yx x

Sistema triangular inferior – resolução com número de operações da ordem de n2


Forma de Newton

= =0 0 0[ ]a y x yLinha 1:

Linha 2: 10 0 11( )a x x ya+ − =−

=−

1 01

1 0

y yax x

−= =−

1 01 0 1

1 0

[ ] [ ][ , ]

y x y xa y x xx x

diferença dividida de ordem 1

00a y=

Resolução do sistema triangular inferior → subs tuição descendente

−

− −

= − −

−−

− − − −

1 0 1 1

2 0

0 0 1 1

2 0 2 1 2 2

0

0 0

1

0 0( ) 0( ) 0

( ) ( )( )

111

00

(

1

)( )

..( ) . )

000

(( )n nn n n nn n n

x x x x a y

x

x x a yx x

x x x x x

a y

x ax x x yx x

110 0 1( )xy x ya + − =



Forma de Newton

Linha 3: 2 0 2 0 2 10 21 2( ) ( )( )x x x x xa aa x y+ − + − − =

Resolução do sistema triangular inferior → subs tuição descendente

−

− −

= − −

−−

− − − −

1 0 1 1

2 0

0 0 1 1

2 0 2 1 2 2

0

0 0

1

0 0( ) 0( ) 0

( ) ( )( )

111

00

(

1

)( )

..( ) . )

000

(( )n nn n n nn n n

x x x x a y

x

x x a yx x

x x x x x

a y

x ax x x yx x

2 0 21 0

2 00

2 101

2( ) ( )( )x x x x xy yyx

a x yx

+ − + − − =−−

−− − −−

=− −

1 02 0 2 0

1 02

2 0 2 1

( )

( )( )

y yy y x xx xa

x x x x

−− − − − −−

=− −

+2 0 1 0 1 0 2 0

1 02

2 0 2 1

1 1( )( ) ( )( )

( )( )

x xy y x x y y x xx xa

x x x x


Forma de Newton

− − − −−

−

− −

=− −

−1 0 2 1

1 02

2

12 0

2 1

0 1 01

0

0 ( )( )

( )( )

(( ) ( )) ( )y y yy y x xx xa

x

y x

x x

x

x

x x

− − −

=− −

−−

−12 1 0 2 1

1 02

2 0 2 1

01 ( )(( )

(

) )(

)( )

y y x xx xa

x x x

y x

x

y x −− −− −

=−

1 02 1

2 1 1 02

2 0

y yy yx x x xa

x x

−= =−

1 2 0 12 0 1 2

2 0

[ , ] [ , ][ , , ]

y x x y x xa y x x xx x


−− − − − −−

=− −

+2 0 1 0 1 0 2 0

1 02

2 0 2 1

1 1( )( ) ( )( )

( )( )

x xy y x x y y x xx xa

x x x x


Forma de NewtonGeneralizando

−−= =−

1 0 10 1

0

[ ,..., ] [ ,..., ][ , ,..., ] n n

nn

ny x x y x xa y x x x

x x

diferença dividida de ordem n, calculada recorrendo a 2 diferenças

divididas de ordem n – 1

Desenvolvendo é possível concluir que o coeficiente an é a diferença dividida de ordem n

ordem n – 1

ordem n

Linha n: 2 00 2 0 2 1 0 11 2 1( ) ( )( ) ... ( )( )...( )n n n n n nx x x x x x x x xa x xa yaa x −+ − + − − + + − − − =

0 0 1

0 1

1 02 1

1 0 2 1 1 00

1 0

1

2 0

( ) ( )( ) ...

... ( )( )...( )

n n n

n n n n nn

y yy yy

a

x x x x x x

x

y x x x xyx x x x

x x x x x y−

+ − + − − +

−− −

+

− −−

− − − =

−−


Forma de NewtonAs diferenças divididas podem ser obtidas por recorrência através de

−+

+ −=−

1 11

[ ,..., ] [ ,..., ][ , ,..., ]ii i

ii

k kk

k

y x x y x xy x x xx x

Para se obter as diferenças divididas é usual construir uma tabela de diferenças divididas

=[ ]i iy x y ordem 0

=−

=−

−= =

−=

−= −=−

−−

= =−

0 0 0

1 00 1

1 0

1 2 0 11 1 1 0 1 2

2 0

0 1 2 32 1

1 2 1 2 3 0 1 22 1

3 0

2 3 1 22 2 2 1 2 3

3 1

2 3

[ ][ ] [ ]

[ , ]

[ , ] [ , ][ ] [ , , ]

[ , , , ][ ] [ ]

[ , ] [ , , ] [ , , ]

[ , ] [ , ][ ] [ , , ]

[ , ]

x y x yy x y x

y x xx x

y x x y x xx y x y y x x x

x xy x x x x

y x y xy x x y x x x y x x x

x xx x

y x x y x xx y x y y x x x

x x

y x x−

=−

=

3 2

3 2

3 3 3

[ ] [ ]

[ ]

y x y xx x

x y x y


Forma de Newton= + − + − − + − − − +

0 1 2

0 1 0 2 0 1 3 0 1 2

( ) ( ) ( )

( ) ( ) ( )( ) ( )( )( ) ...W x W x W x

p x a a x x a x x x x a x x x x x x

Dispondo duma interpolação pn – 1(x) num conjunto de n pontos, se pretendermos

adicionar mais um ponto à interpolação, o novo termo de ordem n que se adiciona à

expressão é um termo que se anula em todos os n nós anteriormente existentes.

Então, o novo polinómio (que interpola os n+1 pontos) corresponde a adicionar um novo

termo ao polinómio anteriormente existente, pn – 1(x)

= + − + − − = + − − − − 0 1( )

0 1 0 2 12 0

( )

1( ) ( ) ( )( ) 1 ( 2) *( 2)( 3)

2W x W x

p x a a x x a x x x x x x x

Exemplo: polinómio que interpola os pontos (x0,y0)=(2,1), (x1,y1)=(3,2), (x2,y2)=(5,1) é:


Forma de Newton

0 1 2

2 2

( ) (

3 3

3 2

) ( )

( )

3

)

0

(

1 2 ( 2)( 3)( 5)( ) ( 2) ( 2)( 3)

11 ( 2) *( 2 ( 2)( 3)( 5) ( 2)( 3))( 3) ( )2

( 5)

W x W x W x

W x W x

x x x

x x x x

p x a a x a x x

x x x p

a

xx xa a

= + − + − − +

= + − −

− − −

− − −− − + = − − −+

O termo W2(x)=(x–2)(x–3)(x–5) anula-se em todos os nós anteriormente existentes, x=x0=2, x=x1=3, x=x2=5, pelo que para o polinómio p3(x) interpolar estes nós, os coeficientes a0, a1 e a2 definidos anteriormente mantêm-se iguais, bastando determinar o novo coeficiente a3

− − = + −2 ( )

3 2 3 ( 2)( 3)( ) ( ) ( 5)W x

xax xp x p x

Se pretendermos adicionar à lista o ponto (x3,y3)=(6,2) resulta:

0 1 2

0 1 2

2

0 1 0 2 0 1

0 1

0

( ) ( ) ( )

( ) ( ) ( )

2

( )

1

3

23 3 ( )( )( )

( 2)(

( ) ( ) ( )( )

( 2) ( 2)( 3) 3)( 5)

p x

W x W x W x

W x W x W x

p x a a x x a x x x x

a a x a x

x x x x x x

x x x

a

ax

= + − + − − +

= + − + − −

−

−+

− −

− −


Interpolação Inversa

nós valores nodais

x0 y0 = f(x0)x1 y1 = f(x1)… …xn yn = f(xn)

= ≈ou seja , ( ) ( )ny f x p x

−

Ω = = ≈1

0 1||

( )

Se ( ) possuir inversa (em inter( , ,..., )) então podemos escrever, ( ) ( )n

f y

nf x x x x x g y p y

nós valores nodais

y0 x0 = g(y0)y1 x1 = g(y1)… …yn xn = g(yn)

ou seja, trocamos o “papel” de x com o de y

Interpolação inversa


Interpolação Inversa

= == = =

≈

Se pretendermos determinar a solução da equação ( ) 0, ou seja se pretendermosencontrar o ponto tal que ( ) 0, então podemos estimar a solução calculando

( ) ( ) e tomar para estimativa

y f xx z f z x z

p y g y = = ≈ = =da solução o valor ( 0) ( 0)z p y g y z

xx0

x2x1

f(x)

z

y1

y2

y0

y x

y0 y2y1

z

x1

x2

x0

≈( ) ( )p y g y

y

0


Teoremas

∈ Ω ∈Ω ∃ ∈Ω = ( )0 1 0 1

1( ) ( ) , , ,..., (nós distintos) , : [ , ,..., ] ( )

!k k

k kf x C x x x f x x x fk

ξ ξ

Teorema:

Demonstração: considere-se a função, Ek(x) = f(x) – pk(x), onde pk é o polinómio de grau ≤ kque interpola os k+1 nós distintos, x0, x1, …, xk.

Notar que para k=1, o teorema corresponde ao teorema do valor intermédio.

+ ΩΩ

−

A função ( ) possui (pelo menos) 1 zeros no intervalo ,

a derivada ' ( ) possui zeros em , a segunda derivada possui 1 zeros e assim sucessivamente.

k

k

E x k

E x kk

ΩA derivada de ordem possui um zero no intervalo .Designando por esse zero resulta,

kξ

= − =( ) ( ) ( )( ) ( ) ( ) 0k k kk kE f pξ ξ ξ =( ) ( )pelo que ( ) ( )k k

kf pξ ξ

= =( )0 1Atendendo que ( ) ! ! [ , ,... ]k

k k kp k a k y x x xξ

=( )0 1resulta ( ) ! [ , ,..., ]k

kf k y x x xξ

xx0

f(x)

x2 x3

p(x)

x1

E(x)


Teoremas

ξ ξ

+

+

∈ Ω ≤ ∈Ω

∀ ∈Ω ∃ ∈Ω = − = ⋅+

10 1

( 1)

( ) ( ) , ( ) (de grau n) que interpola , ,..., (nós distintos)1

, : ( ) ( ) ( ) ( ) ( )( 1)!

nn n

nn n n

f x C p x x x x

x E x f x p x f W xn

Teorema:

Demonstração:

Tendo em atenção o teorema anterior

+ += 1 1( ) ( ) , porque é nó de interpolação de ( )n nf p px x xx

+= ⋅+

( 1)1Como é um ponto qualquer do intervalo resulta, ( ) ( ) ( )

( 1)!n

n nx xE f Wn

xξ

+ = + ⋅1 0 1( ) ( ) [ , ,..., , ] ( )n n n np x p xx f x x x W x

+= = + ⋅0 11( ) ( ) ( ) [ , ,..., , ] ( )n nn nf p p f x xx x Wx x x x

+ = − = − = ⋅0 11( ) ( ) ( ) ( ) ( ) [ , ,..., , ] ( )n n nn n nx x x x x xE f p p p f x x W xx

+= ⋅ = ⋅+

( 1)0 1

1( ) [ , ,..., , ] ( ) ( ) ( )

( 1)!n

n n n nE f x x x W f Wx xn

x xξ

+∈Ω 1 0 1, ( ) interpola os nós , ,..., ,n np x x xx x x


Teoremas+= ⋅ ∈ = Ω

+( 1)

0 11

( ) ( ) ( ) , inter( , ,... )( 1)!

nn n nE x f W x x x x

nξ ξ

+∞∞

≤ ⋅+

( 1)1( ) ( ) ( )

( 1)!n

n nE x f W xn

ξ

Majorando,

Prova-se que,

+−∞ =

≤ × = = −111,...,

!( ) , max max

4n

n i i ii n i

nW x h h h x x

pelo que,

++

∞≤ ×

+

1( 1)( ) ( )

4( 1)

nn

nhE x fn

ξ

xx0 x2 x3x1

W4(x)

x4

h é o espaçamento máximo dos nós


Rigidez dos polinómios→ ∞ = − →Será que quando , ( ) ( ) ( ) 0 ?n nn E x f x p xQuestão:

Resposta: ????

n=4 n=7

( )= +( ) sin exp( 1)f x x



Por vezes os polinómios desenvolvem oscilações, que se podem tornar mais acentuadas com o aumento do grau do polinómio

Resposta: Nem sempre

2

1( )1 25

f xx

=+ ⋅

n=6

n=10Ex: Função de Runge

Nota: Interpolações com nós equidistantes



Por vezes os polinómios desenvolvem oscilações, que se podem tornar mais acentuadas com o aumento do grau do polinómio

Resposta: Nem sempre

→ Não utilizar n elevado

→ outras técnicas: por exemplo,

splines, mínimos quadrados , …

→ Contradição com Teo. de Weierstrass?

→ Não. O teorema de Weierstrass não diz que o polinómio (que aproxima a função) é obtido por interpolação

n=8


Nós de Chebyshev

depende do comportamento

da função

Majorante do erro de interpolação, ( 1)1( ) ( ) ( )( 1)!

nn nE x f W x

nξ+

∞ ∞∞≤ ⋅

+

depende da posição dos nós

majorante do erro da

interpolação

Para nós equidistantes, Wn(x) atingemaiores valores absolutos junto dasextremidades (e menores valores juntodo centro)

Ideia: posicionar os nós de modo aminimizar o máximo absoluto de Wn(x)

=> Juntar os nós na extremidade (econsequentemente afastar no centro)

W6(x) com nós equidistantes


Nós de Chebyshev∞

Demonstra-se que o menor valor para ( ) ocorre se os nós forem posicionados nos

zeros dos polinómios de Chebyshev.nW x

os zeros dos polinómios dPara o intervalo [ 1, 1] e Chebyshev são dados , por,ξ ∈ − +

− =

(2 1)=cos , 1, 2, ... ,

2kk k N

Nπξ Nesta expressão, N representa

o número de pontos

W6(x) com nós equidistantes W6(x) com nós de Chebyshev


Nós de Chebyshev

ξ∈Ω =

∈ − + ∈Para o intervalo genérico [ , ] podemos efectuar uma transformação de coordenadasde [ 1, 1] para [ , ]

x a bx a b

++

+∞ ∞

−≤ ×× +

1( 1)

2 1

( )( )

2 ( 1)!

nn

n n

b aE x fn

− += × + ×1 1( )

2 2x a bξ ξξ

a b x-1 +1 ξ

−∞

− + ≤Demonstra-se que, , utilizando nós de Chebyshev, no intervalo [ ] 2, (1 )1 nnW x

+∞ ∞∞

≤ ⋅+

( 1)1pelo que ( ) ( ) ( )

( 1)!n

n nE x f W xn

ξ +∞ ∞

≤ ×× +

( 1)1( )

2 ( 1)!n

n nE x fn

a expressão do majorante do erro de interpolação utilizaPara um

ndo nós intervalo genérico

de Chebyshev tom [

a a forma, ]a bΩ =

x(ξ)


Nós de Chebyshev

n=6 n=10

2

1( )1 25

f xx

=+ ⋅

Ex: Função de Runge Interpolação com nós de Chebyshev


Nós equidistantes

n=6n=10

2

1( )1 25

f xx

=+ ⋅

Ex: Função de Runge Interpolação com nós equidistantes


( )0 1 0 1

1( ) ( ) , , ,..., , : [ , ,..., ] ( )

!k k

k kf x C x x x f x x x fk

ξ ξ∈ Ω ∈Ω ∃ ∈Ω =Teorema:

xx0 x1ξ

1

10 1

||( ) (

0 1 0

) intervaloque contém

,

11, [ , ] '( ) , inter( , )

o

o

f x f xx x x x

k f x x f x xξ ξ↑

−−

= = ∈ = Ω

→↑

= =1 0

dif. divididaconflue

0 1 0 0

nte

0lim [ , ] [ , ] '( )x x

f x x f x x f x xx0 x1ξ

Diferenças divididas confluentes

0 1 2 0 1 2''( )

2, [ , , ] , inter( , , )2

fk f x x x x x xξ ξ= = ∈ xx1 x2ξx0

xx1 x2ξx0

Assim, caso haja informação das derivadas, podemos inserir essa informação na tabela de diferenças divididas através de diferenças divididas confluentes

→= =

1 2 0

00 1 2 0 0 0,

''( )lim [ , , ] [ , , ]

2x x x

f xf x x x f x x x



( ) ( ){ }0

22 2

1, ,( ) 1 2 ( ) ( ) ( )

n

k k k k k k k

k

np x L x x x L x y x x L x y=

+ = − ⋅ ⋅ − ⋅ ⋅ + − ⋅ ⋅

1 nó

0 1

0 1 0 1

Teorema: O polinómio de grau que nos nós distintos , ,..., interpola os, , ,valores , ,..., e as derivadas , ,..., é:

1

2n s

n

n n

x x x

y y y y y

n

y

+

≤ +

No caso de, em todos os nós, existir informação da função e da(s) derivada(s)


Neste(s) caso(s) há formulas específicas para a interpolação

Para o caso da interpolação de Hermite (função e primeira derivada) em n+1 nós, o erroassociado à interpolação é:

Nota: Lk(x) são os polinómios de Lagrange

( )(2

2 11

2)2

2( )

( ) ( ) ( ) ( ) ,(2 2)!

n

n nnfx f x p x W x

nE ξ ξ+

+

+= − = ⋅ ∈Ω+


SplinesSplines: Troços de polinómios assegurando uma certa continuidade na transição dos troços

Spline linear – troços de rectas com continuidade da função na transição dos troços

xx0

h1

x2 x3

S1(x)

x1

troço 1 troço 2 troço 3

S2(x) S3(x)

h2 h3

Spline quadrático – troços de parábolas com continuidade da função e da 1ª derivada na transição dos troços

Spline cúbico – troços de cúbicas com continuidade da função, da 1ª derivada e da 2ª derivada na transição dos troços

equação duma recta

continuidade na função

Spline linear

1 1,...,, maxi i i ii n

h x x h h− == − =


Spline quadráticoSpline quadrático: Troços de polinómios de grau ≤ 2, que se unem de modo contínuo e com tangente também contínua

xx0

h1

x2 x3

S1(x)

x1


S2(x)

S3(x)

h2 h3

equação duma

parábola

continuidade na função e na tangente 1 1,...,

, maxi i i ii nh x x h h− =

= − =

• Troços de parábolas

• Condições a impor:

→ Interpolar os pontos(=> continuidade na função)

→ Con nuidade na derivada(na transição dos troços)


Spline quadrático

20 1 2( )iS x a a x a x= + +

Reescrevendo a equação da recta na forma de Lagrange,

( )1troço ,i ii x x x−→ ∈

1 2'( ) 2iS x a a x = +x

xi-1hi

Si(x)

xi

troço i

−−

− −= + 11

,( ) i ii i i

i i

x x x xS x m mh h

ou seja, as derivadas à esquerda e à direita do nó xi são iguais e valem mi mi é a derivada do spline no nó xi1

,( )'( ), ( )

i ii i

i i

i

i

S x mS x m

S x m+

= → ==

xx1 x3

S2(x)

x2

troço 2 troço 3

S3(x)Exemplo:

3 2

2

2

2 2 22,(

' ))

(,(

)S x mm

S xS x

m

=

→ =

=

tangente com declive m2

Ou seja, escrevendo a equação da derivada do spline na forma de Lagrange, garantimos desde logo a continuidade na derivada na transição entre os troços

++ +

+ +

− −→ = +11 1

1 1

, ( ) i ii i i

i i

x x x xS x m mh h


Spline quadráticoDesenvolvendo a expressão da derivada, podemos reescrevê-la na forma de Newton,

Primitivando a expressão,

( )11 1

2)2

( i ii i i

i

m mS x m x xxh

C−− −

−= + − +

Impondo a condição do spline passar no ponto da esquerda do intervalo

xxi-1

Si(x)

troço i

yi-1

xi

1 1( ) 1,2, ... ,ii iS i nx y− −= =

( )1 1 121

1 1

02

i ii i

ii i ix x ym mm x C

h− −

=

− − −−−+ − + = 1 1 1i i iy mC x− − − −=

( ) 1 11 121

1( )2

i ii i i

ii i i

m mS x m x x x m xh

y−− − −− −

−= + − + −pelo que,

−−

− −= + 11

,( ) i ii i i

i i

x x x xS x m mh h

− −−−

− + − − = +

1 1 1

1,( )

i

i ii

h

ii

i i

ii

x xx x x xS x m mh h

− −−

− − = − +

1 1

1,( ) 1 i ii i i

i i

x x x xS x m mh h

( )11 1

,( ) i ii i i

i

m mS x m x xh

−− −

− = + −


Spline quadráticoReescrevendo a expressão anterior,

Impondo a condição do spline passar no ponto da direita do intervalo

xxi-1

Si(x)

troço i

yi

xi

( ) 1,2, ... ,ii iS i nx y= =

( ) ( )211 1 1 1( ) 1,2, ... ,

2i i

i i i i ii

m mS x y m x x x x i nh

−− − − −

−= + − + − =

Nesta expressão os mi são incógnitas m0 , m1 , … , mn n+1 incógnitas

( ) ( )211 1 1 12

i i

i ii i i

h h

i i iii

m my m x xx xh

y−− − − −

−+ − + − =

11 1 2

i ii i i i i

m my m h h y−− −

− + + = 1

12i i

i i im m h y y−

−+

= −

11 2 i i

i ii

y ym mh

−−

− + = 1

12 1,2, ... ,i ii i

i

y ym m i nh

−−

− = − =

Esta expressão representa n equações


Spline quadráticoResumindo, a expressão do spline quadrático é

Para determinarmos as incógnitas recorremos à condição

( ) ( )211 1 1 1( ) 1,2, ... ,

2i i

i i i i ii

m mS x y m x x x x i nh

−− − − −

−= + − + − = onde os mi são incógnitas n+1 incógnitas

112 1,2, ... ,i i

i ii

y ym m i nh

−−

−= − = n equações

n equações a n+1 incógnitas necessário fornecer condição suplementar

fornecer o valor de um dos mi vulgar fornecer m0 (uma estimativa de f’(x0))

Resulta o sistema

0

11

fornecer

2 1,2, ... ,i ii i

i

my ym m i n

h−

−

− = − =

0

11

fornecer

2 1,2, ... ,i ii i

i

my ym m i n

h−

−

⇔ − + = =


Spline quadráticoO sistema pode ser colocado na forma matricial (com duas diagonais não nulas)

( )( )( )

( )( )

1 1 0 1

2 2 1 2

3 3 2 3

2 2 3 2

1 1 2 1

0 00 0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0 0

0 0 0 0 0

0 0 0 0

1 1 2

1 1 2

1 2

1 1 2

1 1 2

1

0

0 0 0 2

,

0 1

1

0

n n n n

n n n n

n n

m y

m y y h

m y y h

m y y h

m y y h

m y y h

m y

− − − −

− − − −

− − − = − −

( )1n ny h−

−


Spline cúbicoSpline cúbico: Troços de polinómios de grau ≤ 3, que se unem de modo contínuo, com tangente contínua e com curvatura também contínua

xx0

h1

x2 x3

S1(x)

x1


S2(x)

S3(x)

h2 h3

equação duma cúbica

continuidade na função,

na tangente e na curvatura

1 1,...,, maxi i i ii n

h x x h h− == − =

• Troços de cúbicas

• Condições a impor:

→ Interpolar os pontos(=> continuidade na função)

→ Con nuidade na 1ª derivada(na transição dos troços)

→ Con nuidade na 2ª derivada (na transição dos troços)


Spline cúbico2 3

0 1 2 3( )iS x a a x a x a x= + + +

Reescrevendo a equação da recta na forma de Lagrange,

( )1troço ,i ii x x x−→ ∈

xxi-1

hi

Si(x)

xi

troço i

11

,,( ) i ii i i

i i

x x x xS x M Mh h

−−

− −= +

ou seja, as 2ª derivada à esquerda e à direita do nó xi são iguais e valem Mi Mi é a 2ª derivada do spline no nó xi1

,,( ) ,,( ),, ( )i i

i i

i i

i

i

S x MS x M

S x M+

= → ==

xx1 x3

S2(x)

x2

troço 2 troço 3

S3(x)Exemplo:

3 2

2

2

22

2 2

,,(

,,,,(

( )

))

S

MM

M

xx

x

SS

=

→ =

=

2ª derivada com valor M2

Ou seja, escrevendo a equação da 2ª derivada do spline na forma de Lagrange, garantimos desde logo a continuidade na 2ª derivada na transição entre os troços

11 1

1 1

,, ( ) i ii i i

i i

x x x xS x M Mh h+

+ ++ +

− −→ = +

21 2 3

,( ) 2 3iS x a a x a x = + + 2 3,,( ) 2 6iS x a a x = +


Spline cúbicoPrimitivando 2 vezes a expressão,

O termo α x + β pode ser escrito na forma de Lagrange,

11

,,( ) i ii i i

i i

x x x xS x M Mh h

−−

− −= + ( ) ( )1

21

2

2, )

2( i i

i i ii i

x x x xS x M M

h hα−

−

− − = + +

( ) ( )11

3 3

6( )

6i i

i i ii i

x x x xS x M M

h hxα β−

−

− − + + +=

( ) ( )3 3

111( )

6 6i i

i i ii i

i i

i i

x x x xc dh

x x xh

xS x M M

h h− −

−− −+

− −= + +

1com ; i i

i i

d c cx dxh h

α β −− −= =


Spline cúbico

Impondo as condições de interpolação

xxi-1

Si(x)

troço i

yi-1

xi

1 1( )1,2, ... ,

( ) i

i i

ii

iSi

x yS y

nx

− −== =

pelo que, para o troço i, a equação do spline é,

( ) ( )3 31 1

1( )6 6

i i i ii i i

i i i i

x x x x x x x xS x M M c dh h h h

− −−

− − − −= + + +

yi

113

11

(( )6

i

iii i

i

h

iixx

hxxM M −− −

−

−− +

1

0

13

1)6

i

ii

i

h

ii

i

xxc dh h

xx −−

=

−−−+ +

1

0

1

(

i

i

i

ii

h

xM

y

x−

=

−=

− 3 31

0

) ( )6 6

i

iii

h

i

i

i

ixxxM ch

xh

=

−−−+ +

1

0 i

ii

h

i

i ixx ydh h

=

−

− + =

2

1 1

2

6

6

ii i

ii i

hc y M

hd y M

− −= −

=

−

( ) ( )3 3 2 21 1

1 1 1( ) 1,2, ... ,6 6 6 6i i i i i i

i i i i i i ii i i i

x x x x h x x h x xS x M M y M y M i nh h h h

− −− − −

− − − −= + + − + − =

Nesta expressão os Mi são incógnitas M0 , M1 , … , Mn n+1 incógnitas


Spline cúbicoDerivando a expressão anterior (expressão do spline) resulta,

Impondo a continuidade na primeira derivada

( ) ( )2 2 2 21

1 1 11 1,( )

2 2 6 6i i i i


x x x x h hS x M M y M y Mh h h h

−− − −

− − = − + − − + −

( ) ( ) ( )2 2

1 11 1

,( )2 2 6

i i i i ii i i i i

i i i

x x x x y y hS x M M M Mh h h

− −− −

− − − = − + + − −

xx1 x3

S2(x)

x2

troço 2 troço 3

S3(x)2 2,( )S x− 3 2

,( )S x+

1, ,( ) ( )i ii iS Sx x+

− +=

1

(,( ) ii

iii

xS M

xx −

−= − ( )

2 21

0

11

) ( )2 2 6

i

i i i ii i

h

ii i i

i x y y hM M Mh h h

x=

− −−

− −+ + − −

1

2

11

11

(( ), ( )2

i

iiii i

i

i

h

ii

xS xxM M

xh

x

+

++

++

−−= − +

( )2

1 11

1

0

1

)2 6

i i ii i

i i

y y hM Mh h

+ +

+ +

=

+−+ − −

( ) ( ) ( )2 2

1 1 11 1

1 1 1

, ( )2 2 6

i i i i ii i i i i

i i i


+ + ++ +

+ + +

− − − = − + + − −


Spline cúbico

pelo que a continuidade na derivada resulta,

1, ,( ) ( )i ii iS Sx x+

− +=

( )11

,( )2 6

i i i ii i i i

iix h y y hS M M M

h−

−−

= + − −

( )1 1 111

2

1

, ( )2 6i i i i

i i ii

ii x h y y hS M M Mh+

+ + ++

+

− = − + − −

( ) ( )1 1 1 1

11 12 6 2 6i i i

i i i i i i i i

i ii ii

h y y h h yM M y hh h

M MM M− + + +

++−

− − + − − = − + − −

11

1 1 1

11

2 21,2, ... , 1

6 6 6i i i i i

ii

ii i

i ii

h h h h y y y yM M i nh h

M + + ++

−−

+

+ − − + + = − = −

Esta expressão representa n – 1 equações

1

2

11

11

(( ), ( )2

i

iiii i

i

i

h

ii

xS xxM M

xh

x

+

++

++

−−= − +

( )2

1 11

1

0

1

)2 6

i i ii i

i i

y y hM Mh h

+ +

+ +

=

+−+ − −

1

(,( ) ii

iii

xS M

xx −

−= − ( )

2 21

0

11

) ( )2 2 6

i

i i i ii i

h

ii i i

i x y y hM M Mh h h

x=

− −−

− −+ + − −


Spline cúbicoResumindo, a expressão do spline cúbico é

Para determinarmos as incógnitas recorremos à condição

onde os Mi são incógnitas n+1 incógnitas

n – 1 equações

n – 1 equações a n+1 incógnitas necessário fornecer 2 condições suplementares

( ) ( )3 3 2 21 1

1 1 1( ) 1,2, ... ,6 6 6 6i i i i i i


x x x x h x x h x xS x M M y M y M i nh h h h

− −− − −

− − − −= + + − + − =

1 1 1 11 1

1

2 2 1,2, ... , 16 6 6

i i i i i i i ii i i

i i

h h h h y y y yM M M i nh h

+ + + −− +

+

+ − −+ + = − = −


Spline cúbicoHipótese 1: Spline natural

Na ausência de outra informação considerar 0 0,,( ) 0,,( ) 0n n

S x M

S x M

= =

= =

ou seja, o sistema a resolver é

1

0

1 1 1 11

1

1,2, ... , 16 3 6

0

0


i i

n

h h h h y y y yM M M i nh

M

M

h+ + + −

− ++

=

=

+ − − + + = − = −

No caso de espaçamento uniforme resulta (hi = h ; i = 1, 2, … , n)

+ −− +

− + + + =

=

=

1 11

0

12 2

6 3 6

0

0

i i ii i i

n

h h h y

M

y yh

M

M M M + −− +

− + +

=

+

=

=

1 11 1 2

0

2

0

4 6

0

i i ii i i

n

y y yM M Mh

M

M


Spline cúbicoO sistema pode ser colocado na forma matricial (sistema tridiagonal)

2 1 02

1 3 2 12

2

3

3

2 1 2

1

0

32

0 0 0 0 0 0 00 0 0 0 0

0 0

26

1 4 1 261 4 1

1 4

4 11 4 1 2

0 0 00 0 0 0 0 0

0 0 0 0 0 00

61 4 1

2

0 0 0 00 0 0 0 00 0 0

0

1

0 60 0 10

n

n n n n

n

nn

y y yh

M y y yM hM

MM y y y

hMy

M

M

−

− − − −

−

− +

− +

= − + −

1 2

2

0

n ny yh

− −

+


Spline cúbicoHipótese 2: Spline completo

Caso as 1as derivadas nas extremidades sejam conhecidas considerar 0 0, ,( ), ,( )n n

S x y

S x y

=

=( ) ( ) ( )− −

− −

− − −= − + + − −2 2

1 11 1

,( )2 2 6

i i i i ii i i i i

i i i


−−= − +

10

20 01

1 01

0 1

(( ),( )2

h

xxS Mxx Mx

h( )

=

−+ − − =2

1 0 11 0 0

1 1

0

) ,2 6

y y hM M yh h

− − − = − 1 01 1

0 1 01

,3 6

y yh hM M yh

−

−= − 1

(,( ) nn n

nn

xS M

xx ( )− −

−

=

− −+ + − − =

02 2

1 11

) ( ) ,2 2 6

n

n n n nn n n n

h

n

n

n n

x y y hM M Mx yh h h

−−

− + = − 1

1,

6 3n n n n

n n nn

h h y yM M yh


Spline cúbicoou seja, para spline completo o sistema a resolver é

−−

+ + + −− +

+

+ − − + + = − = −

−+ =

−

−+ = −

1 1 1 11 1

1

1 01 10 1 0

1

11

,

1,2, ... ,

3 6

,6 3

16 3 6

n n n nn n n

n


i i

y yh hM M yh

h h y yM

h h h h y y y yM M M i nh h

M yh

Hipótese 3: spline periódico

Hipótese 3: continuidade da 3ª derivada em x1 e xn–1, …

• • • (consultar referências)

• • • (consultar referências)


Splines cúbicos

Ω = 0 1 0 1(Holladay): [ , ] Nós: , ,..., Valores nodais: , ,...,n n

ba

a b x x x y y yTeorema

Interpretação: De todas as curvas interpoladoras com continuidade até à segundaderivada, o spline cúbico natural é a curva mais direita possível (porque o valor médioda curvatura é mínimo)

[ ]

2

2

De todas as funções com continuidade até à segunda derivada ( ( )) que interpolam os anteriores valores, o spline cúbico natural é a função que minimiza o valor

( ) ''( )

e é uma função única

b

a

f C

J f f x dx

∈ Ω

= .

Nota: J(f) é proporcional à energia de deformação duma viga elástica


Erros da interpolação por splines cúbicos

4

3

4

4

4

4

0 0

1 1

2 2

2

3

2

3

5384

3 1216 24

1 112 3

1

' '

'' ''

''' ''' 12

max mini i

f S C D f h C

f S C D f h C

hf S C D f h Ch

hf S C D f h Ch

h h h h

∞ ∞

∞ ∞

∞ ∞

∞ ∞

− ≤ ⋅ ⋅ =

− ≤ ⋅ ⋅ = +

− ≤ ⋅ ⋅ = + ⋅

− ≤ ⋅ ⋅ = +

= =


Diferença finita regressiva de 2º ordem

−=−

−− −− −=

−−=−

0 0

1 00 1

1 0

1 02 1

2 1 1 01 1 0 1 2

2 0

2 11 2

2 1

2 2

[ , ]

[ , , ]

[ , ]

x yy yy x xx x

y yy yx x x xx y y x x x

x xy yy x xx x

x y

( )0 0 1 0 0 1 2 0 1

0 1 0 1 2 1 0

0 1 2

( ) [ , ] ( ) [ , , ] ( ) ( )

'( ) [ , ] [ , , ] ( ) ( )

''( ) 2 [ , , ]

p x y y x x x x y x x x x x x x

p x y x x y x x x x x x x

p x y x x x

= + ⋅ − + ⋅ − ⋅ −

= + ⋅ − + −

= ⋅

Polinómio interpolador nos nós x0, x1, x2 – tabela diferenças divididas

x0 x1 x2

p(x)



Considerando nós equidistantes, h=x1 – x0, h=x2 – x1

( )

( )

= + ⋅ − + −

−− −− − −= + ⋅ − + −− −

0 1 0 1 2 1 0

1 02 1

1 0 2 1 1 01 0

1 0 2 0

'( ) [ , ] [ , , ] ( ) ( )

'( ) ( ) ( )

p x y x x y x x x x x x x

y yy yy y x x x xp x x x x xx x x x

−− −−− −− − − += ⋅ = ⋅ = ⋅ =−

1 02 1 1 02 1

2 1 1 0 2 1 00 1 2 2

2 0

2''( ) 2 [ , , ] 2 2

2

y yy y y yy yx x x x y y yh hp x y x x x

x x h h

( )−− −−= + ⋅ − + −

1 02 1

1 01 0'( ) ( ) ( )

2

y yy yy y h hp x x x x x

h h

( )− − += + ⋅ − + −1 0 2 1 01 02

2'( ) ( ) ( )

2y y y y yp x x x x x

h h

x0 x1 x2

h h



Particularizando para x=x2, resulta

− − + = + ⋅ − + − 2 2

1 0 2 1 01 02 2

2

2'( ) ( ) ( )

2h h

y y y y yp x xx xh h

x

x

x0 x1 x2

p(x)

A diferença finita designa-se porregressiva porque os nós utilizadossão o ponto onde calculamos a“derivada” e pontos “para trás”

( )− − += + ⋅ +1 0 2 1 02 2

2'( ) 2

2y y y y yp x h h

h h

− += 2 1 02 2

2''( )

y y yp xh

− − += +1 0 2 1 02

2'( ) 3

2y y y y yp x

h h

− += 2 1 02

3 4'( )

2y y yp x

h


Diferença finita regressiva de 2º ordemExpansão em série de Taylor de ordem 2 com resto de ordem 3

ξ− = + × − + × − + × −2 31''( ) '''( )( ) ( ) '( ) ( ) ( ) ( )

2 6f x ff x h f x f x h h h

ξ− = + × − + × − + × −2 32''( ) '''( )( 2 ) ( ) '( ) ( 2 ) ( 2 ) ( 2 )

2 6f x ff x h f x f x h h h

ξ − = − + −2 31''( ) '''( )

( ) ( ) '( )2 6

f x ff x h f x f x h h h (*)

2 32''( ) '''( )( 2 ) ( ) 2 '( ) 4 8

2 6f x ff x h f x f x h h hξ

− = − × + × − × (**)

Efectuando (**) – 4×(*), de modo a anular o termo em f ’’(x), resulta

3 31 2'''( ) '''( )( 2 ) 4 ( ) 3 ( ) 2 '( ) 8 4

6 6f ff x h f x h f x f x h h hξ ξ− − − = − + − +

2 32''( ) '''( )( 2 ) ( ) 2 '( ) 4 8

2 6f x ff x h f x f x h h hξ− = − × + × − ×

ξ− =× × × ×−×− + 2 314 4 4 4''( ) '''( )

( ) ( ) '( )2

46

f x ff x h f x f x h h h×4 (*)

− ×4(**) (*)

(**)


Diferença finita regressiva de 2º ordemDesenvolvendo

3 31 2'''( ) '''( )( 2 ) 4 ( ) 3 ( ) 2 '( ) 8 4

6 6f ff x h f x h f x f x h h hξ ξ− − − = − + − +

3'''( )( 2 ) 4 ( ) 3 ( ) 2 '( ) 4

6ff x h f x h f x f x h hξ

− − − = − + − [ 2 , ]x h xξ ∈ −

ξ− − − + = +

( )

2( 2 ) 4 ( ) 3 ( ) '''( )'( )

2 3hD f x Erro

f x h f x h f x ff x hh

Diferença finita regressiva de 2ª ordem e respectivo erro (de método)

( 2 ) 4 ( ) 3 ( )( )

2hf x h f x h f xD f x

h− − − += 2

2'''( )

'( )3

fE x hξ=


Resumo das principais diferenças finitas

( ) ( )( )h

f x f x hD f x

h− −=

Prim

eira

der

ivad

a

2 ''( ) '''( )E x h f ξ= − ⋅

Prim

eira

ord

em1

1'( ) ''( )

2E x h f ξ= ⋅ ⋅

( ) ( )( )

2hf x h f x h

D f xh

+ − −= 21

1'( ) '''( )

6E x h f ξ= − ⋅ ⋅

( ) ( )( )h

f x h f xD f x

h+ −= 1

1'( ) ''( )

2E x h f ξ= − ⋅ ⋅

Segu

nda

orde

m ( 2 ) 4 ( ) 3 ( )( )

2hf x h f x h f x

D f xh

− + + + −= 22

1'( ) '''( )

3E x h f ξ= ⋅ ⋅

3 ( ) 4 ( ) ( 2 )( )

2hf x f x h f x h

D f xh

− − + −= 22

1'( ) '''( )

3E x h f ξ= ⋅ ⋅

Segu

nda

deriv

ada 2

2

( 2 ) 2 ( ) ( )( )h

f x h f x h f xD f x

h+ − + +=

Segu

nda

orde

m Progressiva

Regressiva 22

( ) 2 ( ) ( 2 )( )h

f x f x h f x hD f x

h− − + −= 2 ''( ) '''( )E x h f ξ= ⋅

22

( ) 2 ( ) ( )( )h

f x h f x f x hD f x

h+ − + −= 2

21

''( ) ''''( )12

E x h f ξ= − ⋅ ⋅Central

Progressiva

Regressiva

Central

Progressiva

Regressiva


Formula de Aitken-Neville

x0 x1 xn xn+1

1 1 0 1 1( ) ( ) ( ) ( ) [ , ,... , ] ( )n n n n n n n np x p x a W x p x y x x x x W x+ + += + ⋅ = + ⋅

( )0,00

1,,, 1

1 00

( )( ) )

( ( )(

) nnnn

n

p x p xx xp x

xp

xx+

+

− ⋅ −= +

−


x0 x1 xn xn+1

p0,n

p1,n

p0,n+1

x2 xn-1



( )000

,20,,3

3 02

1,2( )( )

( )(

( ))

p xx xp x

x x

p xp x

− ⋅ −= +

−


Exemplo: n=2

x0 x1 x2

p0,2(x)

x3

p1,2(x)p0,3(x)



( )1

1,

,,1

, ( )(

( ( ))( ) ) m km

m km k

m km k

m

p p xp x

x xp x

x

x

x+

++

+

− ⋅ −= +

−

Generalizando


x0 x1 xm+k

pm,k

pm+1,k

pm,k+1

xm xm+k+1xm+1



( )1, ,, 1 ,

1

( ) ( ) ( )( ) ( ) m m k m k

m k m km k m

x x p x p xp x p x

x x+

++ +

− ⋅ −= +

−

Se pretendermos efectuar o cálculo do polinómio para x=0 resulta,

( ), 1,, 1 ,

1

(0) (0)(0) (0) m m k m m k

m k m km k m

x p x pp p

x x+

++ +

⋅ − ⋅= +

−

( )1, ,, 1 ,

1

0 0 00 0

( ) ( ) ( )( ) ( ) m m k m k

m k m km k m

x p pp p

x x+

++ +

− ⋅ −= +

−

1 , 1,, 1

1

(0) (0)(0) m k m k m m k

m km k m

x p x pp

x x+ + +

++ +

⋅ − ⋅=

−


Extrapolação de Richardson

Extrapolação de Richardson:

• Sucessão de valores de h: hm, hm+1, …, hm+k

• Calcular diferenças finitas para cada valor de h –> termos de ordem zero dos polinómios

1 , 1,, 1

1

(0) (0)(0) m k m k m m k

m km k m

x p x pp

x x+ + +

++ +

⋅ − ⋅=

−1 , 1,

, 11

m k m k m m km k

m k m

h D h DD

h h+ + +

++ +

⋅ − ⋅=

−

• Usar Aitken-Neville para calcular os valores do polinómio em h=0 (extrapolação)

2 3''( ) '''( )( ) ( ) '( ) ...

2! 3!f x f xf x h f x f x h h h+ = + ⋅ + ⋅ + ⋅ +Expansão em série

2

1

( ) ( ) ''( ) '''( ) ( ) ( )'( ) ... '( )

2! 3!i

i

i

f x h f x f x f x f x h f xf x h h f x C hh h

∞

=

+ − + − = − ⋅ − ⋅ − = + ⋅

1

polinómio na variável

( ) ( )'( ) ( )

k

ii

i

h

f x h f xf x C h p hh

=

+ − + ⋅ =

A diferença finita é o termo de ordem zero do polinómio na variável h


Extrapolação de RichardsonExemplo: f(x)=exp(x), calcular f ’(1) utilizando extrapolação de Richardson

1 0,0 0 1,00,1

1 0

0.2 3.342295 0.4 3.0091750.2 0.4

h D h DD

h h× − × × − ×= =

− −

Nota: Valor exacto, f ’(1)=exp(1)=2.718282Sugestão: considerar h0=0.4, hm=h0/2m

0 ,00( ) ( ) (1.4) (1)

0.4 Diferença finita,0.4

f x h f x f fh Dh

+ − −= → = =

01 ,01

( ) ( ) (1.2) (1)0.2 Diferença finita,

2 0.2h f x h f x f fh D

h+ − −= = → = =

1 , 1,, 1

1

m k m k m m km k

m k m

h D h DD

h h+ + +

++ +

⋅ − ⋅=

−

0,1 2.676055D =Tabela

0,1

0 0,0

1 1,0

h D

h DD

Gráfico

h1 h0

hD0,1

D0,0

D1,0

00, 3.342295D =

01, 3.009175D =


Extrapolação de RichardsonNovo ponto, h2=h0/22=0.1

2 1,0 1 2,01,1

2 1

0.1 3.009175 0.2 2.8588422.708509

0.1 0.2h D h D

Dh h

× − × × − ×= = =− −

Valor exacto,

f ’(1) =2.718282

02 ,02 2

( ) ( ) (1.1) (1)0.1 Diferença finita, 2.858842

2 0.1h f x h f x f fh D

h+ − −= = → = = =

1 , 1,, 1

1

m k m k m m km k

m k m

h D h DD

h h+ + +

++ +

⋅ − ⋅=

−

Tabela0 0,0

0,1

1 1,0

1,1

2 2,

0,2

0

Dh D

h

h

D

DD

D

Gráfico

h1 h0

D0,0

D1,0

h2

D0,1D1,1

D0,2

D2,0

2 0,1 0 1,10,2

2 0

0.1 2.676055 0.4 2.7085092.719327

0.1 0.4h D h D

Dh h

× − × × − ×= = =− −


Extrapolação de RichardsonTabela

= =

= =

= =

=

==

0 0,0

1 1,0

2 2,

0,

,

0

0,1

1

2

1

0.4 3.342295

0.2 3.00912.676055

2.7085092.7193275

0.1 2.858842

7

h DD

Dh D

h DD

Gráfico

h1 h0

D0,0

D1,0

h2

D0,1D1,1

D0,2

D2,0

Valor exacto,

f ’(1) =2.718282

diferenças finitas


Extrapolação de Richardson

A convergência é mais rápida e a formula de Aitken-Neville adaptada resulta

2 4 2

1

( ) ( ) '''( ) '''''( ) ( ) ( )'( ) ... '( )

2 3! 5! 2i

i

i

f x h f x h f x f x f x h f x hf x h h f x C hh h

∞

=

+ − − + − − = − ⋅ − ⋅ − = + ⋅

1 , 1,, 1

1

(0) (0)(0) m k m k m m k

m km k m

x p x pp

x x+ + +

++ +

⋅ − ⋅=

−

2 21 , 1,

, 1 2 21

m k m k m m km k

m k m

h D h DD

h h+ + +

++ +

⋅ − ⋅=

−

Se usarmos diferenças finitas centrais a série de potencias de h só tem potências pares

2 3

2 3

3

''( ) '''( )( ) ( ) '( ) ...

2! 3!''( ) '''( )

( ) ( ) '( ) ...2! 3!

'''( )( ) ( ) 2 '( ) 2 ...

3!

f x f xf x h f x f x h h h

f x f xf x h f x f x h h h

f xf x h f x h f x h h

+ = + ⋅ + ⋅ + ⋅ +

− − = − ⋅ + ⋅ − ⋅ +

+ − − = ⋅ ⋅ + ⋅ ⋅ +

2

2 2

1

polinómio na variável

( ) ( )'( ) ( )2

k

ii

i

h

f x h f x hf x C h p hh

=

+ − − + ⋅ =


Integração Numérica – Regras de Newton-Cotes

primitivaro polinómio

( ) ( ) ( ) ( ) ( ) ( )

b b

n h

a a

I f f x d x p x d x I f

= ≈ =

Aproximar a função integranda por um polinómio interpolador, utilizando para nós deinterpolação os extremos do intervalo e nós igualmente espaçados no interior do intervalo

n=0 (interpolação grau zero) – regras do rectângulo à esquerda, à direita e do ponto médio

( )( ) ( )hI f b a f a= − ⋅

a b

f(x)

Ih(f)

a b

f(x)

Ih(f)

( )( ) ( )hI f b a f b= − ⋅ ( ) ( )2h

a bI f b a f + = − ⋅

a b

f(x)

Ih(f)

(a+b)/2


Integração Numérica – Regras de Newton-Cotesn=1 (interpolação linear) – regra do trapézio

[ ]( ) ( )( ) ( ) ( ) ( )2 2h

f a f b b aI f b a f a f b+ − = − ⋅ = ⋅ +

a b

f(x)

Ih(f)

(a+b)/2

p(x)

a b

f(x)

Ih(f)

p(x)

n=2 (interpolação quadrática) – regra de Simpson

( ) 4 (( ) / 2) ( )( ) ( )6

( ) 4 ( )6 2

hf a f a b f bI f b a

b a a bf a f f b

+ ⋅ + + = − ⋅ − + = ⋅ + ⋅ +


Integração Numérica – Dedução da regra de Simpsonn=2 (interpolação quadrática) – Formula de Lagrange

0 2 0 11 22 0 1 2

0 1 0 2 1 0 1 2 2 0 2 1

2 2

( ) ( ) ( ) ( )( ) ( )( )

( ) ( ) ( ) ( ) ( ) ( )h hh h h h

x x x x x x x xx x x xp x y y yx x x x x x x x x x x x

−− −

− ⋅ − − ⋅ −− ⋅ −= + +− ⋅ − − ⋅ − − ⋅ −

x0IIa

x2IIb

f(x)

Ih(f)

x1

p(x)

h h

0 2 0 11 22 0 1 22 2 2

( ) ( ) ( ) ( )( ) ( )( ) ( )2 2

b b b b

h

a a a a

x x x x x x x xx x x xI f p x dx y dx y dx y dxh h h

− ⋅ − − ⋅ −− ⋅ −= = + +−

0 2 0 11 22 0 1 22 2 2

( ) ( ) ( ) ( )( ) ( )( )2 2

x x x x x x x xx x x xp x y y yh h h

− ⋅ − − ⋅ −− ⋅ −= + +−

0 2 0 11 22 2 2

( ) ( ) ( ) ( )( ) ( )2 2

? ? ?b b b

a a a

x x x x x x x xx x x x dx dx dxh h h

− ⋅ − − ⋅ −− ⋅ − = = =−

= ⋅ + ⋅ + ⋅2 0 0 1 1 2 2( ) ( ) ( ) ( )p x y L x y L x y L x


Integração Numérica – Dedução da regra de SimpsonIntroduzindo a variável auxiliar z = x – x1 (z é a distância a x1)

x0IIa

x2IIb

x1

h h

3 2 321 2

2 2 2 2 2

3 32 2 20 2

2 2 2 2 2

0 12

( ) ( ) ( ) 1 1 1 2( ) ( )

2 2 2 2 3 2 2 3 3

( ) ( ) ( ) ( ) 1 1 1 4 4( ) ( )

3 3 3

( ) ( )2

b h h h

ha h h

b h h h

ha h h

x x x x z z h z z h h hdx dz z zh dzh h h h h

x x x x z h z h z h hdx dz z h dz zhh h h h h

x x x x dxh

+

−− −

+

−− −

− ⋅ − ⋅ −= = − = × − = × =

− ⋅ − + ⋅ −= = − = − × − = × =− − −

− ⋅ − =

3 2 3

22 2 2 2

( ) 1 1 1 2( ) ( )

2 2 2 3 2 2 3 3

b h h h

ha h h

z h z z z h h hdz z zh dzh h h h

+

−− −

+ ⋅ = + = × + = × =

0 1 01 1

2 1 2

x x z x x z hz x x x z x

x x z x x z h− = + − = +

= − = + − = + − = −


Integração Numérica – Dedução da regra de SimpsonPelo que

x0IIa

x2IIb

f(x)

Ih(f)

x1

p(x)

h h

= = =

− ⋅ − − ⋅ −− ⋅ −= = + +−

0 2 0 11 22 0 1 22 2 2

/3 4 /3 /3

( ) ( ) ( ) ( )( ) ( )( ) ( )

2 2

b b b b

h

a a a a

h h h

x x x x x x x xx x x xI f p x dx y dx y dx y dxh h h

2 0 1 24

( ) ( )3 3 3

b

h

a

h h hI f p x dx y y y= = + +

0 1 24 ( )

( )6 6 6

( ) ( ) 4 ( )6 2

h

h

b a b a b aI f y y y

b a a bI f f a f f b

− × − −= + +

− + ⇔ = × + × +

Finalmente, atendendo a que h=(b – a)/2, resulta


Integração Numérica – Regras de Newton-Cotes

0 1

0 1

( ) ( ) ( ) ( ) ( ) [ , ,..., , ] ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) [ , ,..., , ] ( ) ( )

n n n nb b b b b

n h n n n

a a a a a

f x p x E x f x p x f x x x x W x

f x d x p x d x E f x d x p x d x f x x x x W x d x

≈ → = − = ⋅

≈ → = − = ⋅

Para calcularmos o erro associado a cada regra de Newton-Cotes podemos integrar o erro daaproximação efectuada

Em face do valor deste integral é possível deduzir uma expressão específica para cada umadas regras

Trapézio:

Ponto médio:

Simpson:

31''( ) ( )

24hE f b aξ= ⋅ ⋅ −

31''( ) ( )

12hE f b aξ= − ⋅ ⋅ −

51''''( ) ( )

2880hE f b aξ= − ⋅ ⋅ −

ξ ∈[ , ]a b


Tabela com algumas regras de Newton-Cotes

Ponto médio ( )= − × +( ) ( ) ( ) / 2hI f b a f a b 3 (2)1( ) ( )

24hE b a f ξ= ⋅ − ⋅

Trapézio

Simpson

3/8 (de Simpson)

Boole

3 (2)1( ) ( )

12hE b a f ξ= − ⋅ − ⋅

5 (4)1( ) ( )

2880hE b a f ξ= − ⋅ − ⋅

7 (6)1( ) ( )

1935360hE b a f ξ= − ⋅ − ⋅

−= + × =, ( )i i ib aa a i f f a

na0 a1 a2 an•••

a bn+1 pontos

[ ] ( )− × + ×= −= +0 1( ) ( )2

( )2h

b aI f fb a f a b ff

( )−= × + +0 1 2( ) 46h

b aI f f f f

5 (4)1( ) ( )

6480hE b a f ξ= − ⋅ − ⋅( )−= × + + +0 1 2 3( ) 3 38h

b aI f f f f f

( )−= × + + + +0 1 2 3 4( ) 7 32 12 32 790h

b aI f f f f f f

21( ) '( )

2hE b a f ξ= ± ⋅ − ⋅Rectânguloà esquerda, à direita = − × = − ×( ) ( ) ( ) , ( ) ( ) ( )h hI f b a f a I f b a f b

Nota: Algumas formulas de ordem superior exibem pesos negativos, facto considerado indesejável


Integração Numérica – Grau de uma regraUma regra diz-se de grau n se integrar sem erro todos os polinómios de grau ≤n e existirpelo menos um polinómio de grau n+1 que não é integrado exactamente.

ξ= − ⋅ − ⋅31( ) ''( )

12hE b a f

Exemplos:

Da análise (da ordem da derivada) da expressão doerro, constata-se que funções de grau 1 (logo comsegunda derivada nula) são integradas sem erro e quefunções de grau 2 (logo com segunda derivada nãonula) são integradas com erro, logo a regra do trapéziotem grau 1

a b

f(x)

Ih(f)

p(x)

Regra do Trapézio (polinómio interpolador de grau 1)

Atendendo ao grau do polinómio interpolador, a regraintegra (pelo menos) funções lineares.


Integração Numérica – Grau de uma regra

Regra do ponto médio (polinómio interpolador de grau 0)

Contudo, da análise (da ordem da derivada) daexpressão do erro, constata-se que funções de grau 1(logo com segunda derivada nula) são integradas semerro e que funções de grau 2 (logo com segundaderivada não nula) são integradas com erro, logo aregra do ponto médio tem grau 1

Exemplos (cont.):

ξ= ⋅ − ⋅3 (2)1( ) ( )

24hE b a f a b

f(x)

Ih(f)

(a+b)/2

a b

f(x)

Ih(f)

(a+b)/2

Atendendo ao grau do polinómio interpolador, a regraintegra (pelo menos) funções constantes.


Integração Numérica – Grau de uma regra

Regra de Simpson (polinómio interpolador de grau 2)

Contudo, da análise (da ordem da derivada) daexpressão do erro, constata-se que funções de grau 3(logo com quarta derivada nula) são integradas semerro e que funções de grau 4 (logo com quartaderivada não nula) são integradas com erro, logo aregra do ponto médio tem grau 3

Exemplos (cont.):

ξ= − ⋅ − ⋅5 (4)1( ) ( )

2880hE b a fa b

f(x)

Ih(f)

(a+b)/2

p(x)Atendendo ao grau do polinómio interpolador, a regraintegra (pelo menos) funções quadráticas.


Dedução alternativa da regra de Simpson

1 2 3( ) ( ) ( ) ( ) (0) ( ) ( )

h

h

h

I f f x d x A f h A f A f h I f

+

−

= ≈ ⋅ − + ⋅ + ⋅ =

Nota: Devido à linearidade do operador integral, se a regra integrarsem erro os monómios 1, x, x2, ..., xn, então integra sem erro todosos polinómios de grau ≤ n

Exercício: a) Deduzir o valor dos pesos Ai de modo à regra seguinte ter o maior grau possível.

Resolução:

Temos 3 incógnitas (A1, A2, A3)

→ necessitamos de 3 equações

2 20 1 2 0 1 2( ) ... 1 ...n n

n n np x dx a a x a x a x dx a dx a x dx a x dx a x dx= + + + + = + + + +

b) Indicar o grau da regra e a expressão do erro.

1 2 3

( ) ( ) ( )

( ) ( ) (0) ( )

h

h

h

I f f x d x

I f A f h A f A f h

+

−

=

= ⋅ − + ⋅ + ⋅

–h h

f(x)

Ih(f)

0

p(x)



( ) 1f x = →

1 2 3 1 2 3

( ) ( ) ( ) 1 ( ) 2

( ) ( ) (0) ( )

h h

h

h

h h

h

I f f x d x d x x h

I f A f h A f A f h A A A

+ +

+

−

− −

= = = = = ⋅ − + ⋅ + ⋅ = + +

1 2 3 2A A A h + + =

( )f x x= →

2

1 2 3

1 2 3

( ) ( ) ( ) ( ) 02

( ) ( ) (0) ( )( ) 0

h h h

hh h

h

xI f f x d x x d x

I f A f h A f A f hA h A A h

+ + +

−− −

= = = == ⋅ − + ⋅ + ⋅

= ⋅ − + ⋅ + ⋅

1 3 1 30A h A h A A − ⋅ + ⋅ = =

2( )f x x= →

32 3

1 2 32 2 2

1 2 3

2( ) ( ) ( ) ( )3 3

( ) ( ) (0) ( )

( ) 0

h h h

hh h

h

xI f f x d x x d x h

I f A f h A f A f h

A h A A h

+ + +

−− −

= = = == ⋅ − + ⋅ + ⋅= ⋅ − + ⋅ + ⋅

2 2 3

1 3 1 32 23 3

A h A h h A A h ⋅ + ⋅ = + =


Dedução alternativa da regra de SimpsonResulta o sistema de 3 equações lineares (a 3 incógnitas)

1 2 3

1 3

1 3

2

23

A A A h

A A

A A h

+ + =

=

+ =

Solução 1 3

2

26

246

hA A

hA

= = = ×

Ou seja,

3( )f x x= →

( ) ( ) ( )

43

3 3

( ) ( ) ( ) ( ) 04

2( ) 4 06

2 4 0 06

h h h

hh h

h

xI f f x d x x d x

hI f f h f f h

h h h

+ + +

−− −

= = = =

= × − + × +

= × − + × + =

Grau da regra de Simpson

Pelo modo como foi construída, a regra tem (pelo menos) grau 2. Terá grau 3?

( ) ( ) ( )

1 2 3( ) ( ) (0) (

2(

)

) 4 06

h

hhI f f h f

I f A f h A f A f h

f h

= ⋅ − +

= − + × +

⋅ + ⋅

( ) 0 ( )hI f I f = = , pelo que tem grau 3



4( )f x x= →

( ) ( ) ( )

5 54

54 4

2( ) ( ) ( ) ( )5 5

2( ) 4 06

2 44 06 6

h h h

hh h

h

x hI f f x d x x d x

hI f f h f f h

h hh h

+ + +

−− −

= = = =

= − + × +

= + × + =

5 52 4( ) ( )5 6 hh hI f I f = =≠

Terá grau 4?

pelo que não tem grau 4, ouseja a regra de Simpson temgrau 3

Qual a expressão do erro?

A aplicação da regra a um polinómio de grau 3 não origina erro, mas a um polinómio degrau 4 já origina. Então a expressão do erro será do tipo, E = C . f(4)(ξ)

Qual o valor de C?

Se f(x)=x4, então f(4)=24, pelo que E = 24C.5 5

52 4 4Por outro lado, ( ) ( )5 6 15hh hE I f I f h= − = − = −

5 54 1Então, 24 (2 )15 2880

C h C h= − = − 5 (4)1resultando, (2 ) ( )2880

E h f ξ= −


Integração Numérica – Regras de Gauss

nós deinterpola ão

1ç

( ) ( ) ( ) ( ) ( ) ( ) ( )

b b N

n i i h

ia a

I f f x d x p x d x A f x I f↑

=

= ≈ = ⋅ = Em Newton-Cotes os nós de interpolação estão definidos “à partida” (nós equidistantes), oque limita o grau de exactidão da regra de integração

Regras de integração de Newton-Cotes

Regras de integração de Gauss - Nas regras de Gauss a posição dos nós de interpolação éescolhida “do melhor modo possível”

Dispomos de 2N parâmetros (os valores dos pesos Ai e a localização dos pontos xi)

Os pesos e a localização são parâmetros

a definir

1

( ) ( )i

i

N

h i i

i Ax

I f A f x=

= ⋅

→ a regra terá grau 2N – 1


Regra de Gauss com 2 pontos

1

1 1 2 2

1

( ) ( ) ( ) ( ) ( ) ( )hI f f x d x A f x A f x I f

+

−

= ≈ ⋅ + ⋅ =

Nota: Devido à linearidade do operador integral, se a regra integrar sem erro os monómios1, x, x2, ..., xn, então integra sem erro todos os polinómios de grau ≤ n

Exercício: a) Deduzir o valor dos pesos Ai e a localização das abcissas xi de modo à regraseguinte ter o maior grau possível.

Resolução:

Temos 4 incógnitas (A1, A2, x1, x2)


2 20 1 2 0 1 2( ) ... 1 ...n n

n n np x dx a a x a x a x dx a dx a x dx a x dx a x dx= + + + + = + + + +

b) Indicar o grau da regra.

1

1

1 1 2 2

( ) ( ) ( )

( ) ( ) ( )h

I f f x d x

I f A f x A f x

+

−

=

= ⋅ + ⋅



( ) 1f x = →

1 1

1

1

1 1

1 1 2 2 1 2

( ) ( ) ( ) 1 ( ) 2

( ) ( ) ( )h

I f f x d x d x x

I f A f x A f x A A

+ +

+

−

− −

= = = = = ⋅ + ⋅ = +

1 2 2A A + =

( )f x x= →

1 1 12

11 1

1 1 2 2 1 1 2 2

( ) ( ) ( ) ( ) 02

( ) ( ) ( )h

xI f f x d x x d x

I f A f x A f x A x A x

+ + +

−− −

= = = = = ⋅ + ⋅ = ⋅ + ⋅

1 1 2 2 0A x A x ⋅ + ⋅ =

2( )f x x= →

1 1 132

11 1

2 21 1 2 2 1 1 2 2

2( ) ( ) ( ) ( )

3 3

( ) ( ) ( ) ( ) ( )h

xI f f x d x x d x


+ + +

−− −

= = = = = ⋅ + ⋅ = ⋅ + ⋅

2 21 1 2 2

2( ) ( )

3A x A x ⋅ + ⋅ =



3( )f x x= →

1 1 143

11 1

3 31 1 2 2 1 1 2 2

( ) ( ) ( ) ( ) 04

( ) ( ) ( ) ( ) ( )h

xI f f x d x x d x


+ + +

−− −

= = = = = ⋅ + ⋅ = ⋅ + ⋅

3 31 1 2 2( ) ( ) 0A x A x ⋅ + ⋅ =

Resulta o sistema de 4 equações não lineares (a 4 incógnitas)

1 2

1 1 2 2

2 21 1 2 2

3 31 1 2 2

2

0

2( ) ( )

3

( ) ( ) 0

A A

A x A x

A x A x

A x A x

+ =

⋅ + ⋅ =

⋅ + ⋅ = ⋅ + ⋅ =

Solução 1 2

1 2

1

13

A A

x x

= =

− = =

Ou seja, = × − + × += ⋅ + ⋅

1 1 2 2( ) ( ) ( )1 1

( ) 1 13 3hh I fI f A f x A ff x f



4( )f x x= →

1 1 154

11 1

4 4

2( ) ( ) ( ) ( )

5 5

1 1 1 1 1 1 2( ) 1 1

9 9 93 3 3 3h

xI f f x d x x d x

I f f f

+ + +

−− −

= = = =

= × − + × = − + = + =

2 2( ) ( )

5 9 hI f I f = =≠

Grau da regra de Gauss com 2 pontos

Pelo modo como foi construída, a regra tem (pelo menos) grau 3.

Terá grau 4?

pelo que não tem grau 4, ou seja a regrade Gauss com 2 pontos tem grau 3

→ as regras de Gauss com N pontos tem grau 2N – 1


Comparação da regra do trapézio com regra de GaussTrapézio (2 pontos)

[ ]( ) ( ) ( )2h

b aI f f a f b−= ⋅ +

a b

f(x)

Ih(f)

a b

f(x)

Ih(f)

x1 x2

Gauss com 2 pontos1 1 2 2( ) ( ) ( )hI f A f x A f x= ⋅ + ⋅

Para [ , ] [ 1, 1]

1 1( ) 1 1

3 3h

a b

I f f f

= − +

= × − + × +

Grau 1

Grau 3


Regras de Gauss-Legendre

Para Gauss-Legendre os pesos Ai e a localização dos pontos xi encontra-se tabelado para ointervalo [a,b]=[–1,+1].

1

( ) ( ) ( ) ( )

b N

i i h

ia

I f f x dx A f x I f=

= ≈ ⋅ =

( )I

1 1

11 1

(

I

)

( ) ( ) ( ) ( ) ( ) ( ) ( )dx

b N

i i h

a

F

di

I f f x dx f x J d F d A F I f

ξ

ξ

ξ ξ ξ ξ ξ ξ+ +

=− −

= = ⋅ = ≈ ⋅ =

Para utilizarmos a informação das tabelas é necessário efectuar uma mudança de variávelpara o intervalo [–1,+1],

Mudança de variável para o intervalo [–1,+1]

1 1( ) ,

2 2 2dx b ax a b Jd

ξ ξξξ

− + −= × + × = =a b x-1 +1 ξ


Regras de Gauss-Legendre no intervalo [-1,+1]

Nº de pontos, N

1

( ) ( )N

h i i

i

I f A F ξ=

= ⋅

O erro associado às formulas de Gauss-Legendre (com N pontos) é,

Abcissas ξi Pesos Ai

1 0 2

2 1 3± 1

3 ±

0

3 / 58 95 9

η η+= × − × = ∈+ ×

42 1 (2 )

3

( !)( ) ( ) , , [ , ]

(2 1) ((2 )!)N N

h N NNE C b a f C a b

N N

4(3 2 6 / 5) / 7

(3 2 6 / 5) / 7

± −

± +

(18 30) 36

(18 30) 36

+

−

ξ ξ+

−

= 1

1

( ) ( )I f F d


Regras de Gauss – regra de Gauss-LobattoAs regras de Gauss são uma família de regras, à qual a regra de Gauss-Legendre pertence.

↑=

= ≈ ⋅ =escolher a melhor

localização poss v l1

í e

( ) ( ) ( ) ( )

b N

i i h

ia

I f f x dx A f x I f

Existem outras regras de Gauss, pertencentes a esta família

Gauss-Legendre

Gauss-Legendre-Lobatto – regra de Gauss-Legendre que inclui os nós extremos do intervalo

↑=

⋅= ≈ × =+ ⋅ +escolher a melhor

loca1

lização possível

( ) ( ) ( )( ) ( ) ( )

b N

i i h

ia

I f f x dx A f x IA a B f b ff

Os coeficientes A, B, Ai e a posição dos pontos xi são parâmetros a determinar

Nota: Para 2 pontos a regra de Gauss-Lobatto é idêntica à regra do trapézio e para 3 pontosé idêntica à regra de Simpson


[ ] [ ]ξ ξ+

−

⋅ − + + ⋅ − + += ≈ + =1

1

10 ( 1) ( 1 ( ) ( )( ) ( ) ( ) ( )) hI f f x d x IA A f ff ff

Exercício: a) Deduzir o valor dos pesos A0 e A1 e a localização da abcissa ξ de modo à regraseguinte ter o maior grau possível.

Resolução:Temos 3 incógnitas (A0, A1, ξ )


b) Indicar o grau da regra.

[ ] [ ]ξ ξ

+

−

=

= ⋅ − + + ⋅ − +

1

1

0 1

( ) ( ) ( )

( ) ( 1) (1) ( ) ( )h

I f f x d x

I f A f f A f f

( ) 1f x = →[ ] [ ]ξ ξ

+ +

+

−

− −

= = = == ⋅ − + + ⋅ − + == ⋅ + + ⋅ + = +

1 1

1

1

1 1

0 1

0 1 0 1

( ) ( ) ( ) 1 ( ) 2

( ) ( 1) (1) ( ) ( )

(1 1) (1 1) 2 2h

I f f x d x d x x

I f A f f A f f

A A A A

1 2 1A A + =

Regra de Gauss-Lobatto com 4 pontos


2( )f x x= →

( )f x x= →

[ ] [ ]

1 1 12

11 1

0 1

0 1

( ) ( ) ( ) ( ) 02

( ) ( 1) (1) ( ) ( )

( 1 1) ( ) 0h

xI f f x d x x d x

I f A f f A f f

A A

ξ ξξ ξ

+ + +

−− −

= = = == ⋅ − + + ⋅ − + == ⋅ − + + ⋅ − + =

0 0 =

[ ] [ ]

1 1 132

11 1

0 1

2 2 20 1 0 1

2( ) ( ) ( ) ( )

3 3

( ) ( 1) (1) ( ) ( )

(1 1) (( ) ) 2 2h

xI f f x d x x d x

I f A f f A f f

A A A A

ξ ξ

ξ ξ ξ

+ + +

−− −

= = = == ⋅ − + + ⋅ − + == ⋅ + + ⋅ − + = + ⋅

2

1 213

A A ξ + ⋅ =

3( )f x x= →

[ ] [ ]

1 1 143

11 1

0 1

3 30 1

( ) ( ) ( ) ( ) 04

( ) ( 1) (1) ( ) ( )

( 1 1) ( ) 0h

xI f f x d x x d x

I f A f f A f f

A A

ξ ξ

ξ ξ

+ + +

−− −

= = = == ⋅ − + + ⋅ − + == ⋅ − + + ⋅ − + =

0 0 =



4( )f x x= →

[ ] [ ]

1 1 154

11 1

0 1

4 4 40 1 0 1

2( ) ( ) ( ) ( )

5 5

( ) ( 1) (1) ( ) ( )

(1 1) (( ) ) 2 2h

xI f f x d x x d x

I f A f f A f f

A A A A

ξ ξ

ξ ξ ξ

+ + +

−− −

= = = == ⋅ − + + ⋅ − + == ⋅ + + ⋅ − + = + ⋅

4

1 215

A A ξ + ⋅ =


Resulta o sistema de 3 equações não lineares (a 3 incógnitas)

0 1

20 1

40 1

1

13

15

A A

A A

A A

ξ

ξ

+ =

+ ⋅ = + ⋅ =

Solução

0

1

16

56

15

A

A

ξ

= = = ±

Ou seja,

[ ] [ ] [ ]0 1( ) ( 1) ( 1) ( ) ( ) ( 1) ( 1)1 5 1 1

5 56 6hI f A f f A f f f f f fξ ξ = ⋅ − + + + ⋅ − + + = ⋅ − + + + ⋅ + −


5( )f x x= →

[ ]

1 1 165

11 1

5/2 5/2

1 5 1 15 56 6

( ) ( ) ( ) ( ) 06

( ) ( 1) ( 1)

1 5 1 1( 1 1) ( ) 0 0 0

6 6 5 5

h

xI f f x d x x d x

I f f f f f

+ + +

−− −

= = = =

= × − + + + × + =

= × − + + × − + =

−

+ =

0 0, ou seja, ( ) ( )hI f I f = =

Grau da regra de Gauss-Lobatto com 4 pontos

Pelo modo como foi construída, a regra tem (pelo menos) grau 4.

Terá grau 5?

pelo que a regra de Gauss-Lobattocom 4 pontos tem grau 5



6( )f x x= →

[ ]

1 1 176

11 1

6/2 6/2 3

1 5 1 15 5

2( ) ( ) ( ) ( )

7 7

( ) ( ) ( 1) ( 1)

1 5 1 1 1 5 2 1 1 26(1 1) ( ) 2

6 6 5 5 6 6 5 3 7

6

7

6

5 5

h h

xI f f x d x x d x

I f I f f f f f

+ + +

−− −

= = = =

= = × − + + + × + =

= × + + × + = × + × = + =

−

2 26( ) ( )

7 75 hI f I f = =≠

Terá grau 6?

pelo que não tem grau 6, ou seja a regra deGauss-Lobatto com 4 pontos tem grau 5

→ as regras de Gauss-Lobatto com N pontos tem grau 2N – 3


O erro associado às formulas de Gauss-Lobatto (com N pontos) é,3 4

2 1 (2 2)3

( 1) (( 2)!)( ) ( ) , , [ , ]

(2 1) ((2 2)!)N N

h N NN N NE C b a f C a b

N Nη η− − − ⋅ − ⋅ −= × − × = ∈

− × −


Regra do trapézio corrigidaO polinómio interpolador pode interpolar também derivada(s) da função. O caso mais usualé considerar uma interpolação de Hermite utilizando para nós os extremos do intervalo [a,b]

Aproximando a função a integrar por um polinómio de Hermite, com informação da funçãoe da primeira derivada, considerando para nós de interpolação os pontos extremos dointervalo resulta,

5 (4)1( ) ( )

720hE b a f ξ= ⋅ − ⋅

A formula possui dois termos: um termo que tem os valores dafunção, e que é idêntico à regra do trapézio, e um termo que têm osvalores das derivadas. O termo das derivadas pode ser entendidocomo uma correcção ao termo que têm os valores da função. Poressa razão a regra designa-se por regra do trapézio corrigida

Regra do Trapézio corrigida (polinómio interpolador de grau 3)

a b

f(x)

Ih(f)

p(x)[ ] [ ]2

( ) ( ) (( )

( ) ( )12

)2

' 'hb aI f f a b a f a ff b b−+ ⋅ −−= ⋅ +

A regra possui grau 3 e a expressão do erro correspondente é,


Regras compostasUm modo de reduzir o erro cometido no cálculo aproximado do integral é subdividir ointervalo [a,b] em N subintervalos e aplicar as regras “básicas” anteriormente estudadas.

Em termos genéricos a regra do trapézio compostapode ser apresentado como

Ex: Regra do trapézio composta com 3 subintervalos iguais (N=3)

[ ] [ ] [ ]≈ ⋅ + + ⋅ + + ⋅ +1 1 20 32( ) ( ) ( ) ( ) ( ) ( )2 2 2h h hf a f f f f fa a aa a

−

=

= + +

1

1

1

1 1( )( ) )( ()2 2

N

i

h f a f bf aI f h

f(x)

a0IIa

a1 a2 a3IIb

h=(b – a)/Nh=

=

= = + + 31 2

0 1 2

( ) ( ) ( ) ( ) ( )

a ba ab

a a a a a

I f f x dx f x dx f x dx f x dx

= + = + +

1 2( ) ( )1 1( ) (

2( ))

2 hf ah If bf a a ff

= + + +

0 3II II

1 21 1( ) ( ) ( ) ( )2 2

a b

af f aah f a f


Regras compostasPara determinarmos o erro cometido podemos somar a contribuição do erro cometido emcada um dos subintervalos.

Resumindo, o erro da regra do rectângulo composta é

ξ ξ ξ= + + = − ⋅ − ⋅ − ⋅ − ⋅ − ⋅ − ⋅0 1 1 2 2 3

3 3 3[ , ] [ , ] [ , ] [ , ] 1 0 1 2 1 2 3 2 3

1 1 1( ) ''( ) ( ) ''( ) ( ) ''( )

12 12 12a b a a a a a aE E E E a a f a a f a a f

ξ ξ−= − ⋅ ⋅ ∈2( ) ''( ) , [ , ]12h

b aE f h a b

f(x)

a0IIa

a1 a2 a3IIb

h = − − = ⋅

( ) /h b a Nb a N h

ξ ξ ξ= − ⋅ ⋅ − ⋅ ⋅ − ⋅ ⋅3 3 31 2 3

1 1 1''( ) ''( ) ''( )

12 12 12h f h f h f

ξ=

=

= − ⋅ ⋅3

3

1

1''( )

12

N

i

i

h f

ξ−= − ⋅ ⋅2( )''( )

12b a h f

ξ= − ⋅ ⋅ ⋅31''( )

12h N f


Tabela com algumas regras compostas

Ponto médio composta −

=

= × + 1

1

( ) ( 2)N

h i

i

I f h f a h ξ−= ⋅ ⋅ 2''( )24h

b aE f h

Trapézio composta

Simpson composta

Trapéziocorrigidacomposta

ξ−= − ⋅ ⋅ 2''( )12h

b aE f h

ξ−= − ⋅ ⋅(4) 4( )2880hb aE f h

−= = + ×, ib ah a a i h

Na0 a1 a2 an•••

a bN intervalos

−

=

= × + +

1

1

1 1( ) ( ) ( ) ( )

2 2

N

h i

i

I f h f a f a f b

−

−

= =

= × + + × + × +

1

1

1 1

( ) ( ) ( ) 2 ( ) 4 ( 2)6

N N

h i i

i i

hI f f a f b f a f a h

ξ−= ± ⋅ ⋅'( )2h

b aE f hRectângulo

à esquerda, à direita composta

−

= =

= × = × 1

1 1

( ) ( ) , ( ) ( )N N

h i h i

i i

I f h f a I f h f a

[ ]−

=

= × + + + −

1 2

1

1 1( ) ( ) ( ) ( ) '( ) '( )

2 2 12

N

h i

i

hI f h f a f a f b f a f b ξ−= ⋅ ⋅(4) 4( )720hb aE f h


Integração com splines – integração com splines cúbicos• Um modo de obter regras de integração semelhante às compostas é utilizando splines.

• A utilização de splines de grau zero conduz às regras do rectângulo compostas, enquantoa integração com spline de grau 1 conduz à regra do trapézio composta.

• A utilização de splines de grau superior conduz a regras diferentes das regras compostasanteriormente estudadas.

Integração com splines cúbicos

No troço i o spline cúbico é dado por

11

( ) ( ) ( ) ( ) ( )i

i

xb b N

i

ia a x

I f f x dx S x dx S x dx I S

−=

= ≈ = =

3 31

1

2 21

1 1

( ) ( )( )

6 6

6 6

i ii i i

i i

i i i ii i i i

i i

x x x xS x M Mh h

h x x h x xy M y Mh h

−−

−− −

− −= + +

− −+ − + −

S1(x)

x0IIa

x1 x2 x3IIb

hi = xi – xi-1hi

S3(x)S2(x)


Integração com splines – integração com splines cúbicosPrimitivando

1 1

3 3 2 21 1

1 1 1( ) ( )

( )6 6 6 6

i i

i i

x x

i i i i i ii i i i i i i

i i i ix x

x x x x h x x h x xS x dx M M y M y M dxh h h h

− −

− −− − −

− − − −= + + − + −

1

4 4 2 2 2 21 1

1 1 1( ) ( ) ( ) ( )

24 24 6 2 6 2

i

i

x

i i i i i ii i i i i i

i i i i x

x x x x h x x h x xM M y M y Mh h h h

−

− −− − −

− − − −= + − + −

− −

3 3 2 2

1 1 124 24 6 2 6 2i i i i i i

i i i i i ih h h h h hM M y M y M− − −

= + + − + −

( ) ( )3

1 12 24i i

i i i ih hy y M M− −= + +−

Somando a contribuição de todos os troços resulta

( ) ( )3

1 1

1

( ) ( )2 24

b N

i ii i i i

ia

h hI S S x dx y y M M− −

=

= = + +

−

Nota: a expressão tem umaparte idêntica à regra dotrapézio composta mais umtermo correctivo com base nos“momentos” (2as derivadas)


Integração adaptativa

Há 2 variantes da integração adaptativa:

•não iterativa – o número de vezes que se efectuam subdivisões é de apenas uma

•iterativa – o número de vezes que se efectuam subdivisões não é definida à partida (éresultado da verificação do critério do erro)

• Considerar uma subdivisão inicial do intervalo [a,b]

• Distribuir a tolerância disponível pelos troços

• Tendo em conta a expressão teórica do erro,estimar para cada troço a correspondente derivada

• Estimar o erro cometido em cada troço

• No caso do erro exceder a tolerância atribuída aesse troço, então subdividir devidamente o troço

ξ ξ= ⋅ ⋅ ≈( ) )( ) (( ) ( ),pk k ki i i if DfE C h

O método tenderá a colocar mais subintervalos onde a correspondente derivada for maior

aIIa0 a1 a2

bIIa3

ε1 ε2 ε3

ε ii

hb a

ε ε=−

≈ ⋅ ⋅( )ki

piE DC h

Integração adaptativa – método que procura que o resultado obtido tenha um erro inferior auma tolerância ε especificada pelo utilizador


Integração adaptativa não iterativaExemplo: Utilizando a regra do trapézio composta

Para o troço [ai-1, ai] de dimensão hi

31 ''( )12h i iE f hξ= − ⋅ ⋅

Na integração adaptativa não iterativa, geralmente, aderivada é aproximada por uma diferença finita apropriada

( )1 1

2

( ) 2 ( 2) ( )''( )2

i i ii i

i

f a f a h f af Dh

ξ − −− ⋅ + +≈ =I I

A estimativa de erro para o troço i é obtida através de

3112i i iE D h= − ⋅ ⋅I I

f(x)

a0IIa

a1 a2 a3IIb

ai-1 ai-1+h/2 ai

hi

[ ]1

1

( ) ( ) ( )2

N

ih i i

i

hI f f a f a−

=

= +


Integração adaptativa não iterativaSe a estimativa de erro |Ei| for superior à tolerância εi que está disponível para esse troço,então esse troço é subdividido em mi subintervalos de modo ao erro nesse troço passar aser inferior à tolerância disponível.

Erro para o troço i após a subdivisão em mi subintervalos

3112i i i iE m D h = × − ⋅ ⋅

I I ai-1 ai

hi

mi subintervalos

i i ih h m=ih

31

12i

i i ii

hE m Dm

= × − ⋅ ⋅

I I 32

1 112

i

i i ii

E

E D hm

= × − ⋅ ⋅

I I2

1i i

i

E Em

= ×

Pretendemos que, após a subdivisão do troço i, o erro nesse troço seja inferior à tolerânciadisponível para esse troço,

2

1i i i i

i

E Em

ε ε< × < 2 ii

i

Em

ε >

Recuperando a expressão do erro para o intervalo i resulta 3112i i i im D h ε> ⋅ ⋅I I

expressão onde se admitiu que aderivada em cada subintervalo éaproximada por Di’’


Integração adaptativa não iterativaA expressão anterior pode ser reescrita em termos da tolerância “total” ε

ε ε

εεε

= ⋅ ⋅− − > = ⋅ ⋅> ⋅ ⋅ −

3

2

3

112

12112

ii

i i

i i ii

i i i i

hD hb a b am D hh

m D h b a

I I

I I

I Iε

− > ⋅ ×

12i i ib am D hI I

1DI I

m1 subintervalos

2DI I

m2 subintervalos

3DI I

m3 subintervalos

a0 a1 a2 a3

aIIa0 a1 a2

bIIa3troço 1 troço 2 troço 3Subdivisão inicial

Subdivisão final

calcular I1 calcular I2 calcular I3

Ih=I1+I2+I3


Integração adaptativa iterativaComparativamente ao método não iterativo, o algoritmo iterativo descrito em seguidapossui as seguintes diferenças:

• se a estimativa do erro for superior à tolerânciapermitida a esse troço, então o troço é divididoao meio

• o número de vezes que se efectuam subdivisõesnão é definido à partida (é resultado daverificação do critério do erro)

• a estimativa de erro é actualizada para os novostroços

• a estimativa do erro não recorre a diferençasfinitas – em cada troço o erro é estimadorecorrendo a 2 aproximações do integral paraesse troço

ai bitroço em avaliação

subdivisão em 2 troços

2 1( )iE I Iα≈ ⋅ −

2 1( )E I Iα≈ ⋅ −

subdivisão em 2 troços

troço em avaliação


Integração adaptativa iterativaExemplo: Calcular I(f) utilizando a regra do trapézioadaptativa iterativa com um tolerância ε = 10–2

Regra do trapézio – dedução da estimativa do erro

3

1 1 2

( ) 1''( )12 1

b aE f ξ−= − ⋅ ⋅

3

2 2 2

( ) 1''( )12 2

b aE f ξ−= − ⋅ ⋅

1

3

1 1'' ''( )

( ) '' 112

D f

b aE I I Dξ≈

− = − ≈ − ⋅ ⋅

32

2

( ) 1''( ) ''( )

12 12h hb ahN

b a b aE f h E fN

ξ ξ−=

− −= − ⋅ ⋅ = − ⋅ ⋅a b

N subintervalos

a b

1 subintervalo

a bc

2 subintervalos

2

3

2 2'' ''( )

( ) 1''12 4

D f

b aE I I Dξ≈

− = − ≈ − ⋅ ⋅

3

2 1( ) 1

(*) (**) '' 112 4

b aI I D− − − ≈ − ⋅ ⋅ −

2 13''

( ) 1112 4

I IDb a

− ≈

− − ⋅ −

Para 1 subintervalo, N=1

Para 2 subintervalos, N=2

Subtraindo as 2 expressões

(*)

(**)

1

0

exp(5 )( ) ( ) , ( )

90xI f f x dx f x= =


Integração adaptativa iterativa

11 1( ) ( )2 2i iI h f a f b = +

Substituindo a aproximação da derivada na expressão do erro para 2 troços

ai bi

1 subintervalo

h

ai bi

2 subintervalos

h/2

ci

h/2

21 1( ) ( ) ( )

2 2 2i i ihI f a f c f b = + +

3

3

2 13

2 12

3

2

( ) 3( )

( )121

112 43 4

( ) 1 4

''

22

1''

4

D

D

I Ib a

I IE

b aE

b ab a

− ≈ −− ⋅ − ≈ ⋅ ⋅ ⋅− ≈ − ⋅ ⋅

−−−−

( )2 2 113

E I I ≈ ⋅ −

Resumindo, para um troço de dimensão h, o valor da regrado trapézio com 1 e com 2 subintervalos é

( )2 2 113

E I I≈ ⋅ −e o erro (para 2 subintervalos) pode ser estimado por



Opção: divisão inicial em 2 troços ii

hb a

ε ε=−

0IIa

1IIb

tolerância ε = 1x10–2

h=1

1/20 tolerância εi = 5x10–3

h=1/2


h=1/2


Integração adaptativa iterativaTroço [0, 1/2] , h=1/2 , tolerância εi = 5x10–3

1

0

exp(5 )( ) ( ) , ( )


0 1/2

1 subintervalo

0 1/21/4

2 subintervalos

11 1 1(0) (1 2) 0.0366182 2 2

I f f = + =

21 1 1(0) (1 4) (1 2) 0.0280044 2 2

I f f f = + + =

( ) 32 2 1 2

10.002871 3 10

3E I I E −= ⋅ − = − ×Estimativa de erro

3 32 3 10 5 10 iE ε− −= × < × = OK


Integração adaptativa iterativaTroço [1/2, 1] , h=1/2 , tolerância εi = 5x10–3

1

0

exp(5 )( ) ( ) , ( )


1/2 1

1 subintervalo

1/2 13/4

2 subintervalos

11 1 1(1 2) (1) 0.4460992 2 2

I f f = + =

21 1 1

(1 2) (3 4) (1) 0.3411644 2 2

I f f f = + + =

( ) 32 2 1 2

10.034978 35 10


ε− −= × > × =

3 32 35 10 5 10 iE

dividir o troço [1 2 , 1] em dois troços



0IIa

1IIb

tolerância ε = 1x10–2

h=1

13/4 εi = 2.5x10–3

h=1/4

3/41/2 εi = 2.5x10–3

h=1/4


h=1/2

erro |E2|= 3x10–311/2 tolerância εi = 5x10–3

h=1/2

erro |E2|= 35x10–3

ii

hb a

ε ε=−


Integração adaptativa iterativaTroço [1/2, 3/4] , h=1/4 , tolerância εi = 2.5x10–3

1/2 3/4

1 subintervalo

1/2 3/45/8

2 subintervalos

11 1 1

(1 2) (3 4) 0.0759774 2 2

I f f = + =

21 1 1

(1 2) (5 8) (3 4) 0.0696008 2 2

I f f f = + + =

( ) 32 2 1 2

10.002126 2.1 10


3 32 2.1 10 2.5 10 iE ε− −= × < × = OK


Integração adaptativa iterativaTroço [3/4, 1] , h=1/4 , tolerância εi = 2.5x10–3

3/4 1

1 subintervalo

3/4 17/8

2 subintervalos

11 1 1

(3 4) (1) 0.2651864 2 2

I f f = + =

21 1 1

(3 4) (7 8) (1) 0.2429268 2 2

I f f f = + + =

( ) 32 2 1 2

10.007420 7.4 10


ε− −= × > × =

3 32 7.4 10 2.5 10 iE

dividir o troço [3 4 , 1] em dois troços



0 1tolerância ε = 1x10–2

1/20 εi = 5x10–3

|E2|= 3x10–311/2 εi = 5x10–3

|E2|= 35x10–3

3/41/2 εi = 2.5x10–3

h=1/4

|E2|= 2.3x10–313/4 εi = 2.5x10–3

h=1/4

|E2|= 7.4x10–3

17/8 εi = 1.25x10–3

h=1/8

7/83/4 εi = 1.25x10–3

h=1/8


Integração adaptativa iterativaTroço [3/4, 7/8] , h=1/8 , tolerância εi = 1.25x10–3

3/4 7/8

1 subintervalo

3/4 7/813/16

2 subintervalos

11 1 1

(3 4) (7 8) 0.0846958 2 2

I f f = + =

21 1 1

(3 4) (13 16) (7 8) 0.08270816 2 2

I f f f = + + =

( ) 32 2 1 2

10.000662 0.7 10

3E I I E −= ⋅ − = − = ×Estimativa de erro

3 32 0.7 10 1.25 10 iE ε− −= × < × = OK


Integração adaptativa iterativaTroço [7/8, 1] , h=1/8 , tolerância εi = 1.25x10–3

7/8 1

1 subintervalo

7/8 115/16

2 subintervalos

11 1 1

(7 8) (1) 0.1582318 2 2

I f f = + =

21 1 1

(7 8) (15 16) (1) 0.15451916 2 2

I f f f = + + =

( ) 32 2 1 2

10.001237 1.24 10


3 32 1.24 10 1.25 10 iE ε− −= × < × = OK



0 1tolerância ε = 1x10–2

1/20 εi = 5x10–3

|E2|= 3x10–311/2 εi = 5x10–3

|E2|= 35x10–3

3/41/2 εi = 2.5x10–3

|E2|= 2.3x10–313/4 εi = 2.5x10–3

|E2|= 7.4x10–3

17/8 εi = 1.25x10–37/83/4 εi = 1.25x10–3

|E2|= 0.7x10–3 |E2|= 1.24x10–3



0 1I[0, 1]= ?

1/20I[0, 1/2]= 0.028004

11/2

3/41/2

I[1/2, 3/4]= 0.069600 13/4

17/8

I[7/8, 1]= 0.154519

7/83/4

I[3/4, 7/8]= 0.082708

Valor obtido para o integral



→ = 0.334832hI

0 1/2 13/4 7/8

Valor exacto ( )1

5 15 5

0

0

e 1 1( ) e e 1 0.327585

90 5 90 450

xxI f dx= = = − =

×Erro efectivo

[0, 1] [0, 1 2] [1 2, 3 4] [3 4, 7 8] [7 8, 1]

0.028004 0.069600 0.082708 0.154519 0.334832

I I I I I= + + + =

= + + + =

efectivo exacto 0.327585 0.334832 0.007247aproximadoE I I= − = − = −

2 2efectivo 0.7 10 1 10 (tolerância)E ε− −× < × =

Valor obtido para o integral


Método de Romberg

[ ]1 2

1

(4) 4

1 1( ) ( ) ( ) ( ) '( ) '( )2 2 12

( )720

N

h i

ih h

h h

hI f h f a f a f b f a f bI I E

b aE I I f hξ

−

=

= × + + + − → = +− = − = ⋅ ⋅

Para a regra do trapézio corrigida composta

[ ]2 4

1,0

( )

1 2(4) 4

1

1 1( ) ( ) ( ) '( ) '( ) ( )2 2 12 720

h

C h hT

N

h h i

i

h b aI I E h f a f a f b f a f b f hξ−

=

− = + = ⋅ + + + − + ⋅ ⋅

O

2 2Se ( ) é possível demonstrar quenf x C +∈

regrado

trap

2 4 6 2 2 2,0 1

éz

2

i

3

o

( )n nh nI T C h C h C h C h h +

↑

= + + + + + +… O

regrado

tr

2 4,0

apézio

1ou seja, ( )hI T C h h↑

= + +O


Método de Romberg0Considere-se uma sequência = e aplique-se a regra do trapézio

2k k

hh

06

22

, 14 ( ) (*)kk k kkh I T C h C h h+→ += + O

2 46

1 1,0 1 2 ( )2 2 2

k k kk k k

h h hh I T C C h+ + = → = + + +

O

Eliminando o termo h2 do erro da aproximação

4 622

21,0 1

1 ( )4

**1 (4

)k kk k C hI T C h h+ = + ++ O

1,0 ,4 6

204 (**) (*)4

4 1 1 )4 (k k k kI I T h hT C+ + − +

−

− = − O

,1

1,0 4 6,2

0 1 ( )4

44 1

k

kk k

k

T

T TC h hI + −

− =

−+

O

1,0 ,0,1

4,

4 1k k

k

T TT + −

=−

Ou seja, a aproximação do integral com erro de ordem h4 é

4 62,1

1 ( )4 kk kI T C h h− += O


Método de Romberg

2 46

2 1,0 1 2 ( )4 4 4

k k kk k k

h h hh I T C C h+ + = → = + + +

O


4 622

21,0 1

1 ( )4

(*14

)k kk k C hT C h hI + + = ++ O

2,0 1,04 6

23 2

14 (**) (* 1) 44 4

4 ( )k k k kI h hI CT T+ + + − +

− =− − O

1,1

2,0 63

1 42

,044

1 ( )41

k

kk

kk

T

T TI C h h

+

+ + − +−

=−

O

De modo análogo ao anterior, considerando agora hk+1 e hk+2

2 46

1 1,0 1 2 ( )2 2 2

k k kk k k

h h hh I T C C h+ + = → = + + +

O

22,0 12

4 624

1 ( )1 **4

(4

)k kk k C hI T C h h+ = + ++ O


2,0 1,01,1

4,

4 1k k

k

T TT + +

+

−=

−1,14 6

23

1 ( )4k k kC hI T h+ −= +O


Método de Romberg


2 2 61,1 ,

2144 (**) ) (4 )(* kk kI I T T h+ − =− +− O

,2

22,0 1, 602

44 1

( )

kT

k kk

T ThI + +−

=−

+

O

Combinando as expressões de Tk,1 e de Tk+1,1, com o intuito de eliminar o termo h4,


21,1 ,1

,2 2

4,

4 1k k

k

T TT + −

=−,2

6( )k kI T h+= O

164

1, 23

14

( *)( *)k k kI T C h h+= +− O

164

, 2 ((14

) *)k k kI T C h h+= − O


Método de Romberg

O procedimento efectuado pode ser generalizado, de modo a eliminar-se os sucessivostermos de h2m, conseguindo-se assim aproximações com erro de ordem h2m+2.

1, 1 , 1,

4,

4 1

mk m k m

k m m

T TT + − −−

=−,

2 2( )mkk mI T h ++= O

A formula de recorrência para Tk,m surge por vezes escrita na forma

1, 1 , 1, 1, 1 4 1

k m k mk m k m m

T TT T + − −

+ −

−= +

−

O método de Romberg é normalmente aplicado com a regra do trapézio, mas também podeser aplicado com outras regras tais como a regra do ponto médio ou de Simpson (esteúltimo caso requereria uma redefinição da formula de recorrência)

A formula de recorrência poderia ter sido deduzida através da formula de Aitken-Neville(tal como se efectuou para o método de Richardson)


Método de RombergFormula de recorrência

1, 1 , 1,

44 1

mk m k m

k m m

T TT + − −−

=−

Tabela

0,2

1

0 0,0

1 1,0

2 2,0

3 3,0

0,3

,2

0,1

1,1

2,1

TT

h T

h T

h T

h

T

T

T

T

T

Erro de ordem h2

Regra dos trapézios

Erro deordem h4

Erro deordem h6

Erro deordem h8

Formula de recorrência


Método de Romberg

Ex: Aplicar método de Romberg (utilizando a regra do trapézio) ao cálculo do integral2

1

0

xI e dx−=

Regra dos trapézios

Opção: iniciar processo com 2 subintervalos Nota: resolução em precisão simples

1

0 0

1

,01 1 1, , ( ) ( ) ( )

2 2 2 2k

N

k ik

i

hh h T h f a f a f b−

=

= = = × + +

0 0,01 1 1 1 10, 2, , (0) (1) 0.73137002 2 2 2 2

k N h T f f f = = = = × + + =

1 1,01 1 1 1 1 3 11, 4, , (0) (1) 0.74298384 4 2 4 2 4 2

k N h T f f f f f = = = = × + + + + =

2,0212, 8, , 0.74586538

k N h T= = = = =…

3 3,013, 16, , 0.7465842

16k N h T= = = = =…

0 11/2

0 11/2 3/41/4

0 11/2 3/41/4

0 11/2 3/41/4


Método de Romberg

Formula de recorrência 1, 1 , 1,

44 1

mk m k m

k m m

T TT + − −−

=−

Tabela

Para m=1

11,0 ,0

,1 1

44 1k k

k

T TT + −

=−

1,0 0,00,1

4 4 0.7429838 0.7313700 0.74685514 1 4 1T T

T× − × −= = =

− −

2,0 1,01,1

4 4 0.7458653 0.7429838 0.74682584 1 4 1T T

T× − × −= = =

− −

3,0 2,02,1

40.7468238

4 1T T

T× −

= = =−

0,1

1,1

0 0,0

1 1,0

2 2,0

3 3,

2

0

,1

0.7468551

0.7

1 2 0.7313700

1 4 0.7429838

1 8 0.7

468258

0.7468238458653

1 16 0.7465842

h T

h T

h T

h

T

T

T

T

== =

= =

=

=

==

=

=


Método de Romberg


44 1

mk m k m

k m m

T TT + − −−

=−

Tabela

Para m=2

21,1 ,1

,2 2

44 1k k

k

T TT + −

=−

1,1 0,10,2

16 16 0.7468258 0.7468551 0.746823816 1 16 1T T

T× − × −= = =

− −

2,1 1,11,2

16 16 0.7468238 0.7468258 0.746823716 1 16 1T T

T× − × −= = =

− −

0,10 0,0

1 1,0

2 2,0

3 3,0

0

1,1

,2

1,2

2,1

0.7468551

0.7

1 2 0.7313 0.7468238

0.74

700

1 4 0.7429838

1 8 0.7458653

1 16 0.7465842

468258

0.7468

6

238

8237

h T

h T

h T

T

T

T

h

T

T

T

= =

= =

= =

= =

=

=

=

=

=


Método de Romberg


44 1

mk m k m

k m m

T TT + − −−

=−

Tabela

Para m=3

31,1 ,1

,3 3

44 1k k

k

T TT + −

=−

1,2 0,20,3

64 64 0.7468237 0.7468238 0.746823764 1 64 1T T

T× − × −= = =

− −

0 0,0

1 1,0

2 2,0

3 3

0,1

1,1

0,2 0,3

1,

0

2

1

,

2,

1 2 0.7313700

1 4 0.74

0.7468551

0.746825

0.7468238

0.74682378

0.746823

29838

1 8 0.7458653

1 16 0.746584

0.746 7

8

23

2

8T

T

h T

h T

h T

h T

T

T

T

T

= =

= =

= =

=

=

=

=

=

=

=

=


Equações não lineares – processo iterativo

Sucessão iterativa: x0, x1, x2, x3, …

zx

f(x)ySeja f(x) uma função e considere-se a equação f(x)=0.A solução da equação designa-se por raiz da equação ou por zero da função (z)

k ke z x= −

0 1 2 3

0 1 3 4

, , , ,, , , , 0 0

k

k k

x x x x z x ze e e e e z x

→ ⇔ →→ ⇔ = − →

Pretendemos que a sucessão de valores xk convirja para a solução z, ou seja, que oerro ek tenda para zero

zx

f(x)

x0

e0

x1x2x3 x4

e3

e2e1


Equações não lineares - multiplicidadeSeja um zero da função ( ). Se ( ) , então tem multipliciade ssemz f x f x C z m∈Propriedade:

( 1) ( )( ) 0 , '( ) 0 , ''( ) 0 , , ( ) 0 , ( ) 0m mf z f z f z f z f z−= = = = ≠

Ex:

zx

f(x)y

zx

f(x)y

zx

f(x)y

m=1zero simples

m=2zero duplo

m=3zero triplo


Equações não lineares – ordem de convergênciaPara avaliar a rapidez de convergência dum método indicamos a ordem de convergência(e a constante de erro assimptótico)

0 1 1

0 1 1

, , , , ,, , , , , 0

k k

k k

x x x x ze e e e

+

+

→→

k ke z x= −

Se a partir duma certa ordem k,

10 , com 1 ,k

k

em M p

e+< ≤ ≤ < ∞ ≥p então o método tem ordem de convergência p

Se existir o limite , 1lim k

kk

e

ec+

→∞=p , então c é a constante de erro assimptótico

Neste caso é usual escrever-se 1k ke ec+ p


Equações não lineares – ordem de convergênciaExemplo: Admitamos que e7 = 4x10–3 (erro na iteração 7)

Hipótese 1) c = 1/2

( )3 38 7

1Hipótese 1a) 1 4 10 2 10

2pp e c e − −= → ⋅ = × × = ×

1

1k ke ec+ p

( )38 7

1Hipótese 1b) 4 10 8 10

2pp e c e − −= → ⋅ = × × = ×

2 62

Hipótese 2) p = 1

( )3 38 7

1 1Hipótese 2a) c 4 10 2 10

2 2pe c e − −= → ⋅ = × × = ×

1

( )38 7

11 1Hipótese 2b) c 4 10 1 10

4 4pe c e − −= → ⋅ = × × = × 3

e7= 4x10–3 e8= 2x10–3

e7= 4x10–3 e8= 8x10–6

e7= 4x10–3 e8= 2x10–3

e7= 4x10–3 e8= 1x10–3


Método da falsa posição – ordem de convergência 1

( )( ) ( )

[ , ]

( )

k k

k kk k

k k

k k

a f af b f af a bb a

b f b

−=−

Equação da recta secante (ou seja do polinómio interpolador) – tabela dif. divididas

( )( ) ( )( ) ( ) k k

k kk k

f b f ap x f a x ab a

−= + ⋅ −−

xk+1

f(bk)

akbk

f(ak)

f(x)

( )( ) ( ) [ , ]k k k kp x f a f a b x a= + ⋅ −

• intervalo [ak,bk] que contém raiz• recta secante (interpolação linear)• determinar intersecção da secante com eixo dos xx• escolher novo intervalo [ak+1,bk+1] que contenha a raiz



Intersecção da recta secante com o eixo dos xx

1 1( ) 0k kx p x+ +→ = ( )1( ) [ , ] 0k k k k kf a f a b x a+ + ⋅ − =

1( )

[ , ]k

k kk k

f ax af a b+ − = −

( )1[ , ] ( )k k k k kf a b x a f a+ ⋅ − = −

1( )

( ) ( )k

k kk k

k k

f ax a f b f ab a

+ = − −−

A expressão anterior pode escrever-se de outro modo,

xk+1

f(bk)

akbk

f(ak)

f(x)

1( )

[ , ]k

k kk k

f ax af a b+ = −

Explicitando a diferença dividida resulta,

1( ) ( )

( ) ( )k k k k

kk k

f b a f a bxf b f a+

⋅ − ⋅=−



a0b0=b1

f(x)

a1=a2 a3

y

b2=b3

+ = −1( )

[ , ]k

k kk k

f ax af a b


Majorante do erro – critério de paragem

Expansão de f(x) em série de Taylor em torno de z (série de ordem 0 com resto de ordem 1)

II0

( ) ( ) ( ) ( ) , inter( , )'k k kf x f z f z x z xξ ξ= + ⋅ − ∈

Ou seja, na iteração k, o erro é limitado por

( ) ( )'k kf x f z xξ = ⋅ −( )( )'k

k

f xz x

f ξ − =

minorante majoranteprimeira primeiraderivada derivada

( ) ( )| | , ( )'

k

k kk

e

f x f xz x m f M

M mξ

↑ ↑

≤ − ≤ ≤ ≤

( ) ( )k kk

f x f xe

M m≤ ≤


Majorante do erro – critério de paragem

O resultado anterior pode ser utilizado num critério de paragem.Se pretendermos um erro inferior a uma tolerância ε,

( )kk k

f xe e

mε ε< ≤ < ( )kf x m ε < ×

ou seja, paramos de iterar quando o valor dafunção for inferior ao produto da tolerância (ε )com o minorante da primeira derivada (m)

f(x)

xk

f(xk)z

ek

Este critério de paragem pode ser aplicado à generalidade dos métodos que abordam aequação na forma f(x)=0


Métodos da falsa posição modificado e da secante

Se no intervalo [ak,bk] a função não alterara sua concavidade (se for côncava ou se forconvexa) o método da falsa posição tende aimobilizar um dos extremos e a convergirlentamente para a raiz

Possíveis alterações ao método da falsa posição:

y

ak-1bk

f(x)

ak ak+1

• método da falsa posição modificado (algoritmo de Illinois)

• método da secante


Método da falsa posição modificado (algoritmo de Illinois)

→ Se um dos extremos se mantiverem duas iterações seguidas, entãodividir (sucessivamente) o valor dafunção desse extremo por 2

1

( ) ( )( ) ( )

k k k kk

k k

f b a f a bx

f b f a+⋅ − ⋅

=−

12

12

ak-1bk

f(x)

ak xk+1

f(bk)

f(bk)/2

1( ) ( )

( ) ( )k k k k

kk k

f b a f a bxf b f a+

⋅ − ⋅=−


Método da secante

• os pontos utilizados para traçar a secante são as duas últimas estimativas

11

( )[ , ]

kk k

k k

f xx xf x x+

−

= −

1 11

1

( ) ( )( ) ( )

k k k kk

k k

f x x f x xxf x f x

− −+

−

⋅ − ⋅=−

−−

−

−=−

11

1

( ) ( )c/ [ , ] k k

k kk k

f x f xf x xx x

Modo alternativo de escrita

Formula de iteração

• o método pode não convergir

• o intervalo pode não conter a raiz

• se convergir a ordem de convergência é de (1 5) / 2 1.618+ = (para zeros simples)

x–1x3

f(x)

x1 x2

y

x0

x4


Ordem de convergência supralinear (>1)Exercício: a) Demonstre que se {xk} for uma sucessão que converge supralinearmente para z,

então +

→ ∞

−=

−1lim 1k k

kk

x xz x

b) Qual é a importância deste resultado do ponto de vista dum critério de paragens de iterações?

+

→ ∞

−−1lim k k

kk

x xz x

a)

+

→ ∞

− + −=

−1lim k k

kk

x z z xz x

+

→ ∞

− ± −=

−1lim k k

kk

z x z xz x

+

→ ∞

−= ±

−1lim 1 k

kk

z xz x

+

→ ∞

−= ±

−11 lim k

kk

z xz x

Pela definição de ordem de convergência, 1, 0 k

k

ek N m M

e+∀ > < ≤ ≤ < ∞p

1 1k k

k k

e z xM M

e z x+ +−

≤ ⇔ ≤−p pEntão, 1k

kk

z xM z x

z x+−

≤ × −−

1p−


Ordem de convergência supralinear (>1)

b)

Para convergência supralinear (p>1),

1lim lim 0 porque e 1 (logo 1 0)kk kk k

k

z xM z x x z p p

z x+

→ ∞ → ∞

−≤ × − = → > − >

−1p−

Resumindo, + +

→ ∞ → ∞

− −= ± = ± =

− −1 1lim 1 lim 1 0 1k k k

k kk k

x x x zz x z x

ou seja, assimptoticamente (i.e., no limite) 1k k ke x x+ −∼

pelo que |xk+1 – xk| representa uma estimativa do erro |ek|

logo se pretendermos um erro inferior a ε (i.e., |ek|< ε )

poderemos terminar o processo iterativo quando a diferença entre duas iterações seguidasfor inferior a ε , ou seja, quando |xk+1 – xk| < ε

Nota: efectuamos o cálculo de xk+1 para que o erro de xk seja inferior a ε (|ek|=|z – xk|)

11lim 1k k

k k k kkk

x xx x z x e

z x+

+→ ∞

−= → − − =

−∼


Método de Newton

=( )

[ , ] ( )( )

'k k

k k k

k k

x f xf x x f x

x f x

• estimativa inicial x0• recta tangente (interpolação linear de Hermite)

• determinar intersecção da tangente com eixo dos xx

Equação da recta tangente (ou seja do polinómio interpolador de Hermite) – tabela dif. div.

( )= + ⋅ −( ) ( ) ( )'k k kp x f x f x x x

x0

x1

yf(x)

x2


Método de Newton

Intersecção da recta tangente com o eixo dos xx

1 1( ) 0k kx p x+ +→ = ( )+ + ⋅ − =1( ) ( ) 0'k k k kf x f x x x

+ − = −1( )( )'

kk k

k

f xx xf x

( )+ ⋅ − = −1( ) ( )' k k k kf x x x f x

A expressão anterior pode ser escrita na forma,

+ = −1( )( )'

kk k

k

f xx xf x

x0

x1

yf(x)

x2

+ = + = −1( )

,( )'

kk k k k

k

f xx x h hf x

(Se convergir,) A ordem de convergência do método de Newton é de 2 (para zeros simples)


Método de Newton – ordem e condições de convergênciaProb:

b) Aproveite para estabelecer uma condição que garanta a convergência.

a,b)

{ }0Seja [ , ], [ , ], considere a sequência gerada pelo método de Newtonkz a b x a b x∈ ∈

a) Demonstre que, se ( ) 0, [ , ], se a sequência convergir para , então aordem de convergência é de 2 (quadrática).

'f x x a b z≠ ∈

Método de Newton 1( )( )'

kk k

k

f xx xf x+ = −

x0

x1

yf(x)

x2

Expansão de f(z) em série de Taylor em torno de xk(série de ordem 1 com resto de ordem 2)

I

2

I0

( )( ) ( ) ( ) ( ) ( )2

inter( , ) [ , ]

'''k k k k

k

ff z f x f x z x z x

z x a b

ξ

ξ

= + ⋅ − + ⋅ −

∈ ⊂


Para haver convergência o erro têm de diminuir (duma iteração para a outra), pelo que

II0

2( )( ) ( ) ( ) ( ) ( )2'''k k k kff z f x f x z x z xξ= + ⋅ − + ⋅ − 2( )( ) ( ) ( ) ( )

2''' k k k kff x z x f x z xξ

⋅ − = − − ⋅ −

2( ) ( ) ( )( ) 2 ( )

''' '

kk k

k k

f x fz x z xf x f x

ξ − = − − ⋅ −

⋅1

2( ) ( ) ( )( ) 2 ( )

''' '

k

k kk k

x

kf x fz x z xf x f x

ξ

+

= − − ⋅ −⋅

21

( ) ( )2 ( )

'''k k

k

fz x z xf x

ξ+ = − ⋅ −

⋅ 21

21

( ) ( )2 ( )

'''

k ke e

k kk

fz x z xf x

ξ

+

+ − = − ⋅ −⋅

21

( )2 ( )

'''k k

k

fe ef x

ξ+ = − ⋅

⋅

221

12k kMe em+ ≤ ⋅

⋅

Tomando o minorante para a primeira derivada e o majorante para a segunda derivada

1

2

0 ( )

( )

'

''

km f x

f Mξ

< ≤

≤

20

1

12M em

≤⋅

(*)

21

12 kk ke eM em+ ⋅

⋅⇔ ≤ ⋅

Condição de convergência (do método de Newton)

Método de Newton – ordem e condições de convergência


1( )

2 ( )'''k k

k

fe ef x

ξ+ = − ⋅

⋅2

Retomando a equação (*), se existir convergência, então

Método de Newton – ordem e condições de convergência

1 ( )lim lim

2 ( )'''

k

k kkk

e ff xe

ξ+

→ ∞ → ∞ =

⋅21 ( )

lim2 ( )

'''

k

kk

e f zf ze

+

→ ∞ =

⋅2

pelo que, para zeros simples, a ordem de convergência do método de Newton é quadrática


Método de Muller – interpolação quadrática

• utilizar as 3 últimas aproximações para efectuar interpolar quadrática (xk–2, xk–1, xk)

• novo ponto é a intersecção da parábola com o eixo dos xx

2 1

1 2 1

( ) ( ) [ , ] ( )[ , , ] ( ) ( )k k k k

k k k k k

p x f x f x x x xf x x x x x x x

−

− − −

= + ⋅ −+ ⋅ − ⋅ −

1 2

1 2 1 2

c/ [ , , ][ , ] [ , ] [ , ]( )

k k k k

k k k k k k k

k k

A f x x xB f x x f x x f x xC f x

− −

− − − −

== + +=

Polinómio interpolador

22( ) ( ) ( )k k k k kp x A x x B x x C= ⋅ − + ⋅ − +

que pode escrever na forma

xk–1xk–2

f(x)

xk

y

p(x)

xk+1

Nota: Para iniciar o processo, terão de ser fornecidas 3 estimativas iniciais (x–2, x–1, x0)


Método de Muller – interpolação quadrática

• das duas intersecções existentes, escolher a mais próxima da última aproximação(ou seja, o sinal ± é escolhido de modo a maximizar o denominador, para que |xk+1–xk| seja mínimo)

• o método permite obter raízes complexas

Nova aproximação (intersecção com xx)

( )1 1

2 2

2

4k

k k

k k k k

Cx xB B A C

+⋅= −

± − ⋅ ⋅

(Se convergir,) A ordem de convergência do método de Muller é de 1.84 (para zeros simples)

xk–1xk–2

f(x)

xk

y

p(x)

xk+1

que se pode escrever na forma

( )+

− ± − ⋅ ⋅= +

⋅

12 2

1

4

2k k k k

k kk

B B A Cx x

A

1 2

1 2 1 2

c/ [ , , ][ , ] [ , ] [ , ]( )

k k k k

k k k k k k k

k k

A f x x xB f x x f x x f x xC f x

− −

− − − −

== + +=


Interpolação quadrática inversa

• utilizar as 3 últimas aproximações para efectuar interpolar quadrática (xk–2, xk–1, xk)

• novo ponto (xk+1) é o valor do polinómio interpolador inverso em y=0

Nota: Para iniciar o processo, terão de ser fornecidas 3 estimativas iniciais (x–2, x–1, x0)

• Se a função tiver inversa, podemos utilizar interpolação inversa, x= f –1(y) x= g(y)

x= f –1(y) x= g(y)

yk–1yk–2

g(x)

yky

p(x)xk

xk–1

xk+1

xk–2

f(x)

x

y

xk–1xk xk–2


Interpolação quadrática inversa

yk–1yk–2

g(x)

yky

p(x)xk

xk–1

xk+1

xk–2

2 1 1 2 1( ) ( ) [ , ] ( ) [ , , ] ( ) ( )k

k k k k k k k k k

x

p y g y g y y y y g y y y y y y y− − − −= + ⋅ − + ⋅ − ⋅ −

1 11

1 1 1

1 1 21 2

2 2 1 1 2

( ) ( ) 1[ , ]

[ , ]

[ , ] [ , ] 1 1 1[ , , ]

[ , ] [ , ]

k k k kk k

k k k k k k

k k k kk k k

k k k k k k k k

g y g y x xg y yy y y y f x x

g y y g y yg y y yy y y y f x x f x x

− −−

− − −

− − −− −

− − − − −

− −= = =− −

−= = − − −

Polinómio interpolador

Atendendo a que

1 2 1 1 2 1(0) [ , ] [ , , ]k k k k k k k k k kx p x g y y y g y y y y y+ − − − −= = − ⋅ + ⋅ ⋅

A nova aproximação é o valor do polinómio em y=0

(Se convergir,) A ordem de convergência da interpolação quadrática inversa é 1.84

A formula iterativa da interpolação quadrática inversa pode escrever-se

1 21

1

1

2 1[1 1

[ [ , ], ] , ]k

k kk k

k k

k k k k k k

y yy y f x x f x x

yx xf x x

−

− − −+

− −

⋅+ − −

= −


Iterativo ponto fixo – condição e ordem de convergência

1

II II

( ) ( ) ( ) ( ) , inter( , ) [ , ]'

kx

k k k

z

g z g x g z x z x a bξ ξ

+

= + ⋅ − ∈ ⊂

1 ( ) ( )'k kz x g z xξ+ = + ⋅ −1

1 ( ) ( )'k k

k

e e

kz x g z xξ+

+ − = ⋅ − 1 ( )'k ke g eξ+ = ⋅

Considerando ( ) , [ , ]'g M x a bξ ≤ ∈

1 ( )'k

k

e ge

ξ+ =

1k

k

eM

e+ ≤podemos escrever

A condição para o processo ser convergente é M < 1 (ou seja, se g(x) for contractiva).

Além disso constatamos que a ordem de convergência é 1 (linear).

Teorema: Se [ , ], ( )I a b g I I= ⊂ e se g(x) for contractiva em I (i.e., se |g’(x)|<1 em I), então

a solução do problema x=g(x) é única e o processo xk+1=g(xk) converge para z seja qual for a

estimativa inicial 0 [ , ]x I a b∈ =

Expansão de g(z) em série de Taylor em torno de xk (expansão de ordem 0 com resto de ordem 1)



|g’(x)|>1 – processo divergente

x0

y=xy=g(x)

z

x1 =g(x0)

x1

x2 =g(x1)

x2

|g’(x)|<1 – processo convergente

x0

y=xy=g(x)

z

x1 =g(x0)

x1

x2 =g(x1)

x2



x0

y=x

y=g(x)

z

x1 =g(x0)

x1

x2 =g(x1)

x2 x0

y=x

y=g(x)

z

x2 =g(x2)

x2

x1 =g(x0)

x1

|g’(x)|>1 – processo divergente|g’(x)|<1 – processo convergente


Iterativo ponto fixo – majorante do erro, critério de paragemExercício:

1 1 / ( ) 1 , [ , ]1

'k k kMe x x c g x M x a bM+ +≤ − ≤ < ∈

−

1

II II

( ) ( ) ( ) ( ) , inter( , ) [ , ]'

kx

k k k

z


+

= + ⋅ − ∈ ⊂

1 ( ) ( )'k kz x g z xξ+ = + ⋅ −1

1 ( ) ( )'k k

k

e e

kz x g z xξ+

+ − = ⋅ − 1 11 ( ) ( )' kk kkz x g z x x xξ+ + + − = ⋅ −− +

1 1 1( ) ( ) ( ) ( )' 'k k k kz x g z x g x xξ ξ+ + + − = ⋅ − + ⋅ − [ ] 1 11 ( ) ( ) ( ) ( )' 'k k kg z x g x xξ ξ+ + − ⋅ − = ⋅ −

1 11 ( ) ( )' 'k k kg z x g x xξ ξ+ +− − = −

Considerando ( ) 1 , [ , ]'g M x a bξ ≤ < ∈

1 1(1 ) k k kM z x M x x+ + − − ≤ −

1

1 11ke

k k kMz x x xM

+

+ + − ≤ −−

pode ser utilizado num critério de paragem do método ponto fixo.

Se pretendermos |ek+1|< ε, então paramos quando

Mostre que no método iterativo de ponto fixo

1 11k k kMe x xM+ +≤ −

−

1 11

1 k k k kM Mx x x xM M

ε ε+ +−− < − <

−



Iterativo ponto fixo – ordem de convergência

1

II II

( ) ( ) ( ) ( ) , inter( , ) [ , ]'

kx

k k k

z


+

= + ⋅ − ∈ ⊂

1 ( ) ( )'k kz x g z xξ+ = + ⋅ −1

1 ( ) ( )'k k

k

e e

kz x g z xξ+

+ − = ⋅ − 1 ( )'k ke g eξ+ = ⋅ 1 ( )'k

k

e ge

ξ+ =

A ordem de convergência é linear (p=1) e

a constante de erro assimptótico é C=|g’(z)|

Quanto mais pequena for a constante de erro assimptótico, mais rápida é a convergência

Uma forma de melhorar a rapidez de convergência é utilizar aceleração de Aitken

1lim lim ( ) ( )' 'k

k kk

eg g z

eξ+

→∞ →∞ = =

1lim k

kkp

eC

e+

→∞=



Iterativo ponto fixo – aceleração de AitkenTomando o anterior resultado 1 ( ) ( ) , inter( , )'k k k k kz x g z x z xξ ξ+− = ⋅ − ∈

1 ( )'k kz x z x+− ≅ ⋅ −G1Considerando ( ) ( )' ' 'k kg gξ ξ +≅ ≅ G

2 1 1 1 1( ) ( ) , inter( , )'k k k k kz x g z x z xξ ξ+ + + + +− = ⋅ − ∈

2 1( )'k kz x z x+ +− ≅ ⋅ −G

2 1 1( )'k k k kx x x x+ + +− ≅ ⋅ −G

–

2 1

1

' k k

k k

x xx x

+ +

+

−

−≅G

(*)

1 ( )'k kz x z x+− ≅ ⋅ −G ( ) 11 ' 'k kz x x+ ⋅ − ≅ − ⋅G G 1

1''

k kx xz + − ⋅ ≅

−GG

Desenvolvendo a equação (*)

Tendo em conta a aproximação para G’

21 1k k kx x x

z

+ +− ⋅

≅

2 1k k k kx x x x+ +− ⋅ + ⋅

1k kx x+ −

1 2 1

1

k k k k

k k

x x x xx x

+ + +

+

− − +−

+ +

+

+

+

+

+

− ⋅ ≅

−−−−−

2 1

1

2 1

1

1

1

k k

k k

k k

k k

k kx xx xx xx

x x

x

z+ − ⋅≅−

1

1''

k kx xz GG


Iterativo ponto fixo – aceleração de Aitken

Contudo

21 2

2 12k k k

k k k

x x xzx x x

+ +

+ +

− ⋅ ≅

− + −

22 1

2 12k k k

k k k

x x xzx x x

+ +

+ +

⋅ − ≅

− +

12 2 2

2 1 12

2 12 2k k k kk k k k kk kkx x x xxx xxx x x x + ++ + + +− ⋅ + +⋅ −−⋅ − = ⋅

( ) ( )2 22 11

212 2k k kk kk k kk x xxx x x x xx+ ++ += ⋅ + + ⋅−+ −− ⋅

( ) ( )22 1 12k k k k k kx x x x x x+ + += − + ⋅ − −

pelo que2

2 1

2 12k k k

k k k

x x xzx x x

+ +

+ +

⋅ −≅− +

( )21

2 12k k

kk k k

x xz x

x x x+

+ +

− ≅ −

− +

Aceleração de Aitken – gerar uma nova sequência x’k (a partirde sequência obtida por iterativo ponto fixo)

( )21

2 1

'2

k kk k

k k k

x xx x

x x x+

+ +

−= −

− +

A nova sequência x’k tem na mesma ordem de convergência 1,mas possui melhor constante de erro assimptótico C=|g’(z)|2

x’3

x0

x1

x2

x3

x4

x5

x’0x’1x’2

… …

sequ

ênci

a po

r pon

to fi

xo

Aitken


Método de Stefensen

As iterações geradas por ponto fixo são efectuadas sobre a aproximação gerada por Aitken

x0

g(x0)

pont

o fix

o

Aitken

g(g(x0))

x1

g(x1)

g(g(x1))

Aitkenx2

g(x2)

g(g(x2))

Aitkenx3

g(x3)

g(g(x3))

Aitkenx4

…

…

A formula de iteração resultante é[ ]

( )

2

1

( )( ) 2 ( )

k kk k

k k k

g x xx x

g g x g x x+

−= −

− ⋅ +

O método de Steffensen pode ser apresentado como um método iterativo de ponto fixo

[ ]( )

2

1

( )( ) c/ ( )

( ) 2 ( )k k

k k k kk k k

g x xx G x G x x

g g x g x x+

−= = −

− ⋅ +

iteração 1 iteração 2 iteração 3


Método de Stefensen

O método de Steffensen também pode ser obtido através da intersecção da recta y=x com arecta secante que passa por 2 pontos gerados pelo método iterativo de ponto fixo

xk

y=xy=g(x)

z

g(xk)

g(xk)

g(g(xk))

xk+1

Pode mostrar-se que:1) Se g’(z)≠1, então z é um ponto fixo atractivo de G(x) e G’(z)=02) Não é necessário que z seja ponto fixo atractivo de g(x) para haver convergência3) Se existir convergência, esta é supralinear (quadrática)


Determinação de todos os zeros de polinómiosSeja uma equação polinomial pn(x)=0.

grau grau

1

1

1( ) ( ) ( )

n n

n np x x z q x↑ ↑

−

−= − ⋅

Uma alternativa é considerar a forma implícita de qn–1(x), i.e. considerar

o problema pode passar a ser expresso na forma qn–1(x)=0 (deflação).

Por vezes é vantajoso desenvolver esta forma implícita – é o caso do método de Newton

11

( )( ) ( ) nn

p xf x q xx z−= =

−

11

( )( ) ( ) nn

p xf x q xx z−= =

−

II

1( )( )'

k

kk k k k

k

h

f xx x x hf x+ = − = +

II

1( )( )'

k

kk k k k

k

h

q xx x x hq x+ = − = +

( )c/( )'

kk

k

q xhq x

= −

Contudo, o cálculo dos coeficientes de qn–1(x) está sujeito a erros de arredondamento.

Após obtermos um primeiro zero z1 ,


Determinação de todos os zeros de polinómios

( )( )'

q xhq x

= −

1

( )( ) p xq xx z

=−

1 ( )( )'q x

h q x = −

21 1

1

( ) ( )( )

( )

'p x p xx z x z

p xx z

−− −= −

−1

( ) 1( )'p xp x x z

= − +−

ou seja, após conhecermos z1, o processo iterativo do método de Newton é1

1 1

1 ( ) 1 ( ) 1c/( ) ( )' 'k k

kk k k k k

p x p xhh p x x z p x x z

−

= − + = − + − − 1k k kx x h+ = +

É possível generalizar o processo e, após conhecermos z1, z2, z3, … , zm

1k k kx x h+ = +1

1 2

( ) 1 1 1c/( )' k

kk k k k m

p xhp x x z x z x z

−

= − + + + + − − −

Nota: É vantajoso determinar os zeros por ordem crescente do seu valor absoluto(devido à extrema sensibilidade dos zeros dos polinómios de grau elevado a pequenas perturbações nos coeficientes)

12 2

1 1 1

( ) ( ) ( ) ( ) ( )( )( ) ( )

' '' p x x z p x p x p xq xx z x z x z

⋅ − −= = −− − −



Cálculo de z1:

= − + − = − +3 2 2( ) 6 11 6 '( ) 3 12 11p x x x x p x x x

+ = +1k k kx x h

Exemplo: Obter o zeros reais de p(x) = x3 – 6 x2 + 11 x – 6 pelo método de Newton conjugadocom a técnica de deflação. Sugestão: partir sempre da origem, x0 = 0

= − ( )( )'

kk

k

p xhxp

=0 0x 0 0 0( ) ( ) 0,54(54( (0) )0)' 'ph p x xp p= −= − =

= + =1 0 0 0,54(54)x x h = − =1 1 1( ) ( ) 0,3034986652'h p x xp

= + =2 1 1 0,8489532106x x h = − =2 2 2( ) ( ) 0,1257208604'h p x xp

= + =4 3 3 0,9990915481x x h

= + =3 2 2 0,974674071x x h −= − = × 23 3 3( ) ( ) 2,441747706 10'h p x xp

−= − = × 44 4 4( ) ( ) 9,072166289 10'h p x xp

= + =5 4 4 0,9999987647x x h −= − = × 65 5 5( ) ( ) 1,235290422 10'h p x xp

= + =6 5 5 1x x h = = → =6 1( ) (1) 0 1p x p z


Determinação de todos os zeros de polinómiosCálculo de z2:

+ = +1k k kx x h

=0 0x1

0

1

00 0 1

'( ) 1 1,2( )

'(0) 1(0) 0 1

p xhx x z

ppp

−− = − + −

= − + =

−

1 0 0 1,2x x h= + =

2 1 1 1,753846154x x h= + =

4 3 3 1,99847524x x h= + =

3 2 2 1,959397304x x h= + = 23 3,907793615 10h −= ×

34 1,52244232 10h −= ×

5 4 4 1,999997682x x h= + = 65 2,317894627 10h −= ×

6 5 5 2x x h= + = 6 2( ) (2) 0 2p x p z= = → =

2

1

3( )Usar deflação, 11

) 1 6( 6p xq xx

x x xz x

− + −=−

=−

( )'( )

kk

k

q xhq x

= −1

( )( )' k

k

xq xq

−

= −

1

1

( ) 1( )' k

k k

p xp x x z

−

= − + −

1

11

1 1 1

'( ) 10,5538461538

( )p xhx x zp

−

= − + = −

1 1z =

1

22

2 2 1

'( ) 10,2055511499

( )p xhx x zp

−

= − + = −



Cálculo de z3:

+ = +1k k kx x h

=0 0x1

00

0 0 1 2

1'( ) 1 1 3( )

'(0) 1 1(0) 0 1 0 2k

pp xhx x z x zp p

−−

= − + + = = − + + − − − −

1 0 0 3x x h= + = 1 3( ) (3) 0 3p x p z= = → =

1

3

2

2( )Usar deflação, ( )( )(

6 11 6( 1)( 2) )

p xq x x xx

xx z xx z

= −−−

+ −=−−

( )'( )

kk

k

q xhq x

= −1

( )( )' k

k

xq xq

−

= −

1

1 2

( ) 1 1( )' k

k k k

p xp x x z x z

−

= − + + − −

1 21, 2z z= =

1 2 31, 2, 3z z z→ = = =


Propriedades

11 12 13 14 11 12 13 14

31 32 33 34 31 32 33 34

41 4

21 22 23 24 21 22 23 24

21 22 23

2 43 44 41 42

24

43 44

0 0 00 0 00 00 0 0

11

11

a a a a a a a a

a a a a a a a aa a a a a a a

a a a a a a a aa a a am m m m

am

=

+ + + +

( ) 1 1 1AB B A− − −=

Exemplo: Adicionar à linha 3, a linha 2 multiplicada por um factor multiplicativo m

2) Inversa do produto

1) Combinação linear de linhas duma matriz – soma de uma linha com outra linhamultiplicada por um factor multiplicativo


Método de Gauss(1) 1 1 1

1 1 22 1 1

A A= = −

(1) (1) (2)M A A=

21

31

12

1 1 1 1 1 1 1 1 1 1 11 1 1 2 1 1 1 2 2 10 1 2 1 1 0 1 2 1 1 1

0 0 0 00 00

10mm

⇔ = = − − −

− −

−−

−−

(2) (2) (3)M A A=

3212

00

1 1 1 11 1 1 1 1 1 10 1 2 10 1 2 1 2 1

30 1 1 1 1 10 1

0 00 00 0

0 0 00

20m

⇔ = = − − − − − − − − − −

ou seja (1)

(1) (1) (2) (2) (1) (3)

(2) (2) (3) M

A AM A A M M A AM A A

== ==

(3)M A A =


Método de Gauss

MA U=

( )

32

1(2) 110 10 1

0 1

0 00

2

001

0

10m

M− =

A(3) é triangular superior, i.e., (3)

0

0 0

1 1 12 1

32

A U = = −

−

Então 1A M U− =

( ) 11 (2) (1)M M M−− = ( ) ( )1 1(1) (2)M M

− −=

( )21

( )

1

1

3

11 1

1 10 1

0 01

0 00 0

0 12mm

M−

= =

( ) ( )1 11 (1) (2)M M M− −− =

21 21

31 32 31 32

0 00 0 0 0 0 11 1 11

12

2

11 0 1 10 1 0 1 1 1

000 0 0m m

m m m m

= = =


Método de GaussM –1 é triangular superior, i.e.,

Então 1A M U A L U−= ⇔ =

Factores Matriz do final damultiplicativos factorização de Gauss

1 1 1 11 1 11 1 2 11 1 2

1 32 1 1 2 12 2

0 00 0

0 0

↑ ↑

⇔ = − − −

21

31 32

1 1111

1

0 00 0

11

1

2

00

2

mm

L

m

M− = = =

Ou seja, através do método de Gauss podemos obter uma factorização A=LU.

Depois de se obter uma factorização A=LU (que requer um número de operações da ordemde n3) o sistema de equações é resolvido mediante uma substituição descendente seguidaduma substituição ascendente (que requerem um número de operações da ordem de n2)


Método de Gauss – escolha de Pivot

− − =

1

2

3

4

4 1 3 0 20 0 1 1 00 2 2 4 40 1 0 1 0

xxxx

Resolva pelo método de Gauss o sistema de equações

⏐ = ⏐ = −

=

=

−

(1) (1) (2) (2)

32

42

40 1 1 02 2 4 41 0

2 01 0

1 3 0 20001 0

A A

m

b

m

b

??

i) Condensação

Troca de linhas (ou seja de equações)

para que o pivot não seja nulo

⏐ = −

−

(2) (2)

2 2 4 40 1

4 1 3 0 2000

1 01 0 1 0

A b

=

⏐ =

−

−

(2) (2)

42

2 2 4 40 1 1 01 0 1 0

4 1 3 0 2

1

0002

A

m

b ⏐ = − − − −

−

(3) (3)

2 2 4 40 1

4 1 3 0 2000

1 00 1 1 2

A b


Método de Gauss – escolha de Pivot

− ⏐ = − − − − = −43

(3) (3) 4 1 3 0 20 2 2 4 4

10 0 1 01 1 20 01 1

A

m

b − ⏐ = − − −

(4) (4) 4 1 3 0 20 2 2 4 40 0 1 1 00 0 0 2 2

A b

− ⋅ = − =4 42 2 1x x

ii) Substituição ascendente

−

= − − −

1

2

3

4

00 00 0

4 1 3 0 22 2 4 4

10

1 02 2

xxxx

− = = =3 4 3 40 1x x x x− ⋅ + ⋅⋅ + ⋅ + ⋅ = = = −3 4

2 3 4 24 (2 4 )

2 2 4 4 12x xx x x x

− − ⋅⋅ + − ⋅ = = =2 31 2 3 1

2 ( 3 ) 34 3 2

4 2x xx x x x

Em aritmética em ponto flutuante, devido aos arredondamentos, em geral a escolha de pivot é vantajosa mesmo quando o pivot não é nulo


Método de Gauss – importância da escolha de Pivot

b) Resolva agora efectuando escolha de pivot (na mesma com uma mantissa com 4 dígitos)

a) Resolva o sistema pelo método de Gauss utilizando uma mantissa com 4 dígitos, ou seja, simulando os cálculos em FP(10,4,2,A)

= −

1

2

0.0003 1.246 1.2490.4370 2.402 1.968

xx

⏐ = − 2

) (

1

(1 1) 0.0003 1.246 1.2490.4370 2.402 1.968

Am

b

Considere o sistema de equações

⏐ =

(2) (2)

(2) (2)

22 2

0.0003 1.246 1.2490 a b

A b

= − × = − − ×(2) (1) (1)22 22 21 12 2.402 1457 1.246a a m a = − −2.402 1815. 422 = −1817. 402

= − × = − ×(2) (1) (1)2 2 21 1 1.968 1457 1.249b b m b = −1.968 1819.793 = −1818. 032= −1.968 1820

= − −2.402 1815

Nota:solução exacta,x1=10, x2=1

i) Condensaçãoa)

= = → =21 210.4370

1456.6(6) 14570.0003

m m


⏐ = − −

(2) (2) 0.0003 1.246 1.2490 1817 1818

A b

−= = → =−2 2

18171.00065 1.001

1818x x


= − −

1

2

0.0003 1.246 1.2490 1817 1818

xx

− ⋅⋅ + ⋅ = = 21 2 1

1.249 1.2460.0003 1.246 1.249

0.0003xx x x − ×= 1.249 1.246 1.001

0.0003

−=11.249 1.247 246x

0.0003 = 0.020.0003

= 6.666(6)

A solução obtida é x1=6.667, x2=1.001, enquanto a solução exacta é x1=10, x2=1.

Comparando com a solução exacta verifica-se que o valor obtido para x1 possui 33% de erro.

= 6.667




− ⏐ = 2

) )

1

(1 (1 0.4370 2.402 1.9680.0003 1.246 1.249

Am

b − ⏐ =

(2) (2)2

(2) (2)

2 2

0.4370 2.402 1.9680 a b

A b

− −= = × → = ×4 421 21

0.00036.8649886 10 6.865 10

0.4370m m

−= − × = − × × −(2) (1) (1) 422 22 21 12 1.246 6.865 10 2.402a a m a = +1.246 0.001648973

= 1.248

−= − × = − × ×(2) (1) (1) 42 2 21 1 1.249 6.865 10 1.968b b m b = −1.249 0.001351 032

= +1.246 0.001649

b) Uma forma de minimizar os problemas numéricos é efectuar escolha de pivot

⏐ = −

(1) (1) 0.0003 1.246 1.2490.4370 2.402 1.968

A b

i) Condensação

− ⏐ =

(1) (1) 0.4370 2.402 1.9680.0003 1.246 1.249

A b

Troca de linhas (ou seja de equações)

= 1.247649

= −1.249 0.001351

= 1.247649 = 1.248

para que o pivot tenha maior

valor absoluto


− ⏐ =

(2) (2) 0.4370 2.402 1.9680 1.248 1.248

A b

= = → =2 21.248

1 1.0001.248

x x


− =

1

2

0.4370 2.402 1.9680 1.248 1.248

xx

+ ⋅⋅ − ⋅ = = 21 2 1

1.968 2.4020.4370 2.402 1.968

0.4370xx x x + ×= 1.968 2.402 1.000

0.4370

+=11.968 2.402

0.4370x = 4.370

0.4370= 10

A solução obtida é x1=10, x2=1, igual à solução exacta



Tipos de Pivot

−

− − − − − −

=

( ) ( ) ( ) ( ) ( ) ( )( ) 1, 1 1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ) ( )1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ), , , ,

( ) ( ) ( ) ( ), , , ,

,

0

0 0

0 0

0 0

k k k k k kk k k p q n

k k k k kk k k k k p k q k n

k k k kk k k p k q k n

k k k kp k p p p q p n

q k

a a a a a aA

a a a a a

a a a a

a a a a

a

( ) ( ) ( ) ( ), , ,

( ) ( ) ( ) ( ), , , ,0 0

k k k kq p q q q n

k k k kn k n p n q n n

a a a

a a a a

Tipos de pivot: - pivot parcial- pivot total- pivot diagonal

- pivot parcial com patamar

submatrizactiva

Os candidatos a pivot encontram-se na submatriz activa


Pivot parcial( ) ( ) ( ) ( ) ( ) ( )( ) 1,1 1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ) ( )1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ), , , ,

( ) ( ) ( ) ( ), , , ,

,

0

0 0

0 0

0 0





q k

a a a a a aA

a a a a a

a a a a

a a a a

a

−

− − − − − −

=

( ) ( ) ( ) ( ), , ,

( ) ( ) ( ) ( ), , , ,0 0

k k k kq p q q q n


a a a

a a a a

pivot parcial – os candidatos a pivot são os elementos da coluna k da submatriz activa

( ) ( )Escolher: max k kik pki ka a

≥→ =

Trocar linha com linha p k→

linha p

linha k

coluna k


Algoritmo pivot parcial

# inicialização# linha pi

para # para to

vot

# escolha do elemento pivot# para todos as entradas aba

1 até

ixo de

das as colunas a conde 1

para 1 até s

nsar

kk

kk

k n

p kpivot a

k

k n ai

=

= −

==

= +

# 1. escolher pivot

# linha pivot

# troca de linhas# se , então trocar linha com a

e então

se entãopara até

linha # para todas as entradas não nulas das linhas a trocar

ik

ik

a pivotp ipivot a

p

k

p k p kkj k n

>

= =

=

≠=

≠

# para todas as linhas abaixo da linha # factor multiplicativo

# pa

para 1 até /

para 1 a ra todos as entradas nté ão nulas dessa

kj

kj pj

pj

ik ik kk

aux aa aa aux

i k nm a a

j k n

k

= = =

= +=

= +

# 2. condensação

linha

ij ij ik kja a m a

= −


Pivot total( ) ( ) ( ) ( ) ( ) ( )( ) 1,1 1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ) ( )1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ), , , ,

( ) ( ) ( ) ( ), , , ,

,

0

0 0

0 0

0 0





q k

a a a a a aA

a a a a a

a a a a

a a a a

a

−

− − − − − −

=

( ) ( ) ( ) ( ), , ,

( ) ( ) ( ) ( ), , , ,0 0

k k k kq p q q q n


a a a

a a a a

pivot total – os candidatos a pivot são todos os elementos da submatriz activa

( ) ( )

,Escolher: max k k

ik pqi j ka a

≥→ =

Trocar linha com linha e trocar coluna com coluna p k q k→

linha p

linha k

coluna qcoluna k


Pivot diagonal( ) ( ) ( ) ( ) ( ) ( )( ) 1,1 1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ) ( )1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ), , , ,

( ) ( ) ( ) ( ), , , ,

,

0

0 0

0 0

0 0





q k

a a a a a aA

a a a a a

a a a a

a a a a

a

−

− − − − − −

=

( ) ( ) ( ) ( ), , ,

( ) ( ) ( ) ( ), , , ,0 0

k k k kq p q q q n


a a a

a a a a

pivot diagonal – os candidatos a pivot são os elementos da diagonal da submatriz activa

( ) ( )Escolher: max k kii qqi ka a

≥→ =

Trocar linha com linha e trocar coluna com coluna q k q k→

linha q

linha k

coluna qcoluna k

Com pivot diagonal, a simetria duma matriz

simétrica é mantida


Pivot parcial com patamar( ) ( ) ( ) ( ) ( ) ( )( ) 1,1 1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ) ( )1, 1 1, 1, 1, 1,

( ) ( ) ( ) ( ), , , ,

( ) ( ) ( ) ( ), , , ,

,

0

0 0

0 0

0 0





q k

a a a a a aA

a a a a a

a a a a

a a a a

a

−

− − − − − −

=

( ) ( ) ( ) ( ), , ,

( ) ( ) ( ) ( ), , , ,0 0

k k k kq p q q q n


a a a

a a a a

pivot parcial – os candidatos a pivot são os elementos da coluna k da submatriz activa

( ) ( )Escolher: max k kik pki ka a

≥→ =

( ) ( )Trocar linha com linha se: , 0 1k kpk kkp k a aτ τ→ > ≤ ≤

linha p

linha k

coluna k

Com patamar, a troca só é efectuada

se “valer a pena”, i.e., se apk for francamente superior a akk

τ é o valor do patamar


Propriedades

1 0 0 00 1 0 00 0 1 00 0 0 1

I

=

(1) A troca de 2 linhas duma matriz A pode ser traduzida pela pré-multiplicação duma matriz de permutação elementar P(e)

(2) A matriz de permutação elementar P(e) pode ser obtida a partir da matriz identidade à qual é efectuada as correspondentes trocas de linhas

Ex: matriz 4x4, matriz elementar de troca das linhas 1 e 3

( )

0 0 00 1 0 0

0 0 01

1

1

0 0 0

eP

=

Linhas1234

3214

Ex: troca das linhas 1 e 3 duma matriz A (4x4)

11 12 13 14 31 32

21 22 23 24 21 22 23 24( )

31 32 33 34

41 42 43 44 41 42 43 44

11 12 1

3 34

3 1

3

4

0 0 00 1 0 0

0 0 00 0 0 11

1

P e

a a a aa a a a a a a a

A P Aa a a aa a a

a a aa

a a a

a

a

a a a a

= = =


Propriedades

(4) As matrizes de permutação elementar P(e) possuem as seguintes propriedades:

( ) ( )) e ei P P I⋅ = ( )( ) ( ))Te eii P P=

(3) A troca de 2 colunas duma matriz A pode ser traduzida pela pós-multiplicação duma matriz de permutação elementar P(e)

Ex: troca das colunas 2 e 3 duma matriz A (4x4)

11 12 13 14 11 14

21 22 23 24 21 24( )

31 3

13

23

2 33 34 31 34

41 42 43

12

22

3

44 41 44

33

43

2

42

1 0 0 00 0 00 0 0

0 11

1

0 0

P e

aa

a a a a a aa a a a a a

A A Pa a a a a

aa

aa a a a

aaaa aa

= = =

1 0 0 00 1 0 00 0 1 00 0 0 1

I

=

Ex: matriz 4x4, matriz elementar de troca das colunas 2 e 3

( )

1 0 0 00 0 00 0 0

1

11

0 0 0

eP

=

Colunas1 2 3 4

1 3 2 4


Propriedades

21(1) (1)

3

31

1

41 41

21

1 0 0 0 1 0 0 0 1 0 0 00 0 0 1 0 0 0 0

)0 0 0 0 1 0 0 00 0 0 1 0 0 1 0 1

1 11 1

0

mi P

mM

m m

mm

− = = − − −

−−

(5) Troca de 2 linhas duma matriz de factores multiplicativos

21(1)

31

41

1 0 0 01 0 00 1 00 0 1

mM

mm

− = − −

Ex: Seja

e considere-se a troca das linhas 2 e 3

Linhas1234

1324

(1)

1 0 0 00 0 00 0 00

10 1

1

0

P

=


Propriedades

2

(1) (1) (1)

41

1

1

3

1 0 0 01 0 00 1 00 0 1

mP M P

m

m = −

−−

ou seja

pelo que P(1)M(1)P(1) corresponde a trocar os factores multiplicativos das linhas 2 e 3

21(1) (1) (1)

31

41

1 0 0 0 1 0 0 0 1 0 0 00 0 0 1 0 0

10 0 0

)0 0 0 0 1 0 0 0 00 0

1

0 1 01

1 0 0 0

1

0 1

mii P M P

mm

− = − −

41

3

21 21

1 3

41

1

1 0 0 0 1 0 0 0 1 0 0 00 0 0 0 0 0 0

0 0 0 0 0 0 00 0 1 0 0 0

11

1 0 01

1

1 11m

mm

m

m

m

−

=

− −=

− −

−


Condensação de Gauss com pivot parcial

(1)

(1 ) (1) (1)

(2) (1) (1 ) (1) (1) (1)

P

P

A A

A P A

A M A M P A

=

=

= =

(1)

0 0 1 00 1 0 01 0 0 00 0 0 1

P =

Ex: Seja A uma matriz (4x4) e considere-se o processo decondensação de Gauss com troca das linhas 1 e 3 na primeirafase da condensação, troca das linhas 2 e 3 na segunda fase etroca das linhas 3 e 4 na terceira fase da condensação

Fase 1 – troca das linhas 1 e 3 seguida de condensação


(2) (1) (1) (1)

(2 ) (2) (2) (2) (1) (1) (1)

(3) (2) (2 ) (2) (2) (1) (1) (1)

P

P

A M P A

A P A P M P A

A M A M P M P A

=

= =

= =

(1) 21

31

41

1 0 0 01 0 00 1 00 0 1

mM mm

− = − −

(2)

1 0 0 00 0 1 00 1 0 00 0 0 1

P =

(2)

32

42

1 0 0 00 1 0 00 1 00 0 1

M mm

= − −

11 12 13 14

21 22 23 24

31 32 33 34

41 42 43 44

a a a aa a a aA a a a aa a a a

=


Condensação de Gauss com pivot parcialAtendendo a que P(2)P(2) = I, podemos reescrever A(3) na forma

(1) (2)

(21)Corresponde atrocar os factoresmultiplicativos de

devido a ,ou seja, a trocar os

coeficientes

(3) (2) (2) (1) (1) (1)

(

das linhas 2 e 3

3) (2) (2) (1) (1(2) (2) ) (1)

M P

P

A M P M P A

A M P M P PP A

==

Linhas1234

3214

3124

(21) (2) (1)

0 0 1 01 0 0 00 1 0 00 0 0 1

P P P = =


(3) (2) (2) (1) (2) (2) (1) (1)

(3 ) (3) (3) (3) (2) (2) (1) (2) (2) (1) (1)

(4) (3) (3 ) (3) (3) (2) (2) (1) (2) (2) (1) (1)

P

P

A M P M P P P A

A P A P M P M P P P A

A M A M P M P M P P P A

=

= =

= =

(3)

1 0 0 00 1 0 00 0 0 10 0 1 0

P =

(3)

43

1 0 0 00 1 0 00 0 1 00 0 1

Mm

= −


Condensação de Gauss com pivot parcialAtendendo a que P(3)P(3) = I, podemos reescrever A(4) na forma

Linhas1234

3214

3124

(321) (3) (2) (1)

0 0 1 01 0 0 00 0 0 10 1 0 0

P P P P = =

Ou seja, a partir do método de Gauss com pivot parcial, podemos obter uma factorizaçãoPA = LU, onde P é a matriz que considera todas as permutações efectuadas e L é a matrizdos factores multiplicativos tendo em consideração as trocas de linhas efectuadas àposteriori

(2) (3)

Corresponde atrocar os factoresmultiplicativos de

devido a ,ou sej

(4) (3) (3) (2) (2) (1) (2) (2) (

a, a trocar oscoeficientes das

linhas 3 e

1) (1)

(4) (3 (3)) (3) (2) (2) (1) (2)

4

(3

M P

A M P M P M P P P A

A M P P P M PPM

==

(1) (3)

(321)Corresponde atrocar os factoresmultiplicativos de

devido a ,ou seja, a trocar os

coeficientes das linhas

) (2(3) ) (1( ) (3

e 4

1) )

3

P

M P

P PP AP

3142

( )( ) ( )(321)

(4) (3) (3) (2) (3) (3) (2) (1) (2) (3) (3) (2) (1) (1)

P PM

A M P M P P P M P P P P P A

=

=

(4)A M P A =

1 (4)M A P A− = LU P A =


Normas de matrizesNormas de vectores

= + + +1 21 nx x x x

( )= + + +1

2 2 2 21 22 nx x x x norma Euclideana

∞ ==

1, ,max ii n

x x

norma 1

norma do máximo (ou do infinito)

Normas de matrizes (mxn)

≤ ≤=

= 1 11

maxm

i jj ni

A a

∞ ≤ ≤=

= 11

maxn

i ji mj

A a

( )= =

=

12

2

1 1

m n

i jF

i j

A a norma de Frobenius


Número de condição (de matrizes)Admitindo que existem perturbações nos valores da matriz A e do vector b, então resultam perturbações na solução do sistema x

[ ] { } { }⋅ = → ⋅ =A x b A x b ( ) ( )δ δ δ⋅ = ⇔ + ⋅ + = + A x b A A x x b b

(i) Admitir que apenas existem perturbações no 2º membro (e consequentemente na solução)

( )δ δ⋅ + = +A x x b b δ δ ⋅ + ⋅ = +A x A x b bδ δ⋅ =

⋅ =

A x bA x b

⋅ =A x b = ⋅ ≤ ⋅b A x A x ≤b

xA

(*)

δ δ δ δ−⋅ = = ⋅1A x b x A b δ δ δ− − = ⋅ ≤ ⋅1 1x A b A bδ δ− ≤ ⋅1x b

Ax x


Número de condição (de matrizes)

Tendo em atenção (*),

(ii) Perturbações na matriz A (e consequentemente na solução)

δ δ δ− −≤ ⋅ ≤ ⋅1 1x b bA A

bx xA

≤b

xA

δ δ− ≤ ⋅ ⋅

con

1

d A

x bA A

x bδ δ

≤ ⋅cond x b

Ax b

−= ⋅ 1cond A A AO número de condição duma matriz traduz, em termos relativos,a relação entre as perturbações na solução x e as perturbaçõesno segundo membro b.

Analogamente se demonstra que a relação entre as perturbações na solução x e asperturbações da matriz A também dependem do número de condição da matriz

Um número de condição elevado indica que as perturbações do segundo membro sãoampliadas sobre a solução do sistema


Efeito dos erros de arredondamento

(i) Pivot parcial – em certos casos patológicos γ pode ser muito elevado podendo atingir ovalor máximo de 2n – 1. Contudo, estes casos patológicos são raros e na prática a factorizaçãocom pivot parcial é em geral numericamente estável

Na resolução dum sistema (de dimensão n) em ponto flutuante, devido aosarredondamentos, a factorização obtida não é exactamente igual à matriz original

Pode demonstrar-se que os elementos da matriz erro são majorados por:

matrizdos

erros

L U A E↑

⋅ = +

1ije n u γ α≤ ⋅ ⋅ ⋅1

- constante da ordem da unidade de arredondamento - maior elemento (em módulo) de

- factor de crescimento dos coef. de durante factorizaçãoij

uAA

αγ

(ii) Pivot total – o majorante de γ cresce lentamente (com o aumento da dimensão dosistema) não se conhecendo casos para os quais seja superior a n. Logo a utilização de pivottotal é numericamente estável.


Efeito dos erros de arredondamentoIntroduzindo o conceito de resíduo,

Atendendo a que

r b A x= − ⋅

onde o número de condição surge novamentecomo factor de ampliação

( )r b A x A x A x A x x A xδ= − ⋅ = ⋅ − ⋅ = ⋅ − = ⋅

1r A x x A rδ δ −= ⋅ = ⋅ δ − = ⋅1x A r 1x A rδ − ≤ ⋅

(de fácil cálculo após se obter )x

1x rA

x xδ − ≤ ⋅

bA x b A x b x

A⋅ = ⋅ ≥ ≥

Então 1 1x r x rA A

bx x xA

δ δ− −≤ ⋅ ≤ ⋅δ − ≤ ⋅ ⋅

con

1

d A

x rA A

x b

δ≤ ⋅cond

x rA

x bOu seja

Resumindo, o número de condição da matriz desempenha um papel fundamental nos errosexistente na solução do sistema de equações


Matriz em banda0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 00 0 0 0 0

0 0 0 0 0 0 00 0 0 00 0 0 0 0

0 00 0 0 0 0 0 00 0 0 00 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0

XX X

X XX X XX X X

X X X X XX X

X X X XXX X X X

X XX X X X

X X XX X X X

XX X X X X

X XX

XX

XX

XX

XX

XX

X

XX XX

X X XX

A

XXXX X X

=

largura de banda superior: número de diagonais não nulas, acima da diagonal principal

largura de banda inferior: número de diagonais não nulas, abaixo da diagonal principal

largura de banda = largura de banda superior + largura de banda inferior + 1ou seja, largura de banda é o numero total de diagonais não nulas

Documents

matematica computacional