Matemática Numérica I

MATEMÁTICA NUMÉRICA I

TEXTOS DE APOIO

JOSÉ A. FERREIRA E GONÇALO PENA

DEPARTAMENTO DE MATEMÁTICA

UNIVERSIDADE DE COIMBRA

2015/2016

CONTEÚDO

1 Erros. Propagação de Erros. 11.1 Erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Propagação de erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3 Propagação de erros em sistemas lineares . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.3.1 Alguns resultados matriciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3.2 Número de condição de uma matriz . . . . . . . . . . . . . . . . . . . . . . . . 13

1.4 Alguns problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2 Métodos Numéricos para Sistemas de Equações Lineares 172.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2 Métodos Directos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2.1 Método de eliminação de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2.2 Métodos de factorização . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.3 Métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.3.1 Métodos Estacionários . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352.3.2 Métodos não Estacionários . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 462.3.3 Algumas considerações finais . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

2.4 Alguns problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

3 Interpolação polinomial 653.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.2 Polinómio interpolador de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . 673.3 Polinómio interpolador de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 733.4 Polinómio interpolador de Lagrange segmentado . . . . . . . . . . . . . . . . . . . . . 793.5 Polinómio interpolador de Hermite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 843.6 Splines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 933.7 Interpolação de funções definidas em R2 . . . . . . . . . . . . . . . . . . . . . . . . . . 1013.8 Alguns problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

4 Equações não lineares 1134.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1134.2 O método da bissecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1154.3 O método de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1184.4 O método do ponto fixo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1224.5 O método da falsa posição . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1284.6 Alguns problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

CAPÍTULO 1ERROS. PROPAGAÇÃO DE ERROS.

1.1 Erros

Consideremos a expressão clássica de Newton para a energia cinética de um corpo com massa m0

KN =12

m0v2,

em que v denota a velocidade do corpo. A energia cinética KN é, em certas condições, uma aproximaçãopara a energia cinética definida por Einstein

KE = mc2−m0c2,

em que m0 representa a massa do corpo em repouso, c a velocidade da luz e m a massa relativa do corpoem movimento que é dada por

m =m0√1− v2

c2

.

De facto, suponhamos que v c. Comecemos por notar que

KE = m0c2

1√

1− v2

c2

−1

= m0c2

(1√

1+ x−1), com x =−v2

c2 .

O Teorema de Taylor constitui a ferramenta fundamental para estabelecer o resultado pretendido.

1

2 MATEMÁTICA NUMÉRICA I J. A. FERREIRA, G. PENA

Teorema 1.1 (Teorema de Taylor). Seja f : I −→ R uma função com derivadas até à ordem ncontínuas num intervalo aberto I contendo a e que admite derivada de ordem n+1 nesse intervalo.Então, para x ∈ I, existe η(x) ∈ I tal que

f (x) = Pn(x;a)+Rn(x;a), (1.1.1)

em que o polinómio de Taylor e o resto de Lagrange são definidos, respectivamente, por

Pn(x;a) =n

∑j=0

f ( j)(a)j!

(x−a) j, Rn(x;a) =f (n+1)(η(x))(n+1)!

(x−a)n+1.

Pelo Teorema de Taylor, com a = 0 e n = 1, temos

1√1+ x

= 1− 12

x+38

1√(1+η(x))5

x2.

Assim

KE =12

m0v2 +38

m0c2 1√(1+ η(v))5

v4

c4 ,

em que η(v) = η(− v2

c2

)∈(− v2

c2 ,0). Logo KE ' KN e o erro cometido nesta aproximação é dado por

Erro =38

m0c2 1√(1+η(v))5

v4

c4 .

Este erro, que é usualmente designado por erro de truncatura, satisfaz

Erro≤ 38

m01√(

1− v2

c2

)5

v4

c2 .

e constitui o erro cometido ao substituir a expressão da energia cinética de Einstein pela energia cinéticade Newton.

É frequente a substituição de modelos matemáticos complexos por modelos mais simples sendo, nestecaso, necessário indicar o erro que se comete nesta substituição. Um outro exemplo significativo dasubstituição anterior ocorre no contexto dos modelos diferenciais.

Consideremos o seguinte problema diferencial de condição inicial

u′(t) = f (t,u(t)), t ∈ (0,T ],u(0) = u0.

(1.1.2)

O conjunto de problemas do tipo anterior que admitem uma solução explícita é bastante reduzidoe a única possibilidade de descrever, pelo menos de modo aproximado, a dinâmica de tais sistemas, érecorrer a modelos discretos. Seja ti, i = 0, . . . ,N uma partição uniforme de [0,T ] com t0 = 0, tN = T eti− ti−1 = h. Pelo Teorema de Taylor temos

u′(ti) =1h

(u(ti+1)−u(ti)

)− h

2u′′(ξi), ξi ∈ (ti, ti+1). (1.1.3)

A relação anterior permite-nos substituir o problema de condição inicial (1.1.2) pelo modelo discretoseguinte

ui+1 = ui +h f (ti,ui), i = 0, . . . ,N−1,u0 dado.

(1.1.4)

J. A. Ferreira, G. Pena Matemática Numérica I 3

Este modelo define recursivamente ui, i = 1, . . . ,N, tendo-se u(ti) ' ui, i = 1, . . . ,N. Vejamos seguida-mente qual o erro que se comete ao substituir o modelo contínuo (1.1.2) pelo modelo discreto (1.1.4).

Consideremos, em (1.1.4), ui substituído por u(ti). Obtemos

1h

(u(ti+1)−u(ti)

)= f (ti,u(ti))+Erro, (1.1.5)

em que Erro representa o erro de truncatura.Combinando as igualdades (1.1.3) e (1.1.5), segue-se que

u′(ti)+12

hu′′(ξi) = f (ti,u(ti))+Erro,

e uma vez que u satisfaz (1.1.2), concluímos

Erro =12

hu′′(ξi).

Atendendo a que u′(t) = f (t,u(t)), então, com f = f (x,y), vem

u′′(t) =ddt

f (t,u(t)) =∂ f∂x

(t,u(t))+∂ f∂y

(t,u(t))u′(t),

ou ainda,

u′′(t) =∂ f∂x

(t,u(t))+∂ f∂y

(t,u(t)) f (t,u(t)).

Deste modo, admitindo que ∣∣∣∣∂ f∂x

∣∣∣∣≤M,

∣∣∣∣∂ f∂y

∣∣∣∣≤M e | f | ≤M,

deduzimos|Erro| ≤ 1

2hM(1+M).

A estimativa anterior permite concluir que o erro de truncatura tende para zero quando o diâmetro dapartição em [0,T ] tende para zero.

Retomemos a lei de Newton para a energia cinética de um corpo em movimento. Nesta lei m0representa a massa do corpo em repouso que é determinada utilizando um instrumento de medida.Deste modo, o valor tomado para m0 é, em geral, apenas aproximado m0 estando afetado de um errodesconhecido. À quantidade

∆m0 = m0−m0

chamamos erro absoluto de m0. É claro que tal quantidade só nos dá informação significativa quandocomparado com a quantidade exacta m0. Surge deste modo o conceito de erro relativo

δm0 =∆m0

m0.

Observamos que

δm0 =∆m0

m0 +∆m0,

em que, em geral, ∆m0 é pequeno quando comparado com m0. É conveniente, em termos práticos,considerar |∆m0|

|m0| como indicador do erro relativo |δm0|, isto é,

|δm0| '|∆m0||m0|

.


Uma questão natural que se coloca diz respeito à qualidade do resultado obtido quando utilizamos, porexemplo, a expressão da energia cinética de Newton com um valor aproximado para a massa m0 e para avelocidade v e são conhecidos majorantes ε1 e ε2 dos seus erros absolutos ∆m0 e ∆v, respetivamente.

Notamos que se tem

∆KN = KN(m0,v)−KN(m0,v).

Atendendo a que vale a seguinte representação

∆KN = KN(m0,v)−KN(m0,v)+KN(m0,v)−KN(m0,v),

obtemos

∆KN =∂KN

∂m0(m∗0,v)∆m0 +

∂KN

∂v(m0,v∗)∆v,

em que m∗0 ∈ (m0,m0) ou m∗0 ∈ (m0,m0) e v∗ ∈ (v,v) ou v∗ ∈ (v,v).Como m∗0 ∈ [m0− ε1,m0 + ε1] e v∗ ∈ [v− ε2,v+ ε2], deduzimos para |∆KN | a seguinte estimativa

|∆KN | ≤maxΩ

∣∣∣∣∂KN

∂m0

∣∣∣∣ε1 +maxΩ

∣∣∣∣∂KN

∂v

∣∣∣∣ε2,

em que Ω = [m0− ε1,m0 + ε1]× [v− ε2,v+ ε2].

Podemos facilmente deduzir uma estimativa para o erro relativo δKN , análoga à anterior, sendoconhecidos majorantes δi, i = 1,2, dos erros relativos δm0, δv, respetivamente. Pode provar-se que vale aseguinte estimativa

|δKN | ≤maxΩ

∣∣∣∣m0

KN

∣∣∣∣maxΩ

∣∣∣∣∂KN

∂m0

∣∣∣∣δ1 +maxΩ

∣∣∣∣v

KN

∣∣∣∣maxΩ

∣∣∣∣∂KN

∂v

∣∣∣∣δ2, (1.1.6)

em que, como anteriormente, Ω = [m0− ε1,m0 + ε1]× [v− ε2,v+ ε2], em que εi, i = 1,2, são majorantesdos erros absolutos |∆m0| e |∆v|.

Introduzimos até ao momento alguns conceitos envolvendo alguns tipos de erros bem como a suapropagação. Os conceitos básicos na teoria de erros são os conceitos de erro absoluto e erro relativo.Estes erros são usualmente designados por erros de observação pois estão associados à medição de umadeterminada quantidade. Um outro tipo de erros surge quando efetuamos arredondamentos. De factoconsideremos um determinado número

x = aα1 . . .αtαt+1 . . . ,

em que a ∈ N0, αi ∈ 0,1, . . . ,9. Admitamos que efetuamos um arredondamento no dígito αt+1 seguidode truncatura. Se αt+1 ≤ 4, então obtemos

x = a.α1 . . .αt .

Por outro lado, se αt+1 ≥ 5, então

x = a.α1 . . .(αt +1).

Em qualquer uma das duas situações anteriores,

|∆x| ≤ 5×10−(t+1).


1.2 Propagação de erros

Consideremos y = g(x1, . . . ,xn) em que g : Rn −→R admite derivadas parciais contínuas num determi-nado conjunto que especificamos seguidamente. Seja xi ' xi, i = 1, . . . ,n e admitamos que são conhecidosmajorantes εi, i = 1, . . . ,n, para os erros absolutos ∆xi, i = 1, . . . ,n. Determinamos seguidamente umaestimativa para o erro absoluto ∆y. Observamos que apenas vamos estender a expressão obtida na secçãoanterior no contexto da energia cinética de Newton.

Tem-se sucessivamente

∆y = y− y

= g(x1, . . . ,xn)−g(x1, . . . ,xn)

=n

∑i=1

∂g∂xi

(Xi)∆xi

em que Xi ∈Ω =n

∏j=1

[x j− ε j,x j + ε j], i = 1, . . . ,n. Logo

|∆y| ≤n

∑i=1

maxΩ

∣∣∣∣∂g∂xi

∣∣∣∣εi.

No que diz respeito ao erro relativo, admitamos agora que são conhecidos majorantesδi, i = 1, . . . ,n, dos erros relativos δxi, i = 1, . . . ,n. Atendendo a que se tem sucessivamente

δy =∆yy

=1

g(x1, . . . ,xn)

(g(x1, . . . ,xn)−g(x1, . . . ,xn))

=1

g(x1, . . . ,xn)

n

∑i=1

xi∂g∂xi

(Xi)δxi,

deduzimos

|δy| ≤n

∑i=1

maxΩ

∣∣∣∣xi

g(x1, . . . ,xn)

∣∣∣∣maxΩ

∣∣∣∣∂g∂xi

∣∣∣∣δi,

em que, como anteriormente, Ω =n

∏j=1

[x j− ε j,x j + ε j] e |∆x j| ≤ ε j, j = 1, . . . ,n.

1.3 Propagação de erros em sistemas lineares

A propagação do tipo de erros introduzidos é relevante no contexto dos sistema lineares. De factoconsideremos o sistema

Ax = b

e seja b um vetor cujas componentes foram obtidas por medição. Neste caso o erro absoluto do vetor dasaproximações ∆b é um vetor. Do mesmo modo, sendo x tal que

Ax = b,

o erro absoluto desta solução aproximada é o vetor ∆x cujas componentes são os erros absolutos dascomponentes xi. A fim de “medir” o vetor ∆x necessitamos de introduzir o conceito de norma e estudaralguns resultados matriciais.


1.3.1 Alguns resultados matriciais

Definição 1.1. Seja V um espaço vectorial real. A aplicação ‖·‖ : V −→ R+0 tal que

• ‖v‖ ≥ 0, ∀v ∈V, ‖v‖= 0 =⇒ v = 0,

• ‖αv‖= |α|‖v‖, ∀v ∈V,∀α ∈ R,

• ‖v+w‖ ≤ ‖v‖+‖w‖, ∀v,w ∈V,

é dita norma em V .

Exemplo 1.1. As aplicações

• ‖x‖∞ = maxi=1,...,n |xi|,

• ‖x‖1 = ∑ni=1 |xi|,

• ‖x‖2 =

(∑

ni=1 x2

i

)1/2

,

são normas em Rn.

Definição 1.2. Num espaço vectorial V , duas normas ‖·‖ e ‖·‖v dizem-se equivalentes se existemconstantes positivas K1 e K2 tais que

K1‖w‖ ≤ ‖w‖v ≤ K2‖w‖, ∀w ∈V.

Observamos que as normas ‖·‖1,‖·‖2 e ‖·‖∞ são equivalentes em Rn. Certas normas de uso frequentesão induzidas por produtos internos, cuja definição relembramos de seguida.

Definição 1.3. Seja V um espaço vectorial real ou complexo. A aplicação (·, ·) : V ×V −→K, comK= R ou K= C tal que

• (v,v)≥ 0, ∀v ∈V, (v,v) = 0 =⇒ v = 0,

• (v,w) = (w,v), ∀v,w ∈V,

• (αv,w) = α(v,w),(v+w,r) = (v,r)+(w,r), ∀v,w,r ∈V, ∀α ∈K,

é dita produto interno em V.

No espaço vectorial real Rn, a aplicação

(x,y) = xty =n

∑i=1

xiyi, x,y ∈ Rn,

é um produto interno. No espaço vectorial complexo Cn, a aplicação

(x,y) = xty =n

∑i=1

xiyi, x,y ∈ Cn,

é um produto interno.


Seja V um espaço vectorial com produto interno (·, ·). A aplicação

‖v‖=√(v,v), v ∈V,

é uma norma.

Exemplo 1.2. Consideremos agora o espaço das matrizes de ordem n com entradas em R que denotamospor Mn(R). As aplicações seguintes:

‖A‖1 = maxj=1,...,n

n

∑i=1|ai j|, A = [ai j]

ni, j=1 ∈Mn(R), (1.3.1)

‖A‖∞ = maxi=1,...,n

n

∑j=1|ai j|, A = [ai j]

ni, j=1 ∈Mn(R), (1.3.2)

‖A‖F =

( n

∑i, j=1

a2i j

) 12

, A = [ai j]ni, j=1 ∈Mn(R), (1.3.3)

‖A‖= maxi, j=1,...,n

|ai j|, A = [ai j]ni, j=1 ∈Mn(R), (1.3.4)

são normas em Mn(R). A norma ‖·‖F é dita norma de Frobenius e admite a representação

‖A‖F =√

tr(AAt),

em que tr(AAt) representa o traço de AAt .

Definição 1.4. Uma norma ‖·‖m em Mn(R) diz-se compatível ou consistente com uma normavectorial ‖·‖v em Rn se

‖Ax‖v ≤ ‖A‖m‖x‖v, ∀A ∈Mn(R),∀x ∈ Rn.

Observamos que a norma de Frobenius é compatível com a norma ‖·‖2. De facto tem-se sucessiva-mente

‖Ax‖22 =

n

∑i=1

( n

∑j=1

ai jx j

)2

≤n

∑i=1

n

∑j=1

a2i j

n

∑j=1

x2j

=

( n

∑i=1

n

∑j=1

a2i j

)‖x‖2

2

= ‖A‖2F‖x‖2

2.

Definição 1.5. Se a norma matricial ‖·‖ em Mn(R) satisfaz

‖AB‖ ≤ ‖A‖‖B‖, ∀A,B ∈Mn(R),

então diz-se que ‖·‖ é submultiplicativa.


Em Mn(R), as normas ‖·‖1 e ‖·‖∞ são submultiplicativas. Contudo a norma (1.3.4) não tem estapropriedade. De facto, sejam A = [ai j]

ni, j=1 e B = [bi j]

ni, j=1 definidas por

ai j = 1, bi j = 1, i, j = 1, . . . ,n.

Temos ‖A‖= ‖B‖= 1 e ‖AB‖= n.Introduzimos seguidamente o conceito de norma matricial induzida por uma norma vectorial ou

subordinada a uma norma vectorial.

Definição 1.6. A aplicação

‖A‖= sup06=x∈Rn

‖Ax‖v

‖x‖v, A ∈Mn(R),

em que ‖·‖v é uma norma vectorial em Rn, é uma norma em Mn(R) que é designada norma matricialinduzida pela norma vectorial ‖·‖v (ou subordinada a uma norma vectorial ‖·‖v).

De facto, é manifesto que‖A‖ ≥ 0, ∀A ∈Mn(R),

e se A ∈Mn(R) é tal que ‖A‖= 0, então ‖Ax‖v = 0, para todo x em Rn. Logo

Ax = 0, ∀x ∈ Rn.

Assim, escolhendo x de modo adequado, concluímos A = 0.Verifica-se facilmente que se tem

‖αA‖= |α|‖A‖, ∀α ∈ R, ∀A ∈Mn(R),

e

‖A+B‖= sup06=x∈Rn

‖(A+B)x‖v

‖x‖v

≤ sup06=x∈Rn

‖Ax‖v

‖x‖v+ sup

06=x∈Rn

‖Bx‖v

‖x‖v

= ‖A‖+‖B‖.

A representação‖A‖= sup

‖x‖v=1‖Ax‖v, A ∈Mn(R),

para a norma induzida é equivalente à introduzida anteriormente.Observamos que a norma induzida é submultiplicativa. De facto

‖AB‖= sup‖x‖=1‖ABx‖v

≤ sup‖x‖v=1

‖A‖‖Bx‖v

= ‖A‖ sup‖x‖v=1

‖Bx‖v

= ‖A‖‖B‖.


Algumas das normas matriciais introduzidas anteriormente são normas induzidas por normas vetoriais.De facto, a norma matricial ‖·‖∞ é induzida pela correspondente norma vectorial. Para provar esteresultado notamos que se tem sucessivamente

‖A‖∞ = sup‖x‖∞=1

maxi

∣∣∣∣n

∑j=1

ai jx j

∣∣∣∣

≤ sup‖x‖∞=1

maxi

n

∑j=1|ai j|‖x‖∞

= maxi

n

∑j=1|ai j|.

Por outro lado,

‖A‖∞ ≥maxi

∣∣∣∣n

∑j=1

ai jx j

∣∣∣∣, x ∈ Rn, ‖x‖∞ = 1.

De facto, fixemos, para i = 1,2, . . . ,n, o vetor x definido por x j =ai j

|ai j|se ai j 6= 0 e x j = 0 se ai j = 0.

Então

‖A‖∞ ≥n

∑j=1|ai j|,

para i = 1, . . . ,n, e portanto

‖A‖∞ ≥maxi

n

∑j=1|ai j|,

o que conclui a prova da seguinte igualdade

‖A‖∞ = maxi

n

∑j=1|ai j|.

No que diz respeito à norma matricial ‖·‖1, começamos por notar que se tem

‖A‖1 = sup‖x‖1=1

n

∑i=1

∣∣∣∣n

∑j=1

ai jx j

∣∣∣∣

≤ sup‖x‖1=1

n

∑j=1

n

∑i=1|ai j||x j|

= maxj

n

∑i=1|ai j|.

Por outro lado, notemos que

‖A‖1 ≥n

∑i=1

∣∣∣∣n

∑j=1

ai jx j

∣∣∣∣, ‖x‖1 = 1.

Fixemos, para `= 1, . . . ,n, o vetor x tal que x` = 1 e x j = 0, j 6= `. Então

n

∑i=1

∣∣∣∣n

∑j=1

ai jx j

∣∣∣∣=n

∑i=1|ai`|,

Logo‖A‖1 ≥max

j∑i=1|ai j|,


o que conclui a prova da seguinte igualdade

‖A‖1 = maxj

∑i=1|ai j|.

Com o objetivo de estabelecer uma caracterização para a norma ‖A‖2 recordamos o conceito de raioespectral de uma matriz A ∈Mn(R): ρ(A) = maxi |λi|, em que λi denota um valor próprio de A.

Observamos que vale o seguinte resultado relativamente a matrizes reais simétricas (At = A).

Teorema 1.2. Se A ∈Mn(R) é real simétrica, então os valores próprios de A são reais e existeuma matriz ortogonal Q (Q−1 = Qt) tal que

QtAQ = D,

em que D é a matriz diagonal dos valores próprios de A.

Recordamos ainda que se A é simétrica e positiva semi-definida, isto é,

xtAx≥ 0, x ∈ Rn,

então os valores próprios de A são não negativos.

O teorema seguinte estabelece a caracterização de ‖A‖2.

Teorema 1.3. Se A ∈Mn(R), então

‖A‖2 =√

ρ(AtA) =√

ρ(AAt).

Se A é simétrica, então‖A‖2 = ρ(A).

Demonstração. Consideremos a matriz AtA. Esta matriz é real simétrica e é positiva semi--definida pois

xt(AtA)x = (Ax)t(Ax) = ‖Ax‖22 ≥ 0, ∀x ∈ Rn.

Logo os seus valores próprios são reais não negativos e existe uma matriz ortogonal Q tal que

QtAtAQ = D,

em que D é a matriz diagonal dos valores próprios µi de AtA.


Tem-se sucessivamente

‖A‖2 = sup‖x‖=1‖Ax‖2

= sup‖x‖=1

√xtAtAx

= sup‖x‖=1

√xtQDQtx

= sup‖x‖=1

√(Qtx)tD(Qtx)

≤ sup‖x‖=1

√max

iµi‖Qtx‖2

2

=√

maxi

µi.

Por outro lado,‖A‖2 ≥ ‖Ax‖2, ‖x‖2 = 1.

Notamos que sendo x vetor próprio de AtA associado ao valor próprio µ , temos

‖Ax‖2 =√

xtAtAx =√

µxtx =√

µ.

Esta desigualdade permite concluir‖A‖2 ≥

√ρ(AtA)

o que finaliza a prova da igualdade‖A‖2 =

√ρ(AtA).

Suponhamos agora que A é real simétrica. Então AtA = A2 e portanto

ρ(AtA) = ρ(A2) = (ρ(A))2.

A finalizar esta breve digressão pelas normas matriciais, estabelecemos seguidamente a existênciade inversa da matriz I−A quando ‖A‖< 1 e a norma ‖·‖ é submultiplicativa. Para tal, consideremos noespaço Mn(R) sucessões

A(m) =[a(m)

i j

]n

i, j=1, m ∈ N.

Atendendo a que Mn(R) é um espaço normado, introduzimos no que se segue o conceito de limite deA(m), m ∈ N.

Definição 1.7. Seja A(m) =[a(m)

i j

]n

i, j=1, m ∈ N, uma sucessão em Mn(R) e A ∈Mn(R). Se

limm

∥∥A(m)−A∥∥= 0,

então dizemos que A(m) converge para A e escrevemos limm

A(m) = A.

Exemplo 1.3. A sucessão

A(m) =

[1

1+m2 msin(

1m

)

1 m2

1+m2

], m ∈ N,


converge, relativamente à norma ‖·‖∞, para a matriz[

0 11 1

].

Consideremos agora a sucessãom

∑j=0

A j, m ∈ N e suponhamos que ‖A‖< 1. Temos sucessivamente

(I−A)m

∑j=0

A j =m

∑j=0

A j−m+1

∑j=1

A j

= I−Am+1,

em que Am+1→ 0. Logo

limm

((I−A)

m

∑j=0

A j)= I.

De igual modo se prova que

limm

m

∑j=0

A j(I−A) = I.

Das duas últimas igualdades deduzimos que a inversa de I−A existe e

(I−A)−1 = limm

m

∑j=0

A j =∞

∑j=0

A j. (1.3.5)

Vejamos agora como determinar um enquadramento para∥∥(I−A)−1

∥∥. De (1.3.5) vem

∥∥(I−A)−1∥∥≤∞

∑j=0‖A‖ j

=1

1−‖A‖ .

Se ‖I‖= 1, então tem-se

1 = ‖I‖=∥∥(I−A)−1(I−A)

∥∥≤∥∥(I−A)−1∥∥(1+‖A‖),

deduzimos1

1+‖A‖ ≤∥∥(I−A)−1∥∥,

o que finaliza a prova do teorema seguinte:

Teorema 1.4. Seja ‖·‖ uma norma submultiplicativa tal que ‖I‖ = 1. Se ‖A‖ < 1 então existe(I−A)−1,

(I−A)−1 =∞

∑j=0

A j,

e1

1+‖A‖ ≤∥∥(I−A)−1∥∥≤ 1

1−‖A‖ .


1.3.2 Número de condição de uma matriz

Seja A não singular e consideremos o sistema Ax = b. Seja x a solução do sistema

Ax = b

em que b é o vetor das aproximações dos termos independentes. Seja ∆b = b− b o erro absoluto e∆x = x− x o vetor dos erros absolutos das componentes xi, i = 1, . . . ,n.

Exemplo 1.4. A título ilustrativo consideremos o seguinte exemplo:

A =

[0.55 0.423

0.484 0.372

], b =

[0.1270.112

].

A solução do sistema Ax = b é x = (1,−1). Consideremos agora

b =

[0.127070.11228

].

Então x = (1.7,−1.91). Observamos que uma pequena perturbação em b induz uma perturbação signifi-cativa em x tendo-se

∆b =

[−0.00007−0.00028

], ‖∆b‖∞ = 0.00028,

∆x =[−0.70.91

], ‖∆x‖∞ = 0.91,

e ainda‖∆b‖∞

‖b‖∞

= 0.00220472,‖∆x‖∞

‖x‖∞

= 0.91.

Este exemplo mostra que uma pequena perturbação no vetor dos termos independentes pode induzirperturbações significativas na solução do sistema.

De que modo podemos, a priori, garantir que pequenas perturbações nos dados induzem pequenasperturbações na solução do sistema?

De Ax = b, vemA(x−∆x) = b−∆b.

Uma vez que Ax = b, vem A∆x = ∆b e portanto ∆x = A−1∆b. Por outro lado, na hipótese das normasmatricial e vectorial serem compatíveis, temos

‖b‖ ≤ ‖A‖‖x‖.

Logo

‖∆x‖‖x‖ ≤

∥∥A−1∥∥‖∆b‖‖x‖

≤∥∥A−1∥∥‖A‖‖∆b‖

‖b‖ .

A partir da desigualdade anterior concluímos que a quantidade

cond(A) = ‖A−1‖‖A‖


tem um papel de relevo no comportamento de ‖∆x‖‖x‖ relativamente a ‖∆b‖

‖b‖ . De facto, se cond(A) é pequeno,

então ‖∆x‖‖x‖ é pequeno desde que ‖∆b‖

‖b‖ seja pequeno. Mais ainda

∀ε > 0 ∃η > 0 : ∀∆b tal que‖∆b‖‖b‖ < η =⇒ ‖∆x‖

‖x‖ < ε,

pois é suficiente tomar η < εcond(A) .

Observamos que se cond(A) é pequeno, então mesmo para perturbações do vetor dos termos inde-pendentes com significado, a perturbação na solução pode tornar-se arbitrariamente pequena. Neste casodizemos que o sistema é bem condicionado. Se pequenas perturbações dos termos independentes induzemgrandes perturbações no vetor solução, então o sistema Ax = b diz-se mal condiconado. A cond(A)chamamos número de condição da matriz A.

Observamos que a matriz A referida no exemplo introdutório tem a seguinte inversa

A−1 =

[−2818.18 3204.553.666.67 −4166.67

]

e portanto cond1(A) = 7621.83.O número de condição de uma matriz é, em certos casos, facilmente determinado a partir dos valores

próprios da matriz A. Comecemos por notar que se A é não singular, então λ = 0 não é valor próprio deA. Mais ainda, os inversos dos valores próprios de A são os valores próprios de A−1. Se A é simétrica,pelo Teorema 1.3,

‖A‖2 = maxi|λi| e

∥∥A−1∥∥2 =

1mini |λi|

,

o que permite estabelecer o seguinte resultado:

Teorema 1.5. Se A é real simétrica, então

cond2(A) =maxi |λi|mini |λi|

.

Consideremos agora x uma aproximação para a solução de Ax = b. O vetor dos resíduos

r = b−Ax

e o número de condição de A, cond(A), é usado seguidamente para estimar o erro‖∆x‖‖x‖ . De facto, tem-se

∆x = x− x

= A−1b−A−1Ax

= A−1(b−Ax)

= A−1r.

Logo‖∆x‖‖x‖ ≤ cond(A)

‖r‖‖b‖ .

Estudámos o efeito dos erro nos dados e dos resíduos no erro da solução do sistema. É de salientarque a matriz do sistema pode estar também afectada de erros, isto é, é conhecida a matriz A com um erro


∆A = A−A. Neste caso, o erro é uma matriz cujas entradas são os erros absolutos das entradas de A. Se xé a solução do sistema Ax = b e x é tal que Ax = b, então

(A−∆A)(x−∆x) = b−∆b

e portanto(A−∆A)∆x = ∆b−∆Ax,

ou ainda, supondo que A é não singular,

(I−A−1∆A)∆x = A−1(∆b−∆Ax). (1.3.6)

Suponhamos agora que∥∥A−1∆A

∥∥< 1 e que a norma ‖·‖ é submultiplicativa. Pelo Teorema 1.4, existe(I−A−1∆A)−1 e, de (1.3.6), deduzimos

∆x = (I−A−1∆A)−1A−1(∆b−∆Ax),

ou ainda‖∆x‖ ≤

∥∥(I−A−1∆A)−1∥∥∥∥A−1∥∥(‖∆b‖+‖∆A‖‖x‖). (1.3.7)

Utilizando a estimativa para ‖(I−A−1∆A)−1‖ estabelecida no Teorema 1.4 concluímos

‖∆x‖‖x‖ ≤

∥∥A−1∥∥

1−∥∥A−1

∥∥‖∆A‖

(‖∆b‖‖x‖ +‖∆A‖

). (1.3.8)

Usando em (1.3.8) a desiguldade ‖x‖ ≥ ‖A‖−1‖b‖ obtemos finalmente

‖∆x‖‖x‖ ≤

cond(A)1− cond(A)‖δA‖(‖δb‖+‖δA‖) (1.3.9)

em que

‖δA‖= ‖∆A‖‖A‖ , ‖δb‖= ‖∆b‖

‖b‖ .

Provámos o resultado seguinte:

Teorema 1.6. Seja A uma matriz não singular e ‖·‖ uma norma submultiplicativa. Se ∆A é talque ∥∥A−1∥∥‖∆A‖< 1,

então o erro relativo ‖δx‖ satisfaz (1.3.9).

1.4 Alguns problemas

1. Determine o número de termos que deve considerar em ∑∞j=1(−1) j+1 x j

j para determinar umaaproximação para f (x) = ln(1+ x), x ∈ [0,1], com erro inferior a 10−3.

2. Calcule ln(2) com erro inferior a 10−2.

3. Estabeleça a estimativa (1.1.6).


4. Indique uma estimativa para os erros absoluto e relativo das seguintes quantidades

x = 8.923, y = 0.79, z = 1.6542,

sabendo que foram obtidas por arredondamento.

5. Considere w = ln(

x4y3

z2

). Sejam x = 1.24,y = 3.174,z = 4.1 obtidos por conveniente arredonda-

mento e truncatura. Determine uma estimativa para o erro que se comete ao substituir w porw = ln

(x4y3

z2

).

6. Considere que x = 0.037 tem 3 casas decimais correctas.

(a) Determine estimativas para os erros absoluto e relativo.(b) Seja f (x) =

√1− x. Determine estimativas para os erros absoluto e relativo.

7. Determine uma aproximação polinomial para f (x) = sinh(x) em [0,1], com erro inferior a 10−2.

8. Mostre que as aplicações

‖x‖∞ = maxi=1,...,n

|xi|, ‖x‖1 =n

∑i=1|xi|, ‖x‖2 =

( n

∑i=1

x2i

)1/2

,

são normas em Rn.

9. Mostre que as normas indicadas na alínea anterior são equivalentes.

10. Mostre que as aplicações (1.3.1), (1.3.2), (1.3.3) e (1.3.4) são normas em Mn(R).

11. Mostre que a matriz

0 1 11 0 −1−1 −1 0

é semelhante a uma matriz diagonal D = (1,0,−1) e determine a matriz dos vetores próprios. Estaúltima matriz é ortogonal?

12. Mostre que se A é não singular, então λ = 0 não é valor próprio de A.

13. Mostre que se A é não singular, então os valores próprios de A−1 são os inversos dos valores própriosde A.

14. Considere a matriz

A =

0.1 −0.2 00 0.5 −0.1

0.8 0.02 0.1

.

Sem determinar explicitamente I−A, mostre que (I−A)−1 existe e determine um majorante para‖(I−A)−1‖∞.

15. Seja R = I−CA em que ‖R‖< 1 e I,C,A ∈Mn(R). Mostre que A e C são não singulares.

16. Mostre que se A é não singular, então

cond2(A)≥maxi |λi|mini |λi|

.

17. Considere o sistema Ax = b e seja x tal que Ax = b. Mostre que

1cond(A)

‖∆b‖‖b‖ ≤

‖∆x‖‖x‖ ≤ cond(A)

‖∆b‖‖b‖ .

CAPÍTULO 2MÉTODOS NUMÉRICOS PARA SISTEMAS DE EQUAÇÕES LINEARES

2.1 Introdução

Sistemas de equações lineares surgem em várias áreas das ciências e da engenharia. A utilização demétodos eficientes para a determinação da sua solução é um requisito essencial aquando da resolução deproblemas aplicados. O objectivo deste capítulo é o estudo de métodos numéricos para a resolução destetipo de sistemas. Os métodos que serão objecto de estudo dividem-se em duas grandes classes: métodosdirectos e métodos iterativos.

Designamos por métodos directos os métodos que permitem obter a solução do sistema num númerofinito de operações. Os métodos deste grupo baseiam-se fundamentalmente no processo de eliminação deGauss. Os métodos numéricos que geram uma sucessão de aproximações para a solução do sistema linearpertencem à classe dos métodos iterativos. Esta classe divide-se em dois grandes grupos: os métodosestacionários e os métodos não estacionários.

A escolha de um método directo ou iterativo está intrinsecamente ligada à estrutura da matriz dosistema. Assim os primeiros métodos são privilegiados para sistemas com matrizes densas, isto é, matrizescom poucas entradas nulas. Por outro lado, para sistemas em que a matriz associada é de grandesdimensões e com poucas entradas não nulas é conveniente a utilização de métodos iterativos.

2.2 Métodos Directos

2.2.1 Método de eliminação de Gauss

Consideremos o sistema Ax= b em que A∈Mn(R), b∈Rn e A é não singular. Esta última propriedadeé equivalente às seguintes:

• existe A−1,

• A tem característica n,

17


• det(A) 6= 0,

• as n colunas de A são linearmente independentes,

• as n linhas de A são linearmente independentes,

• o sistema homogéneo Ax = 0 tem apenas a solução nula.

O primeiro método que vamos considerar é o método de eliminação de Gauss. Este método é compostode duas fases: a fase descendente e a fase ascendente. A fase descendente é baseada em operações pivotaisque consistem em substituir cada linha pela soma, elemento a elemento, dela própria com um múltiploescalar de outra. Por aplicações sucessivas de operações deste tipo, na primeira fase o sistema inicial ésubstituído por um número finito de sistemas equivalentes obtendo-se, no fim da fase descendente, umsistema triangular superior. Na segunda fase (ascendente), a solução do sistema triangular obtido no fimda primeira fase é obtida por substituição de baixo para cima.

Fase descendente:

Consideremos o sistema

a11x1 + a12x2 + . . . + a1nxn = b1a21x1 + a22x2 + . . . + a2nxn = b2a31x1 + a32x2 + . . . + a3nxn = b3

......

. . ....

...an1x1 + an2x2 + . . . + annxn = bn

.

Como a matriz é não singular, por uma troca de linhas conveniente, podemos assumir que a11 6= 0.Substituam-se a linhas Li, i = 2, . . . ,n, por

Li−ai1

a11L1, i = 2, . . . ,n.

Obtemos o sistema equivalente

a11x1 + a12x2 + . . . + a1nxn = b1

a(2)22 x2 + . . . + a(2)2n xn = b(2)2

a(2)32 x2 + . . . + a(2)3n xn = b(2)3...

. . ....

...a(2)n2 x2 + . . . + a(2)nn xn = b(2)n

,

em que

a(2)i j = ai j−ai1

a11a1 j, i, j = 2, . . . ,n, b(2)i = bi−

ai1

a11b1.

Se a(2)22 6= 0 (notemos que podemos reduzir-nos sempre a esta situação mediante uma troca apropriada delinhas), consideremos sobre o último sistema, a seguinte operação

Li = Li−a(2)i2

a(2)22

L2, i = 3, . . . ,n.

Obtemos o sistema linear equivalente


a11x1 + a12x2 + a13x3 + . . . + a1nxn = b1

a(2)22 x2 + a(2)23 x3 + . . . + a(2)2n xn = b(2)2

a(3)33 x3 + . . . + a(3)3n xn = b(3)3...

. . ....

...a(3)n3 x3 + . . . + a(3)nn xn = b(3)n

,

em que

a(3)i j = a(2)i j −a(2)i2

a(2)22

a(2)2 j , i, j = 3, . . . ,n, b(3)i = b(2)i −a(2)i2a22

b(2)2 .

Repetindo o processo anterior obtemos finalmente o seguinte sistema triangular

b11x1 + b12x2 + . . . + b1nxn = c1b22x2 + . . . + b2nxn = c2

. . ....

...bnnxn = cn

. (2.2.1)

Uma forma simples de implementar a fase descendente da eliminação de Gauss (assumindoque não é necessária a troca de linhas no processo) está descrita no seguinte algoritmo:

Algoritmo 2.1 Fase descendente da eliminação de Gauss

function [A,b] = faseDescendente(A,b)[n,m]=size(A);if (n 6= m), error('Apenas para matrizes quadradas'); endfor k=1:n−1

if (A(k,k) ==0), error('Pivot nulo'); endfor i=(k+1):n

multiplicador = A(i,k)/A(k,k);A(i,1:n) = A(i,1:n) − multiplicador*A(k,1:n);b(i) = b(i) − multiplicador*b(k);

endendreturn

Notamos que tanto a matriz triangular superior como o vector resultantes da fase descen-dente são armazenados nas variáveis A e b. Um exercício a cargo do leitor é o de reformularo algoritmo anterior de forma a que apenas sejam alteradas, durante a fase descendente, asentradas da diagonal e acima da diagonal da matriz A, reduzindo assim o número de operaçõesaritméticas usadas.

Fase ascendente

A solução do sistema (2.2.1) é obtida por substituição de ”baixo para cima”:

xn =cn

bnn, xi =

1bii

(ci−

n

∑`=i+1

bi`x`

), i = n−1,n−2, . . . ,1.


A fase descendente da eliminação de Gauss pode ser implementada através do seguintealgoritmo:

Algoritmo 2.2 Substituição ascendente

function [x] = substituicaoAscendente(A,b)[n,m]=size(A);if (n 6= m), error('Apenas para matrizes quadradas'); endx = zeros(n,1);x(end) = b(n)/A(n,n);for i=n−1:−1:1x(i) = (b(i) − sum(A(i,i+1:n)'.*x(i+1:n)))/A(i,i);

endendreturn

As operações realizadas sobre o sistema Ax = b podem ser organizadas num quadro chamado matrizampliada

a11 a12 a13 . . . a1n b1a21 a22 a23 . . . a2n b2a31 a32 a33 . . . a3n b3

......

.... . .

......

an1 an2 an3 . . . ann bn

.

De facto, considerando sobre o quadro anterior as seguintes operações

a(2)i j = ai j−ai1

a11a1 j, i, j = 2, . . . ,n, b(2)i = bi−

ai1

a11bi, i = 2, . . . ,n,

obtemos

a11 a12 a13 . . . a1n b1

0 a(2)22 a(2)23 . . . a(2)2n b(2)2

0 a(2)32 a(2)33 . . . a(2)3n b(2)3...

......

. . ....

...0 a(2)n2 a(2)n3 . . . a(2)nn b(2)n

.

Considerando agora sobre o último quadro as operações

a(3)i j = a(2)i j −a(2)i2

a(2)22

a(2)2 j , i, j = 3, . . . ,n, b(3)i = b(2)i −a(2)i2

a(2)22

b(2)2 , i = 3, . . . ,n,

deduzimos

a11 a12 a13 . . . a1n b1

0 a(2)22 a(2)23 . . . a(2)2n b(2)2

0 0 a(3)33 . . . a(3)3n b(3)3...

......

. . ....

...0 0 a(3)n3 . . . a(3)nn b(3)n

.


O procedimento anterior conduz-nos ao seguinte quadro

a11 a12 a13 . . . a1n b1

0 a(2)22 a(2)23 . . . a(2)2n b(2)2

0 0 a(32)33 . . . a(3)3n b(3)3

......

.... . .

......

0 0 0 . . . a(n)nn b(n)n

que está associado ao sistema triangular

a11x1 + a12x2 + a13x3 + . . . + a1nxn = b1

a(2)22 x2 + a(2)23 x3 + . . . + a(2)2n xn = b(2)2

a(3)33 x3 + . . . + a(3)3n xn = b(3)3. . .

......

a(n)nn xn = b(n)n

A solução, tal como anteriormente, obtém-se por substituição, de baixo para cima. É de salientar quesendo A não singular, o sistema é possível e determinado tendo solução única.

O método de Gauss-Jordan designa outro método directo em que ao fim da primeira fase do métodode eliminação de Gauss se considera uma nova fase de eliminação mas agora de baixo para cima.Consideremos sobre o último sistema as seguintes operações pivotais

Li = Li−a(i)in

a(n)nn

Ln, i = n−1, . . . ,1.

onde assumimos que a(1)1 j = a1 j. Obtemos

a11x1 + a12x2 + a13x3 + . . . + a1,nxn = b1

a(2)22 x2 + a(2)23 x3 + . . . + a2,nxn = b(2)2

a(3)33 x3 + . . . + a(3)3,n−1xn = b(3)3. . . . . .

......

a(n−1)n−1,n−1xn−1 + a(n−1)

n−1,nxn = b(n−1)n−1

a(n)nn xn = b(n)n

O processo repete-se e conduz-nos a um sistema diagonal

d1x1 = c1d2x2 = c2

d3x3 = c3. . .

...dnxn = cn

cuja solução é x =( c1

d1, . . . , cn

dn

).

Observamos que o processo de eliminação (ascendente) acabado de descrever pode ser implementadonas mesmas linhas do Algoritmo 2.1.


Se A é singular com característica r < n, a primeira fase do método de eliminação de Gauss conduz-nosa um sistema triangular inferior do seguinte tipo

a11x1 + a12x2 + a13x3 + . . . + . . . + a1nxn = b1

a(2)22 x2 + a(2)23 x3 + . . . + . . . + a(2)2n xn = b(2)2

a(3)33 x3 + . . . + . . . + a(3)3n xn = b(3)3. . .

......

a(r)rr xr + . . . + a(r)rn xn = b(r)r

obtendo-se a solução deste sistema por substituição.Na descrição do processo de eliminação de Gauss, assumimos que o candidato a pivot é não nulo.

Se num determinado passo da primeira fase, o candidato a pivot, app, é nulo, então devemos procurar aprimeira linha Li, i = p+1, . . . ,n em que aip 6= 0. Seja L` a linha escolhida. Então devemos trocar a linhaLp com a linha L`, isto é,

Lp = L`, L` = Lp,

e continuar com o processo de eliminação de Gauss. Observamos que poderíamos procurar o candidato apivot entre as colunas Ci, i = p+1, . . . ,n.

A definição do pivot, em cada passo do processo, é fundamental para manter controlada a propagaçãodos erros de arredondamento. De facto, se o candidato a pivot é muito pequeno (em valor absoluto) emrelação às outras entradas da matriz, obtemos uma matriz resultante com entradas muito grandes. A fimde evitar tal situação, devemos escolher em cada passo da eliminação de Gauss, a entrada dominante (emmódulo) e efectuar trocas de linhas (ou colunas) de modo a que tal entrada seja o pivot.

Introduzimos a seguinte routina auxiliar, responsável pela troca de linhas de uma matriz:

Algoritmo 2.3 Troca das linhas i e j

function [A] = trocaLinhas(A,i,j)x=A(i,:); A(i,:) = A(j,:); A(j,:) = x;return

Para combinar o algoritmo anterior com o algoritmo da fase descendente, basta substituira quinta linha do Algoritmo 2.1 pelo Algoritmo 2.4 e obtemos uma implementação da fasedescendente com escolha parcial de pivot.

Algoritmo 2.4 Escolha parcial de pivot

% encontrar indice correspondente ao novo pivot[C,I] = max(abs(A(k:end,k)));indice = k+I(1)−1;

A = trocaLinhas(A,k,indice); b = trocaLinhas(b,k,indice);

O exemplo seguinte ilustra a eficiência da escolha parcial de pivot.

Exemplo 2.1. No sistema x1 + 200x2 = 100x1 + x2 = 1

,


cuja solução é x =(100

199 ,99

199

)=(0.502 . . . ,0.497 . . .

), o candidato a pivot na primeira fase do processo é

a11 = 1. Se consideramos em cada operação apenas números com 3 casas decimais, obtemos x2 = 0.6 ex1 = 0.497 sendo significativo o erro na primeira componente. Por outro lado, a escolha parcial de pivotconduz-nos ao seguinte sistema

200x2 + x1 = 100x2 + x1 = 1

,

e, efectuando o processo de eliminação de Gauss considerando em cada operação 3 dígitos, obtemos

200x2 + x1 = 1000.995x1 = 0.5

,

e finalmente x = (0.503,0.497). De facto, a escolha do elemento máximo (em módulo) em de cada linha(coluna) reduz a propagação dos erros de arredondamento.

Observamos, por último, que o MATLAB possui um comando para resolver sistemas linearesutilizando métodos directos. De facto, para obter a solução x de Ax = b, basta escrever x=A\b.Esta operação analisa a estrutura da matriz A (simétrica, diagonal, tridiagonal, etc...) e escolhede uma lista de algoritmos implementados, o que acha mais adequado para determinar a soluçãodo sistema (mais pormenores podem ser encontrados em [6]). Realçamos também, que emborado ponto de vista algébrico esta operação é equivalente a inv(A)*b, o cálculo (numérico) dainversa de uma matriz é uma operação computacionalmente muito mais exigente que a soluçãodo sistema linear Ax = b, e portanto, deve ser uma prática a evitar.

2.2.2 Métodos de factorização

A aplicação do método de eliminação de Gauss, ou a sua variante Gauss-Jordan, ao sistema linearAx = b altera a matriz do sistema bem como o vector dos termos independentes. Assim, a resolução deum novo sistema em que há apenas a alteração do vector dos termos independentes requer a aplicaçãocompleta das duas fases do método.

Admitamos que é conhecida uma factorização

A = LR

em que L e R são matrizes triangular inferior e superior, respectivamente. Assim o sistema inicial éequivalente a

LRx = b.

Seja y = Rx. A solução do sistema obtém-se considerando os dois sistemas triangulares

Ly = b

Rx = y,

ou seja

y1 =b1

`11, yi =

1ìi

(bi−

i−1

∑j=1

ì jy j

), i = 2, . . . ,n,

e

xn =yn

rnn, xi =

1rii

(yi−

n

∑j=i+1

ri jx j

), i = n−1, . . . ,1.


Vejamos como determinar matrizes L e R nestas condições. Um dos processos mais imediatosconsiste em considerar o processo de eliminação de Gauss em termos matriciais. Observamos que o passoelementar

Li = Li−ai`

a``L`

sobre a matriz A traduz-se pelo seguinte produto matricial

EiÀ

em que Ei` denota a seguinte matriz elementar

Ei` =

1 0 0 . . . . . . . . . 00 1 0 . . . . . . . . . 00 . . . . . . . . . . . . . . . 00 0 − ai`

a``. . . 1 . . . 0

0 . . . . . . . . . . . . . . . 00 . . . . . . . . . . . . 1 00 . . . . . . . . . . . . . . . 1

.

Assim a eliminação dos elementos abaixo da posição (1,1) no primeiro passo elementar que consiste emconsiderar sobre A as operações pivotais

Li = Li−ai1

a11L1, i = 2, . . . ,n,

é equivalente ao produto

1 0 0 . . . . . . 0−a21

a111 0 . . . . . . 0

−a31

a110 1 . . . . . . 0

. . . . . . . . . . . . . . . 0−an1

a11. . . . . . . . . . . . 1

A.

Denotemos por E1 a primeira matriz do produto anterior. O segundo passo elementar é feito sobre amatriz resultante do produto anterior, digamos sobre

[a(2)i j

]em que a(2)i1 = 0, i = 2, . . . ,n. Este passo é

equivalente a multiplicar a matriz[a(2)i j

]por

E2 =

1 0 0 . . . . . . 00 1 0 . . . . . . 0

0 −a(2)32

a(2)22

1 0 . . . 0

......

.... . . 0

0 −a(2)n2

a(2)22

0 . . . . . . 1

.

Assim, se[a(3)i j

]denota a matriz resultante do segundo passo elementar, a(3)i1 = 0, i = 2, . . . ,n, e

a(3)i2 = 0, i = 3, . . . ,n, esta matriz satisfaz

E2E1A =[a(3)i j

].


Repetindo este processo, no fim da fase descendente de eliminação de Gauss realizada sobre a matriz A,obtemos finalmente

En−1 · · ·E2E1A = R

em que R denota a matriz triangular superior obtida no fim da eliminação de Gauss. Seja T o produtoEn−1 . . .E2E1. Esta matriz é triangular inferior e não singular pois Ei é não singular sendo a sua inversa amatriz que se obtém de Ei substituindo as suas entradas não diagonais pelos seus simétricos. Assim

L = T−1 = E−11 E−1

2 · · ·E−1n−1

e portantoA = LR.

Observamos que a matriz L é relativamente fácil de construir: é uma matriz triangular inferior comelementos diagonais iguais a 1 e cujos elementos sob a diagonal são os simétricos dos coeficientes usadosnas operações pivotais da eliminação, cada um deles na respectiva posição que figura na correspondentematriz elementar Ei.

O cálculo das matrizes L e R pode ser implementado no seguinte algoritmo (adaptação doAlgoritmo 2.1):

Algoritmo 2.5 Factorização LR

function [L,A] = factorizacaoLR(A)[n,m]=size(A);L = eye(n,m);if (n 6= m), error('Apenas para matrizes quadradas'); endfor k=1:n−1

if (A(k,k) ==0), error('Pivot nulo. E necessaria troca de linhas'); endfor i=(k+1):n

L(i,k) = A(i,k)/A(k,k);A(i,1:n) = A(i,1:n) − L(i,k)*A(k,1:n);

endendreturn

Suponhamos que num determinado passo da fase descendente da eliminação de Gauss, o candidatoa pivot é nulo, digamos a(p)

pp = 0, isto é, o passo elementar p− 1 conduz-nos, em termos matriciais, àigualdade

Ep−1Ep−2 · · ·E1A = A(p)

para uma determinada matriz A(p). Tal como referido anteriormente, procuramos na coluna p e entre aslinhas p+1, . . . ,n, a primeira linha que apresenta uma componente não nula. Seja a(p)

`p 6= 0 a entrada coma propriedade anterior. O passo seguinte consiste em trocar as linhas p e `. Por P p denotamos a matriz depermutação que se obtém da matriz identidade de ordem n trocando as linhas p e `. O produto de P p poruma outra matriz obtém-se desta última matriz trocando as linhas p e `. Então

P pEp−1Ep−2 · · ·E1A = P pA(p)

em que a componente (p, p) de P pA(p) é não nula. Sobre esta última matriz efectuamos mais um passoelementar de eliminação de Gauss. Repetido o processo somos conduzidos a

En−1 · · ·Pkm · · ·EpP pEp−1Ep−2 · · ·E1A = R


em que R é uma matriz triangular superior.Uma vez que as matrizes de permutação e as matrizes elementares são invertíveis, da última igualdade

obtemos

A = E−11 · · ·E−1

p−2P p · · ·Pkm · · ·E−1n−1R.

Seja P o produto das matrizes de permutação que surgem na fase de eliminação. A fim de obter afactorização com as propriedades referidas devemos multiplicar a última igualdade por P obtendo-se

PA = PE−11 · · ·E−1

p−2P p · · ·Pkm · · ·E−1n−1R,

isto é,

PA = LR com L = PE−11 · · ·E−1

p−2P p · · ·Pkm · · ·E−1n−1.

Provámos o seguinte resultado:

Teorema 2.1. Se A é uma matriz quadrada de ordem n não singular, então existem matrizes deordem n P, L e R, tais que P é uma matriz de permutação, L é uma matriz triangular inferior comelementos diagonais iguais a 1 e R é uma matriz triangular superior, que verificam PA = LR. SeP = I, então a factorização é única.

Demonstração. Provemos a unicidade. Suponhamos que

A = L1R1, A = L2R2.

Então

L−12 L1 = R2R−1

1 .

Uma vez que L1 e L−12 são matrizes triangulares inferiores com elementos diagonais iguais à unidade, vem

L−12 L1 = I, R2R−1

1 = I.

Logo

L2 = L1,R2 = R1.

Uma implementação possível do cálculo da factorização apresentada, pode ser baseada noprocesso da eliminação de Gauss descrito e na escolha parcial de pivot:


Algoritmo 2.6 Factorização PA=LR com escolha parcial de pivot

function [P,L,A] = factorizacaoPALR(A)[n,m]=size(A);L = eye(n,m);P = eye(n,m);if (n 6= m), error('Apenas para matrizes quadradas'); endfor k=1:n−1

[C,I] = max(abs(A(k:end,k)));indice = k+I(1)−1;A = trocaLinhas(A,k,indice); P = trocaLinhas(P,k,indice);L(:,1:k−1) = trocaLinhas(L(:,1:k−1),k,indice);

for i=(k+1):nL(i,k) = A(i,k)/A(k,k);A(i,1:n) = A(i,1:n) − L(i,k)*A(k,1:n);

endendreturn

A partir da factorização A = LR podemos facilmente obter

A = LDR, (2.2.2)

em que L foi definida anteriormente, D é a matriz diagonal cujas entradas são os pivots e R é a matriztriangular superior com elementos diagonais igual à unidade

R =

1 r12r11

r13r11

. . . r1nr11

0 1 r23r22

. . . r2nr22

. . . . . . . . . . . . . . .0 0 0 . . . 1

.

A construção da solução do sistema Ax = b utilizando a factorização (2.2.2) é feita seguindo os passosseguintes:

1. Determinar y tal que Ly = b;

2. Determinar z tal que Dz = y;

3. Determinar x tal que Rx = z.

No caso particular de A ser simétrica e positiva definida e não efectuarmos troca de linhas na fase deeliminação de Gauss, A admite a seguinte factorização

A = LLt ,

em que L = LD e D é uma matriz diagonal obtida de D calculando a raiz quadrado dos elementos diagonaisde D. De facto, como A admite a factorização A = LR, At admite a factorização At = RtLt . Da hipóteseAt = A, concluimos que

A = LR = RtLt .


Reescrevendo a igualdade anterior usando (2.2.2), obtemos

L(DR) = Rt(DLt).

Atendendo à unicidade da factorização, deduzimos

L = Rt

o que prova o pretendido.Utilizando a factorização A = LLt a construção da solução do sistema Ax = b é feita seguindo os

passos seguintes

1. Determinar a solução de Ly = b;

2. Determinar a solução de Ltx = y.

Vamos seguidamente utilizar na resolução do sistema Ax = b outros tipos de factorizações: a facto-rização QR e a factorização em valores singulares. Estas factorizações são introduzidas de modo a terem cada passo da resolução do sistema Ax = b, sistemas auxiliares bem condicionados (com matrizesassociadas bem condicionadas). De facto, os métodos de factorização anteriores, baseados na eliminaçãode Gauss, embora computacionalmente simples, não nos permitem garantir que tenhamos sistemas bemcondicionados pois não controlamos, no caso quadrado, os números de condição das matrizes envolvidas.

A decomposição QR, em que Q é uma matriz ortogonal, é obtida utilizando a conhecida classe dematrizes ditas de Householder. Estas matrizes têm a seguinte estrutura particular

H = I−2wwt ,

em que w é um vector unitário (relativamente à norma ‖·‖2). Observamos que as matrizes de Householdersão simétricas pois Ht = H e são ortogonais. De facto

HH = (I−2wwt)(I−2wwt) = I−4wwt +4wwtwwt .

Atendendo a que w é unitário, isto é ‖w‖2 = 1, vem

HH = I.

Vejamos como, utilizando produtos à esquerda da matriz A, construímos uma matriz triangular superior.Seja a1 a primeira coluna da matriz A. Então a primeira coluna b1 do produto de HA é dada por HAe1, ouseja

b1 = Ha1.

Pretende-se que esta coluna tenha apenas a entrada (1,1) não nula. Façamos a escolha de w de modo aque o resultado Ha1 satisfaça este requisito, isto é, Ha1 = βe1.

Observamos queHa1 = a1−2wwta1

e tomemosw =

v‖v‖2

, v = a1−βe1, β = sgn(a11)‖a1‖2.

Notamos que(a1−βe1)

ta1 = ‖a1‖22−βa11,

e‖v‖2

2 = (a1−βe1)t(a1−βe1) = ‖a1‖2

2−2βa11 +β 2.


Atendendo à escolha de β deduzimos

(a1−βe1)ta1 = ‖a1‖2

(‖a1‖2− sgn(a11)a11

),

e‖v‖2

2 = 2‖a1‖2(‖a1‖2− sgn(a11)a11

),

em que ‖a1‖2− sgn(a11)a11 > 0 pois caso fosse nulo teríamos a1 com a forma pretendida.Vejamos que Ha1 tem a forma pretendida. De facto

Ha1 = a1−2‖v‖2

2vvta1

= a1−2(a1−βe1)(a1−βe1)

ta1

2‖a1‖2(‖a1‖2− sgn(a11)a11

)

= βe1.

O processo anterior é utilizado seguidamente na determinação de uma factorização QR de uma matrizA do tipo m×n, m≤ n (a generalização para o caso m > n é semelhante). Vimos que se a1 é a primeiracoluna de A então a matriz de Householder de ordem m

Hm = I− 2‖v1‖2

2v1vt

1, v1 = a1−βe1, β = sgn(a11)‖a1‖2,

é tal que

HmA =

a(1)11 a(1)12 a(1)13 . . . a(1)1n

0 a(1)22 a(1)23 . . . a(1)2n

0 a(1)32 a(1)33 . . . a(1)3n

0...

.... . .

...0 a(1)m2 a(1)m3 . . . a(1)mn

Fixemos agora o vector a(1)2 =[a(1)22 a(1)32 · · ·a

(1)m2

]t. A matriz de Householder de ordem m−1

Hm−1 = I− 2‖v2‖2

2v2vt

2, v2 = a(1)2 −βe1, β = sgn(a(1)22

)∥∥a(1)2

∥∥2,

é tal que

Hm−1a(1)2 = βe1,

em que, neste caso, e1 representa o primeiro vector da base canónica de Rm−1. Logo

Hm−1

a(1)22 a(1)23 . . . a(1)2n

a(1)32 a(1)33 . . . a(1)3n...

.... . .

...a(1)m2 a(1)m3 . . . a(1)mn

=

a(2)22 a(2)23 . . . a(2)2n

0 a(2)33 . . . a(2)3n...

.... . .

...0 a(2)m3 . . . a(2)mn

e portanto, denotando S(2) a seguinte matriz por blocos

S(2) =[

1 0t

0 Hm−1

],


temos

S(2)HmA =

a(1)11 a(1)12 a(1)13 . . . a(1)1n

0 a(2)22 a(2)23 . . . a(2)2n

0 0 a(2)33 . . . a(2)3n

0...

.... . .

...0 0 a(2)m3 . . . a(2)mn

.

O processo repete-se. Fixamos a(2)3 =[a(2)33 · · ·a

(2)m3

]t. A matriz de Householder de ordem m−2

Hm−2 = I− 2‖v3‖2

2v3vt

3, v3 = a(2)3 −βe1, β = sgn(a(2)33

)∥∥a(2)3

∥∥2,

é tal queHm−2a(2)3 = βe1,

em que, neste caso, e1 representa o primeiro vector da base canónica de Rm−2. Logo

Hm−2

a(2)33 a(2)34 . . . a(2)3n

a(2)43 a(2)44 . . . a(2)4n...

.... . .

...a(2)m3 a(2)m4 . . . a(2)mn

=

a(3)33 a(3)34 . . . a(3)3n

0 a(3)44 . . . a(3)4n...

.... . .

...0 a(3)m4 . . . a(3)mn

e portanto, denotando S(3) a seguinte matriz por blocos

S(3) =[

I2 0t

0 Hm−2

],

temos

S(3)S(2)HmA =

a(1)11 a(1)12 a(1)13 a(1)14 . . . a(1)1n

0 a(2)22 a(2)23 a(2)24 . . . a(2)2n

0 0 a(3)33 a(3)34 . . . a(3)3n

0 0 0 a(3)44 . . . a(3)4n

0 0 0 a(3)54 . . . a(3)5n

0...

......

. . ....

0 0 0 a(3)m4 . . . a(3)mn

.

No fim da aplicação do processo descrito obtemos

S(m−1) · · ·S(2)S(1)A = R,

em que S(1) = Hm e R é uma matriz triangular superior do tipo m×n.Recordemos que S( j), j = 2, . . . ,m−1,Hm são matrizes ortogonais e tais que S( j)S( j) = I. Logo

A =m−1

∏j=1

S( j)R,

isto éA = QR,

em que Q = ∏m−1j=1 S( j).



Teorema 2.2. Se A é uma matriz real do tipo m×n então existem matrizes Q e R, respectivamente,ortogonal e triangular superior, tais que

A = QR.

O procedimento descrito para calcular a factorização QR é facilmente traduzido paraMATLAB no seguinte algoritmo:

Algoritmo 2.7 Algoritmo para calcular a A = QR usando matrizes de Householder.

function [Q,R] = qrHouseholder(A)[m,n] = size(A);Q = eye(m);for k=1:min(m−1,n)ak = A(k:end,k);betak= sign(ak(1));ek = zeros(size(ak)); ek(1) = 1;vk = ak + sign(ak(1))*norm(ak)*ek;Hk = eye(m−k+1) − 2*vk*vk'/(vk'*vk);Sk = [eye(k−1) zeros(k−1,m−k+1); zeros(m−k+1,k−1) Hk];A = Sk*A;Q = Q*Sk;

endR = A;

Observamos que existem outras forma de calcular a factorização QR introduzida nestasecção, nomeadamente, utilizando o processo de ortogonalização de Gram-Schmidt. No entanto,o Algoritmo 2.7 não demonstra uma performance razoável para matrizes de grandes dimensões.Outras implementações da factorização QR podem ser encontradas em [5]. O leitor podecomparar a performance destas implementações com a routina do MATLAB para calcular estafactorização, [Q,R]=qr(A).

O processo de construção de uma solução do sistema Ax = b recorrendo à factorização QR é análogoao descrito quando foi usada a factorização LR. Assim a solução do sistema anterior é obtida resolvendoos dois sistemas seguintes:

• Qy = b,

• Rx = y.

Relativamente ao primeiro sistema temos cond2(Q) = 1 e consequentemente

‖δy‖2 ≤ ‖δb‖2.

A mesma conclusão não vale para o sistema triangular.


Vejamos agora outro tipo de factorização matricial que quando usada na construção da solução de umsistema induz sistemas em que, no caso quadrado, os números de condição das matrizes é igual à unidaderelativamente à norma ‖·‖2.

Comecemos por considerar o caso quadrado não singular. Seja A uma matriz de ordem n não singular.A matriz AtA é uma matriz quadrada simétrica e positiva semidefinida. Logo os seus valores próprios sãonão negativos. Seja λ = 0 um valor próprio de AtA e x 6= 0 o vector próprio associado. Então

0 = xtAtAtx = ‖Ax‖22,

ou seja, Ax = 0 com x 6= 0. Este facto permite-nos concluir que A é singular o que contraria a hipótese.Assim os valores próprios de AtA são positivos.

Atendendo a que AtA é real simétrica, existe uma matriz ortogonal V de ordem n tal que

V tAtAV = D,

em que D é a matriz dos valores próprios de AtA.Seja S a matriz diagonal cujas entradas não nulas são

√λi. Da última decomposição vem

V tAtAV S−1 = S.

Tomemos U = AV S−1. A matriz U é facilmente determinada a partir da matriz V

u j =1√λ j

Av j, j = 1, . . . ,n.

A matriz U é ortogonal pois tem-se sucessivamente

U tU = S−1V tAtAV S−1

= S−1DS−1 = I.

As matrizes U,V e S permitem estabelecer a factorização seguinte

A =USV t . (2.2.3)

Introduzimos algumas designações: os escalares√

λ j, j = 1, . . . ,n, são chamados valores singulares damatriz A, os vectores v j, j = 1, . . . ,n, são chamados vectores singulares de A à direita enquanto que osvectores u j, j = 1, . . . ,n, são chamados vectores singulares de A à esquerda.


Teorema 2.3. Se A é não singular, então A admite a factorização em valores singulares (2.2.3)em que U é a matriz dos vectores singulares de A à esquerda, V é a matriz dos vectores singularesde A à direita e S é a matriz dos valores singulares de A.

A aplicação da factorização em valores singulares na resolução de sistemas de equações linearesconduz-nos a sistemas em que o número de condição da matrizes desses sistemas é igual à unidade. Defacto, a solução do sistema Ax = b, em que A admite a factorização (2.2.3), é feita seguindo os passosseguintes:

• Determinar a solução de Uz = b,

• Determinar a solução de Sy = z,

• Determinar a solução de V tx = y.

É de salientar que cond2(U) = cond2(V t) = 1.


No que segue generalizamos o Teorema 2.3 a matrizes não quadradas. Provemos o seguinte teorema:

Teorema 2.4. Seja A uma matriz do tipo m× n, m ≤ n, com característica r. Então existemnúmeros não negativos

λ1 ≥ λ2 ≥ ·· · ≥ λr > 0 = λr+1 = · · ·= λn = 0, (2.2.4)

e matrizes ortogonais U de ordem m e V de ordem n tais que

A =USV t . (2.2.5)

em que S é a matriz diagonal do tipo m×n, s j j =√

λ j, j = 1, . . . ,r, e com as restantes entradasnulas.

Demonstração. Seja A uma matriz do tipo m×n, m ≤ n, com característica r. Consideremos a matrizquadrada de ordem n, AtA. Esta matriz é real simétrica positiva semidefinida. Sejam V = [v1 · · ·vn] amatriz ortogonal dos vectores próprios e sejam λ j ≥ 0, j = 1, . . . ,n, os correspondentes valores próprios.Provemos que A e AtA têm o mesmo núcleo, isto é,

N (AtA) = N (A).

Seja x∈N (AtA). Então xtAtAx = 0, e portanto ‖Ax‖2 = 0. Desta última igualdade vem Ax = 0 e portantoN (AtA)⊆N (A). A outra inclusão é imediata.

Atendendo a que dim(N (A)) = n− r, tem-se dim(N (AtA)) = n− r e portanto car(AtA) = r. Estefacto permite concluir que AtA tem r valores próprios positivos que podemos supor ordenados segundo(2.2.4). Definamos os seguintes r vectores u j, j = 1, . . . ,r, como

u j =1√λ j

Av j, j = 1, . . . ,r.

Os vectores anteriores são ortonormados. De facto, tem-se sucessivamente

(ui,u j) =1√λiλ j

(Avi,Av j)

=

√λ j

λi(vi,v j)

=

√λ j

λiδi j.

Se r < m, existem vectores ur+1, · · · ,um em Rm tais que U = [u1 · · ·um] é uma matriz ortogonal. É claroA,U,S e V verificam a igualdade (2.2.5).

Os vectores u j e v j, j = 1, . . . ,r, são chamados vectores singulares de A à esquerda e direita, respecti-vamente. Os escalares

√λ j, j = 1, . . . ,r, são chamados valores singulares de A. A factorização (2.2.5),

designada decomposição em valores singulares de A, pode ser usada para calcular Ax, x ∈ Rn. De facto,uma vez que v j, j = 1, . . . ,n é uma base ortonormada de Rn, x admite a decomposição

x =n

∑j=1

(x,v j)v j.


Logo

Ax =r

∑j=1

√λ j(x,v j)u j.

A decomposição em valores singulares permite-nos averiguar da possibilidade ou impossibilidade deum sistema bem como determinar a solução, caso exista.

Teorema 2.5. Seja A uma matriz do tipo m×n, característica r e com a decomposição em valoressingulares (2.2.5). Então o sistema Ax = y é possível se e só se (y,z) = 0 para todo z ∈N (At).

Caso Ax = y seja possível, uma solução admite a seguinte representação

x =n

∑j=1

1√λ j

(y,u j)v j. (2.2.6)

Demonstração.

1. Condição necessária: Suponhamos que Ax = y é possível. Para z ∈N (At), temos

(y,z) = (Ax,z) = (x,Atz) = 0.

Se y é ortogonal a N (At), vem

y =r

∑j=1

(y,u j)u j. (2.2.7)

Se x = ∑nj=1 γ jv j, então

y = Ax =n

∑j=1

γ j

√λ ju j,

e portanto, atendendo à unicidade da representação (2.2.7), concluímos

γ j

√λ j = (y,u j), j = 1, . . . ,r, e γ j = 0, j > r,

ou seja, vale a representação (2.2.6).

2. Condição suficiente: Suponhamos que y é ortogonal a N (At). Então y admite a representação(2.2.7). Seja x dado por (2.2.6). Atendendo a que se tem sucessivamente

Ax =r

∑j=1

1√λ j

(y,u j)Av j

=r

∑j=1

(y,u j)u j = y,

concluímos a demonstração.

O resultado anterior permite introduzir o conceito de pseudo-inversa ou inversa de Moore-Penrose deuma matriz. Consideremos o seguinte operador A+ : Rm −→ Rn definido por

A+y =r

∑j=1

1√λ j

(y,u j)v j.


O operador A+ ou a sua matriz

V[

S−1r 00 0

]U t ,

que também representamos por A+, é designado operador pseudo-inverso ou inverso de Moore-Penrose.

A implementação do cálculo da decomposição em valores singulares de uma matriz é umprocedimento relativamente complexo e que sai fora do âmbito deste curso. Mais detalhessobre esta construção e a sua concretização em MATLAB em [5]. Finalmente, observamos quea decomposição em valores singulares pode ser calculada em MATLAB atrav1és da instrução[U,S,V]=svd(A).

2.3 Métodos Iterativos

Os métodos numéricos estudados na secção anterior são designados métodos directos e permitemobter a solução do sistema Ax = b num número finito de operações. Nesta secção pretendemos estudarmétodos que permitem construir uma sucessão de vectores que deverá convergir para a solução do sistemae que são chamados métodos iterativos. Os métodos iterativos admitem a expressão geral

x(m+1) = B(m)x(m)+ c(m), m ∈ N0, (2.3.1)

sendo portanto necessário uma aproximação x(0) dita aproximação inicial. Se B(m) e c(m) são indepen-dentes de m, então o método iterativo diz-se estacionário. Caso contrário, diz-se não estacionário. Nasecção seguinte iremos estudar alguns métodos estacionários. Alguns métodos não estacionários serãoestudados na secção 2.3.2.

2.3.1 Métodos Estacionários

Consideremos o sistema Ax = b em que A ∈Mn(R), e os elementos diagonais de A são não nulos.Os métodos que consideramos seguidamente são estabelecidos a apartir da decomposição

A = D−L−U, (2.3.2)

em que D é a matriz diagonal de entradas aii 6= 0, i = 1, . . . ,n, L é a matriz triangular inferior, comelementos diagonais nulos e ì j =−ai j, i > j, e U é a matriz triangular superior com elementos diagonaisnulos e ui j = −ai j, i < j. Observamos que em MATLAB, as matrizes D, L e U definidas podem sercalculadas usando as instruções diag, tril e triu. De facto, D=diag(diag(A)), L=−tril(A,−1)e U=−triu(A,1).

Método de Jacobi

Consideremos em Ax = b a decomposição (2.3.2). Obtemos o sistema equivalente

x = D−1(L+U)x+D−1b. (2.3.3)

Seja x(0) uma aproximação inicial e seja x(1) o vector que se obtém quando substituímos, no segundomembro de (2.3.3), x por x(0), isto é, seja x(1) definido por

x(1) = D−1(L+U)x(0)+D−1b.


Considerando agora x(1) obtemos x(2)

x(2) = D−1(L+U)x(1)+D−1b.

O método introduzido é designado método de Jacobi e é tem a seguinte expressão

x(m+1) = D−1(L+U)x(m)+D−1b. (2.3.4)

As componentes do vector x(m+1) são facilmente obtidas a partir das componentes de x(m)

x(m+1)i =−

n

∑i6= j=1

ai j

aiix(m)

j +bi

aii, i = 1, . . . ,n.

O método de Jacobi transcreve-se para a linguagem de programação do MATLAB de formaimediata, como vemos pelo algoritmo seguinte:

Algoritmo 2.8 Método de Jacobi (com controlo do número máximo de iterações).

function [x]=jacobi(A,b,x0,nmax)D = diag(diag(A)); L = −tril(A,−1); U = −triu(A,1);

iter=0; x=x0;for i=1:nmax

x = D\((L+U)*x + b);end

Observamos a introdução da variável de entrada nmax, utilizada para limitar o número deiterações do método. Veremos mais tarde outros critérios que permitam terminar o ciclo dométodo. Estes critérios geralmente envolvem a monitorização quando é que uma determinadaquantidade (por exemplo, a norma do resíduo r(m) = Ax(m)− b) atinge uma tolerância decontrolo predefinida.

O método (2.3.4) pode ser ainda escrito na forma

x(m+1) = (I−D−1A)x(m)+D−1b,

o que permite uma implementação mais simples e eficiente pois não há necessidade de calcularas matrizes U e L.

Método de Gauss-Seidel

Se em vez da igualdade (2.3.3) considerarmos

x = (D−L)−1Ux+(D−L)−1b. (2.3.5)

somos conduzidos ao conhecido método de Gauss-Seidel

x(m+1) = (D−L)−1Ux(m)+(D−L)−1b. (2.3.6)

Se substituirmos a relação anterior por

(D−L)x(m+1) =Ux(m)+b,


obtemos explicitamente as componentes do vector da iteração m+1

x(m+1)1 =−

n

∑j=2

a1 j

a11x(m)

j +b1

a11,

x(m+1)i =−

i−1

∑j=1

ai j

aiix(m+1)

j −n

∑j=i+1

ai j

aiix(m)

j +bi

aii, i = 2, . . . ,n

O método de Gauss-Seidel pode ser obtido a partir do método de Jacobi se na definição da componente iiteração m+1 substituirmos as componentes j, para j < i, por x(m+1)

j .

A implementação do método de Gauss-Seidel segue as mesmas linhas do Algoritmo 2.8:

Algoritmo 2.9 Método de Gauss-Seidel (com controlo do número máximo de iterações).

function [x]=gaussSeidel(A,b,x0,nmax)D=diag(diag(A)); L = −tril(A,−1); U = −triu(A,1);

iter=0; x=x0;for i=1:nmaxx = (D−L)\(U*x + b);

end

Embora tenhamos usado o comando \ do MATLAB , o sistema com matriz D− L étriangular inferior e portanto, poderíamos ter usado também os algoritmos estudados na secção2.2.1.

Exemplo 2.2. Ilustramos seguidamente o comportamento dos dois métodos anteriores considerando osistema não singular

4 1 1−1 3 01 0 2

x1x2x3

=

623

,

cuja solução é o vector x = [1 1 1]t . Na Figura 2.1 ilustramos o comportamento do erro∥∥e(m)

∥∥1 =

∥∥x− x(m)∥∥

1, m = 0, . . . ,10,

quando consideramos os métodos de Jacobi (2.3.4) e Gauss-Seidel (2.3.6) com a aproximação inicialnula. Observamos que o erro diminui mais acentuadamente quando consideramos o segundo método.

Métodos de relaxação

O método de Jacobi (2.3.4) pode ser reescrito na forma equivalente

x(m+1) = x(m)+D−1(b−Ax(m)). (2.3.7)

O termo D−1(b−Ax(m)) pode ser interpretado como uma correcção à m-ésima aproximação. Se consi-derarmos que esta correcção é do tipo ωD−1(b−Ax(m)), com ω > 0, obtemos o método de Jacobi comrelaxação

x(m+1) = x(m)+ωD−1(b−Ax(m)), (2.3.8)


0 2 4 6 8 1010−13

10−10

10−7

10−4

10−1

número de iterações

∥ ∥ e(m

)∥ ∥1

JacobiGauss-Seidel

Figura 2.1: Comportamento do erro ‖e(m)‖1 associado à aproximação x(m) definida pelos métodos deJacobi e Gauss-Seidel.

ou ainda, componente a componente,

x(m+1)i = x(m)

i +ω1aii

(bi−

n

∑j=1

ai jx(m)j

), i = 1, . . . ,n. (2.3.9)

O método anterior admite a representação equivalente

x(m+1) = D−1((1−ω)D+ω(L+U)

)x(m)+ωD−1b, m = 0, . . . . (2.3.10)

O método de Gauss-Seidel (2.3.6) uma representação análoga a (2.3.7)

x(m+1) = x(m)+D−1(b+Lx(m+1)− (D−U)x(m)).

Se D−1(b+Lx(m+1)− (U −D)x(m)

)for interpretado como uma correcção de x(m), então considerando

a correcção anterior com peso ω > 0, obtemos o método de Gauss-Seidel com relaxação ou método desobre-relaxação sucessiva

x(m+1) = x(m)+ωD−1(b+Lx(m+1)− (D−U)x(m)), (2.3.11)

ou ainda, componente a componente,

x(m+1)i = x(m)

i +ω1aii

(bi−

i−1

∑j=1

ai jx(m+1)j −

n

∑j=i

ai jx(m)j

), i = 1, . . . ,n. (2.3.12)

O método de Gauss-Seidel com relaxação admite a representação equivalente

x(m+1) = (D−ωL)−1((1−ω)D+ωU)x(m)+ω(D−ωL)−1b, m = 0, . . . . (2.3.13)

Exemplo 2.3. A título ilustrativo, consideremos o sistema[

2 1−1 3

][xy

]=

[12

],


0 2 4 6 8 10

10−4

10−3

10−2

10−1

100

número de iterações

∥ ∥ e(m

)∥ ∥1

ω = 12

ω = 1

Figura 2.2: Comportamento do erro∥∥e(m)

∥∥1 relativo à aproximação x(m) definida pelo método (2.3.11)

para diferentes valores de ω .

cuja solução é o vector x =(1

7 ,57

)t . Na Figura 2.2 ilustramos o comportamento do erro∥∥e(m)

∥∥1,

m = 0, . . . ,10, relativamente à escolha do parâmetro ω, quando consideramos o método (2.3.11) com aaproximação inicial nula. O aumento do parâmetro de relaxação induz uma diminuição do erro

∥∥e(m)∥∥

1.

Observamos que do ponto de vista computacional, ambos os métodos de Jacobi e deGauss-Seidel com relaxação podem ser facilmente implementados seguindo as linhas dosAlgoritmos 2.8 e 2.9.

Convergência

O teorema que apresentamos seguidamente estabelece uma condição necessária e suficiente para aconvergência do método iterativo (2.3.1) quando B(m) = B.

Começamos por notar que o erro na iteração m é definido por

e(m) = x− x(m).

O método iterativo estacionário x(m+1) = Bx(m)+ c, m ∈ N0, diz-se convergente se

limm→∞

∥∥e(m)∥∥= 0.

Observamos que se existir convergência relativamente a uma norma, então, atendendo a que em Rn

todas as normas são equivalentes, concluímos convergência relativamente a qualquer outra norma.

Teorema 2.6. O método iterativo

x(m+1) = Bx(m)+ c, m ∈ N0,


é convergente para todo c ∈ Rn, se e só se ρ(B)< 1.

Demonstração.

• Condição suficiente: Se ρ(B)< 1, então existe em Rn, uma norma induzida ‖·‖ tal que ‖B‖< 1 (ademonstração pode ser encontrada em [5]). Provemos que a sucessão

(x(m)

)mN é de Cauchy em

Rn. Temos, para m≥ p,∥∥x(m)− x(p)

∥∥≤m−1

∑j=p

∥∥x( j+1)− x( j)∥∥,

e ∥∥x( j+1)− x( j)∥∥≤ ‖B‖ j

∥∥x(1)− x(0)∥∥. (2.3.14)

Logo∥∥x(m)− x(p)

∥∥≤∥∥x(1)− x(0)

∥∥m−1

∑j=p‖B‖ j,

≤∥∥x(1)− x(0)

∥∥‖B‖p∞

∑j=0‖B‖ j

≤∥∥x(1)− x(0)

∥∥ ‖B‖p

1−‖B‖ .

Tomando na desigualdade anterior p→+∞, concluímos

limp→∞

∥∥x(m)− x(p)∥∥= 0,

isto é,(x(m)

)mN é uma sucessão de Cauchy e portanto é convergente.

• Condição necessária: Se ρ(B)≥ 1, então tomemos c = x(0) = x, em que x é o vector próprio asso-ciado ao valor próprio λ tal que ρ(B) = |λ |. Uma vez que x(m) admite, neste caso, a representação

x(m) =

( m

∑j=0

λ j)

x,

concluímos que x(m) é divergente.

Suponhamos que o método iterativo converge e seja x(m)→ x. Atendendo à definição da sucessão x(m),concluímos

x = Bx+ c.

Seja m > p. Como vimos anteriormente, se ‖B‖< 1, então

∥∥x(m)− x(p)∥∥≤

∥∥x(1)− x(0)∥∥ ‖B‖

p

1−‖B‖ .

Tomando na desigualdade anterior m→ ∞, obtemos

∥∥x− x(p)∥∥≤

∥∥x(1)− x(0)∥∥ ‖B‖

p

1−‖B‖ . (2.3.15)

O resultado seguinte resume as conclusões anteriores.


Teorema 2.7. Se A = [ai j]ni, j=1 é tal que

q∞ = maxi=1,...,n

n

∑i 6= j=1

∣∣∣∣ai j

aii

∣∣∣∣< 1 ou q1 = maxj=1,...,n

n

∑j 6=i=1

∣∣∣∣ai j

aii

∣∣∣∣< 1,

então, para qualquer b ∈Rn, a sucessão(x(m)

)m∈N definida pelo método de Jacobi, qualquer que

seja a aproximação inicial x(0), converge para a solução do sistema Ax = b e vale a estimativa

∥∥x− x(p)∥∥≤

∥∥x(1)− x(0)∥∥ qp

µ

1−qµ, µ = 1,∞. (2.3.16)

Demonstração. A prova deste resultado é imediata pois a matriz do método de Jacobi B = D−1(L+U)verifica ‖B‖µ = qµ , µ = 1,∞.

Exemplo 2.4. Retomemos o sistema do Exemplo 2.2. A matriz do método de Jacobi, para este sistema, édada por

B =

0 −14 −1

413 0 0−1

2 0 0

,

e ‖B‖∞ = 12 . Logo a sucessão gerada pelo método de Jacobi é convergente e verifica-se para o erro a

estimativa (2.3.15). Uma vez que x(1) = (3/2,2/3,3/2)t , concluímos

∥∥x− x(p)∥∥≤ 3

2p , p≥ 0.

Se pretendermos determinar uma aproximação para a solução tal que o erro é inferior a uma determinadaquantidade ε, então devemos calcular p tal que

32p < ε.

Utilizando o Teorema 2.7 podemos estabelecer uma nova condição suficiente para a convergência dométodo de Jacobi. Para o efeito introduzimos o conceito de matriz diagonal dominante.

Definição 2.1. Se A = [ai j]ni, j=1 é tal que

|aii| ≥n

∑i 6= j=1

|ai j|, i = 1, . . . ,n, (2.3.17)

então A diz-se diagonal dominante por linhas. Se, em (2.3.17), todas as desigualdades são estritas,então A diz-se estritamente diagonal dominante por linhas. Se

|a j j| ≥n

∑j 6=i=1

|ai j|, j = 1, . . . ,n, (2.3.18)

então A diz-se diagonal dominante por colunas. Se, em (2.3.18), todas as desigualdades são estritas,então A diz-se estritamente diagonal dominante por colunas.


Como corolário do Teorema 2.7, vale o seguinte resultado:

Corolário 2.1. Se A é estritamente diagonal dominante por linhas, então o método de Jacobiconverge.

Consideramos seguidamente uma condição suficiente para a convergência do método de Gauss-Seidel.Para concluir a convergência deste método, relativamente à norma ‖·‖∞, temos que estabelecer umacondição suficiente que nos permita concluir

∥∥(D−L)−1U∥∥

∞= sup‖z‖∞=1

∥∥(D−L)−1Uz∥∥

∞< 1. (2.3.19)

Seja x = (D−L)−1Uz. Então

• |x1| ≤ ∑nj=2

∣∣∣∣a1 ja11

∣∣∣∣,

• |x2| ≤∣∣∣∣

a21a22

∣∣∣∣|x1|+∑nj=3

∣∣∣∣a2 ja22

∣∣∣∣,

...

• |xn| ≤ ∑n−1j=1

∣∣∣∣an jann

∣∣∣∣|x j|.

Sejam p j, j = 1, . . . ,n, definidos recursivamente por

• p1 = ∑nj=2

∣∣∣∣a1 ja11

∣∣∣∣,

• pi = ∑i−1j=1

∣∣∣∣ai jaii

∣∣∣∣p j +∑nj=i+1

∣∣∣∣ai jaii

∣∣∣∣, i = 2, . . . ,n−1,

• pn = ∑n−1j=1

∣∣∣∣an jann

∣∣∣∣p j.

Se

p = maxi=1,...,n

pi < 1, (2.3.20)

então ‖x‖∞ < 1, e portanto, de (2.3.19), concluímos

∥∥(D−L)−1U∥∥

∞< 1.


O teorema seguinte resume o resultado demonstrado.

Teorema 2.8. Se A = [ai j] é tal que vale a condição (2.3.20), então o método de Gauss-Seidelconverge qualquer que seja a aproximação inicial x(0) e tem-se

∥∥x− x(m)∥∥≤

∥∥x(1)− x(0)∥∥ pm

1− p. (2.3.21)

Como corolário do teorema anterior, vale o seguinte resultado.

Corolário 2.2. Se A é estritamente diagonal dominante por linhas, então o método de Gauss-Seidel é convergente qualquer que seja a aproximação inicial x(0).

Observamos que a matriz tridiagonal

A =

2 −1−1 2 −1

. . . . . . . . .−1 2 −1

−1 2

(2.3.22)

não é estritamente diagonal dominante mas satisfaz as hipóteses do Teorema 2.8.

No caso de uma matriz com um número de linhas conhecido, a determinação do valor p doteorema anterior pode ser feita com o auxílio do seguinte programa:

N=100;A= 2*diag(ones(1,N)) − diag(ones(1,N−1),1) − diag(ones(1,N−1),−1);

p = zeros(1,N);p(1) = sum(abs(A(1,2:end)/A(1,1)));for i=2:N−1p(i) = sum( abs( A(i,1:i−1)/A(i,i) ).*p(1:i−1) ) + ...

sum(abs(A(i,i+1:N)/A(i,i)));endp(N) = sum(abs((A(N,1:N−1).*p(1:N−1))/A(N,N)));

maximo = max(p)

O valor na variável maximo, que no caso de considerarmos a matriz (2.3.22) com 100linhas, é igual a 0.998046875. Este valor indica que embora o método seja convergente, paragarantirmos uma certa precisão na solução de acordo com a estimativa (2.3.21), temos deefectuar um número considerável de iterações do método de Gauss-Seidel, uma vez p estámuito próximo de 1.


A determinação de uma aproximação para a solução do sistema Ax = b com uma determinada precisãopode ser feita utilizando as estimativas para o erro anteriormente estabelecidas.

Observação 2.1 Observamos que a velocidade de convergência de um método iterativo está associada aoraio espectral da sua matriz de iteração. De facto, consideremos dois métodos iterativos convergentes

x(m+1)1 = B1x(m)

1 + c1, m = 0, . . . ,

ex(m+1)

2 = B2x(m)2 + c2, m = 0, . . . ,

e sejam e(m)1 = x∗− x(m)

1 , e e(m)2 = x∗− x(m)

2 os correspondentes erros absolutos. Notamos que se tem

e(m+1)1 = B1e(m)

1 , e(m+1)2 = B2e(m)

2 , m = 0, . . . .

Sejam λi, i = 1, . . . ,n, os valores próprios de B1 e βi, i = 1, . . . ,n, os valores próprios de B2. Sejamvi, i = 1, . . . ,n e ui, i = 1, . . . ,n, os correspondentes vectores próprios que admitimos, sem perda degeneralidade, constituirem conjuntos ortonormados. Atendendo a

e(0)1 =n

∑j=1

ξ jv j, e(0)2 =n

∑j=1

γ ju j,

deduzimos

e(m+1)1 =

n

∑j=1

ξ jλ m+1j v j, e(m+1)

2 =n

∑j=1

γ jβ m+1j u j,

ou ainda ∥∥e(m+1)1

∥∥22 =

n

∑j=1

ξ 2j |λ j|2(m+1),

∥∥e(m+1)2

∥∥22 =

n

∑j=1

γ2j |β j|2(m+1).

As duas últimas igualdades permitem estabelecer, na hipótese de γk 6= 0 e |βk|= ρ(B2),

∥∥e(m+1)1

∥∥22∥∥e(m+1)

2

∥∥22

≤ρ(B1)

2(m+1)∑

nj=1 ξ 2

j

ρ(B2)2(m+1)γ2k

.

Assim, uma vez que ρ(B1)< ρ(B2), vem

limm→+∞

∥∥e(m+1)1

∥∥22∥∥e(m+1)

2

∥∥22

= 0,

e portanto, a partir de certa ordem, ∥∥e(m+1)1

∥∥2 ≤

∥∥e(m+1)2

∥∥2.

A desigualdade anterior, mostra-nos que o erro induzido por um método iterativo cuja matriz de iteraçãotem menor raio espectral é inferior, pelo menos a partir de certa ordem, ao induzido por um método commatriz de iteração com maior raio espectral.

Exemplo 2.5. Retomemos o Exemplo 2.2 e apliquemos o método de Jacobi com relaxação com parâmetroω > 0. Na Figura 2.3, traçamos a dependência do raio espectral da matriz B(ω) = (1−ω)I+ωD−1(L+U) de iteração deste método relativamente a ω .

Como ρ(B(1

2

))> ρ(B(1)), o método de Jacobi com relaxação converge mais rapidamente para a

solução do sistema no caso ω = 1 do que no caso ω = 12 . Esta conclusão está também de acordo com os

resultados da Figura 2.2.


0 0.5 1 1.5 2

0.2

0.4

0.6

0.8

1

1.2

1.4

ω

ρ(B(ω

))

Figura 2.3: Gráfico de ρ(B(ω)) para a matriz do Exemplo 2.2, para ω ∈ [0,2].

O Exemplo anterior motiva para a importância da escolha do parâmetro de relaxação nos métodos deJacobi e Gauss-Seidel. É claro que o parâmetro de relaxação deve ser escolhido de modo a que o raioespectral das matrizes de iteração dos métodos anteriores tenham raio espectral inferior à unidade. Maisainda, dentro do conjunto de parâmetros de relaxação admissíveis, deve ser escolhido o parâmetro quetorne mí nimo tal raio espectral.

Teorema 2.9. Se a matriz de iteração B = D−1(L+U) do método de Jacobi tem apenas valorespróprios reais inferiores a 1, então o raio espectral da matriz B(ω) = (1−ω)I +ωD−1(L+U)de iteração do método de Jacobi com relaxação é mínimo em ωopt definido por

ωopt =2

2−λmax−λmin,

e

ρ(B(ωopt)) =λmax−λmin

2−λmax−λmin,

em que λmax e λmin são os valores próprios máximo e mínimo de B, e, se λmax +λmin 6= 0, entãoρ(B(ωopt))< ρ(BJ).

Demonstração. Se λ é valor próprio de B e x o seu vector próprio correspondente, então 1−ω +ωλ évalor próprio de B(ω) e x é o vector próprio associado. Reciprocamente, se β é valor próprio de B(ω) e xo vector próprio correspondente, então β−1+ω

ω é valor próprio de B e x o vector próprio correspondente.Logo os valores próprios de B(w) são reais e

βmax = 1−ω +ωλmax, βmin = 1−ω +ωλmin.

Se o maior e menor valores anteriores são iguais em módulo, então o raio espectral toma o seu valormínimo. Se ωopt é tal

1−ωopt +ωoptλmin =−1+ωopt −ωoptλmax,


então, para ω = ωopt , ρ(B(ω)) toma o seu valor mínimo.

Exemplo 2.6. Notamos que a matriz do Exemplo 2.2 está nas condições do teorema anterior pois osvalores próprios da matriz de iteração do método de Jacobi são

λ1 = 0, λ2 =

√6

12e λ3 =−

√6

12.

Logo, o parâmetro de relaxação que minimiza o ρ(B(ω)) é ωopt = 1, o que também é confirmado pelográfico da Figura 2.3.

No resultado seguinte estabelecemos uma condição necessária para a convergência do método deGauss-Seidel com relaxação.

Teorema 2.10. Se o método de Gauss-Seidel com relaxação converge, então ω ∈ (0,2).

Demonstração. Seja B(ω) a matriz de iteração do método de Gauss-Seidel com relaxação. Notamos que

det(B(ω)) = ∏i

λi,

em que λi, i = 1, . . . ,n, representam os valores próprios de B(ω). Mas det(B(ω)) = (1−ω)n, e portanto,

(1−ω)n = ∏λi.

Se o método é convergente, então ρ(B(ω))< 1, o que implica

|1−ω|< 1.

Estabelecemos finalmente, sem demonstração, o teorema de convergência do método de Gauss-Seidelcom relaxação:

Teorema 2.11 ([3]). Se A é simétrica e positiva definida, então o método de Gauss-Seidel comrelaxação converge qualquer que seja a aproximação inicial x(0) e qualquer que seja ω ∈ (0,2).

2.3.2 Métodos não Estacionários

Método da descida mais rápida

Consideremos o sistema Ax = b em que A ∈Mn(R). A determinação da solução x∗ do problemaanterior pode ser feita determinando a solução x∗ ∈ Rn de um problema de minimização, isto é, x∗ ∈ Rn

tal queΦ(x∗) = min

x∈RnΦ(x), (2.3.23)

em que

Φ : Rn −→ R, Φ(x) =12

xtAx− xtb.

Seja (x,Φ(x)) um ponto da superfície z = Φ(x). No ponto (x,Φ(x)) o vector gradiente ∇Φ(x) definea direcção e sentido segundo o qual Φ apresenta, localmente, o maior crescimento. Logo −∇Φ(x) definea direcção e o sentido segundo o qual ocorre, localmente, a maior decrescimento de Φ.


Exemplo 2.7. Seja A= 100I2 e b= 0. Então Φ(x1,x2) = 100(x21+x2

2) e ∇Φ(x,y) = 100(x1,x2). No ponto(1,1) o vector (100,100) define a direcção e sentido em que Φ cresce e (−100,−100) define a direcçãoe sentido de maior decrescimento de Φ. Este comportamento é ilustrado na figura 2.4.

−2 −1 0 1 2 −2

0

20

500

−2 −1 0 1 2−2

−1

0

1

2

Figura 2.4: Superfície z = 100(x2 + y2).

A equivalência entre os dois problemas pode ser demonstrada quando A é simétrica e positiva definida.

Teorema 2.12. Seja A ∈Mn(R) simétrica e positiva definida. x∗ ∈ Rn é solução do sistemaAx = b se e só se x∗ satisfaz (2.3.23).

Demonstração. Se x∗ satisfaz (2.3.23), então ∇Φ(x∗) = 0. Mas ∇Φ(x) = 12(A

t +A)x− b, e, como A ésimétrica, concluímos Ax∗ = b.

Suponhamos agora que Ax∗ = b. Tem-se sucessivamente, com z = x∗− y,

Φ(x∗) = Φ(y+ z)

= Φ(y)+12

ztAy+12

ytAz− ztb+12

ztAz

= Φ(y)+ zt(Ay−b)+12

ztAz

= Φ(y)− 12

ztAz.

Uma vez que A é positiva definida, concluímos

Φ(x∗)≤Φ(y), ∀y ∈ Rn

sendo a igualdade verificada apenas quando y = x∗.

Consideramos agora o seguinte esquema iterativo

x(m+1) = x(m)+α(m)d(m), m = 0,1, . . . , (2.3.24)

em que d(m) e α(m) serão definidos de modo conveniente. Notamos que o vector α(m)d(m) pode servisto como uma correcção à iteração x(m). Assim, com base nas considerações anteriores sobre −∇Φ(x)


representar a direcção de maior decrescimento da função Φ no ponto x, consideremos as direcções d(m)

comod(m) =−∇Φ

(x(m)

).

Observamos que no caso de A ser simétrica, então

∇Φ(x(m)

)= Ax(m)−b,

e portanto

d(m) = r(m), (2.3.25)

em que r(m) é o vector dos resíduos da aproximação x(m).O coeficiente α(m) é escolhido de modo a que Φ

(x(m+1)

)tenha, neste parâmetro, o valor mínimo, isto é,

de modo a ser válida a seguinte igualdade

ddα(m)

(Φ(x(m)+α(m)d(m))

)= 0.

Uma vez que a equação anterior é equivalente à equação

∇Φ(x(m)+α(m)r(m)

)tr(m) = 0,

ou ainda, a(Ax(m)+α(m)Ar(m)−b)tr(m) = 0,

obtemos

α(m) =r(m)tr(m)

r(m)tAr(m). (2.3.26)

O método de descida (também designado método da descida mais rápida) é definido pela seguintesequência de passos: dado x(0) ∈ Rn, determinar, para m = 0, . . . ,

1. d(m) = r(m) = b−Ax(m),

2. α(m) =r(m)tr(m)

r(m)tAr(m),

3. x(m+1) = x(m)+α(m)r(m).

Dado um número máximo de iterações a realizar pelo método da descida, a solução aproxi-mada do sistema linear Ax = b pode ser determinada pelo seguinte algoritmo em MATLAB:

Algoritmo 2.10 Método da descida (com controlo do número máximo de iterações).

function [x]=metodoDescida(A,b,x0,nmax)

x = x0;for i=1:nmaxr = b − A*x;alpha = (r'*r)/(r'*A*r);x = x + alpha*r;

end

return


O método da descida mais rápida é um método não estacionário. De facto, o esquema iterativo (2.3.24)admite a representação

x(m+1) = (I−α(m)A)x(m)+α(m)b,

com α(m) definido por (2.3.26).No resultado seguinte estabelecemos a convergência do método apresentado. Se A é simétrica e

positiva definida, então os valores próprios de A são reais positivos. Sejam λmax e λmin os valores própriosde A, máximo e mínimo, respectivamente. São válidas as desigualdades

1≤ (y,Ay)(y,A−1y)(y,y)(y,y)

≤ (λmax +λmin)2

4λmaxλmin, (2.3.27)

para todo y ∈ Rn−0, e que são usualmente designadas desigualdades de Kantorovich. A convergênciado método é estabelecida relativamente à norma ‖·‖A e que é definida por ‖x‖A =

√xtAx.

Teorema 2.13. Se A ∈Mn(R) é simétrica definida positiva, então os erros e(m+1) = x∗− x(m+1),e(m) = x∗−x(m), em que x(m+1),x(m) são definidos por (2.3.24) com d(m) e α(m) dados por (2.3.25)e (2.3.26), respectivamente, satisfazem

∥∥e(m+1)∥∥2A ≤

(λmax−λmin

λmax +λmin

)2∥∥e(m)∥∥2

A, m = 0, . . . . (2.3.28)

Demonstração. Os erros e(m+1) e e(m) verificam a igualdade seguinte

e(m+1) = e(m)−α(m)r(m).

Tem-se, sucessivamente,

∥∥e(m+1)∥∥2A =

(e(m+1),Ae(m+1))

=(e(m)−α(m)r(m),A

(e(m)−α(m)r(m)

))

=(e(m),Ae(m)

)−2α(m)

(e(m),Ar(m)

)+α(m)2(r(m),Ar(m)

)

=∥∥e(m)

∥∥2A−2α(m)

(r(m),r(m)

)+α(m)2(r(m),Ar(m)

)

=∥∥e(m)

∥∥2A−2

(r(m),r(m)

)(r(m),Ar(m)

)(r(m),r(m))+

( (r(m),r(m)

)(r(m),Ar(m)

))2(

r(m),Ar(m))

=∥∥e(m)

∥∥2A−

(r(m),r(m)

)2

(r(m),Ar(m)

)

=∥∥e(m)

∥∥2A

(1−

(r(m),r(m)

)2

(r(m),Ar(m)

)(e(m),Ae(m)

))

Atendendo a que se tem (e(m),Ae(m)

)=(e(m),r(m)

)=(A−1r(m),r(m)

),

obtemos

1≥(r(m),r(m)

)2

(r(m),Ar(m)

)(e(m),Ae(m)

) =(r(m),r(m)

)2

(r(m),Ar(m)

)(A−1r(m),r(m)

) ≥ 4λmaxλmin

(λmax +λmin)2 ,


e portanto∥∥e(m+1)∥∥2

A ≤∥∥e(m)

∥∥2A

(1− 4λmaxλmin

(λmax +λmin)2

)

=∥∥e(m)

∥∥2A(λmax−λmin)

2

(λmax +λmin)2

o que finaliza a demonstração.

Corolário 2.3. Se A ∈Mn(R) é simétrica definida positiva, então os erros e(m+1) = x∗− x(m+1),e(m) = x∗− x(m), em que x(m+1), x(m) são definidos por (2.3.24) com d(m) e α(m) dados por(2.3.25) e (2.3.26), respectivamente, satisfazem

∥∥e(m)∥∥2

A ≤(

cond2(A)−1cond2(A)+1

)2m∥∥e(0)∥∥2

A, m = 1, . . . . (2.3.29)

Nas condições do Teorema 2.13, concluímos a convergência do método da descida mais rápidadefinido por (2.3.24), (2.3.25) e (2.3.26).

Uma questão natural a que precisamos de dar resposta diz respeito à determinação do número mínimode iterações que devemos considerar de modo a que o erro seja inferior a uma quantidade definida a priori.A estimativa estabelecida no último resultado não permite responder a esta questão pois o erro inicial é,em geral, desconhecido. Estabelecemos no que se segue, uma estimativa para

∥∥e(m)∥∥

A que é independentede∥∥e(0)

∥∥. Seja

ξ =

(cond2(A)−1cond2(A)+1

)2

.

Observamos que, para qualquer ε 6= 0, tem-se

∥∥e(0)∥∥2

A ≤(∥∥e(1)

∥∥A +∥∥x(1)− x(0)

∥∥A

)2

≤ (1+ ε2)∥∥e(1)

∥∥2A +

(1+

1ε2

)∥∥x(1)− x(0)∥∥2

A.

De (2.3.28) vem ∥∥e(1)∥∥2

A ≤ ξ∥∥e(0)

∥∥2A,

e portanto (1− (1+ ε2)ξ

)∥∥e(1)∥∥2

A ≤ ξ(1+

1ε2

)∥∥x(1)− x(0)∥∥2

A.

Se fixarmos ε tal que

ε2 =12

1−ξξ

,

então∥∥e(1)

∥∥2A ≤

2ξ1−ξ

(1+

2ξ1−ξ

)∥∥x(1)− x(0)∥∥2

A.

Atendendo a que se tem ∥∥e(m+1)∥∥2A ≤ ξ m

∥∥e(1)∥∥2

A,


concluímos ∥∥e(m+1)∥∥2A ≤ ξ m 2ξ

1−ξ

(1+

2ξ1−ξ

)∥∥x(1)− x(0)∥∥2

A. (2.3.30)

Esta desigualdade pode ser utilizada para determinar o número mínimo de iterações que devemosconsiderar para que o erro da iteração m seja inferior a uma dada tolerância Tol, isto é,

∥∥e(m)∥∥

A ≤ Tol.

É suficiente tomar m tal que

ξ m 21−ξ

(1+

2ξ1−ξ

)∥∥x(1)− x(0)∥∥2

A ≤ Tol2.

Observamos que se conhecermos o valor ξ , então podemos reformular o algoritmo 2.10para parar quando o número de iterações atinge um certo valor nmax ou o erro

∥∥e(m)∥∥ na

iteração a ser calculada atinge uma determinada tolerância Tol:

Algoritmo 2.11 Método da descida (com controlo do número máximo de iterações e tolerânciapré-definida).

function [x,niter,estimErro] = metodoDescida(A,b,x0,xi,nmax,Tol)

x = x0;r = b − A*x;alpha = (r'*r)/(r'*A*r);x = x + alpha*r;aux = (2/(1−xi))*(1 + 2*xi/(1−xi)) * (x−x0)'*A*(x−x0);erro = sqrt(xi * aux);niter = 1;

while ( (niter < nmax) && (erro > Tol) )niter = niter + 1;r = b − A*x;alpha = (r'*r)/(r'*A*r);x = x + alpha*r;erro = sqrt(xi^niter * aux);

end

majoranteErro = erro;return

Em MATLAB , o valor de ξ pode ser calculado usando a função cond para matrizes dedimensões pequenas, ou uma sua estimativa, condest:

xi = ((cond(A)−1)/(cond(A)+1))^2;

Método dos gradientes conjugados

No método da descida mais rápida, as direcções d(m) são o vector resíduo que coincide com−∇Φ(x(m)). Seja A simétrica e positiva definida. Seguidamente consideramos direcções A-ortogonais,


isto é, direcções ortogonais relativamente ao produto interno induzido por A

(u,v)A = (u,Av).

Consideremos novamente o esquema iterativo

x(m+1) = x(m)+α(m)d(m), m = 0, . . . , (2.3.31)

onde agora as direcções d( j), j = 0, . . . , são escolhidas A-ortogonais. Seja α(m) tal como no caso dométodo da descida mais rápida, isto é, α(m) é escolhido de modo que Φ(x(m)+α(m)d(m)) seja mínimo.Este parâmetro é especificado determinando α(m) tal que

ddα(m)

(Φ(x(m)+α(m)d(m))

)= 0.

A equação anterior é equivalente a

∇Φ(x(m)+α(m)d(m))td(m) = 0,

e, uma vez que∇Φ(x(m)+α(m)d(m)) = A(x(m)+α(m)d(m))−b,

obtemos

α(m) =

(r(m),d(m)

)(d(m),d(m)

)A

. (2.3.32)

A determinação da nova direcção d(m+1) é feita considerando d(0) = r(0) e

d(m+1) = r(m+1)+β (m)d(m), m = 1, . . . , (2.3.33)

em que β (m) é calculado de modo a que d(m+1) e d(m) sejam A-ortogonais. Assim obtemos

β (m) =−(r(m+1),d(m))A(d(m),d(m)

)A

. (2.3.34)

Seja r(m+1) = b−Ax(m+1). Atendendo a que x(m+1) é definido por (2.3.31), obtemos para r(m+1) aexpressão

r(m+1) = r(m)−α(m)Ad(m). (2.3.35)

Provemos o seguinte resultado:

Teorema 2.14. Se a sucessão x( j) é definida por (2.3.31) com

α( j) =

(r( j),r( j)

)(d( j),d( j)

)A

, (2.3.36)

a sucessão das direcções d( j), m = 0,1,2 . . . , é definida por (2.3.33) e os coeficientes β ( j) sãodefinidos por (2.3.34), então r( j), j = 0,1,2 . . . ,m são ortogonais e d( j), j = 0,1,2, . . . ,msão A-ortogonais, em que m é o menor valor tal que r(m+1) = 0.


Demonstração. É fácil verificar que(r(1),r(0)

)= 0 e

(d(1),d(0)

)A = 0. Suponhamos que r( j), j =

0, . . . ,k são ortogonais e d( j), j = 0, . . . ,k são A-ortogonais. Provemos que r( j), j = 0, . . . ,k+1 sãoortogonais e d( j), j = 0, . . . ,k+1 são A-ortogonais.

Comecemos por provar que(r(k+1),r(k)

)= 0. Atendendo a que r(k+1) = r(k)−α(k)Ad(k),

(r(k+1),r(k)

)=(r(k),r(k)

)−α(k)

(Ad(k),r(k)

),

e portanto, utilizando a expressão (2.3.36), concluímos o pretendido.Provemos agora que

(r(k+1),r( j)

)= 0, j = 0, . . . ,k. É fácil mostrar que

(r(k+1),r(0)

)= 0. De facto,

(r(k+1),r(0)

)=(r(k),r(0)

)−α(k)

(Ad(k),d(0)),

e, pela hipótese de indução(r(k+1),r(0)

)= 0.

Consideremos agora j = 1, . . . ,k− 1. Atendendo à hipótese de indução e ao facto de d(k) ser A-ortogonal a d( j), tem-se sucessivamente

(r(k+1),r( j))=

(r(k),r( j))−α(k)

(Ad(k),r( j))

=−α(k)(Ad(k),r( j))

=−α(k)(Ad(k),d( j)−β ( j−1)d( j−1))

=−α(k)(Ad(k),d( j))+α(k)β ( j−1)

(Ad(k),d( j−1))

= 0.

No que diz respeito à A-ortogonalidade do conjunto d( j), j = 0, . . . ,k+1, comecemos por notarque se tem (

d(k+1),Ad(k))=(r(k+1),Ad(k))+β (k)

(d(k),Ad(k))

e, atendendo à definição (2.3.34), vem de imediato(d(k+1),Ad(k)

)= 0.

Consideremos agora(d(k+1),Ad( j)

), j = 0, . . . ,k−1. Atendendo a que d(k) e d( j) são A-ortogonais,

tem-se sucessivamente(d(k+1),Ad( j))=

(r(k+1),Ad( j))+β (k)

(d(k),Ad( j))

=(r(k+1),Ad( j)).

Conjugando a última representação com (2.3.35), obtemos

(d(k+1),Ad( j))= 1

α( j)

((r(k+1),r( j+1))−

(r(k+1),r( j))),

e, atendendo à ortogonalidade do conjunto r( j), j = 0, . . . ,k+1, concluímos(d(k+1),Ad( j))= 0.

No Teorema 2.15 provamos que r(k) é ortogonal a d(k−1). Como consequência deste resultado,estabelecemos a equivalência das representações (2.3.32) e (2.3.36). De facto, tem-se sucessivamente

(r(k),d(k))=

(r(k),r(k)−β (k−1)d(k−1))

=(r(k),r(k)

)−β (k−1)

(r(k),d(k−1))

=(r(k),r(k)

),

o que prova a equivalência das duas representações.


Teorema 2.15. O resíduo r(k+1) é ortogonal às direcções d( j), j = 0, . . . ,k.

Demonstração. É manifesto que r(1) é ortogonal a d(0) pois esta direcção coincide com o resíduo inicial.Admitamos agora que r(`) é ortogonal às direcções d( j), j = 0, . . . , `−1, e provemos que r(`+1) é ortogonalàs direcções d( j), j = 0, . . . , `.

Comecemos por provar que r(`+1) é ortogonal a d(`). Tem-se(r(`+1),d(`)

)=(r(`),d(`)

)−α(`)

(Ad(`),d(`)

),

ou ainda, atendendo à definição (2.3.36) para α(`) e à hipótese de indução,(r(`+1),d(`)

)=(r(`)−α(`)Ad(`),d(`)

),

=(r(`),d(`)

)−∥∥r(`)

∥∥22.

Atendendo a que d(`) = r(`)+β (`−1)d(`−1), concluímos finalmente que(r(`+1),d(`)

)= 0.

Provemos agora que (r(`+1),d( j))= 0, j = 0, . . . , `−1.

Tem-se (r(`+1),d( j))=

(r(`),d( j))−α(`)

(Ad(`),d( j)).

Atendendo a que, pela hipótese de indução,(r(`),d( j)

)= 0, e uma vez que d(`) e d( j) são A-ortogonais,(

Ad(`),d( j))= 0, concluímos que r(`+1) e d( j), j = 0, . . . , `−1, são ortogonais.

Notemos que β (m) pode ser calculado a partir apenas dos resíduos r(m+1) e r(m). De facto,

r(m+1)− r(m) =−α(m)Ad(m),

e portantoα(m)

(r(m+1),Ad(m)

)=(r(m),r(m+1))−

∥∥r(m+1)∥∥22.

Uma vez que r(m) e r(m+1) são ortogonais, vem

(r(m+1),Ad(m)

)=−

∥∥r(m+1)∥∥2

2α(m)

.

Substituindo em (2.3.34) obtemos

β (m) =

∥∥r(m+1)∥∥2

2

α(m)∥∥d(m)

∥∥2A

.

Conjugando agora a identidade anterior com (2.3.32), obtemos

β (m) =

∥∥r(m+1)∥∥2

2∥∥r(m)∥∥2

2

. (2.3.37)

O método dos gradientes conjugados resume-se na seguinte sequência de passos: dado x(0),r(0) = b−Ax(0), seja d(0) = r(0). Para m = 0, . . . , calcular


1. α(m) utilizando (2.3.32);

2. x(m+1) utilizando (2.3.31);

3. r(m+1) = b−Ax(m+1);

4. β (m) utilizando (2.3.34) ou (2.3.37);

5. d(m+1) utilizando (2.3.33).

A transcrição do algoritmo anterior para MATLAB pode ser feita na seguinte função:

Algoritmo 2.12 Método dos gradientes conjugados (com controlo do número máximo deiterações.

function [x,res] = gradientesConjugados(A,b,x0,nmax)r = b − A*x0; d = r; x = x0;for i=1:nmax

alpha = (r'*d)/(d'*A*d); x = x + alpha*d; r = b − A*x;beta = −(r'*A*d)/(d'*A*d); d = r + beta*d;

endres = norm(r);

Notemos que a função anterior devolve a solução aproximada do sistema linear pelo métododos gradientes conjugados, assim como a norma-2 do resíduo.

O método dos gradientes conjugados conduz-nos à solução do sistema Ax = b na n-ésima iteração.

Teorema 2.16. Se x(m), m = 0, . . . ,n, definida pelo método dos gradientes conjugados, então x(n)

é solução de Ax = b.

Demonstração. Seja e(0) o erro inicial. Atendendo a que d( j), j = 0, . . . ,n− 1, são A-ortogonais, d( j),j = 0, . . . ,n−1, é uma base the Rn. Sejam δ j, j = 0, . . . ,n−1, tais que

e(0) =n−1

∑j=0

δ jd( j).

Atendendo a que se tem

d(i)Ae(0) =n−1

∑j=0

δ jd(i)Ad( j),

para i = 0, . . . ,n−1, em que

d(i)Ad( j) =(d(i),d( j))

A = 0, para i 6= j,

obtemos

δi =

(d(i),e(0)

)A∥∥d(i)

∥∥2A

.


No que segue obtemos para δi outra representação. Atendendo a que(d(i),d(`)

)A = 0, para `= 0, . . . , i−1,

tem-se

δi =1

∥∥d(i)∥∥2

A

(d(i),A

(e(0)−

i−1

∑`=0

α(`)d(`)

))

=1

∥∥d(i)∥∥2

A

(d(i),Ae(i)

)

pois x(i) = x(0)+∑i−1`=0 α(`)d(`). Notamos que Ae(i) = r(i), e portanto, de (2.3.32), vem

δi = α(i), i = 0, . . . ,n−1.

Logoe(n) = x∗− x(n)

= x∗− x(0)−n−1

∑`=0

α(`)d(`)

= e(0)−n−1

∑`=0

δ`d(`)

= 0,

isto é x(n) = x∗.

O resultado anterior permite concluir que a n-ésima aproximação definida pelo o método dos gradi-entes conjugados (2.3.31), (2.3.32), (2.3.34), (2.3.33), coincide com a solução do sistema Ax = b quepretendemos resolver. No entanto, atendendo aos erros de arredondamento, o método gera uma sucessãode aproximações que converge para a solução de Ax = b mais rapidamente do que a sucessão gerada pelométodo da descida mais rápida.

Teorema 2.17. Se A ∈Mn(R) é simétrica e definida positiva e x(m), m = 0, . . . , é definida pelométodo dos gradientes conjugados (2.3.31), (2.3.32), (2.3.34), (2.3.33). então o erro e(m+1) =x∗− x(m+1) verifica

∥∥e(m+1)∥∥A ≤

2(√

cond2(A)+1√cond2(A)−1

)m+1

+

(√cond2(A)−1√cond2(A)+1

)m+1

∥∥e(0)∥∥

A, m = 0, . . . . (2.3.38)

Na Figura 2.5 ilustramos o comportamento dos factores

mdr(cond2(A),m) =

(cond2(A)−1cond2(A)+1

)m+1

e

mgc(cond2(A),m)) =2

(√cond2(A)+1√cond2(A)−1

)m+1

+

(√cond2(A)−1√cond2(A)+1

)m+1 ,

para m = 10,20. Observamos que, para estes valores, mdr(cond2(A),m)> mgc(cond2(A),m).


0 20 40 60 80 100

0

0.2

0.4

0.6

cond2(A)

mdr(cond2(A),10)mgc(cond2(A),10)

0 20 40 60 80 100

0

0.1

0.2

0.3

0.4

cond2(A)

mdr(cond2(A),20)mgc(cond2(A),20)

Figura 2.5: Comportamento dos coeficientes mdr(cond2(A),m) e mgc(cond2(A),m) para m = 10 (es-querda) e m = 20 (direita).

2.3.3 Algumas considerações finais

A questão natural que se coloca quando implementamos um dos métodos iterativos estudados, dizrespeito à paragem, isto é: como podemos definir a paragem do processo iterativo? Há vários critérios deparagem para os processos iterativos.

1. O critério de paragem pode ser estabelecido considerando um número máximo de iterações previa-mente definido. Este foi um dos critérios usados em todos os algoritmos deste capítulo.

2. Um outro critério de paragem está associado ao cálculo do resíduo. Poderemos estabelecer que oprocesso iterativo pára quando o resíduo de tal iteração for inferior a uma tolerância especificada apriori, isto é, o processo iterativo termina na iteração x(m) tal que

∥∥r(m)∥∥=

∥∥b−Ax(m)∥∥< Tol,

relativamente a uma determinada norma ‖·‖.

3. Poderemos pretender determinar uma approximação para a solução de Ax = b com erro inferior auma tolerância especificada. Neste caso teremos de determinar o número de iterações m tal que

∥∥e(m)∥∥< Tol.

É manifesto que tal só pode ser feito quando conhecemos majorantes para o erro∥∥e(m)

∥∥ que dependede m mas é independente do erro inicial (que é desconhecido). Por exemplo, para os métodosestacionários em que a matriz de iteração B satisfaz ‖B‖< 1, foi estabelecido

∥∥e(m)∥∥≤ ‖B‖m

1−‖B‖∥∥x(1)− x(0)

∥∥.

Assim, considerando a desigualdade

‖B‖m

1−‖B‖∥∥x(1)− x(0)

∥∥< Tol,

obtemos o número de iterações que devemos considerar. Este tipo de critério foi usado no Algoritmo2.11.


4. Um outro critério de paragem está associado ao comportamento dos termos consecutivos. Assimpodemos estabelecer parar o processo iterativo quando

∥∥x(m+1)− x(m)∥∥< Tol.

Se considerarmos um sistema em que os vectores do processo iterativo têm norma muito pequena,então o processo iterativo pára podendo eventualmente ocorrer

∥∥e(m)∥∥ 0.

5. Um outro critério de paragem do tipo do anterior consiste em parar o processo iterativo∥∥x(m+1)− x(m)

∥∥∥∥x(m)

∥∥ < Tol.

Este critério não vai permitir a ocorrência da situação anormal que pode ocorrer quando usamos ocritério da alí nea anterior.


1. Determine, utilizando o método de eliminação de Gauss, a solução dos seguintes sistemas

(a)

x + y + z = 12x − y + z = 1x − y + z = 1

(b)

2x − y + z + t = 14x + y − z = −1−x + y = 0

2. Determine, para as matrizes associadas aos sistemas da alínea anterior, a factorização LR.

3. Indique como pode obter a factorização QR de uma matriz A utilizando o método de Gram-Schmidt.

4. Utilizando transformações de Householder, determine a factorização QR da matriz

A =

1 6 −12 1 −12 2 6

.

5. Considere A ∈Mn(R) e QR a sua factorização em que Q é ortogonal e R é triangular superior.Mostre que

(a) ‖Q‖1 ≤√

n,

(b) cond1(Q)≤ n,

(c) n−1cond1(R)≤ cond1(A)≤ ncond1(R),

(d) cond2(A) = cond2(R).

6. Recorrendo à factorização QR da matriz A ∈Mn(R), prove que

|det(A)| ≤n

∏j=1‖a j‖2.


7. Determine a factorização em valores singulares da matriz

A =

1 10 11 0

.

Calcule A+.

8. Considere o sistema Ax = b, em que

A =

1 0 2−1 1 −21 1 2

, b =

a−aa

, a ∈ R.

Determine, utilizando a decomposição em valores singulares de A, a solução de Ax = b.

9. Seja A ∈Mn(R) não singular. Mostre que

‖A‖2 = s1,∥∥A−1∥∥

2 =1sn, cond2(A) =

s1

sn,

em que si =√

λi, s1 ≥ s2 ≥ ·· · ≥ sn são os valores singulares de A.

10. Mostre que ‖A‖F =(∑

ni=1 s2

i)1/2

, em que si, i = 1, . . . ,n, são os valores singulares de A.

11. Seja A ∈Mm×n(R) com colunas linearmente independentes. Demonstre as propriedades seguintes:

(a) Se m = n e A é não singular, então A−1 = A+;

(b) (At)+ = (A+)t ;

(c) Se A é simétrica, então a sua pseudo-inversa também é simétrica;

(d) AA+A = A;

(e) A+AA+ = A+;

(f) (AA+)t = (AA+);

(g) (A+A)t = A+A.

12. Seja x∗ a solução dos mínimos quadrados do sistema impossível Ax = b. Estabeleça as equaçõesnormais AtAx∗ = Atb.

13. Diga, utilizando a factorização QR a decomposição em valores singulares, como determinar asolução dos mínimos quadrados de Ax = b.


A =

−8 1 11 −5 11 1 −4

, b =

116−7

.

Seja x(0) o vector nulo. Determine a segunda iteração para a solução do sistema dado

(a) com o método de Jacobi,

(b) com o método de Gauss-Seidel.

Conclua que os dois métodos geram sucessões que convergem para a solução de Ax = b.


15. Considere os sistema

x1 + 5x2 + 2x4 = 8x1 + x2 + 3x3 = 54x1 − x2 + x3 − x4 = 32x1 − x2 + x3 + 8x4 = 10.

(a) A partir da matriz do sistema pode concluir que o método de Jacobi gera uma sucessãoconvergente qualquer que seja a aproximação inicial?

(b) Reescreva o sistema dado de modo a poder aplicar o referido método e determine a segundaiteração partindo da solução nula.


A =

0 6 3α 3 21 1 8

, b =

444

.

Seja x(0) o vector nulo.

(a) Poderá aplicar o método de Gauss-Seidel ao sistema anterior? Em caso negativo reescreva-ode modo conveniente.

(b) Determine um intervalo de variação de α para os quais o método de Gauss-Seidel gera umasucessão convergente.

17. O sistema linear

2x1 − x2 + x3 = −12x1 + 2x2 + 2x3 = 4−x1 − x2 + 2x3 = −5

tem solução (1,2,−1).

(a) Mostre que a matriz de iteração do método de Jacobi tem raio espectral√

52 . O que pode

concluir da aplicabilidade deste método?

(b) Mostre que a matriz de iteração do método de Gauss-Seidel tem raio espectral 12 . O que pode

concluir da aplicabilidade deste método?

18. O sistema linear

x1 + 2x2 − 2x3 = 7x1 + x2 + x3 = 2

2x1 + 2x2 + x3 = 5

tem solução (1,2,−1).

(a) Mostre que a matriz de iteração do método de Jacobi tem raio espectral nulo. O que podeconcluir da aplicabilidade deste método?

(b) Mostre que a matriz de iteração do método de Gauss-Seidel tem raio espectral 2. O que podeconcluir da aplicabilidade deste método?

19. Considere o sistema Ax = b com

5 1 −10 3 12 2 4

, b =

549

.


(a) Mostre que os métodos de Jacobi e de Gauss-Seidel geram sucessões convergentes para asolução de Ax = b.

(b) Sejam BJ e BGS as matrizes de iteração dos métodos de Jacobi e de Gauss-Seidel, respectiva-mente. Utilizando os dois métodos anteriores, determine o número mínimo de iterações quedeve considerar para que

∥∥e(m)J

∥∥1 < 10−3,

∥∥e(m)GS

∥∥1 < 10−3,

em que e(m)J e e(m)

GS representam os erros das aproximações definidas pelos dois métodosconsiderados.

20. Mostre que o método de Gauss-Seidel é convergente quando aplicado ao sistema Ax = b com

A =

2 −1−1 2 −1

. . . . . . . . .−1 2 −1

−1 2

.

21. Considere o método estacionário x(m+1) = Bx(m)+ c, m = 0, . . . , e seja h(m+1) = (B− I)x(m)+ c.Neste caso tem-se x(m+1) = x(m)+h(m+1), m = 0, . . . . Mostre que as correcções h(m), m = 0, . . . ,satisfazem

h(m+1) = Bh(m).

22. Considere o método estacionário x(m+1) = Bx(m)+c, m = 0, . . . , convergente e seja γ(m) = x(m+1)−x(m). Mostre que relativamente a alguma norma ‖.‖ se tem

∥∥e(m)∥∥≤ 1

1−‖B‖∥∥γ(m)

∥∥.

23. Considere o seguinte método iterativo

x(m+1) =

[ω −1

4+

ω2

1 1

]x(m)+ c(ω), m = 0, . . . .

Determine a variação de ω > 0 para a qual o método anterior é convergente.

24. Mostre que o método de Gauss-Seidel com relaxação é convergente quando aplicado ao sistemaAx = b com

A =

4 −1 −1−1 4 −1−1 −1 4

para ω ∈ (0,2) e para qualquer aproximação inicial.

25. Considere o sistema [2 11 2

][x1x2

]=

[01

]

que tem a solução (−1/3,2/3).

(a) Determine x(3) utilizando o método da descida mais rápida com x(0) = (0,0).


(b) Seja e(m) o erro da aproximação x(m) e seja A da matriz do sistema anterior. Defina K(λ1,λ2)tal que seja válida a estimativa

∥∥e(m)∥∥2

A ≤(

λ1−λ2

λ1 +λ2

)2(m−1)

K(λ1,λ2)∥∥x(1)− x(0)

∥∥2A,

em que λ1 > λ2 são os valores próprios de A.

(c) Determine o número mínimo de iterações tal que∥∥e(m)

∥∥2A ≤ 10−3.

26. Mostre que a sucessão dos resíduos r(m),m = 0, . . . , determinados pelo método da descida maisrápida para o sistema Ax = b, verifica

r(m+1) = (I−α(m)A)r(m), m = 0, . . . .

27. Seja x(m),m = 0, . . . , a sucessão de aproximações para a solução x∗ do sistema Ax = b definida pelométodo da descida mais rápida. Mostre que se o erro e( j), para algum j ∈ N0, é vector próprio deA, então x( j+1) = x∗.

28. Considere o sistema Ax = b, em que A é simétrica e positiva definida. Seja x(m), m = 0, . . . , asucessão definida pelo método da descida mais rápida. Mostre que

(a) e(m+1) é ortogonal a Ar(m);

(b) r(m+1) é ortogonal a r(m).

29. Seja A uma matriz real simétrica e positiva definida. Mostre que se d( j), j = 0, . . . ,m, são A-ortogonais, então são linearmente independentes.

30. Prove que os vectores próprios de uma matriz A real simétrica e positiva definida são A-ortogonais.

31. Considere o sistema Ax = b com

A =

2 −1 0−1 2 −10 −1 2

, b =

101

.

que tem por solução (1,1,1). Seja x(0) = (0,0,0).

(a) Prove que a matriz do sistema anterior é simétrica e positiva definida.

(b) Seja x(m) definida pelo método da descida mais rápida. Determine uma estimativa para∥∥e(m)

∥∥A

(c) Determine o número mínimo de iterações que deve considerar para que∥∥e(m)

∥∥A < 10−2

(d) Utilizando o método dos gradientes conjugados, determine x(2).

(e) Calculde r(2). O que pode concluir?

32. (a) Considere o sistema Ax = b, em que A é não singular, e o método iterativo

x(m+1) = x(m)+α(m)d(m), m = 0, . . . ,

em que a direção de descida d(m) é o resíduo r(m). Seja Φ(x) = ‖b−Ax‖22. Mostre que o valor

mínimo de Φ(x(m)+α(m)d(m)) ocorre para α(m) definido por

α(m) =

(r(m),Ar(m)

)∥∥Ar(m)

∥∥22

.


(b) Considere o sistema [2 −1−1 2

][x1x2

]=

[31

],

que tem por solução (3,1).Determine x(3), com x(0) = (0,0), utilizando o método apresentado na alínea anterior, ométodo da descida mais inclinada e o método dos gradientes conjugados.

(c) Comente os resultados obtidos.


CAPÍTULO 3INTERPOLAÇÃO POLINOMIAL

3.1 Introdução

Dada uma função f definida num certo intervalo [a,b], a questão central que consideramos nestecapítulo é a construção de uma função polinomial que, num sentido a especificar, a aproxima. A existênciadesta aproximação para funções contínuas é garantida pelo teorema clássico de Weierstrass:

Teorema 3.1. Se f : [a,b]−→ R é contínua, então, para ε > 0, existe um polinómio P definidoem [a,b] tal que

‖ f −P‖∞ < ε, (3.1.1)

em que ‖ f −P‖∞ = maxx∈[a,b]

| f (x)−P(x)|

O Teorema de Weierstrass garante a existência de uma função polinomial que aproxima a funçãodada uniformemente no intervalo. A utilização de funções polinomiais na substituição de certas funçõesjustifica-se pelas boas propriedades de diferenciabilidade e integrabilidade que estas apresentam. De facto,suponhamos que pretendemos calcular ∫ 1

0e−x2

dx.

A função integranda não é primitivável, e portanto não é possível efectuar o cálculo anterior sem substituira função f (x) = e−x2

,x ∈ [0,1]. Uma substituição natural consiste em tomar

e−x2=

∞

∑j=0

(−1) jx2 j

j!, x ∈ R,

65


e obter ∫ 1

0e−x2

dx =∞

∑j=0

(−1) j

(2 j+1) j!,

obtendo-se o integral pretendido como a soma de uma série. No entanto, poderíamos pretender determinarum valor aproximado para o mesmo integral com um erro inferior a uma tolerância ε , definida a priori.Assim, a substituição da função integranda por uma função polinomial P cuja existência é estabelecida noTeorema de Weierstrass permite obter

∫ 1

0e−x2

dx'∫ 1

0P(x)dx,

em que ∣∣∣∣∫ 1

0e−x2

dx−∫ 1

0P(x)dx

∣∣∣∣< ε.

No caso de termos funções bastante regulares, por exemplo, se as sucessivas derivadas de f existem esão limitadas, então a construção de um polinómio nas condições do Teorema de Weierstrass pode serfeita considerando o seguinte procedimento:

1. fixar em [a,b] uma partição uniforme de espaçamento h = b−an+1 em que a = x0, b = xN , e obter o

conjunto de pontos(xi, f (xi)), i = 0, . . . ,n.

2. determinar em seguida o polinómio Pn, de grau inferior ou igual a n, tal que

Pn(xi) = f (xi), i = 0, . . . ,n.

Neste caso, o número de pontos fixados (ou seja, o espaçamento) é induzido pelo majorante dasderivadas da função e deverá ser tal que

‖ f −Pn‖∞ < ε.

O último problema que especificamos é usualmente conhecido por problema de interpolação e temsignificado num contexto mais geral do que o exposto e que especificamos: dado (xi,yi), i = 0, . . . ,n,um conjunto de valores em que xi, i = 0, . . . ,n, estão por ordem crescente, pretendemos determinar umafunção g numa determinada classe de funções P e definida num domínio adequando tal que

g(xi) = yi, i = 0, . . . ,n.

No caso particular de yi, i = 0, . . . , serem a imagem de xi, i = 0, . . . ,n, por uma função desconhecidaf , então, o problema de interpolação permite determinar g tal que

g(xi) = f (xi), i = 0, . . . ,n,

e, deste modo, determinar, para x ∈ [x0,xn], x 6= xi,

f (xi)' g(x).

A questão que permanece diz respeito à estimação do erro

| f (x)−g(x)|.

Na resolução do problema anterior, colocam-se três questões:

1. Como escolher a classe P?


2. Como determinar g ∈P?

3. Como quantificar o erro f (x)−g(x)?

Relativamente à primeira questão, iremos considerar a classe dos polinómios de grau menor ou igual an, Pn, e restringimo-nos apenas à interpolação polinomial. No que diz respeito à segunda questão, iremosdeterminar P ∈Pn tal que

P(xi) = f (xi), i = 0, . . . ,n. (3.1.2)

Este polinómio é chamado polinómio interpolador de f nos pontos xi, i = 0, . . . ,n, e a estes pontoschamamos nodos de interpolação. O estudo do erro é feito aquando da construção do polinómiointerpolador.

3.2 Polinómio interpolador de Lagrange

Consideremos uma função f definida num intervalo [a,b] e, neste intervalo, a partição

a = x0 < x1 < · · ·< xn−1 < xn = b. (3.2.1)

Pretendemos construir um polinómio P que verifique P(xi) = f (xi), i = 0, . . . ,n.Consideremos as funções

ì(x) =n

∏j=0, j 6=i

x− x j

xi− x j, i = 0, . . . ,n. (3.2.2)

As funções anteriores são chamadas polinómios de Lagrange e têm um papel importante na construção dopolinómio interpolador. Notemos que os polinómios de Lagrange têm grau n e verificam

ì(x j) =

1, i = j,0, i 6= j,

e portanto é fácil concluir que ì, i = 0, . . . ,n, constituem uma base do espaço dos polinómios de graumenor ou igual a n, Pn. A existência e unicidade do polinómio interpolador é estabelecida no teoremaseguinte.

Teorema 3.2. Dado o conjunto de pontos (xi, f (xi)), i = 0, . . . ,n, de abcissas distintas, o polinó-mio

Pn(x) =n

∑i=0

f (xi)ì(x)

é o único polinómio de grau menor ou igual a n que é interpolador de f nos pontos xi, i = 0, . . . ,n.

Demonstração. Atendendo às propriedades dos polinómios de Lagrange, concluímos que Pn tem graumenor ou igual a n e verifica Pn(xi) = f (xi), i = 0, . . . ,n.

Provemos seguidamente a unicidade. Seja Q outro polinómio, de grau menor ou igual a n, que verifica

Q(xi) = f (xi), i = 0, . . . ,n.

Denotemos por R a função polinomial

R(x) = Pn(x)−Q(x).

Esta função tem grau menor ou igual a n e tem por zeros xi, i = 0, . . . ,n. Logo R(x) = 0, ∀x ∈ [a,b], poisum polinómio de grau menor ou igual a n tem quando muito n zeros.


Exemplo 3.1. Consideremos os pontos(0, 1

2

),(1

2 ,1),(1, 1

4

)e determinemos o polinómio que passa nos

pontos dados.Temos

`0(x) = 2(

x− 12

)(x−1), `1(x) =−4x(x−1), `2(x) = 2x

(x− 1

2

),

e portanto

P2(x) =(

x− 12

)(x−1)−4x(x−1)+

12

x(

x− 12

).

Observação 3.1 No teorema anterior é estabelecida a unicidade do polinómio interpolador. No entanto,atendendo a que na sua construção utilizámos os polinómios de Lagrange, o polinómio interpolador Pn é,por abuso de linguagem, chamado polinómio interpolador de Lagrange. Como veremos posteriormente,outras designações serão consideradas para o mesmo polinómio e que estão relacionadas com o processousado na sua construção.

A implementação do cálculo do polinómio interpolador de Lagrange, pode ser feita re-correndo a duas funções: uma que calcule os polinómios de Lagrange (3.2.2) e outra queimplemente a fórmula do Teorema 3.2:

Algoritmo 3.1 Polinómios de Lagrange

function [y] = polinomioLagrange(x,i,xPts)

y=ones(size(x));for k=1:length(xPts)

if ( k 6= i )y = y.*(x−xPts(k))/(xPts(i)−xPts(k));

endendreturn

Algoritmo 3.2 Polinómio interpolador de Lagrange

function [y] = polinomioInterpoladorLagrange(x,xPts,yPts)

y = zeros(size(x));for k=1:length(xPts)y = y + yPts(k)*polinomioLagrange(x,k,xPts);

endreturn

Observamos que a variável x, que contem as coordenadas dos pontos onde pretendemos cal-cular o polinómio interpolador de Lagrange, pode ser um vector ou uma matriz de coordenadas.

Vejamos seguidamente como estimar o erro que se comete ao aproximar f por Pn.


Teorema 3.3. Seja f : [a,b] −→ R tal que f ∈ Cn([a,b]) e existe f (n+1) em (a,b). Em [a,b]consideremos a partição (3.2.1), e seja Pn o polinómio interpolador de Lagrange de f nos n+1pontos xi, i = 0, . . . ,n. Se x ∈ [a,b], então existe η(x) ∈ (a,b) tal que

f (x)−Pn(x) =f (n+1)(η(x))(n+1)!

n

∏j=0

(x− x j). (3.2.3)

Demonstração. Seja x ∈ [a,b]. Se x = xi então teorema está demonstrado. Se x 6= xi, i = 0, . . . ,n,consideremos a função

g(t) = f (t)−Pn(t)− [ f (x)−Pn(x)]n

∏j=0

t− x j

x− x j, t ∈ [a,b].

A função g verificag(x) = 0, g(xi) = 0, i = 0, . . . ,n,

e portanto, g tem pelo menos n+2 zeros. Assim, g(n+1) tem pelo menos um zero em (a,b), isto é, existeη(x) ∈ (a,b) tal que g(n+1)(η(x)) = 0. Temos então

0 = f (n+1)(η(x))− [ f (x)−Pn(x)]1

∏nj=0(x− x j)

( dn+1

dtn+1

n

∏j=0

(t− x j))

t=η(x).

Da igualdade anterior vem

0= f (n+1)(η(x))− [ f (x)−Pn(x)](n+1)!

∏nj=0(x− x j)

,

o que prova o pretendido.

Atendendo ao resultado anterior podemos determinar uma estimativa para o erro cometido ao substituira função f pelo seu polinómio interpolador Pn num ponto do intervalo [a,b]. Se existe M > 0 tal que

∣∣ f (n+1)(x)∣∣≤M, ∀x ∈ [a,b],

então

| f (x)−Pn(x)| ≤M

(n+1)!

n

∏j=0|x− x j|. (3.2.4)

Exemplo 3.2. Consideremos a função f (x) = sin(π

2 x), x ∈ [0,1] e os valores

xi 0 14

12

34 1

f (xi) 0 sin(π

8

) √2

2 sin(3π

8

)1

Com base na tabela anterior, determinemos os polinómios P2 e P4. Estes polinómios admitem a represen-tação

P2(x) =

√2

2`1(x)+ `2(x),


0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

fP2pontos de interpolação

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

fP4pontos de interpolação

Figura 3.1: Gráficos de f (x) = sin(π

2 x)

e P2, à esquerda, e P4, à direita, com os respectivos pontos deinterpolação usados.

em que

`1(x) =x(x−1)12

(12 −1

) , `2(x) =x(x− 1

2

)

1(1− 1

2

) ,

e

P4(x) = sin(π

8

)`1(x)+

√2

2`2(x)+ sin

(3π8

)`3(x)+ `4(x),

em que

`1(x) =x(x− 1

2

)(x− 3

4

)(x−1)

34

(34 − 1

2

)(34 −1

) , `2(x) =x(x− 1

4

)(x− 3

4

)(x−1)

12

(12 − 1

4

)(12 − 3

4

)(12 −1

) ,

`3(x) =x(x− 1

4

)(x− 1

2

)(x−1)

34(

34 − 1

4)(3

4 − 12

)(34 −1

) , `4(x) =x(x− 1

4

)(x− 1

2

)(x− 3

4

)

1(1− 1

4

)(1− 1

2

)(1− 3

4

) .

Na Figura 3.1 ilustramos o comportamento de f e dos seus polinómios interpoladores P2 e P4. Observamosque o aumento do grau induz um aumento de precisão.

Exemplo 3.3. Consideremos a função f (x) = ln(x), x ∈ [1,2] e os valores

xi 1 1.5 2

f (xi) 0 ln(1.5) ln(2)

Determinemos o polinómio interpolador de Lagrange de f nos pontos dados na tabela. Este polinómio,para os dados anteriores, admite a representação

P2(x) =−4ln(1.5)(x−1)(x−2)+2ln(2)(x−1)(x−1.5).

Determinemos uma aproximação para ln(9

8

). Utilizando o polinómio anterior obtemos

ln(9

8

)' P2

(98

)= 0.11241.

A estimativa (3.2.4) permite estabelecer∣∣∣∣ln(9

8

)−P2

(98

)∣∣∣∣≤23!

2183 ≤ 0.01368.


Podemos obter uma estimativa para o erro independente de x ∈ [a,b]. De facto, atendendo a que

|x− x j| ≤ b−a, j = 0, . . . ,n, x ∈ [a,b],

vem, de (3.2.4),

‖ f −Pn‖∞ ≤M

(n+1)!(b−a)n+1,x ∈ [a,b]. (3.2.5)

A estimativa do erro calculada usando esta desigualdade é, de um modo geral, muito grosseira. Estadesigualdade pode ser utilizada para determinar o polinómio interpolador que tenha, em qualquer pontodo intervalo da partição, um erro inferior a uma quantidade estabelecida ε . De facto, basta determinar ntal que

M(n+1)!

(b−a)n+1 < ε (3.2.6)

e, em seguida, considerar para o valor determinado n uma partição uniforme de espaçamento h=(b−a)/n.

Exemplo 3.4. Consideremos novamente f (x) = ln(x), x ∈ [1,2] e determinemos qual deve ser o menorgrau do polinómio interpolador para que o erro cometido ao aproximar f por esse polinómio seja inferiora 10−2.

Atendendo a que

f (n)(x) = (−1)n+1 (n−1)!xn

para n ∈ N, de (3.2.6) vem

| ln(x)−Pn(x)| ≤1

(n+1), x ∈ [1,2].

Assim, o valor de n pretendido é 100.

Vejamos seguidamente como determinar uma estimativa para o erro em função de

h = maxi=1,...,n

(xi− xi−1).

Provemos que

maxx∈[a,b]

∣∣∣∣n

∏j=0

(x− x j)

∣∣∣∣≤hn+1n!

4. (3.2.7)

Para n = 1 temosg(x) = |(x− x0)(x− x1)|=−(x− x0)(x− x1), x ∈ [a,b],

e esta função tem um máximo em x∗ = x0+x12 cujo valor é

g(

x0 + x1

2

)=

∣∣∣∣(

x0 + x1

2− x0

)(x0 + x1

2− x1

)∣∣∣∣≤h2

4.

Suponhamos a veracidade da desigualdade (3.2.7) para n e provemos a sua veracidade para n+1, isto é,

maxx∈[a,b]

∣∣∣∣n+1

∏j=0

(x− x j)

∣∣∣∣≤hn+2(n+1)!

4, (3.2.8)

em que x0 = a e xn+1 = b. Notemos que se x ∈ [a,b], então x ∈ [a,xn] ou x ∈ [xn,b]. Consideremosx ∈ [a,xn]. Atendendo a ∣∣∣∣

n+1

∏j=0

(x− x j)

∣∣∣∣=∣∣∣∣

n

∏j=0

(x− x j)

∣∣∣∣|x−b|, (3.2.9)


∣∣∣∣n

∏j=0

(x− x j)

∣∣∣∣≤hn+1n!

4

e|x−b| ≤ (n+1)h, x ∈ [a,b],

concluímos o pretendido.Quando x ∈ [xn,b], em vez de (3.2.9), devemos considerar

∣∣∣∣∣n+1

∏j=0

(x− x j)

∣∣∣∣∣= |(x−a)||n+1

∏j=1

(x− x j) |,

e a demostração de (3.2.8) é análoga à anterior.Provámos o seguinte corolário:

Corolário 3.1. Seja f : [a,b]−→ R tal que f ∈Cn([a,b]), existe f (n+1) em (a,b) e existe M queverifica ∣∣ f (n+1)(x)

∣∣≤M, x ∈ [a,b].

Em [a,b] consideremos a partição (3.2.1), e seja Pn o polinómio interpolador de Lagrange de fnos n+1 pontos xi, i = 0, . . . ,n. Então

‖ f −Pn‖∞ ≤M

4(n+1)hn+1, x ∈ [a,b], (3.2.10)

em que h = maxi=1,...,n

(xi− xi−1).

Particularizemos este corolário ao caso em que a partição é uniforme. Estabelecemos facilmente oseguinte resultado.

Corolário 3.2. Seja f : [a,b]−→ R tal que f ∈Cn([a,b]), existe f (n+1) em (a,b) e existe M queverifica ∣∣ f (n+1)(x)

∣∣≤M, x ∈ [a,b].

Em [a,b] consideremos a partição uniforme (3.2.1), em que xi− xi−1 =b−a

n , i = 1, . . . ,n, e sejaPn o polinómio interpolador de Lagrange de f nos n+1 pontos xi, i = 0, . . . ,n. Então

‖ f −Pn‖ ≤M

4(n+1)nn+1 (b−a)n+1,x ∈ [a,b]. (3.2.11)

O corolário anterior permite calcular

1. o erro cometido ao substituir f pelo seu polinómio interpolador em qualquer ponto de [a,b],

2. o polinómio interpolador de Lagrange que aproxima f em qualquer ponto do intervalo [a,b] comum erro inferior a ε.

No exemplo anterior considerámos a desigualdade (3.2.5) e determinámos n = 100. Usando (3.2.11)para resolver o problema apresentado naquele exemplo obtemos n = 3.


3.3 Polinómio interpolador de Newton

A expressão do polinómio interpolador foi obtida considerando os polinómio de Lagrange ì(x),i = 0, . . . ,n, definidos sobre a partição xi, i = 0, . . . ,n. A determinação do polinómio interpolador deuma função usando as funções ì exige um grande esforço computacional e não permite obter o polinómiointerpolador de grau n a partir do conhecimento do polinómio interpolador de grau n− 1. Estes doisfactores levam-nos a considerar outro modo de obter o polinómio interpolador.

Consideremos uma função f definida num intervalo [a,b] onde consideramos a partição (3.2.1).Seja Pn o polinómio interpolador de Lagrange para f relativamente a esta partição, ou seja, nos pontosxi, i = 0, . . . ,n, considerados em [a,b]. Este polinómio tem grau menor ou igual a n e, portanto, pertenceao espaço Pn - espaço dos polinómmios de grau menor ou igual a n. Neste espaço consideremos asfunções

φ0(x) = 1, φ j(x) =j−1

∏i=0

(x− xi), j = 1, . . . ,n.

O conjunto φ j(x), j = 0, . . . ,n constitui uma base do espaço Pn. Logo existem constantes a j,j = 0, . . . ,n, tais que o polinómio interpolador de Lagrange é dado por

Pn(x) =n

∑j=0

a jφ j(x).

Determinemos os coeficientes da combinação linear anterior. Atendendo a que Pn(x0) = f (x0) vema0 = f (x0). Considerando esta constante em Pn(x) obtemos

Pn(x) = f (x0)+n

∑j=1

a jφ j(x).

Atendendo a que Pn(x1) = f (x1), vem

a1 =f (x1)− f (x0)

x1− x0

que denotamos por f [x0,x1] e chamamos diferença dividida de primeira ordem.Substituindo em Pn(x) e atendendo a que Pn(x2) = f (x2), deduzimos

a2 =f (x2)− f (x0)− f [x0,x1](x2− x0)

(x2− x0)(x2− x1)=

f [x2,x1]− f [x1,x0]

x2− x0

que denotamos por f [x0,x1,x2].Podemos deste modo obter um processo recursivo para a determinação dos coeficientes do polinómio

Pn.Seja então

f [xi,xi+1] =f (xi+1)− f (xi)

xi+1− xi.

A diferença dividida de ordem k ≥ 2 é definida em função das diferenças divididas de ordem k−1 por

f [xi,xi+1, . . . ,xi+k+1] =f [xi+1, . . . ,xi+k+1]− f [xi, . . . ,xi+k]

xi+k+1− xi, k ∈ N.

Atendendo às considerações anteriores vem

a j = f [x0, . . . ,x j], j = 1, . . . ,n.


Substituindo na expressão de Pn(x), obtemos

Pn(x) = f (x0)+n

∑j=1

f [x0, . . . ,x j]j−1

∏i=0

(x− xi). (3.3.1)

O polinómio interpolador determinado é designado polinómio interpolador de Newton.Notemos que, uma vez determinado Pn e pretendermos obter Pn+1, basta fazer

Pn+1(x) = Pn(x)+ f [x0, . . . ,xn+1]n

∏i=0

(x− xi).

A determinação dos coeficientes do polinómio interpolador de Newton pode ser feita organizando oscálculos numa tabela - tabela de diferenças divididas. A tabela seguinte apresenta os referidos coeficientespara o caso n = 3.

xi f (xi) f [·, ·] f [·, ·, ·] f [·, ·, ·, ·]x0 f (x0)

f (x1)− f (x0)x1−x0

x1 f (x1)f [x1,x2]− f [x0,x1]

x2−x0f (x2)− f (x1)

x2−x1f [x0, ., .,x3]

x2 f (x2)f [x2,x3]− f [x1,x2]

x3−x1f (x3)− f (x2)

x3−x2

x3 f (x3)

Exemplo 3.5. Consideremos uma função da qual se conhecem os valores dados no seguinte quadro

xi 1 54

32 2 5

2

f (xi) −1 0 2 3 1

Determinemos, utilizando diferenças divididas, o polinómio interpolador de f nos pontos xi, i = 0,1,2,3.

Atendendo a que com os dados anteriores obtemos a tabela de diferenças divididas

xi f (xi) f [·, ·] f [·, ·, ·] f [·, ·, ·, ·] f [·, ·, ·, ·, ·]1 −1

454 0 8

8 −1632 2 −8 176

152 8

52 3 −6

−452 1

o polinómio interpolador pretendido tem a seguinte expressão

P4(x) =−1+4(x−1)+8(x−1)(

x− 54

)−16(x−1)

(x− 5

4

)(x− 3

2

)+

+17615

(x−1)(

x− 54

)(x− 3

2

)(x−2).


Novamente, à semelhança do que foi feito para o polinómio interpolador de Lagrange,podemos dividir o cálculo do polinómio interpolador de Newton em duas partes: o cálculo asdiferenças divididas intervenientes na fórmula e o cálculo do polinómio segundo a fórmula(3.3.1). As diferenças divididas podem ser calculadas usando a seguinte função

Algoritmo 3.3 Cálculo das diferenças divididas f [x0, . . . ,xn]

function [diferencas] = diferencasDivididas(xPts,yPts)

A=zeros(length(xPts),length(xPts));A(:,1) = yPts;

for col=2:length(xPts)for linha=1:length(xPts)−(col−1)A(linha,col) = (A(linha+1,col−1) − ...

A(linha,col−1))/(xPts(linha+col−1) − xPts(linha));end

end

diferencas = A(1,:);return

Uma vez determinados os coeficientes a j, j = 0 . . . ,n necessários para calcular o polinómiointerpolador de Newton, a determinação do valor deste polinómio num vector (ou matriz) decoordenadas x é feito com o seguinte algoritmo:

Algoritmo 3.4 Cálculo do polinómio interpolador de Newton

function [y] = polinomioInterpoladorNewton(x,xPts,yPts)

diferencas = diferencasDivididas(xPts,yPts);y = diferencas(1)*ones(size(x));aux=ones(size(x));for k=2:length(xPts)

aux = aux.*(x−xPts(k−1));y = y + diferencas(k)*aux;

endreturn

Observamos que a função diferencasDivididas poderia ser modificada para guardartodas as diferenças divididas nela calculadas. Caso pretendessemos acrescentar um ponto deinterpolação, a nova função poderia devolver os coeficientes para o polinómio interpoladorde Newton que contabiliza o novo ponto de interpolação, sem recalcular todas as entradas damatriz A.

A diferença dividida de ordem n e a derivada de ordem n+1 estão intrinsecamente relacionadas comoestabelecemos no resultado seguinte.


Teorema 3.4. Em [a,b] consideremos a partição (3.2.1). Se f : [a,b]−→ R tal que f ∈Cn([a,b])e existe f (n+1) em (a,b), então existe η ∈ (a,b) tal que

f [x0, . . . ,xn]n! = f (n)(η). (3.3.2)

Demonstração. Seja Pn o polinómio interpolador de f em xi, i = 0, . . . ,n. Atendendo a que e(x) =f (x)−Pn(x) tem pelo menos n+1 zeros em [a,b], então e(n) tem pelo menos um zero em (a,b), isto é,existe η ∈ (a,b) tal que se tem sucessivamente

0 = e(n)(η)

= f (n)(η)−P(n)n (η)

= f (n)(η)− f [x0, . . . ,xn]n!

que prova o pretendido.

Seja Pn o polinómio interpolador de f em xi, i = 0, . . . ,n, em que f (xi) < f (xi+1), i = 0, . . . ,n−1.Seja y∗ ∈ ( f (x0), f (xn)). Pretendemos determinar x∗ ∈ [a,b] tal que y∗ = f (x∗). O problema anterior,usualmente designado por problema inverso, pode ser resolvido, pelo menos de modo aproximado,utilizando a teoria de interpolação. De facto, pretendemos x∗ tal que x∗ = f−1(y∗). Seja Gn o polinómiointerpolador para os seguintes dados

(yi, f−1(yi)

), i = 0, . . . ,n.

Então Gn(y∗) constitui uma aproximação para f−1(y∗), ou seja para x∗.

Exemplo 3.6. Consideremos a seguinte tabela

xi 0 14

12

34 1

f (xi) 0 sin(π

8

) √2

2 sin(3π

8

)1

que diz respeito a f (x) = sin(π

2 x), x ∈ [0,1]. Pretendemos determinar x∗ ∈ (0,1) tal que f (x∗) = 0.95.

A tabela seguinte

yi 0 sin(π

8

) √2

2 sin(3π

8

)1

xi 0 14

12

34 1

permite construir o polinómio

P4(y) =14

y(y−

√2

2

)(y− sin

(3π8

))(y−1)

sin(π

8

)(sin(π

8

)−√

22

)(sin(π

8

)− sin

(3π8

))(sin(π

8

)−1)

+12

y(y− sin

(π8

))(y− sin

(3π8

))(y−1)

√2

2

(√22 − sin

(π8

))(√22 − sin

(3π8

))(√22 −1

)

+34

y(y− sin

(π8

))(y−

√2

2

)(y−1)

sin(3π

8

)(sin(3π

8

)− sin

(π8

))(sin(3π

8

)−√

22

)(sin(3π

8

)−1)

+y(y− sin

(π8

))(y−

√2

2

)(y− sin

(3π8

))(1− sin

(π8

))(1−

√2

2

)(1− sin

(3π8

)) .

(3.3.3)


Atendendo a que P4(0.95) = 0.8211, obtemos x∗ ' 0.8211. Notamos que f (0.8211) = 0.960774. NaFigura 3.2 ilustramos o comportamento de P4 bem como da função inversa

f−1(y) =2π

arcsin(y), y ∈ [0,1]. (3.3.4)

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1 f−1

P4pontos de interpolação

Figura 3.2: Gráficos de P4, definido por (3.3.3), e de f−1(y), definida por (3.3.4).

Vimos dois processos diferentes de determinar o polinómio interpolador:

1. utilização dos polinómios de Lagrange,

2. utilização das diferenças divididas.

No entanto, quando os pontos da partição são igualmente distanciados, isto é, xi− xi−1 = h, i = 1, . . . ,n, aconstrução do polinómio interpolador pode ser feita usando um processo mais simples e mais eficienterecorrendo às chamadas diferenças progressivas.

Definição 3.1. Seja h > 0 uma quantidade fixa e f uma função definida em x e x+h. Chamamosdiferença progressiva de primeira ordem de f em x a

∆ f (x) = f (x+h)− f (x).

Se f é definida em x+ ih, i = 0, . . . ,n, a ∆n−1(∆ f (x)) chamamos diferença progressiva de ordemn≥ 2 de f em x e esta diferença progressiva é denotada por ∆n f (x).

A diferença progressiva de primeira ordem de f em x está relacionada com a diferença divididaf [x,x+h]. De facto, temos

f [x,x+h] =∆ f (x)

h.

Mais geralmente, vale a igualdade

f [x0,x1, . . . ,xn] =∆n f (x0)

n!hn ,


em que xi = xi−1 +h, i = 1, . . . ,n.Atendendo à igualdade anterior, o polinómio interpolador de f em xi, i = 0, . . . ,n, xi = xi−1 +h, i =

1, . . . ,n, pode ser escrito na forma

Pn(x) = f (x0)+n

∑j=1

∆ j f (x0)

h j j!

j−1

∏i=0

(x− xi).

A finalizar, observamos que o cálculo das diferenças progressivas pode ser organizado num quadro análogoao construído para as diferenças divididas.

Relativamente à implementação do cálculo do polinómio anteriormente definido com dife-renças progressivas, observamos que os Algoritmos 3.3 e 3.4 podem ser facilmente adaptados aeste caso.

Para determinar o polinómio interpolador para uma função f definida num intervalo [a,b] coloca-se aquestão de como escolher a partição no intervalo. Será que podemos escolher a partição de uma formaarbitrária ou a escolha dos pontos da partição tem influência no erro cometido?

Consideremos a função de Runge

f (x) =1

1+25x2 , x ∈ [−1,1], (3.3.5)

e as partições

−1,−34,−1

2,−1

4,0,

14,12,34,1,

e

−1,−78,−5

8,−1

2,0,

12,58,78,1.

−1 −0.5 0 0.5 1

−1

−0.5

0

0.5

1

f P8 pontos de interpolação

−1 −0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1

Figura 3.3: Gráficos da função de Runge (3.3.5) e P8 (esquerda) e gráfico do erro E8 (direita).


Sejam P8 e P∗8 os polinómios interpoladores de f para os dois conjuntos de pontos anteriores edenotemos por E8 e E∗8 os erros seguintes

E8(x) = | f (x)−P8(x)|, x ∈ [−1,1], E∗8 (x) = | f (x)−P∗8 (x)|, x ∈ [−1,1].

Os comportamentos de P8,E8 e de P∗8 ,E∗8 são ilustrados nas Figuras 3.3 e 3.4, respectivamente.

−1 −0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1

f P∗8 pontos de interpolação

−1 −0.5 0 0.5 1

0

5 ·10−2

0.1

0.15

0.2

0.25

Figura 3.4: Gráficos da função de Runge (3.3.5) e P∗8 (esquerda) e gráfico do erro E∗8 (direita)

Podemos facilmente verificar que o máximo do valor absoluto do erro é inferior quando os pontos dapartição são considerados próximos das extremidades do intervalo. Assim, para este caso, a escolha dospontos influência o erro.

3.4 Polinómio interpolador de Lagrange segmentado

Consideremos uma função f definida num intervalo [a,b] onde está definida uma partição uniformecom n+1 pontos (3.2.1). Seja Pn o polinómio interpolador de Lagrange. Provámos que se existe umaconstante positiva M tal que ∣∣ f (n+1)(x)

∣∣≤M, ∀x ∈ [a,b],

então

| f (x)−Pn(x)| ≤M(b−a)n+1

4(n+1)n(n+1) , x ∈ [a,b].

Atendendo a que

limn→+∞

(b−a)n+1

4(n+1)n(n+1) = 0,

concluímos que| f (x)−Pn(x)| → 0, n→+∞, x ∈ [a,b].

Assim, para certas funções, o aumento do grau do polinómio implica um aumento de precisão.No entanto existem funções para as quais não podemos concluir que a um aumento do grau do

polinómio interpolador corresponda um aumento da proximidade entre o polinómio interpolador e afunção interpolada, isto é, uma diminuição do erro | f (x)−Pn(x)| em [a,b]. Um exemplo de uma função


que ilustra a situação anterior é a conhecida função de Runge (3.3.5). De facto, consideremos para afunção (3.3.5) o polinómio interpolador de Lagrange de grau 4 - P4(x) - para a partição

−1,−12,0,

12,1,

e sejaE4(x) = | f (x)−P4(x)|, x ∈ [−1,1].

Seja ainda P8 o polinómio interpolador de Lagrange para a função de Runge (3.3.5) e para a partição

−1,−34,−1

2,−1

4,0,

14,12,34,1

eE8(x) = | f (x)−P8(x)|, x ∈ [−1,1].

Na Figura 3.5 ilustramos o comportamento das funções E4 e E8.

−1 −0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1

E4 E8

Figura 3.5: Gráficos das funções E4 e E8.

Observamos que se temmax

x∈[−1,1]E8(x)> max

x∈[−1,1]E4(x).

A construção de um polinómio interpolador que “aproxime” a função de Runge (3.3.5) com erroinferior ao dos polinómios P4 e P8 é conseguida com polinómios de grau inferior a 4. Assim, considerandoos pontos −1,−1

2 ,0, e 0, 12 ,1, e a decomposição de [−1,1] nos subintervalos [−1,0] e [0,1], podemos

facilmente obter os polinómios interpoladores de Lagrange de grau 2, s(1) e s(2) para a função de Runge(3.3.5) e para os intervalos definidos. Seja

s(x) =

s(1)(x) x ∈ [−1,0]s(2)(x) x ∈ ]0,1].

(3.4.1)

Esta função é um “polinómio” de grau 2, interpolador da função de Runge (3.3.5) nos pontos −1, −12 , 0,

12 e 1.


O polinómio s definido por (3.4.1) é um polinómio interpolador da função f para a partição consideradamas diferente do polinómio de Lagrange anteriormente estudado. Para realçar a diferença existente entreestes polinómios vamos adoptar a designação polinómio interpolador de Lagrange segmentado. Na Figura3.6 ilustramos o comportamento do polinómio s e da função de Runge (3.3.5). A comparação entre ocomportamento do erro do polinómio interpolador segmentado de grau 2 e do polinómio interpolador degrau 4 para os pontos −1,−1

2 ,0,12 ,1, é apresentada na Figura 3.7.

−1 −0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1

f s pontos de interpolação

Figura 3.6: Gráficos da função de Runge (3.3.5) e de s.

−1 −0.5 0 0.5 1

0

0.1

0.2

0.3

0.4

E4 Es2

Figura 3.7: Gráficos do erro associado ao polinómio interpolador segmentado de grau 2, Es2, e ao polinómio

interpolador de grau 4, E4, para os pontos −1,−12 ,0,

12 ,1.


Definição 3.2. Seja f uma função definida num intervalo [a,b] onde consideramos a partição (3.2.1).Se n é divisível por m, então chamamos polinómio interpolador de Lagrange segmentado, de graum, ao polinómio

s(x) =

s(1)(x) x ∈ [x0,xm],

s(2)(x) x ∈ ]xm,x2m],...

s( j)(x) x ∈ ]x( j−1)m,x jm],...

s(q)(x) x ∈ ]xn−m,xn],

em que n = mq, e, para i = 1, . . . ,q, s(i) é o polinómio interpolador de Lagrange de grau m para f em[x(i−1)m,xim] considerando a partição

x(i−1)m < x(i−1)m+1 < · · ·< xim.

Para construir cada um dos ramos do “polinómio” s podemos usar um dos processos introduzidos:utilização dos polinómios de Lagrange ou a utilização das diferenças divididas. Neste caso podemosdefinir, de modo análogo ao anterior, o polinómio interpolador de Newton segmentado.

Exemplo 3.7. Determinemos o polinómio interpolador segmentado quadrático para uma função fdefinida no intervalo [a,b] para a partição

a = x0 < x1 < x2 < · · ·< x6 = b.

Temos, utilizando as funções de Lagrange,

s(x) =

f (a) (x−x1)(x−x2)(x0−x1)(x0−x2)

+ f (x1)(x−x0)(x−x2)(x1−x0)(x1−x2)

+ f (x2)(x−x0)(x−x1)(x2−x0)(x2−x1)

x ∈ [a,x2],

f (x2)(x−x3)(x−x4)(x2−x3)(x2−x4)

+ f (x3)(x−x2)(x−x4)(x3−x2)(x3−x4)

+ f (x4)(x−x2)(x−x3)(x4−x2)(x4−x3)

x ∈ ]x2,x4],

f (x4)(x−x5)(x−x6)(x4−x5)(x4−x6)

+ f (x5)(x−x4)(x−x6)(x5−x4)(x5−x6)

+ f (b) (x−x4)(x−x5)(b−x4)(b−x5)

x ∈ ]x4,b],

e, utilizando as diferenças divididas,

s(x) =

f (a)+ f [a,x1](x−a)+ f [a,x1,x2](x−a)(x− x1) x ∈ [a,x2],

f (x2)+ f [x2,x3](x− x2)+ f [x2,x3,x4](x− x2)(x− x3) x ∈ ]x2,x4],

f (x4)+ f [x4,x5](x− x4)+ f [x4,x5,x6](x− x4)(x− x5) x ∈ ]x4,b].

Seja f : [a,b]−→R uma função da qual se conhece f (xi), i= 0, . . . ,n, em que xi, i= 0, . . . ,n satsifaz(3.2.1). Seja s o polinómio interpolador de Lagrange segmentado de grau m, isto é, para i = 1, . . . ,q, ex ∈ [x(i−1)m,xim], s(x) é dado por

s(x) = s(i)(x),

em que s(i) denota o polinómio interpolador de Lagrange de grau m nos pontos x j, j = (i−1)m, . . . , im.Se x ∈ [x(i−1)m,xim], então

f (x)− s(x) = f (x)− s(i)(x)

=f (m+1)(η(x))(m+1)!

im

∏j=(i−1)m

(x− x j).


Assim, se existir M > 0 tal que

| f (m+1)(x)| ≤M, x ∈ [x(i−1)m,xim], i = 1, . . . ,q

deduzimos

| f (x)− s(x)| ≤ M(m+1)!

im

∏j=(i−1)m

|(x− x j)|.

A partir da desigualdade anterior concluímos várias estimativas para o erro do polinómio interpolador deLagrange segmentado de grau m:

1. ‖ f − s‖∞ ≤M

(m+1)!

(max

j=1,...,q(x jm− x( j−1)m)

)m

,

2. ‖ f − s‖∞ ≤M

4(m+1)hm+1, h = max

i=1,...,n(xi− xi−1),

3. ‖ f − s‖∞ ≤M

4(m+1)nm+1 (b−a)m+1, se a partição é uniforme.

Exemplo 3.8. Consideremos a função f (x)= sin(πx), x∈ [12 ,

32 ], e, neste intervalo, a partição 1

2 ,34 ,1,

54 ,

32 .

Determinemos o polinómio interpolador segmentado quadrático de f nos pontos dados, uma aproximaçãopara f (7

8) e uma estimativa para o erro da aproximação.

O polinómio pretendido é

s(x) =

s(1)(x), x ∈ [12 ,1],

s(2)(x), x ∈ ]1, 32 ],

em que

s(1)(x) = sin(π

2

) (x− 34

)(x−1)(1

2 − 34

)(12 −1

) + sin(3π

4

) (x− 12

)(x−1)(3

4 − 12

)(34 −1

) + sin(π) (x− 1

2

)(x− 3

4

)(1− 1

2

)(1− 3

4

) , x ∈[

12,1]

e

s(2)(x) = sin(π) (x− 5

4)(x− 32)

(1− 54)(1− 3

2)+ sin

(5π4

) (x−1)(x− 32)

(54 −1)(5

4 − 32)

+ sin(3π

2

) (x−1)(x− 54)

(32 −1)(3

2 − 54), x ∈

[1,

32

].

Entãof(7

8

)' s(1)

(78

)= 8(7

8− 3

4

)(78−1)−8√

2(7

8− 1

2

)(78−1).

Podemos agora determinar uma estimativa para o erro da aproximação obtida. Usando (3.2.4)obtemos

| f (x)− s(1)(x)| ≤ π3

3!

∣∣∣∣(

x− 12

)(x− 3

4

)(x−1)

∣∣∣∣, x ∈[

12,1],

e, para x = 78 , vem ∣∣∣∣ f

(78

)− s(7

8

)∣∣∣∣≤π3

210 . (3.4.2)

Usando agora (3.2.5), obtemos

| f (x)− s(1))(x)| ≤ π3

768, x ∈

[12,1]. (3.4.3)

Observamos que a estimativa para o erro (3.4.3) é mais “grosseira” do que a estimativa (3.4.2), eainda que a estimativa (3.4.3) é válida para x ∈ [1, 3

2 ] pois que | f ′′′(x)| ≤ π3, x ∈ [1, 32 ].


3.5 Polinómio interpolador de Hermite

Em muitas questões práticas é conhecida uma função num conjunto de pontos e ainda o coeficienteangular da tangente ao gráfico da função nesses pontos. Nestas circunstâncias, para aproximar a funçãof , determinamos um polinómio que coincida com f nos pontos dados e cujo gráfico tenha nos nodos deinterpolação tangente com o coeficiente angular dado. Mais especificamente, coloca-se o problema:

Seja f definida no intervalo [a,b]. Suponha que são conhecidos f e f ′ nos pontos dapartição (3.2.1). Pretendemos determinar um polinómio H(x) tal que

H(xi) = f (xi), H ′(xi) = f ′(xi), i = 0, . . . ,n. (3.5.1)

O polinómio a determinar, interpolador para f e para a sua derivada, deve verificar as 2n+2 condições(3.5.1). Logo, se existir tal polinómio, este deve ter grau inferior ou igual a 2n+ 1. No resultado queapresentamos seguidamente estabelecemos a existência e unicidade do polinómio que verifica (3.5.1).

Teorema 3.5. Seja f uma função definida num intervalo [a,b] onde consideramos a partição(3.2.1). Dados f (xi), f ′(xi), i = 0, . . . ,n, existe um único polinómio que denotamos por H2n+1, degrau menor ou igual a 2n+1, que verifica (3.5.1)

Demonstração. Provemos a unicidade. Sejam F e G dois polinómios que verificam as condições (3.5.1)e seja R o polinómio seguinte

R(x) = F(x)−G(x), x ∈ [a,b].

Este polinómio tem grau inferior ou igual a 2n+1 e verifica a

R(xi) = R′(xi) = 0, i = 0, . . . ,n.

Então R tem pelo menos 2n+2 zeros reais. Logo R é o polinómio nulo. Concluímos que se existir umpolinómio que verifica (3.5.1) então este polinómio é único.

Provemos a existência. Atendendo às condições consideradas, o polinómio, a existir, pertence aoespaço P2n+1 dos polinómios de grau menor ou igual a 2n+1. Consideremos os polinómios de Lagrange

ì(x) = ∏j 6=i

x− x j

xi− x j, i = 0, . . . ,n.

Definamos as seguintes 2n+2 funções polinomiais hi e hi

hi(x) = [1−2`′i(xi)(x− xi)]ì(x)2, hi(x) = (x− xi)ì(x)2, i = 0, . . . ,n.

As funções anteriores são polinómios de grau menor ou igual a 2n+1 e satisfazem as condições:

hi(x j) = δi j, i, j = 0, . . . ,n, h′i(x j) = 0, i, j = 0, . . . ,n,hi(x j) = 0, i, j = 0, . . . ,n, h

′i(x j) = δi j, i, j = 0, . . . ,n.

Atendendo às condições anteriores, é fácil demonstrar que hi, hi, i = 0, . . . ,n é uma base de P2n+1.Assim, a existir o polinómio H2n+1, este admite a representação

H2n+1(x) =n

∑i=0

Cihi(x)+n

∑i=0

Cihi(x).


com Ci,Ci ∈ R,ß = 0, . . . ,n.Observamos que provar a existência de H2n+1 é equivalente a provar a existência de Ci,Ci, i = 0, . . . ,n,

tais que H2n+1 verifica (3.5.1). Atendendo a que H2n+1(xi) = f (xi), obtemos Ci = f (xi), i = 0, . . . ,n.Por outro lado, de H ′2n+1(xi) = f ′(xi), vem Ci = f ′(xi), i = 0, . . . ,n. Finalmente, obtemos a expressão deH2n+1

H2n+1(x) =n

∑i=0

f (xi)hi(x)+n

∑i=0

f ′(xi)hi(x). (3.5.2)

O polinómio anterior é chamado polinómio interpolador de Hermite de f nos pontos xi, i = 0, . . . ,n.

Exemplo 3.9. Consideremos a tabela

xi 0 12 1

f (xi) 0√

22 1

f ′(xi)π2

π2

√2

2 0

O polinómio interpolador de Hermite, para os dados anteriores, admite a representação

H5(x) = 2π(x−1)2(x− 1

2

)2x+8√

2(x−1)2x2+

+4√

2π(x−1)2(x− 1

2

)x2 +4(1−6(x−1))

(x− 1

2

)2x2.

Atendendo a que a tabela anterior diz respeito à função f (x) = sin(πx2 ), na Figura 3.8 ilustramos o

comportamento do polinómio de Hermite de grau 5 e da sua derivada.

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8

1

f H5 pontos de interpolação

0 0.2 0.4 0.6 0.8 1

0

0.5

1

1.5

f ′ H ′5 pontos de interpolação

Figura 3.8: Gráficos de f (x) = sin(πx2 ) e H5 (esquerda) e das suas derivadas (direita).

Vamos seguidamente utilizar diferenças divididas para construir o polinómio interpolador de Hermite.Na construção do polinómio interpolador de Hermite de grau 2n + 1 consideramos as funções

hi(x),hi(x), i = 0, . . . ,n, que constituem uma base de P2n+1. O polinómio anterior pode ser construídoutilizando outra base para o espaço referido.


Consideremos as funções

ψ0(x) = 1,

ψ1(x) = (x− x0),

ψ2(x) = (x− x0)2,

ψ3(x) = (x− x0)2(x− x1),

...

ψ2n+1(x) = (x− x0)2 · · ·(x− xn−1)

2(x− xn).

(3.5.3)

Estas funções são linearmente independentes. Logo H2n+1 admite a representação

H2n+1(x) =2n+1

∑j=0

A jψ j(x).

Determinemos as constantes Ai, i = 0, . . . ,2n+1. Uma vez que H2n+1(x0) = f (x0), obtemos

A0 = f (x0).

Por outro lado, atendendo a que H ′2n+1(x0) = f ′(x0), deduzimos

A1 = f ′(x0)

que representamos porf [x0,x0].

Substituindo as constantes determinadas em H2n+1 e considerando as condições para x = x1, vem

A2 = f [x0,x0,x1]

eA3 = f [x0,x0,x1,x1].

Repetindo o processo anterior, obtemos os restantes coeficientes Ai, i = 4, . . . ,2n+ 1 e, deste modo,estabelecemos para H2n+1(x) a expressão

H2n+1(x) = f (x0)+ f ′(x0)(x− x0)+ f [x0,x0,x1](x− x0)2

+ f [x0,x0,x1,x1](x− x0)2(x− x1)

+ · · ·++ f [x0,x0, . . . ,xn,xn](x− x0)

2 · · ·(x− xn−1)2(x− xn).

(3.5.4)

O polinómio interpolador de Hermite H2n+1 pode ser determinado utilizando a tabela das diferençasdivididas como veremos seguidamente. Consideremos a mudança de variável

y0 = x0, y1 = x0, y2 = x1, y3 = x1, . . . ,y2n = xn, y2n+1 = xn,

e seja g tal queg(y0) = f (x0),

g(y1) = f (x0),

g(y2) = f (x1),

g(y3) = f (x1),

...

g(y2n) = f (xn),

g(y2n+1) = f (xn).


Podemos então definir o polinómio

H(x) = g(y0)+g[y0,y1](x− x0)+g[y0,y1,y2](x− x0)2

+g[y0,y1,y2,y3](x− x0)2(x− x1)+

+ · · ·++g[y0,y1, . . . ,y2n+1](x− x0)

2(x− x1)2 · · ·(x− xn−1)

2(x− xn)

ondeg[y0,y1] = f [x0,x0],

g[y0,y1,y2] = f [x0,x0,x1],

...

g[y0,y1, . . . ,y2n,y2n+1] = f [x0,x0, . . . ,xn,xn],

(3.5.5)

As constantes (3.5.5) podem ser determinadas utilizando a tabela das diferenças divididas para g.

Exemplo 3.10. Determinemos o polinómio de Hermite de menor grau que é interpolador de f (x) = e−x2,

x ∈ [−1,1].Consideremos (3.5.4) com n = 1 e a tabela seguinte

xi f (xi) f ′(xi)

−1 e−1 2e−1

1 e−1 −2e−1

Com os dados anteriores construimos a seguinte tabela de diferenças divididas

xi f (xi) f [·, ·] f [·, ·, ·] f [·, ·, ·, ·]x0 f (x0)

f ′(x0)x0 f (x0) f [x0,x0,x1]

f [x0,x1] f [x0,x0,x1,x1]x1 f (x1) f [x0,x1,x1]

f ′(x1)x1 f (x1)

que, para o caso particular deste exemplo, toma a forma

xi f (xi) f [·, ·] f [·, ·, ·] f [·, ·, ·, ·]−1 e−1

2e−1

−1 e−1 −e−1

0 01 e−1 −e−1

−2e−1

1 e−1

e portantoH3(x) = e−1 +2e−1(x+1)− e−1(x+1)2, x ∈ [−1,1].


A construção de um polinómio interpolador H que verifique outras restrições de regularidade para asderivadas, além de H ′(xi) = f ′(xi), i = 0, . . . ,n, é feita de modo análogo à construção do polinómio deHermite usando diferenças divididas. O exemplo seguinte ilustra o procedimento a seguir para um casoespecífico.

Exemplo 3.11. Seja f uma função definida num intervalo [a,b] e suponha dada a tabela

xi x0 x1

f (xi) f (x0) f (x1)f ′(xi) f ′(x0) f ′(x1)f ′′(xi) f ′′(x0)

Determinar o polinómio interpolador H(x) que verifique

H(xi) = f (xi), H ′(xi) = f ′(xi), i = 0,1, H ′′(x0) = f ′′(x0).

O polinómio H(x), a existir, vai pertencer ao espaço P4. Consideremos neste espaço as funçõesseguintes

ψ0(x) = 1,

ψ1(x) = (x− x0),

ψ2(x) = (x− x0)2,

ψ3(x) = (x− x0)3,

ψ4(x) = (x− x0)3(x− x1).

(3.5.6)

O conjunto ψ, i = 0, . . . ,4 é uma base de P4 e portanto

H(x) =4

∑j=0

A jψ(x).

Determinemos as constantes A j, j = 0, . . . ,4, utilizando as condições impostas a este polinómio: Aten-dendo a que H(x0) = f (x0), vem A0 = f (x0). Por outro lado, atendendo a que H ′(x0) = f ′(x0), obtemosA1 = f [x0,x0], que, tal como anteriormente, denotamos por f [x0,x0]. Mais ainda, atendendo a queH ′′(x0) = f ′′(x0), vem A2 =

f ′′(x0)2 que representamos por f [x0,x0,x0]. De H(x1) = f (x1), obtemos

A3 = f [x0,x0,x0,x1].

Considerando agora H ′(x1) = f ′(x1), vem

A4 =f ′(x1)− f [x0,x0]−2 f [x0,x0,x0](x1− x0)−3 f [x0,x0,x0,x1](x1− x0)

2

(x1− x0)3 ,

ou seja,A4 = f [x0,x0,x0,x1,x1].

Finalmente estabelecemos

H(x) = f (x0)+ f [x0,x0](x− x0)+ f [x0,x0,x0](x− x0)2

+ f [x0,x0,x0,x1](x− x0)3 + f [x0,x0,x0,x1,x1](x− x0)

3(x− x1).

Observamos que a implementação do cálculo do polinómio de Hermite se pode fazer combase nos Algoritmos 3.3 e 3.4.

Seguidamente vamos estudar o erro cometido ao aproximar uma função pelo seu polinómio de Hermitede grau 2n+1.


Teorema 3.6. Seja f uma função com derivadas até à ordem 2n+1 contínuas em [a,b] e comderivada de ordem 2n+ 2 em (a,b). Em [a,b] consideremos a partição (3.2.1) e seja H2n+1 opolinómio interpolador de Hermite para a função f e relativamente à partição anterior. Sex ∈ [a,b], então existe ε(x) ∈ (a,b) tal que

f (x)−H2n+1(x) =f (2n+2)(ε(x))(2n+2)!

n

∏j=0

(x− x j)2. (3.5.7)

Demonstração. Seja x ∈ [a,b]. Se x = xi, i = 0, . . . ,n, então o teorema está provado. Suponhamos agoraque x 6= xi, i = 0, . . . ,n e consideremos a função

g(t) = f (t)−H2n+1(t)− [ f (x)−H2n+1(x)]n

∏j=0

(t− x j)2

(x− x j)2 , t ∈ [a,b].

A função g tem o valor zero em xi, i = 0, . . . ,n, (zeros duplos) e em x (zero simples). Por aplicaçãosucessiva do teorema de Rolle, concluímos que g(2n+2) tem pelo menos um zero em (a,b), isto é, existeε(x) ∈ (a,b) tal que

0 = g(2n+2)(ε(x)) = f (2n+2)(ε(x))− [ f (x)−H2n+1(x)](2n+2)!

∏nj=0(x− x j)2 .

Esta última igualdade finaliza a demonstração.

A partir o resultado anterior, estabelecemos seguidamente estimativas para ‖ f −H2n+1‖∞ :

1. Se | f (2n+2)(x)| ≤M,x ∈ [a,b], de (3.5.7), vem

| f (x)−H2n+1(x)| ≤M

(2n+2)!

n

∏j=0

(x− x j)2.

2. Atendendo a que para x ∈ [a,b], |x− xi| ≤ (b−a), então, de (3.5.7), deduzimos

‖ f −H2n+1‖∞ ≤M

(2n+2)!(b−a)2n.

3. Provámos anteriormente a desigualdade∣∣∣∣

n

∏j=0

(x− x j)

∣∣∣∣≤hn+1n!

4,

em que h = maxi=1,...,n

(xi− xi−1). Então

‖ f −H2n+1‖∞ ≤Mh2(n+1)(n!)2

16(2n+2).

Observamos que dependendo do comportamento de | f (2n+2)| podemos, ou não, concluir que o aumentodo grau do polinómio interpolador de Hermite implique uma diminuição do erro cometido ao aproximara função por este polinómio. Uma forma de minimizar o erro consiste na utilização de polinómiosinterpoladores de Hermite de baixo grau.


Definição 3.3. Seja f uma função definida num intervalo [a,b] onde consideramos a partiçãoa = x0 < x1 < · · · < xn−1 < xn = b. Sejam m,q ∈ N tais que mq = n. O polinómio s(x) definidopor (3.5.8)

s(x) =

H(1)2m+1(x), x ∈ [x0,xm]

H(2)2m+1(x), x ∈ ]xm,x2m]

...

H( j)2m+1(x), x ∈ ]x( j−1)m,x jm],

...

H(q)2m+1(x), x ∈ ]xn−m,xn],

(3.5.8)

em que H( j)2m+1 denota o polinómio interpolador de Hermite de grau 2m+1 para os pontos

x( j−1)m,x( j−1)m+1, . . . ,x jm, j = 1, . . . ,q,

é designado polinómio interpolador de Hermite segmentado de grau 2m+1.

Exemplo 3.12. Seja f : [−1,1] −→ R a função de Runge (3.3.5). O comportamento do polinómiointerpolador de Hermite de grau 9, H9, para f e para a sua derivada f ′ nos pontos

−1,−12,0,

12,1,

e da sua derivada são ilustrados na Figura 3.9. Na Figura 3.10 ilustramos o comportamento dos errosdos polinómios anteriores.

−1 −0.5 0 0.5 10

0.2

0.4

0.6

0.8

1

f H9 pontos de interpolação

−1 −0.5 0 0.5 1

−2

0

2

f ′ H ′9 pontos de interpolação

Figura 3.9: Gráficos da função de Runge e de H9 (esquerda) e das suas derivadas (direita).


−1 −0.5 0 0.5 1

0

5 ·10−2

0.1

0.15

0.2E9

−1 −0.5 0 0.5 1

0

0.5

1

1.5

| f −H ′9|

Figura 3.10: Gráficos dos erros | f (x)−H9(x)|, x ∈ [−1,1], (esquerda) e de | f ′(x)−H ′9(x)|, x ∈ [−1,1](direita).

Consideremos agora os conjuntos de pontos

−1,−12,0 e 0,

12,1,

e seja s o polinómio interpolador de Hermite segmentado de grau 5, definido por

s5(x) =

1− 10029 x2− 10000

841 x2(x+ 1

2

)− 151250

10933 x2(x+ 1

2

)2

−2156250142129 x2

(x+ 1

2

)2(x+1), x ∈ [−1,0].

1− 10029 x2 + 10000

841 x2(x− 1

2

)− 151250

10933 x2(x− 1

2

)2

+2156250142129 x2

(x− 1

2

)2(x−1), x ∈ ]0,1].

O comportamento do polinómio segmentado anterior e da sua derivada está ilustrado na Figura 3.11. NaFigura 3.12 ilustramos o comportamento dos erros deste polinómio e da sua derivada.

−1 −0.5 0 0.5 10

0.2

0.4

0.6

0.8

1

f s5 pontos de interpolação

−1 −0.5 0 0.5 1

−2

0

2

f ′ s′5 pontos de interpolação

Figura 3.11: Gráficos da função de Runge (3.3.5) e de s5 (esquerda) e das suas derivadas (direita).


−1 −0.5 0 0.5 1

0

5 ·10−2

0.1

| f − s5|

−1 −0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1| f ′− s′5|

Figura 3.12: Gráficos dos erros | f (x)− s5(x)|, x ∈ [−1,1], (esquerda) e de | f ′(x)− s′5(x)|, x ∈ [−1,1](direita).

As Figuras 3.10 e 3.12 ilustram a redução do erro do polinómio de Hermite quando se considera opolinómio de Hermite segmentado de grau 5.

Seja f : [a,b]−→R uma função da qual se conhece f (xi), f ′(xi), i= 0, . . . ,n, em que xi, i= 0, . . . ,nsatisfaz (3.2.1). Seja s o polinómio interpolador de Hermite segmentado de grau 2m+ 1 definido por(3.5.8). Se x ∈ [a,b], então x ∈ [x(i−1)m,xim], para algum i ∈ 1, . . . ,q. Logo

s(x) = s(i)(x),

em que s(i) é o polinómio interpolador de Hermite de grau 2m+1 nos pontos x j, j = (i−1)m, . . . , im.Assim

f (x)− s(x) = f (x)− s(i)(x)

=f (2m+2)(η(x))(2m+2)!

im

∏j=(i−1)m

(x− x j)2.

Assim, se existir M > 0 tal que

| f (2m+2)(x)| ≤M, x ∈ [x(i−1)m,xim], i = 1, . . . ,q

deduzimos

| f (x)− s(x)| ≤ M(2m+2)!

im

∏j=(i−1)m

(x− x j)2.

A partir da desigualdade anterior concluímos várias estimativas para o erro do polinómio interpolador deHermite segmentado de grau 2m+1 segmentado:

1. ‖ f − s‖∞ ≤ M(2m+2)!

(max j=1,...,q(x jm− x( j−1)m)

)2m

,

2. ‖ f − s‖∞ ≤ Mh2m+2(m!)2

16(2m+2)! , h = maxi=1,...,n(xi− xi−1),

3. ‖ f − s‖∞ ≤ M(b−a)2m+2(m!)2

16(2m+2)!nm , se a partição é uniforme.


3.6 Splines

Consideremos para a função de Runge (3.3.5) e o polinómio de Lagrange segmentado quadrático sL,2para os pontos −1,−1

2 ,0,12 ,1. Denotemos por sH,3 o polinómio de Hermite segmentado cúbico para a

função de Runge (3.3.5) nos pontos −1,0,1. O comportamento das derivadas s′L,2 e s′H,3 é ilustrado naFigura 3.13. Na Figura 3.14 ilustramos o comportamento de s′′H,3 e s′′′H,3. Para o exemplo de Runge e paraos polinómios anteriores, podemos mostrar que não existe s′L,2(0), s′H,3 e s′′H,3 são funções contínuas em[−1,1] e não existe s′′′H,3(0).

−1 −0.5 0 0.5 1

−2

0

2

f ′

s′L,2s′H,3

Figura 3.13: Comportamento das derivadas de primeira ordem dos polinómios de Lagrange segmentadoquadrático e de Hermite segmentado cúbico para a função (3.3.5).

−1 −0.5 0 0.5 1

−40

−20

0

f ′′

s′′H,3

−1 −0.5 0 0.5 1

−500

0

500f ′′′

s′′′H,3

Figura 3.14: Comportamento das derivadas do polinómio de Hermite segmentado cúbico de segundaordem (esquerda) e de terceira ordem (direita) para a função (3.3.5).


Consideremos agora a função

f (x) =

0, x ∈ [−1,0]x2, x ∈ ]0,1]

, (3.6.1)

e seja sH o polinómio de Hermite segmentado cúbico para esta função nos pontos

−1,−12,0,

12,1.

Na Figura 3.15 ilustramos o comportamento de sH e s′H . O comportammento de s′′H é ilustrado na Figura3.16. Observamos que esta derivada não existe em x =−1

2 e x = 12 .

−1 −0.5 0 0.5 10

0.2

0.4

0.6

0.8

1 fsH

−1 −0.5 0 0.5 1

−2

0

2

f ′

s′H

Figura 3.15: Comportamento do polinómio de Hermite segmentado cúbico (esquerda) e da sua derivada(direita) para a função (3.6.1).

−1 −0.5 0 0.5 1

−40

−20

0

20

f ′′

s′′H

Figura 3.16: Comportamento da segunda derivada do polinómio de Hermite segmentado cúbico para afunção (3.6.1).


A questão que pretendemos resolver seguidamente diz respeito à construção de um polinómio in-terpolador segmentado com segunda derivada contínua. Com o objectivo de introduzir o processo deconstrução de tal polinómio, começamos por introduzir o conceito de spline interpolador de grau 2.

Definição 3.4. Seja f uma função definida num intervalo [a,b] onde consideramos a partição a =x0 < x1 < · · ·< xn−1 < xn = b. O polinómio segmentado

s(x) =

s(1)(x), x ∈ [x0,x1],

s(2)(x), x ∈ ]x1,x2],...s( j)(x), x ∈ ]x j−1,x j],...s(n)(x), x ∈ ]xn−1,xn],

(3.6.2)

em que s( j) é um polinómio quadrático em ]x j−1,x j], que verifica

s(x j) = f (x j), j = 0, . . . ,n,s é contínua em xi, i = 1, . . . ,n−1,s′ é contínua em xi, i = 1, . . . ,n−1,s′(x0) = f ′(x0)ou s′(xn) = f ′(xn),

é designado spline interpolador quadrático.

Observamos que, por definição, a construção do spline interpolador quadrático requer a determinação de3n parâmetros e este número é igual ao número de condições que são exigidas ao spline s.

Consideremos para s(i)′a seguinte expressão

s(i)′(xi) = Mi

xi− xhi

+Mi+1x− xi−1

hi. (3.6.3)

A definição anterior induz de imediato a continuidade de s′ em xi. De facto, tem-se sucessivamente

limx→x−i

s′(x) = limx→x−i

s(i)′(x)

= limx→x−i

Mixi− x

hi+Mi+1

x− xi−1

hi

= Mi+1

= limx→x+i

Mi+1xi+1− x

hi+1+Mi+2

x− xi

hi+1

= limx→x+i

s(i+1)′(x)

= limx→x+i

s′(x).

De (3.6.3), vem

s(i)(x) = Mi(x− xi)

2

2hi+Mi+1

(x− xi−1)2

2hi+C.

Uma vez que s(xi) = s(i)(xi), para i = 1, . . . ,n, deduzimos

C = f (xi)−Mi+1hi

2.


Deste modo temos

s(i)(x) = Mi(x− xi)

2

2hi+Mi+1

(x− xi−1)2

2hi+ f (xi)−

Mi+1hi

2. (3.6.4)

Para i = 0, temos também s(1)(x0) = f (x0), e portanto

M2 =2h1

(f (x0)− f (x1)

)−M1. (3.6.5)

Das duas relações anteriores obtemos

s(1)(x) = M1(x− x1)

2

2h1+

(2h1

(f (x0)− f (x1)

)−M1

)(x− xi−1)

2

2hi+

2h1

(f (x0)− f (x1)

)−M1. (3.6.6)

Atendendo a que s é contínua em xi, i = 1, . . . ,n−1, vem sucessivamente

limx→x−i

s(x) = limx→x−i

s(i)(x)

= limx→x−i

Mi(x− xi)

2

2hi+Mi+1

(x− xi−1)2

2hi+ f (xi)−

Mi+1hi

2

= f (xi)

(3.6.7)

e

limx→x+i

s(x) = limx→x+i

s(i+1)(x)

= limx→x+i

Mi+1(x− xi+1)

2

2hi+1+Mi+2

(x− xi)2

2hi+1+ f (xi+1)−

Mi+2hi+1

2

= Mihi+1

2+ f (xi+1)−

Mi+2hi+1

2.

(3.6.8)

Conjugando (3.6.7) e (3.6.8) obtemos

Mi+1hi+1

2+ f (xi+1)−

Mi+2hi+1

2= f (xi), i = 1, . . . ,n−1. (3.6.9)

De (3.6.9), a partir de Mi+1 calculamos Mi+2, i = 1, . . . ,n−1, iniciando-se o processo com M2. Observa-mos que, de (3.6.5), dado M1, calculamos M2. Assim, a determinação de todos os coeficientes é feita apartir de (3.6.9), (3.6.5) desde que seja dado M1.

A condição s′(x0) = f ′(x0) induz uma expressão para M1 que depende apenas dos dados do nossoproblema. De facto, considerando s(1)

′(x0), em que s(1)

′(x) é definida por (3.6.3), obtemos

M1 =− f ′(x0). (3.6.10)

Passemos agora a construção da solução do problema inicial: a determinação de um polinómiosegmentado cúbico, interpolador, contínuo e com derivada de primeira e segunda ordens contínuas nosnodos da partição.


Definição 3.5. Seja f uma função definida num intervalo [a,b] onde consideramos a partição a =x0 < x1 < · · ·< xn−1 < xn = b. O polinómio segmentado

s(x) =

s(1)(x), x ∈ [x0,x1],

s(2)(x), x ∈ ]x1,x2],...s( j)(x), x ∈ ]x j−1,x j],...s(n)(x), x ∈ ]xn−1,xn],

(3.6.11)

em que s( j) é um polinómio cúbico em ]x j−1,x j], que verifica

s(x j) = f (x j), j = 0, . . . ,n,s′′ é contínua em xi, i = 1, . . . ,n−1,s′ é contínua em xi, i = 1, . . . ,n−1,s é contínua em xi, i = 1, . . . ,n−1,

(3.6.12)

es′(x0) = f ′(x0)e s′(xn) = f ′(xn) (3.6.13)

ous′′(x0) = 0e s′′(xn) = 0 (3.6.14)

é designado spline interpolador cúbico.Se (3.6.13), então o spline s é designado spline completo. Por outro lado, se (3.6.14), o spline s é

designado spline natural. Neste caso as condições (3.6.14) surgem na literatura substituidas por

s′′(x0) = f ′′(x0)e s′′(xn) = f ′′(xn). (3.6.15)

Observamos que, por definição, a construção o spline interpolador cúbico requer a determinação de 4nparâmetros e este número é igual ao número de condições que são exigidas ao spline s.

Consideremos para s(i)′′

a seguinte expressão

s(i)′′(xi) = Mi

xi− xhi

+Mi+1x− xi−1

hi. (3.6.16)

A definição anterior induz de imediato a continuidade de s′′ em xi. De facto, tem-se sucessivamente

limx→x−i

s′′(x) = limx→x−i

s(i)′′(x)

= limx→x−i

Mixi− x

hi+Mi+1

x− xi−1

hi

= Mi+1

= limx→x+i

Mi+1xi+1− x

hi+1+Mi+2

x− xi

hi+1

= limx→x+i

s(i+1)′′(x)

= limx→x+i

s′′(x).


De (3.6.16) vem

s(i)(x) =−Mi(x− xi)

3

6hi+Mi+1

(x− xi−1)3

6hi−Ci

x− xi

hi+Di

x− xi−1

hi.

O spline cúbico é interpolador. Logo s(xi) = s(i)(xi) = f (xi). Assim

Di = f (xi)−Mi+1h2

i

6,

e portanto


3

6hi+Mi+1

(x− xi−1)3

6hi−Ci

x− xi

hi−(

Mi+1h2

i

6− f (xi)

)x− xi−1

hi.

O spline é uma função contínua nos nodos da partição. Logo

limx→x+i−1

s(x) = s(xi−1)

= f (xi−1)

= limx→x+i−1

s(i)(x)

=Ci +Mih2

i

6,

e portanto

Ci = f (xi−1)−Mih2

i

6.

Obtemos, deste modo, para i = 1, . . . ,n, a seguinte expressão


3

6hi+Mi+1

(x− xi−1)3

6hi

+

(f (xi−1)−Mi

h2i

6

)−x+ xi

hi+

(f (xi)−Mi+1

h2i

6

)x− xi−1

hi.

(3.6.17)

Para determinar completamente o spline interpolador cúbico, temos que determinar os n+1 coeficientesM1, . . . ,Mn+1. Consideramos agora a condição sobre s′(x). Esta derivada é contínua nos pontos da partição.Atendendo a que se tem

limx→x−i

s′(x) = limx→x−i

s(i)′(x)

= Mihi

6+Mi+1

hi

3+

f (xi)− f (xi−1)

hi

elim

x→x+is′(x) = lim

x→x+is(i+1)′(x)

=−Mi+1hi+1

3−Mi+2

hi+1

6+

f (xi+1)− f (xi)

hi+1

concluímos

Mihi

6+Mi+1

hi +hi+1

3+Mi+2

hi+1

6=

f (xi+1)− f (xi)

hi+1− f (xi)− f (xi−1)

hi, i = 1, . . . ,n−1. (3.6.18)

Temos que estabelecer as duas condições restantes. Consideramos no que segue as condições (3.6.13)para o spline completo.


De (3.6.17), para i = 0, N e s′(x0) = f ′(x0), s′(xn) = f ′(xn), vem

M1h1

3+M2

h1

6=− f ′(x0)+

f (x1)− f (x0)

h1(3.6.19)

e

Mnhn

6+Mn+1

hn

3= f ′(xn)−

f (xn)− f (xn−1)

hn. (3.6.20)

Conjugando (3.6.18), (3.6.19) e (3.6.20) obtemos o seguinte sistema linear (tridiagonal) AM = B emque M = [M1 M2 . . . ,Mn+1]

t , e

A =

h13

h16

h16

h1+h23

h26

. . . . . . . . .hn−1

6hn−1+hn

3hn6

hn6

hn3

,

B =

− f ′(x0)+f (x1)− f (x0)

h1f (x2)− f (x1)

h2− f (x1)− f (x0)

h1...

f (xn)− f (xn−1)hn

− f (xn−1)− f (xn−2)hn−1

f ′(x0)− f (xn)− f (xn−1)hn

.

Se a partição é uniforme, então

A = h

13

16 0 . . . 0 0 0

16

23

16 . . . 0 0 0

......

......

......

...0 0 0 . . . 1

623

16

0 0 0 . . . 0 16

13

.

Atendendo a que A é não singular pois λ = 0 não é valor próprio de A 1, concluímos que existe um e umsó spline cúbico natural que é interpolador de f em xi, i = 0, . . . ,n.

O teorema seguinte permite-nos concluir a unicidade do spline interpolador cúbico completo e natural.

Teorema 3.7. Seja f : [a,b]−→ R, f ∈C2[a,b] e, em [a,b], consideremos a partição

a = x0 < x1 < · · ·< xn = b.

Seja s o spline interpolador cúbico completo, então

∫ b

a(s′′(x))2dx≤

∫ b

a( f ′′(x))2 dx. (3.6.21)

1Teorema de Gershgorin: Se λ ∈ C é valor próprio de A = [ai j]_i, j = 1ninMn(R), então λ ∈ ∪ni=1BRi(aii), em que

Ri = ∑i 6= j|ai j|.


Demonstração. Comecemos por notar que se tem

∫ b

a( f ′′(x))2dx =

∫ b

a(s′′(x))2dx+

∫ b

a(e′′(x))2dx−2

∫ b

as′′(x)( f ′′(x)− s′′(x))dx,

em que

e(x) = f (x)− s(x).

Para o último termo do segundo membro da identidade anterior vale a representação

∫ b

as′′(x)( f ′′(x)− s′′(x))dx = ( f ′− s′)s′′|ba−

∫ b

as′′′(x)( f ′(x)− s′(x))dx,

e, atendendo a que s é completo, vem

∫ b

as′′(x)( f ′′(x)− s′′(x))dx =−

∫ b

as′′′(x)( f ′(x)− s′(x))dx.

Uma vez que s denota o spline interpolador cúbico, s′′′(x) = Ki em (xi,xi+1), e portanto

∫ b

as′′′(x)( f ′(x)− s′(x))dx =

n

∑j=0

∫ xi+1

xi

s′′′(x)( f ′(x)− s′(x))dx

=n

∑j=0

Ki

∫ xi+1

xi

( f ′(x)− s′(x))dx

=n

∑j=0

Ki( f (x)− s(x))|xi+1xi

= 0

pois s é interpolador de f nos nodos da partição. Logo

∫ b

a( f ′′(x))2dx =

∫ b

a(s′′(x))2dx+

∫ b

a(e′′(x))2dx.

A identidade anterior permite concluir (3.6.21).

Sejam s1 e s2 dois splines interpoladores cúbicos de f nos nodos xi, i = 0, . . . ,n. Então s1− s2 é umspline interpolador cúbico da função nula. Pelo Teorema 3.7 vem

∫ b

a((s1− s2)

′′(x))2 dx = 0,

e, atendendo a que a função integranda é contínua, obtemos

(s1− s2)′′(x) = 0, x ∈ [a,b].

Da igualdade anterior vem s1(x)− s2(x) = Ax+B,x ∈ [a,b] e, uma vez que s1 e s2 são interpoladores def , concluímos a igualdade dos splines interpoladores cúbicos completos anteriores.


O comportamento do erro do spline interpolador cúbico é estabelecido no teorema seguinte queapresentamos sem demonstração.

Teorema 3.8. Se f ∈C4[a,b], então

‖ f − s‖∞ ≤5

384

∥∥ f (4)∥∥

∞h4,

em que h = maxi=1,...,n hi, hi = xi− xi−1, i = 1, . . . ,n.

3.7 Interpolação de funções definidas em R2

Consideremos a superfície S definida por z = f (x,y) e sejam (xi,y j,zi j), i = 0, . . . ,n, j = 0, . . . ,m,pontos S . Pretendemos averiguar a existência de uma superfície polinomial

zp =n

∑i=0

m

∑j=0

ai jxiy j

tal quezp(xi,y j) = f (xi,y j), i = 0, . . . ,n, j = 0, . . . ,m.

Comecemos por definir o problema de interpolação bidimensional. Seja Ω = [a,b]× [c,d] e f : Ω−→R. Em Ω introduzimos a seguinte partição rectangular

(xi,y j) : i = 0, . . . ,n, j = 0, . . . ,m, x0 = a < x1 < · · ·< xn = b, y0 = c < y1 < · · ·< ym = d. (3.7.1)

Por Pn,m denotemos o seguinte polinómio, de grau menor ou igual a n em x e de grau menor ou igual am em y,

Pn,m(x,y) =n

∑i=0

m

∑j=0

ai jxiy j.

Pretendemos construir Pn,m que verifica

Pn,m(xi,y j) = f (xi,y j) := fi j, i = 0, . . . ,n, j = 0, . . . ,m. (3.7.2)

Consideremos as funções de Lagrange unidimensionais anteriormente introduzidas

ì(x) =n

∏i6= j=0

x− x j

xi− x j, i = 0, . . . ,n,

e

` j(y) =n

∏j 6=i=0

y− yi

y j− yi, j = 0, . . . ,m,

e sejaφi j(x,y) = ì(x)` j(y), i = 0, . . . ,n, j = 0, . . . ,m. (3.7.3)

No teorema seguinte estabelecemos a existência e unicidade do polinómio Pn,m que satisfaz (3.7.2).


Teorema 3.9. Seja f : Ω−→R e Ω = [a,b]× [c,d] onde consideramos a rede rectangular (3.7.1).O polinómio

Pn,m(x,y) =n

∑i=0

m

∑j=0

f (xi,y j)φi j(x,y),

é o único polinómio, de grau menor ou igual a n em x e de grau menor ou igual a m em y, quesatisfaz (3.7.2).

Demonstração. Atendendo à definição de φi j, é fácil mostrar que Pn,m satisfaz (3.7.2).Provemos a unicidade. Sejam

Qn,m(x,y) =n

∑i=0

m

∑j=0

ai jxiy j

e

Rn,m(x,y) =n

∑i=0

m

∑j=0

ai jxiy j,

dois polinómios que verificam (3.7.2). Então

E(x,y) = Qn,m(x,y)−Rn,m(x,y) =n

∑i=0

m

∑j=0

ci jxiy j, ci j = ai j−bi j,

é tal queE(xk,yp) = 0, k = 0, . . . ,n, p = 0, . . . ,m.

Mas, de

0 =n

∑i=0

m

∑j=0

ci jxiky j

p

=n

∑i=0

( m

∑j=0

ci jy jp)xi

k

=n

∑i=0

γixik, k = 0, . . . ,n γi =

m

∑j=0

ci jy jp

deduzimosγi, i = 0, . . . ,n,

e portanto, para i = 0, . . . ,n,0 = γi

=m

∑j=0

ci jy jp, p = 0, . . . ,m,

que permite concluir0 = ci j

= ai j−bi j

O polinómio Pn,m é designado polinómio interpolador de Lagrange de f nos nodos (xi,y j),i = 0, . . . ,n, j = 0, . . . ,m.


Exemplo 3.13. Seja f : [0,1]× [0,1] −→ R definida por f (x,y) = ex+y. O polinómio interpolador deLagrange nos pontos

y j/xi 0 1

0 1 e1 e e2

tem a seguinte representação

P1,1(x,y) = (1− x)+ ex)(1− y)+(e(1− x)+ e2x)y,(x,y) ∈ [0,1]× [0,1].

Na Figura 3.17 ilustramos o comportamento de f e do seu polinómio interpolador P1,1. O comporta-mento do erro | f −P1,1| é ilustrado na Figura 3.18.

0 0.2 0.4 0.6 0.8 1 0

0.5

1

5

0 0.2 0.4 0.6 0.8 1 0

0.5

1

5

Figura 3.17: Gráficos de f (x,y) = ex+y (esquerda) e do polinómio interpolador P1,1 (direita).

0 0.2 0.4 0.6 0.8 1 0

0.5

10

0.5

Figura 3.18: Gráfico de |ex+y−P1,1(x,y)|, (x,y) ∈ [0,1]× [0,1].

Estudamos no que segue o comportamento do erro f −Pn,m. Por Cp(Ω) denotamos o espaço dasfunções contínuas com derivadas parciais até à ordem p contínuas em Ω.


Teorema 3.10. Seja Pn,m o polinómio interpolador de f : Ω −→ R, Ω = [a,b]× [c,d], na rederectangular (3.7.1). Se f ∈Cn+m+2(Ω), então

| f (x,y)−Pn,m(x,y)| ≤1

(n+1)!

∥∥∥∥∂ n+1 f∂xn+1

∥∥∥∥∞

n

∏i=0|x− xi|

+1

(n+1)!(m+1)!

∥∥∥∥∂ n+m+2 f

∂xn+1∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|

n

∏i=0|x− xi|

+1

(m+1)!

∥∥∥∥∂ m+1 f∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|.

(3.7.4)

Demonstração. O erro f −Pn,m admite a seguinte decomposição

f (x,y)−Pn,m(x,y) = E1(x,y)+E2(x,y), (3.7.5)

em que

E1(x,y) = f (x,y)−n

∑i=0

f (xi,y)ì(x) e E2(x,y) =n

∑i=0

(f (xi,y)−

m

∑j=0

fi j` j(y))ì(x).

Para E1 tem-se o seguinte majorante

|E1(x,y)| ≤1

(n+1)!

∥∥∥∥∂ n+1 f∂xn+1

∥∥∥∥∞

n

∏i=0|x− xi|. (3.7.6)

No que diz respeito a E2, observamos que se tem

E2(x,y) =n

∑i=0

R(xi,y)ì(x),

em que

R(x,y) = f (x,y)−m

∑j=0

f (x,y j)` j(y).

Com esta identificação, temos para E2 a decomposição

E2(x,y) =n

∑i=0

R(xi,y)ì(x)−R(x,y)+R(x,y),

em que ∣∣∣∣n

∑i=0

R(xi,y)ì(x)−R(x,y)∣∣∣∣≤

1(n+1)!

∥∥∥∥∂ n+1R∂xn+1

∥∥∥∥∞

n

∏i=0|x− xi|

e∂ n+1R∂xn+1 =

∂ n+1

∂xm+1

(f (x,y)−

m

∑j=0

f (x,y j)` j(y))=

1m+1

∂ n+m+2 f∂xn+1∂ym+1 (x,η)

m

∏j=0

(y− y j),

para algum η ∈ (y0,ym), e portanto∣∣∣∣

n

∑i=0

R(xi,y)ì(x)−R(x,y)∣∣∣∣≤

1(n+1)!(m+1)!

∥∥∥∥∂ n+m+2 f

∂xn+1∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|

n

∏i=0|x− xi|


Para R tem-se a seguinte estimativa

|R(x,y)| ≤ 1(m+1)!

∥∥∥∥∂ m+1 f∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|.

Logo para E2 é válida a seguinte estimativa

|E2(x,y)| ≤1

(n+1)!(m+1)!

∥∥∥∥∂ n+m+2 f

∂xn+1∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|

n

∏i=0|x− xi|+

1(m+1)!

∥∥∥∥∂ m+1 f∂ym+1

∥∥∥∥∞

m

∏j=0|y− y j|.

(3.7.7)Finalmente, conjugando (3.7.5) com (3.7.6) e (3.7.7), obtemos (3.7.4).

Como corolário do Teorema 3.10 é válido o seguinte resultado.

Corolário 3.3. Seja Pn,m o polinómio interpolador de f : Ω −→ R, Ω = [a,b]× [c,d], na rederectangular (3.7.1). Se f ∈Cn+m+2(Ω), então

‖ f −Pn,m‖∞ ≤1

4(n+1)

∥∥∥∥∂ n+1 f∂xn+1

∥∥∥∥∞

hn+1x

+1

16(n+1)(m+1)

∥∥∥∥∂ n+m+2 f

∂xn+1∂ym+1

∥∥∥∥∞

hn+1x hm+1

y

+1

4(m+1)

∥∥∥∥∂ m+1 f∂ym+1

∥∥∥∥∞

hm+1y ,

(3.7.8)

em que hx = maxi=1,...,n xi− xi−1, hy = max j=1,...,m y j− y j−1.

O conceito de polinómio interpolador segmentado pode também ser facilmente introduzido para ocaso bidimensional.

Definição 3.6. Seja f : Ω−→R, Ω = [a,b]× [c,d] e neste conjunto consideremos a rede rectangular(3.7.1) e em que n = pxqx, m = pyqy. O polinómio s definido por

s(x,y) =

s(1,1)(x,y), (x,y) ∈ [x0,xpx ]× [y0,ypy ],

s(2,1)(x,y), (x,y) ∈ ]xpx ,x2px ]× [y0,ypy ],...

s(qx,1)(x,y), (x,y) ∈ ]xn−px ,xn]× [y0,ypy ],...

s(qx,qy)(x,y), (x,y) ∈ ]xn−px ,xn]× [ym−py ,xm],

em que s(i, j)(x,y) é o polinómio interpolador de Lagrange, de grau px em x e py em y, de f nos nodos

x(i−1)px < x(i−1)px+1 < · · ·< xipx , y( j−1)py < y( j−1)py+1 < · · ·< y jpy ,

é chamado polinómio interpolador de Lagrande segmentado, de grau px em x e py em y, de f na rederectangular (3.7.1).


Exemplo 3.14. Seja f : [0,1]× [0,1] −→ R definida por f (x,y) = ex+y. O polinómio interpolador deLagrange segmentado linear nos pontos

y j/xi 0 12 1

0 1 e12 e

12 e

12 e e

32

1 e e32 e2

tem a seguinte representação

s(x,y) =

4((

x− 12

)− e

12 x)(

y− 12

)

+4(−e

12(x− 1

2

)+ ex

)y, (x,y) ∈ [0, 1

2 ]× [0, 12 ](

−2e12(x−1

)+2e

(x− 1

2

))(−2)

(y− 1

2

)

+(−2e(x−1)+2e

32(x− 1

2

))2y, (x,y) ∈ ]1

2 ,1]× [0, 12 ]

4(e

12(x− 1

2

)− ex

)(y−1)

+4(−e(x− 1

2

)+ e

32 x)(

y− 12

), (x,y) ∈ [0, 1

2 ]× ]12 ,1](

−2e(x−1)+2e32(x− 1

2

))(−2)(y−1)

+(−2e

32 (x−1)+2e2

(x− 1

2

))2(y− 1

2

), (x,y) ∈ ]1

2 ,1]× ]12 ,1]

(3.7.9)

Na Figura 3.19 ilustramos o comportamento do erro f (x,y)− s(x,y) para (x,y) ∈ [0, 12 ]× [0, 1

2 ].

00.2

0.4 00.2

0.40

5 ·10−2

0.1

Figura 3.19: Gráfico de |ex+y− s(x,y)|,(x,y) ∈ [0, 12 ]× [0, 1

2 ].

Das Figuras 3.18 e 3.19 concluímos que há uma redução significativa no erro ao substituirmos opolinómio interpolador P1,1 por s definido por (3.7.9) na região [0, 1

2 ]× [0, 12 ]

O comportamento do erros do polinómio interpolador segmentado pode ser estudado a partir doTeorema 3.10. De facto, se (x,y) ∈ ]x(i−1)px ,xipx ]× ]y( j−1)py ,y jpy ], então

| f (x,y)− s(x,y)|= | f (x,y)− s(i, j)(x,y)|


e, atendendo a que se tem

| f (x,y)− s(i, j)(x,y)| ≤ 1(px +1)!

∥∥∥∥∂ px+1 f∂xpx+1

∥∥∥∥∞

ipx

∏`=(i−1)px

|x− x`|

+1

(px +1)!(py +1)!

∥∥∥∥∂ px+py+2 f

∂xpx+1∂ypy+1

∥∥∥∥∞

jpy

∏k=( j−1)py

|y− yk|ipx

∏`=(i−1)px

|x− x`|

+1

(py +1)!

∥∥∥∥∂ py+1 f∂ypy+1

∥∥∥∥∞

jpy

∏k=( j−1)py

|y− yk|,

concluímos

| f (x,y)− s(x,y)| ≤ 1(px +1)!

∥∥∥∥∂ px+1 f∂xpx+1

∥∥∥∥∞

ipx

∏`=(i−1)px

|x− x`|

+1

(px +1)!(py +1)!

∥∥∥∥∂ px+py+2 f

∂xpx+1∂ypy+1

∥∥∥∥∞

jpy

∏k=( j−1)py

|y− yk|ipx

∏`=(i−1)px

|x− x`|

+1

(py +1)!

∥∥∥∥∂ py+1 f∂ypy+1

∥∥∥∥∞

jpy

∏k=( j−1)py

|y− yk|.

Mais ainda

‖ f − s‖∞ ≤1

4(px +1)

∥∥∥∥∂ px+1 f∂xpx+1

∥∥∥∥∞

hpx+1x

+1

14(px +1)(py +1)

∥∥∥∥∂ px+py+2 f

∂xpx+1∂ypy+1

∥∥∥∥∞

hpx+1x hpy+1

y

+1

4(py +1)

∥∥∥∥∂ py+1 f∂ypy+1

∥∥∥∥∞

hpy+1y .

A finalizar esta breve digressão pela interpolação bidimensional, realçamos o facto de existirem parafunções de duas variáveis, extensões naturais dos conceitos de polinómio de Hermite e de spline ([4]).

É comum, em problemas de aproximação de funções, considerar a questão da determinação de umpolinómio interpolador a partir do conhecimento de uma função em pontos de uma rede não rectangular,isto, é, num conjunto de pontos que não estão associados a uma partição rectangular do domínio. Emmuitas aplicações é considerada uma partição do domínio em triângulos sendo os nodos de interpolaçãoos vértices destes conjuntos. É manifesto que a partição em triângulos não poderá ser qualquer.

Definição 3.7. Seja Ω um domínio poligonal de R2. Uma partição em triângulos TH = Ti, i ∈ Iem que Ti denota o triângulo aberto, I denota um conjunto de indices, diz-se admissível para Ω se sãoválidas as condições seguintes:

1. se i, j ∈ I, i 6= j, então Ti∩Tj = /0,

2. se i, j ∈ I, i 6= j, então T i ∩ T j = /0 ou T i ∩ T j = o vértice comum aos dois triângulos ouT i∩T j = lado comum aos dois triângulos,

3. ∪i∈IT i = Ω.

Consideremos então uma triângulação TH admissível para Ω, denotemos por VH = (xi,y j), i ∈I , j ∈J o conjunto de vértices dos triângulos desta triangulação e suponhamos que são conhecidos


f (xi,y j), i ∈I , j ∈J . A construção de um polinómio interpolador s tal que, em cada triângulo T i devertices (xi,y j),(xp,yq),(xr,yt), é um polinómio de grau um em cada uma das variáveis, resume-se àdeterminação dos coeficientes a,b,c tais que

s(x,y) = a+bx+ y,(x,y) ∈ T i.

Uma vez que se pretende que s seja interpolador, os coeficientes anteriores deverão verificar

a+bxi + y j = f (xi,y j)

a+bxp + yq = f (xp,yq)

a+bxr + yt = f (xr,yt).

(3.7.10)

O sistema (3.7.10) é considerado para cada triângulo T i obtendo-se, deste modo, uma função po-linomial segmentada, de grau um em cada variável, interpoladora de f nos vértices dos triângulos deTH .


1. Na seguinte tabela são dados diferentes valores para o peso específico p da água a diferentestemperaturas t (em graus Celsius):

t 0 1 2 3

p 0.999871 0.999928 0.999969 0.999991

Usando interpolação linear, quadrática e cúbica, determine uma aproximação para p quando t = 4oCusando a fórmula interpoladora de Lagrange e de Newton. Compare os resultados obtidos sabendoque o valor exacto 1.000000.

2. Pretende-se construir uma tabela para a função f (x) = ex, com x ∈ [0,1]. Considere o valor dee com erro absoluto inferior a 0.5× 10−6 e uma partição com pontos igualmente distanciados.Determine o diâmetro da partição a considerar de modo que o polinómio interpolador de Lagrangepermita obter uma aproximação para f com um erro inferior a 10−6.

3. Considere a função f (x) = ln(x+1), x ∈ [1,3]. Determine o polinómio interpolador de Lagrangeque aproxima f em [1,3] com um erro inferior a 10−2.

4. Determine o polinómio interpolador de Lagrange segmentado linear para a função f tal que

xi 0 0.25 0.5 0.75 1

p 1 1.25 1 1.5 1

5. Calcule o polinómio interpolador de Lagrange segmentado quadrático para a função do problemaanterior.

6. Determine, utilizando interpolação segmentada linear com dois subintervalos para f (x) = cosx,x ∈ [0,π], uma aproxime cos(20o) e indique um majorante para o erro cometido.


7. Durante a sedimentação da reacção de saponificação entre quantidades equimolares de sódio eacetato de etilo, a concentração c (em gr/`) de cada reagente varia, com o tempo t (em min), deacordo com a equação

1c=

1c0

+ kt,

onde c0 é a concentração inicial e k (gr/(` ·min)) é a constante de reacção. Foram obtidos osseguinte resultados

1/c 24.7 32.4 38.4 45.0 52.3 65.6 87.6 102 154 192

t 1 2 3 4 5 7 10 12 20 25

(a) Obtenha uma estimativa para a concentração inicial.

(b) Obtenha uma estimativa para a concentração ao fim de 15min.

8. Estabeleça, utilizando o polinómio interpolador de Lagrange linear, a seguinte fórmula de integração

∫ b

af (x)dx' b−a

2(

f (a)+ f (b)),

que é designada fórmula dos trapézios. Suponha que | f ′′(x)| ≤M,x ∈ [a,b]. Mostre que∣∣∣∣∫ b

af (x)dx− b−a

2(

f (a)+ f (b))∣∣∣∣≤

M12

(b−a)3.

9. Estabeleça, utilizando o polinómio interpolador de Lagrange linear em [a,b], a seguinte fórmula dederivação numérica

f ′(a)' 1b−a

(f (b)− f (a)

),

que é designada fórmula progressiva. Mostre que

f ′(a)− 1b−a

(f (b)− f (a)

)=−b−a

2f ′′(η), η ∈ (a,b).

10. Estabeleça, utilizando o polinómio interpolador de Lagrange quadrático em [a,b], para os pontosa = x0 < x1 < x2 = b com h = x1−a = b− x1, a seguintes fórmula de derivação numérica

f ′(x1)'12h

(f (b)− f (a)

),

que é designada fórmula centrada. Mostre que

f ′(x1)−12h

(f (b)− f (a)

)=−h2

2f ′′(η), η ∈ (a,b).

11. Mostre que as funções polinomiais definidas por (3.5.3) são linearmente independentes.

12. Mostre que o conjunto ψi, i = 1, . . . ,4, em que ψi, i = 1, . . . ,4, são definidas por (3.5.6), é umabase de P4.

13. Considere a função f (x) = 3xex− e2x. Determine uma aproximação para f (1.03) usando o polinó-mio interpolador de Hermite para os pontos x0 = 0 e x1 = 1.05. Determine uma estimativa para oerro cometido.


14. Considere f (x) = ex. Determine o polinómio interpolador de Hermite de grau 5, H5, usando ospontos x0 = 0,x1 = 1 e x2 = 2. Compare H5(0.25) com f (0.25) e com P2(0.25), em que P2 é opolinómio interpolador de Lagrange de grau 2, para a função f , nos mesmos pontos.

15. Determine o polinómio de grau mínimo que faça a concordância entre a recta

y =−2+12(8− x),

no ponto (8,−2), e a circunferência

(x−1)2 +(y+2)2 = 1,

no ponto (1,−1). Nota: Duas curvas dizem-se concordantes de tiverem a mesma tangente no pontode união.

16. Determine o polinómio interpolador de uma função f para o suporte

f (0) =−1, f ′(0) =−2, f (1) = 0, f ′(1) = 10, f ′′(1) = 40.

17. Considere f (x) = ln(1+ x), x ∈ [0,1]. Determine o polinómio de Hermite H, de menor grau, talque

‖ f −H‖∞ < 2×10−3.

18. Considere f (x) = e−x2, x ∈ [−1,1]. Determine o polinómio interpolador de Hermite segmentado s

para f de menor grau. Calcule s(12) e estime o erro |e−1/4− s(1

2)|.

19. Diga se a função seguinte é um spline ou não

f (x) =

x2 +1, x ∈ [0,1],2x, x ∈ ]1,2],x2−2x+4, x ∈ ]2,3].

20. Diga se existem parâmetros a,b,c, tais que

f (x) =

x2, x ∈ [−1,0],ax3 +bx2 + cx, x ∈ ]0,1],2− x, x ∈ ]1,2],

é um spline cúbico.

21. Estabeleça o sistema que permite calcular o spline interpolador cúbico natural.

22. Calcule o spline cúbico natural para os seguintes dados

xi 0 1 2

f 1 2 1f ′ 0 1 0

23. Considere a seguinte partição uniforme do intervalo [a,b], a = x0 < x1 < · · · < xn = b, comh = xi− xi−1, e os pontos auxiliares

x−3 = a−3h, x−2 = a−2h, x−1 = a−h, xn+1 = b+h, xn+2 = b+2h, xn+3 = b+3h.


Para i =−1, . . . ,n+1, denotemos por Bi(x) o seguinte polinómio segmentado

Bi(x) =1h3

(x− xi−2)3, x ∈ [xi−2,xi−1],

h3 +3h2(x− xi−1)+3h(x− xi−1)2−3(x− xi−1)

3, x ∈ ]xi−1,xi],

h3 +3h2(xi+1− x)+3h(xi+1− x)2−3(xi+1− x)3, x ∈ ]xi,xi+1],

(xi+2− x)3, x ∈ ]xi+1,xi+2],0 nos restantes casos,

que é designado spline-B.

(a) Mostre que Bi ∈C2[a,b] e

· · · xi−2 xi−1 xi xi+1 xi+2 · · ·

Bi 0 0 1 4 1 0 0B′i 0 0 3

h 0 −3h 0 0

B′′i 0 0 6h2 −12

h26h2 0 0

(b) Seja s(x) = ∑n+1i=−1CiBi(x). Mostre que se s é o spline interpolador cúbico completo, então

C = [C−1 . . .Cn+1]t verifica

BC = F,

em que F = [ f ′(x0) f (x0) f (x1) . . . f (xn) f ′(xn)]t e

B =

−3h 0 3

h . . . 0 0 01 4 1 . . . 0 0 0...

......

......

......

0 0 0 . . . 4 1 00 0 0 . . . 1 4 10 0 0 . . . −3

h 0 3h

24. Considere f (x,y) = ln(1+ x+ y),(x,y) ∈ [0,1]× [0,1].

(a) Determine o polinómio interpolador P linear em cada uma das variáveis.

(b) Calcule P(14 ,

14) e determine | f (1

4 ,14)−P(1

4 ,14)|.

(c) Construa o polinómio interpolador segmentado linear s considerando a partição 0, 12 ,1.

(d) Calcule s(14 ,

14) e determine | f (1

4 ,14)− s(1

4 ,14)|.

(e) Construa um polinómio interpolador Pn,m para f tal que

‖ f −Pn,m‖∞ ≤ 10−2.

25. Considere o triângulo Ω e a triangulação T representados na Figura 3.20 e seja f : Ω−→ R2 talque

f (0,0) = 1, f(1

2,0)= f(

0,12

)=

12, f (1,0) = f (0,1) = f

(12,12

)= 0.

Construa o polinómio interpolador segmentado s para f nos pontos dados, de grau um em cada umadas variáveis.


(0,1)

(0,0.5)

(0,0)

(0.5,0.5)

(1,0)(0.5,0)

Figura 3.20: Triangulação do triângulo rectângulo de vértices (0,0),(1,0) e (0,1).

CAPÍTULO 4EQUAÇÕES NÃO LINEARES

4.1 Introdução

A determinação de zeros de equações lineares matriciais foi considerado no Capítulo 2. Neste casoconsiderámos métodos de duas classes: directos e iterativos. Na última classe de métodos obtivémos umasucessão de vectores de aproximações para o vector solução da equação matricial. No presente capítulo onosso objectivo é o estudo de métodos para a determinação da solução de equações não lineares.

Já nos deparámos com equações não lineares polinomiais do tipo

a0 +a1x+ . . .anxn = 0, an 6= 0.

Para estas equações não lineares particulares existem métodos com características próprias. No entantoos métodos que consideraremos neste capítulos poderão ser utilizados neste contexto. No entanto a suaaplicação é mais geral. Observamos que para algumas equações polinomiais existem fómulas explícitaspara os seus zeros.

Consideremos a equação f (x) = 0, f : [a,b]−→ R, em que [a,b] é tal que existe apenas uma raiz x∗

desta equação em [a,b]. Observamos que se f ′ tem sempre o mesmo sinal [a,b] e f (a) f (b)< 0, entãoexiste de facto apenas um zero de f em [a,b].

O nosso objectivo é estudar métodos que permitem construir uma sucessão de aproximações xn, n∈N,para x∗ e tal que

limn→∞

xn = x∗.

Assim sendo, os métodos que iremos estudar pertencem à classe dos métodos iterativos e a definição dasucessão de aproximações requer uma aproximação inicial. Mais ainda, como veremos, a aplicação dosdiversos métodos requer a definição de um intervalo contendo uma só raiz da equação que se pretenderesolver. Este é determinado a partir da representação gráfica da função f determinando um intervalocontendo a abcissa do ponto de intersecção da curva y = f (x) com o eixo das abcissas. Este métodode localização é usualmente designado método gráfico. Observamos que localizar as raízes de uma

113


determinada equação consiste em determinar intervalos contendo uma raiz da equação em estudo e, emseguida, provar analiticamente que de facto em cada intervalo indicado existe uma e uma só raiz.

Consideremos a equação

sin(x)x− x = 0. (4.1.1)

Na Figura 4.1 apresentamos o gráfico de f . Assim, a raiz da equação anterior está no intervalo [π8 ,

π2 ]. De

facto f é contínua neste intervalo e f (π8 ) f (π

2 )< 0. A unicidade resulta do facto de f ′(x)< 0 no intervaloanterior.

0 1 2 3 4 5 6

−6

−4

−2

0

y = sin(x)x − x

y = 0

Figura 4.1: Gráfico de f (x) = sin(x)x − x.

Observamos que caso não seja fácil a determinação gráfica de y = f (x), é por vezes possível reescrevera equação f (x) = 0 na forma equivalente

f1(x) = f2(x),

em que f1 e f2 têm uma representação gráfica mais simples. Assim o zero de f verifica f1(x∗) = f2(x∗), eportanto é a abcissa do ponto de inetrsecção das duas curvas

y = f1(x), y = f2(x).

Consideremos a equação (4.1.1). Esta equação admite a representação equivalente

sin(x) = x2.

Consideremos as duas curvas

y = sin(x), y = x2

que se intersectam apenas num ponto para x ∈ [π4 ,

π2 ]. Na Figura 4.2 apresentamos as duas curvas anterio-

res.


0 1 2 3 4 5 6

−1

−0.5

0

0.5

1 y = sin(x)

y = x2

Figura 4.2: Gráficos de f1(x) = sin(x) e de f2(x) = x2.

Identificado o intervalo contendo apenas a raiz da equação f (x) = 0 que pretendemos aproximar,podemos considerar um dos métodos que estudamos nas secções seguintes. Por exemplo, consideremos asucessão

xn+1 = xn−sin(xn)− x2

n

cos(xn)− sin(xn)xn− xn

, n = 0,1, . . . ,

em que x0 =π2 . Na tabela seguinte apresentamos alguns elementos da sucessão anterior.

xi f (xi)

1.5707963267948966192 −0.93417655442730.9060367009005804134 −0.03734662835940.8768125208725520179 −0.00010964231330.8767262161570021592 −9.67958207316 ·10−10

0.8767262153950624460 −7.54473754399 ·10−20

Um conceito relevante que determina a escolha do método iterativo é o conceito de ordem.

Definição 4.1. Consideremos f : [a,b]−→ R com uma único zero x∗ em [a,b]. O método iterativoxn+1 = g(xn), n ∈ N0, diz-se com ordem p ∈ R+, se existe n0 ∈ N0, tal que en = x∗− xn, paran ∈ N0, n≥ n0, verifica

|en+1| ≤C|en|p, (4.1.2)

onde C denota uma constante positiva, independente de n, e p é a maior constante positiva com apropriedade (4.1.2).

4.2 O método da bissecção

Consideremos o intervalo I0 = [a,b] onde f tem apenas um zero x∗. Sem perda de generalidade,suponhamos que f (a)> 0 e f (b)< 0. Consideremos a seguinte sequência de passos:

1. Seja

x1 =a+b

2.


Se f (x1) < 0, então x∗ ∈ I1 = [a,x1] = [a1,b1]. Caso contrário, x∗ ∈ I1 = [x1,b] = [a1,b1]. Emqualquer um dos casos, |I1|= 1

2 |I0|.2. Seja

x2 =a1 +b1

2.

Se f (x2) < 0, então x∗ ∈ I2 = [a1,x2] = [a2,b2]. Caso contrário x∗ ∈ I2 = [x2,b1] = [a2,b2]. Emqualquer um dos caso |I2|= 1

2 |I1|.3. Seja

x3 =a2 +b2

2.

4. O processo repete-se.

a bx1

I0

I1

I2

Figura 4.3: Ilustração gráfica da sucessão (xn)n∈N e dos intervalos In associados ao método da bissecção.

Obtemos deste modo a sucessão (xn)n∈N dos pontos médios dos intervalos, obtidos por bissecção, eque constituem aproximações para x∗. O processo descrito anteriormente é usualmente designado métododo ponto médio.

Teorema 4.1. Seja f : [a,b] −→ R tal que f é contínua e monótona em [a,b] e f (a) f (b) < 0.Seja (xn)n∈N a sucessão definida seguidamente:

1. I0 = [a0,b0] = [a,b];

2. Se In−1 = [an−1,bn−1] é tal que f (an−1) f (bn−1)< 0, então xn =an−1+bn−1

2 .

Então|In|=

12n |I0| (4.2.1)

elimn→∞

xn = x∗. (4.2.2)

Demonstração. Se f (an−1) f (xn)< 0, então x∗ ∈ In = [an−1,xn] = [an,bn] e x∗ ∈ In = [xn,bn−1] = [an,bn]caso contrário. Assim

|In|=12|In−1|=

12n |I0|.


Notemos que

|xn− x∗| ≤ |In|=12n |I0|

e consequentemente vale (4.2.2).

Se pretendermos determinar uma aproximação, utilizando o método do ponto médio, para o únicozero de f que existe em [a,b] com erro inferior a uma tolerância ε , então o número de iterações n deveverificar

|xn− x∗|< ε.

Atendendo a que se tem

|xn− x∗| ≤ 12n |I0|,

é suficiente tomar n a verificar

n >ln( |I0|

ε)

ln2. (4.2.3)

O algoritmo anterior traduz-se na seguinte rotina de MATLAB

Algoritmo 4.1 Método da bissecção com controlo de número máximo de iterações e erroabsoluto

function [zero,res,niter,erro]=bisseccao(f,a,b,tol,nmax)x=[a (a+b)/2 b]; fx = feval(f,x); niter=0; I=(b−a)/2;

while I≥tol && niter≤nmaxniter=niter+1;if fx(1)*fx(2)<0

x(3)=x(2); x(2)=x(1)+(x(3)−x(1))*0.5;fx=feval(f,x); I=(x(3)−x(1))*0.5;

elseif fx(2)*fx(3)<0x(1)=x(2); x(2)=x(1)+(x(3)−x(1))*0.5;fx = feval(f,x); I=(x(3)−x(1))*0.5;

elsex(2)=x(find(fx==0)); I=0;

endend

zero=x(2); x=x(2); res = feval(f,x); erro=I;return

Observamos que na função anterior é necessária uma implementação inline da função f.

Exemplo 4.1. Consideremos a equação

xcos(x) = 2sin(x).

Esta equação é equivalente atan(x) =

x2.

Considerando as abcissas dos pontos de intersecção das curvas y = tan(x) e 2y = x obtemos as raízes daequação dada. Na Figura 4.4 apresentamos as curvas anteriores.


−3π2

−π2 0 π

23π2

−6

−4

−2

0

2

4

6

Figura 4.4: Gráficos de f1(x) = tan(x) (traço contínuo) e de f2(x) = x2 (traço descontínuo).

Provemos que, em [π, 75 π], existe apenas uma raiz da equação que pretendemos resolver. Comecemos

por notar que se

f (x) = tan(x)− x2,

então f (π) < 0, f (7π5 ) > 0 e f é contínua no intervalo anterior. Uma vez que f ′ > 0 concluímos que

existe de facto apenas um zero de f neste intervalo.A sucessão de aproximações definidas pelo método da bissecção é definida seguidamente:

iteração sinal intervalo

x1 =6π5 − x∗ ∈ [6π

5 , 7π5 ] = [a1,b1]

x2 =13π10 − x∗ ∈ [13π

10 , 7π5 ] = [a2,b2]

x3 =27π20 − x∗ ∈ [27π

20 , 7π5 ] = [a3,b3]

x4 =11π

8 + x∗ ∈ [27π20 , 11π

8 ] = [a4,b4]

x5 =109π

80 + x∗ ∈ [27π20 , 109π

80 ] = [a5,b5]

x6 =217π160 − x∗ ∈ [217π

160 , 109π80 ] = [a6,b6]

Notamos que x7 =8764

π e f (x7)'−0.02097889.

4.3 O método de Newton

Seja f : [a,b]−→ R tal que

1. f (a) f (b)< 0 ( f (a)< 0 e f (b)> 0 sem perda de generalidade),

2. f ′ tem sinal constante em (a,b) ( f ′ > 0 sem perda de generalidade),

3. f ′′ tem sinal constante em (a,b) ( f ′′ > 0 sem perda de generalidade).


x0 x1x2

Figura 4.5: Ilustração gráfica da construção da sucessão (xn)n∈N do método de Newton.

Consideremos a sucessão (xn)n∈N definida da seguinte forma (ver Figura 4.5):

1. Dado (x0,y0), x0 = b, y0 = f (b), a recta tangente à curva y = f (x) que tem a seguinte equação

y− y0 = f ′(x0)(x− x0).

Se x1 denota a abcissa do ponto de intersecção da recta anterior com o eixo das abcissas, então

x1 = x0−f (x0)

f ′(x0).

2. Consideremos a recta tangente à curva y = f (x), no ponto (x1, f (x1)), de equação

y− f (x1) = f ′(x1)(x− x1).


x2 = x1−f (x1)

f ′(x1).

3. O processo repete-se...

4. Consideremos a recta tangente à curva y = f (x) no ponto (xn, f (xn)) de equação

y− f (xn) = f ′(xn)(x− xn).

Se xn+1 denota a abcissa do ponto de intersecção da recta anterior com o eixo das abcissas, então

xn+1 = xn−f (xn)

f ′(xn). (4.3.1)

Obtivémos, deste modo, uma sucessão (xn)n∈N definida recursivamente por (4.3.1). O método descrito édesignado método de Newton ou método de Newton-Raphson.


A função seguinte implementa o método de Newton. Notemos que além da implementaçãoinline de f, é necessário fornecer uma implementação do mesmo tipo para a derivada f'.

Algoritmo 4.2 Método de Newton com controlo de número máximo de iterações e erro absoluto

function [zero,res,niter,erro]=newton(f,df,x0,tol,nmax)x=x0;fx=feval(f,x); dfx=feval(df,x);niter=0; er=tol*10;

while er≥tol && niter≤nmaxniter=niter+1;diff=−fx/dfx;xold=x; x=x+diff;diff=abs (diff);fx=feval(f,x); dfx=feval(df,x);er=abs(x−xold);

endzero=x; res=fx; erro=er;return

Exemplo 4.2. Consideremos a equação

ln(1+ x)− x+1 = 0.

O método gráfico permite facilmente localizar as raízes da equação anterior reescrevendo-a na formaequivalente

ln(1+ x) = x−1.

Na Figura 4.6 ilustramos o comportamento das funções f1(x) = ln(1+ x) e f2(x) = x−1.

−1 0 1 2 3 4

−2

0

2

y = ln(1+ x)y = x−1

Figura 4.6: Gráficos de f1(x) = ln(1+ x) e de f2(x) = x−1.

Provemos que as raízes da equação dada estão nos intervalos [− 910 ,− 1

10 ] e [32 ,

52 ]. Observamos que

a função f (x) = ln(1+ x)− x+1 é contínua nos intervalos indicados. No que diz respeito ao primeirointervalo, temos


1. f(− 9

10

)=−0.402585, f

(− 1

10

)= 0.994639,

2. f ′(x) = 11+x −1 > 0, x ∈

(− 9

10 ,− 110

).

Logo existe uma raiz única em(− 9

10 ,− 110

).

No que diz respeito ao segundo intervalo, temos

1. f(3

2

)' 0.416291, f

(52

)'−0.247237,

2. f ′(x) = 11+x −1 < 0, x ∈

(32 ,

52

).

Logo existe uma raiz única em [32 ,

52 ].

Na tabela seguinte apresentamos alguns elementos da sucessão anterior com x0 =− 910 .

xi f (xi)

−0.85526832300066159066 −0.0776054313105−0.84213564414917271192 −0.00388347780766−0.84140765842997510135 −0.0000106002220641−0.84140566045191836530 −7.93566566427 ·10−11

−0.84140566043696063785 −4.44760152368 ·10−21

Estudamos seguidamente as propriedades da sucessão definida pelo método de Newton.

Teorema 4.2. Seja f : [a,b]−→ R, f ∈C2[a,b], tal que

1. f (a) f (b)< 0,

2. f ′ tem sinal constante em (a,b),

3. f ′′ tem sinal constante em (a,b).

Se (xn)n∈N é a sucessão definida por (4.3.1) com x0 tal que

f (x0) f ′′(x0)< 0, (4.3.2)

então (xn)n∈N converge para a única raiz x∗, em (a,b), da equação f (x) = 0, e o erro en = x∗−xn

verifica

|en+1| ≤M2m|en|2, (4.3.3)

onde0 < m≤ | f ′(x)|, | f ′′(x)| ≤M, x ∈ [a,b],

Demonstração. Suponhamos, sem perda de generalidade que, em (a,b), f ′ > 0, f ′′ > 0 e seja x0 = b.Provemos que (xn)n∈N é uma sucessão limitada inferiormente e monótona decrescente.

1. É claro que x∗ ≤ x0. Suponhamos que x∗ ≤ xn e provemos que x∗ ≤ xn+1. Notemos que se temsucessivamente

0 = f (x∗)

= f (xn +(x∗− xn))

= f (xn)+ f ′(xn)(x∗− xn)+f ′′(ηn)

2(x∗− xn)

2,


em que ηn ∈ (x∗,xn) pois x∗ ≤ xn. É fácil agora provar a igualdade

x∗ = xn−f (xn)

f ′(xn)+

f ′′(ηn)

2 f ′(xn)(x∗− xn)

2, (4.3.4)

e, atendendo à definição (4.3.1), concluímos que (xn)n∈N verifica

x∗− xn+1 =−f ′′(ηn)

2 f ′(xn)(x∗− xn)

2.

Uma vez que f ′ > 0, f ′′ > 0, estabelecemos

x∗− xn+1 ≤ 0, ∀n ∈ N0.

Logo a sucessão (xn)n∈N é limitada.

2. Observamos agora que (xn)n∈N é monótona decrescente. De facto, x1 ≤ x0. Suponhamos quexn ≤ xn−1. De (4.3.1), tem-se sucessivamente

f (xn) = f (xn−1 +(xn− xn−1))

= f (xn−1)+ f ′(xn−1)(xn− xn−1)+f ′′(ηn)

2(xn− xn−1)

2

=f ′′(ηn)

2(xn− xn−1)

2,

com ηn ∈ (xn,xn−1), e portanto f (xn)≥ 0. Assim, uma vez que

xn+1 = xn−f (xn)

f ′(xn),

vem xn+1 ≤ xn.

Provámos que (xn)n∈N é uma sucessão limitada inferiormente e monótona decrescente. Logo éconvergente. Seja x o limite desta sucessão. Atendendo a que f e f ′ são contínuas, de (4.3.1), vem

f (x) = 0,

e, atendendo à unicidade do zero de f em (a,b), concluímos que x = x∗.A estimativa (4.3.3) para o erro deduz-se facilmente de (4.3.4).

Observamos que o Teorema 4.2 é uma condição suficiente para a convergência do método de Newton.Atendendo a este facto, a sucessão gerada pelo Newton pode convergir não sendo válidas algumas dassuas hipóteses. Por outro lado, a estimativa (4.3.3) estabelece que o método de Newton apresenta ordemde convergência igual a 2 e a velocidade da convergência depende do comportamento de f ′ e f ′′.

4.4 O método do ponto fixo

A construção dos métodos estudados para a determinação de uma sucessão de aproximações parao zero x∗ de f , num determinado intervalo, consideram na definição da sucessão a expressão de f . Noentanto, por vezes é conveniente reescrever a equação f (x) = 0 numa equação equivalente

x = g(x).

Neste caso, o zero de f verifica g(x∗) = x∗.


Definição 4.2. Seja g : I −→ R. Se p ∈ I é tal que g(p) = p, então p diz-se ponto fixo de g.

Os pontos fixos de g, num determinado intervalo, são as abcissas dos pontos de interseção da curvay = g(x), x ∈ I, com a recta y = x. A determinação dos zeros de f passa então pela determinação dospontos fixos de g.

Teorema 4.3. Se g : [a,b] −→ R, é contínua e g([a,b]) ⊆ [a,b], então g tem um ponto fixo em[a,b]. Se

|g′(x)| ≤ K < 1,∀x ∈ [a,b], (4.4.1)

então o ponto fixo é único.

Demonstração.Existência: Se g(a) = a ou g(b) = b, então g tem um ponto fixo. Suponhamos agora que g(a)> a e

g(b)< b. Seja h(x) = g(x)− x, x ∈ [a,b]. Atendendo a que h(a)h(b)< 0 e h é contínua, então h tem umzero neste intervalo. Seja p ∈ (a,b) tal que h(p) = 0. Logo p = g(p).

Unicidade: Se p e q são pontos fixos de g, então g(p) = p e g(q) = q. Temos sucessivamente

|p−q|= |g(p)−g(q)|= |g′(η)||p−q|≤ K|p−q|,

em que η pertence ao intervalo de extremos p e q. Logo

(1−K)|p−q| ≤ 0,

e portanto p = q.

Exemplo 4.3. A equaçãox3 +4x2−10 = 0

tem uma única raiz em [1,2]. De facto, f (x) = x3 + 4x2− 10 é tal que f (1) f (2) < 0, f é contínua eportanto existe x∗ ∈ [1,2] tal que f (x∗) = 0. A raiz anterior é única pois f ′(x) = 3x2 +8x > 0, x ∈ (1,2).

Existem várias equações do tipo x = g(x) equivalentes a f (x) = 0 :

1. g(x) = x3 +4x2 + x−10. Uma vez que g(1) =−4 vem que g([1,2]) 6⊆ [1,2].

2. g(x) = 12

√10− x3. Observamos que g(2) =

√2

2 , e portanto g([1,2]) 6⊆ [1,2].

3. g(x) =√

10x+4 . Para esta expresssão temos

(a) g(1) =√

2, g(2) =√

53 ,

(b) g′(x) =−√

102(x+4)3/2 , x ∈ (1,2), e portanto g é decrescente em [1,2].


Logo

g([1,2]) =

[√53,√

2

]⊂ [1,2].

Concluímos que g tem um ponto fixo x∗ em [1,2]. Para estabelecer a unicidade de x∗ provemos que|g′(x)|< 1, x ∈ [1,2]. Seja h(x) = |g′(x)|, x ∈ [1,2]. Esta função é decrescente em [1,2] pois que

h′(x) =−3√

104(x+4)5/2 < 0, x ∈ [1,2].

Logo

|g′(x)|= h(x)≤ h(1) =

√2

10< 1, x ∈ [1,2].

Consideramos seguidamente o problema da construção de um método iterativo para determinar umaaproximação para o único ponto fixo de g : [a,b]−→R. Para x0 ∈ [a,b], consideremos a sucessão (xn)n∈Ndefinida por

xn+1 = g(xn), n = 0,1 . . . . (4.4.2)

O método iterativo definido por (4.4.2) é designado método do ponto fixo. Na Figura 4.7 mostramos umainterpretação gráfica da sucessão gerada pelo método anterior.

x0 x1x2 x3

g(x0)

g(x1)

g(x2)

g(x3)

Figura 4.7: Ilustração gráfica da construção da sucessão (xn)n∈N do método do ponto fixo.

Exemplo 4.4. Retomemos o Exemplo 4.3 em que considerámos a equação

x3 +4x2−10 = 0

com uma única raiz em [1,2]. Reescrevemos, no intervalo anterior, a equação em estudo na formaequivalente

x =

√10

x+4,

em que g(x) =√

10x+4 tem um ponto fixo único em [1,2].

Consideremos agora a sucessão (xn)n∈N definida por

xn+1 =

√10

xn +4, n = 0,1, . . . . (4.4.3)


Na tabela seguinte apresentamos alguns elementos desta sucessão quando x0 = 1 bem como os corres-pondentes valores de | f |.

n xn f (xn)

1 1.4142135623730950488 0.828427124746190097602 1.3590402174294351594 0.101904633624161938713 1.3660182195292663164 0.0130209922238245255394 1.3651297414694172065 0.00165574924180382598955 1.3652427711324894042 0.00021067461276573120107...

......

21 1.3652300134140968453 7.8058660544728304262 ·10−18

Se considerarmos x0 = 2, obtemos x21 ' 1.36523 e | f (x21)| ' 1.201 ·10−17. Por outro lado, com x0 = 1.5obtemos x21 ' 1.36523 e | f (x21)| ' 0. Neste caso, várias aproximações iniciais geram sucessões em queo valor correspondente da função f parecem convergir para zero.

Os resultados apresentados no exemplo anterior pemitem-nos intuir que a sucessão (xn)n∈N definidapor (4.4.3) converge para a raiz da equação pretendida. No resultado seguinte estabelecemos as condiçõesque nos permitem concluir a convergência do método do ponto fixo.

Teorema 4.4. Seja g : [a,b]−→ [a,b] tal que

|g′(x)| ≤ k < 1, ∀x ∈ [a,b]. (4.4.4)

Então, qualquer que seja a aproximação inicial x0 ∈ [a,b], a sucessão gerada pelo métodoiterativo (4.4.2) converge para o único ponto fixo de g em [a,b].

Demonstração. O Teorema 4.3 permite concluir a existência e unicidade do ponto fixo x∗ ∈ [a,b].Sejam xn+1 e xn definidos por (4.4.2) e x0 ∈ [a,b]. Temos, sucessivamente,

xn+1− x∗ = g(xn)−g(x∗)

= g′(ηn)(xn− x∗),

em que ηn pertence ao intervalo de extremos x∗ e xn.Considerando agora (4.4.4) obtemos

|xn+1− x∗| ≤ k|xn− x∗|, (4.4.5)

e portanto|xn− x∗| ≤ kn|x0− x∗|. (4.4.6)

Atendendo a que k < 1, de (4.4.6), vem

limn→+∞

|xn+1− x∗|= 0.


Corolário 4.1. Nas hipóteses do Teorema 4.4 tem-se

|en| ≤ kn maxx0−a,b− x0 (4.4.7)

e|en| ≤

kn

1− k|x1− x0|. (4.4.8)

Demonstração. A prova de (4.4.7) é imediata. Provemos agora (4.4.8). Consideremos x j+1 e x j definidospor (4.4.2). É fácil demonstrar a desigualdade

|x j+1− x j| ≤ k j|x1− x0|.

Para |xn+m− xn|, temos, sucessivamente,

|xn+m− xn| ≤m+n−1

∑j=n|x j+1− x j|

≤m+n−1

∑j=n

k j|x1− x0|

≤ kn|x1− x0|∞

∑j=0

k j

≤ kn

1− k|x1− x0|

sendo a última igualdade obtida atendendo ao facto de k < 1.A desigualdade

|xn+m− xn| ≤kn

1− k|x1− x0|

permite finalmente concluir (4.4.6).

Notamos que as desiguladades (4.4.7) e (4.4.8) podem ser utilizadas para determinar o número deiterações do método do ponto fixo que devemos considerar para que o erro da aproximação determinadaseja inferior a uma tolerância ε, fixada a priori. Assim, de (4.4.7) obtemos

n >ln(

εmaxb−x0,x0−a

)

lnk,

e, de (4.4.8), vem

n >ln(

ε|x1−x0|

)

lnk.

Observamos também que a velocidade de convergência do método do ponto fixo definido por umafunção g é determinada pelo majorante de |g′|. Consideremos dois métodos do ponto fixo, definidos pelasfunções g1 e g1, para determinar uma aproximação para a raiz de uma determinada equação. Se

max[a,b]|g′1|< max

[a,b]|g′2|,

então devemos escolher o método definido pela função g1.


O Teorema 4.4 é uma condição suficiente para a convergência do método do ponto fixo. Na provadeste teorema estabelecemos a desigualdade (4.4.5) que pode ser reescrita na forma equivalente

|en+1| ≤ k|en|.

Esta desigualdade estabelece que, nas condições do Teorema 4.4, o método do ponto fixo apresentaconvergência linear. No resultado que apresentamos seguidamente estabelecemos que, sob determinadascondições, o método do ponto fixo tem convergência quadrática.

Teorema 4.5. Seja g : [a,b] −→ R com um único ponto fixo x∗ em [a,b]. Se g′(x∗) = 0 e g′′ élimitada num intervalo contendo o ponto fixo, então existe δ > 0 tal que, qualquer que sejax0 ∈ I = [x∗−δ ,x∗+δ ], método do ponto fixo (4.4.2) tem convergência quadrática.

Demonstração. Atendendo a que g′(x∗) = 0, fixado k < 1, existe δ > 0 tal que

|g′(x)| ≤ k, ∀x ∈ I.

Mostremos que g(I)⊆ I. Para o efeito, notamos que, para x ∈ I, se tem sucessivamente

|g(x)− x∗|= |g(x)−g(x∗)|= |g′(η)||x∗− x|≤ k|x∗− x|≤ δ ,

em que η ∈ I.Uma vez que g(I)⊆ I e |g′(x)|< 1 para x ∈ I, pelo Teorema 4.4, a sucessão definida pelo método do

ponto fixo converge para x∗. Provemos agora que tal convergência é quadrática. Atendendo a que se tem

xn+1 = g(xn)

= g(x∗)+g′(x∗)(xn− x∗)+12

g′′(η)(xn− x∗)2

=12

g′′(η)(xn− x∗)2,

com η no intervalo definido por xn e x∗, vem

|en+1|=∣∣∣∣12

g′′(η)

∣∣∣∣|en|2.

Finalmente, atendendo a que g′′ é limitada em algum intervalo contendo x∗, concluímos

|en+1| ≤M|en|2.

Observamos por último que o método do ponto fixo se pode implementar de forma análogaao método de Newton (ver Algoritmo 4.2).


4.5 O método da falsa posição

A aplicação do método de Newton para determinar uma aproximação para a solução x∗ de f (x) = 0requer a existência da primeira derivada. Estudamos seguidamente um método iterativo que segue amesma filosofia do método de Newton mas em que a recta tangente é substituída pela recta secante.Suponhamos que f : [a,b]−→ R é tal que f (a) f (b)< 0. Suponhamos que f (a)< 0, f (b)> 0.

a x1 x2 b

Figura 4.8: Construção gráfica da sucessão (xn)n∈N do método da corda falsa..

1. Consideremos x0 = a e a recta secante à curva y = f (x) nos pontos (a, f (a)),(b, f (b)) de equação

y = b+ f [x0,b](x− x0).


x1 = b− f (b)f [x0,b]

.

Comparamos em seguida f (x1) com f (a) e f (b) e definimos I1 = [a1,b1] tal que x∗ ∈ I1.

2. Consideremos a recta secante à curva y = f (x) nos ponto (a1, f (a1)),(b1, f (b1)) de equação

y = f (b1)+ f [a1,b1](x−b1).


x2 = b1−f (b1)

f [a1,b1].

Comparamos em seguida f (x2) com f (a1) e f (b1) e determinamos I2 = [a2,b2] tal que x∗ ∈ I2.

3. O processo repete-se.

4. Seja In = [an,bn] tal que x∗ ∈ In. Consideremos a recta secante à curva y = f (x) no ponto(an, f (an)),(bn, f (bn)) de equação

y = f (bn)+ f [an,bn](x−bn).

Seja xn+1 a abcissa do ponto de intersecção da recta anterior com o eixo das abcissas, então

xn+1 = bn−f (bn)

f [an,bn]. (4.5.1)


O método apresentado é designado método da falsa posição ou método da corda falsa. Em cadaiteração do método, determinamos a abcissa do ponto de intersecção da recta secante ao gráfico de f empontos determinados e, em seguida, determinamos o novo intervalo comparando o sinal da função nestanova aproximação com os extremos considerados de modo a definir o novo intervalo que nos vai permitircalcular a nova aproximação, ver Figura 4.8.

Estudamos seguidamente as propriedades da sucessão definida por (4.5.1) para funções convexas ouconcavas.

Definição 4.3. Seja f : [a,b]−→ R. Se

f (θx1 +(1−θ)x2)≤ θ f (x1)+(1−θ) f (x2), ∀x1,x2 ∈ [a,b], ∀θ ∈ [0,1],

então f diz-se convexa em [a,b]. Se a desigualdade é estrita, então f diz-se estritamente convexa em[a,b].

Sef (θx1 +(1−θ)x2)≥ θ f (x1)+(1−θ) f (x2), ∀x1,x2 ∈ [a,b],∀θ ∈ [0,1],

então f diz-se côncava em [a,b]. Se a desigualdade é estrita, então f diz-se estritamente côncava em[a,b].

Teorema 4.6. Seja f : [a,b]−→ R convexa ou côncava tal que f (a) f (b)< 0. Então a sucessãogerada pelo método da corda falsa converge e existe n0 ∈ N0 tal que

|en+1| ≤M|en|, n≥ n0, (4.5.2)

com M < 1.

Demonstração. Seja f uma função convexa e admitamos que f (a)< 0, f (b)> 0. Seja (xn)n∈N a sucessãodefinida pelo método da corda falsa. Esta sucessão é crescente e limitada superiormente e portanto éconvergente. Seja x este limite. Provemos que x = x∗. Para o efeito suponhamos que f (x) 6= 0. Entãoexistem ε > 0 e n0 ∈ N0 tais que

| f (xn)| ≥ ε, n≥ n0. (4.5.3)

Da definição (4.5.1) com bn = b, vem

xn+1−b =− f (b)f [xn,b]

=−(b− xn)f (b)

f (b)− f (xn)

e portanto

|xn+1−b| ≤ MM+ ε

|xn−b|,

e ainda

|xn+k−b| ≤(

MM+ ε

)k

|xn−b|.

Logo, tomando o limite quando k→ ∞, vem

x = b.


Mas a sucessão (xn)n∈N é construida de modo a ter-se f (xn) f (b)≤ 0. Consequentemente

0≤ limn→∞

f (xn) f (b) = f (b)2 > 0,

o que constitui uma contradição.Para mostrar que a convergência do método da corda falsa é linear, comecemos por observar que a

recta secante definida a partir de (an, f (an)),(b, f (b)) tem a seguinte equação

y = f (an)+ f [an,b](x−an)

ou sejay = f (xn)+ f [xn,b](x− xn)

e portanto a aproximação xn+1 satisfaz

0 = f (xn)+ f [xn,b](xn+1− xn). (4.5.4)

Por outro lado temos

0 = f (xn)+ f [xn,b](x∗− xn)+ f [x∗,xn,b](x∗−b)(x∗− xn). (4.5.5)

De facto, uma vez que se tem

f [x∗,xn,b](x∗−b)(x∗− xn) =− f [xn,b](x∗− xn)+ f [x∗,xn](x∗− xn)

=− f (xn)− f [xn,b](x∗− xn)

concluímos a igualdade (4.5.5). Conjungando (4.5.4) e (4.5.5), obtemos

f [xn,b](x∗− xn+1) =− f [x∗,xn,b](x∗−b)(x∗− xn),

que nos permite estabelecer

limn→∞

|en+1||en|

= limn→∞

∣∣∣∣f [x∗,xn,b](x∗−b)

f [xn,b]

∣∣∣∣

= limn→∞

∣∣∣∣1−f [x∗,xn]

f [xn,b]

∣∣∣∣

=

∣∣∣∣1−f ′(x∗)f ′(ξ )

∣∣∣∣

em que ξ ∈ (x∗,b). Atendendo a que f é convexa deduzimos 0 6= f ′(x∗)f ′(ξ ) < 1.

Uma variante ao método anterior é chamado método da secante. Consideremos as aproximaçãoes x0 ex1 coincidindo com os extremos do intervalo [a,b] onde f tem apenas um zero. A partir da recta secante àcurva y = f (x) e do ponto de intersecção desta recta com o eixo das abcissas, obtemos

x2 = x1−f (x1)

f [x0,x1].

Considerando agora a secante à curva y = f (x) definida pelos pontos (x1, f (x1)) e (x2, f (x2)) somosfacilmente conduzidos a

x3 = x2−f (x2)

f [x1,x2].

O processo repete-se e a partir das aproximações xn−1 e xn obtemos

xn+1 = xn−f (xn)

f [xn−1,xn]. (4.5.6)

Observamos que a única diferença entre o método da secante e o da falsa posição é a determinaçãodos pontos que definem a secante. É claro que não havendo a escolha criteriosa dos pontos que definem asecante pode surgir facilmente uma situação de divergência.



1. (a) Localize a raiz da equação ln(x−1) = 1x−1 .

(b) Aplique e o método da bisseção duas vezes e indique o majorante do erro com que vemafectada a aproximação calculada.

2. (a) Localize as raízes de x2 +4| ln |x||−4 = 0.

(b) Determine a terceira iteração utilizando o método de Newton.

3. O método de Newton para determinar a raiz da equação f (x) = 0 pode ser indentificado comoum método do ponto fixo. Utilizando esta identificação, prove que o método de Newton temefectivamente ordem de convergência quadrática.

4. Utilizando o método de Newton, indique um método iterativo para calcular o inverso de um númeropositivo β .

5. Considere a equação f (x) = 0 que rescrevemos na forma equivalente x = x+α f (x) com α 6= 0.Indique como deve escolher o parêmtro α de modo a aplicar o método do ponto fixo.

6. (a) Mostre que o método iterativo

xn+1 =c

xp−1n

, n = 0, . . . ,

não é “útil” para encontrar as raiz c1p , p≥ 2, em que c > 0.

(b) Mostre que o método de Newton permite obter uma sucessão convergente para c1p .

7. (a) Localize a raiz positiva da equação cosh(x)−2+2x = 0.

(b) Determine a segunda iteração utilizando ao método da corda falsa.


BIBLIOGRAFIA

[1] R. L. Burden and J. D. Faires. Numerical Analysis. Brooks/Cole Publishing Company, Pacific Grove,1997.

[2] W. Gautschi. Numerical Analysis: An Introduction. Birkhäuser, Berlin, 1987.

[3] R. Kress. Numerical Analysis. Springer, New York, 1997.

[4] P. M. Prenter. Splines and Variational Methods. John Wiley & Sonsr, New York, 1976.

[5] A. Quarteroni, R. Sacco, and F. Saleri. Numerical Mathematics. Springer, New York, 2000.

[6] A. Quarteroni and F. Saleri. Scientific Computing with MATLAB and Octave. Springer, New York,2007.

133

Documents

Matemática Numérica I