Embed Size (px)
Citation preview
1 Representacao de Numeros e Erros
1.1 Sistemas de Vırgula Flutuante
Para podermos fazer calculos e necessario antes de mais escolhermos umsistema para representar os numeros com que trabalhamos. Supondo quevamos trabalhar com numeros reais, os sistemas habitualmente utilizadospara os representar sao os sistemas de vırgula (ponto) flutuante. As-sim, vamos comecar por definir matematicamente estes sistemas.
Seja β um numero natural, diferente de 1, a que chamaremos base dosistema. A base e o numero de dıgitos diferentes que usamos para rep-resentar os numeros. A base mais corrente e a decimal, em que se usamdez dıgitos (os algarismos); quando usamos esta base, temos β = 10. Noentanto, se atentarmos a forma como os numeros sao representados inter-namente nos computadores e outros sistemas de calculo, verificaremos quea base aı utilizada e a binaria, ou seja, β = 2, ja que, por razoes tecnicas,e conveniente trabalhar apenas com dois sımbolos diferentes: 0 e 1. Nessecaso, cada sımbolo representado designa-se por bit. Uma vez escolhida abase, qualquer elemento x do sistema de vırgula flutuante representa-se naforma
x = σ × 0.a1a2a3...an × βt (1.1)
onde σ representa o sinal (σ = ±1), ai sao dıgitos da base considerada et e um numero inteiro. A sucessao de dıgitos a1a2a3...an, onde a1 6= 0,designa-se mantissa. Assim, alem da base, qualquer sistema de vırgula flu-tuante caracteriza-se pelo comprimento da mantissa, isto e, o numeron de dıgitos que a compoem. Finalmente, um sistema de vırgula flutu-ante caracteriza-se pelos limites inferior e superior do expoente t, que repre-sentaremos respectivamente por t1e t2. Chegamos assim a seguinte definicao.
Definicao 1. Sistema de vırgula flutuante com base β e n dıgitos namantissa:
V F (β, n, t1, t2) = {x ∈ IR : x = σ×0.a1a2a3...an×βt, σ = ±1, a1 6= 0, t1 ≤ t ≤ t2, t ∈ Z}∪{0}.
De acordo com esta definicao, como e natural, o numero 0 pertence aqualquer sistema VF , embora formalmente ele nao possa ser representado naforma (1.1), ja que o primeiro dıgito da mantissa, por definicao, e diferentede zero.
Exemplo 1.1. Calculadora em que os numeros sao representados nabase decimal, com 12 dıgitos na mantissa e com o expoente entre -99 e 99.Neste caso, o sistema utilizado e V F (10, 12,−99, 99).
1
Exemplo 1.2. Computador em que os numeros sao representados nabase binaria, sendo reservados 56 bits para a mantissa e 8 bits para o ex-poente. Dos 8 bits do expoente 7 sao reservados ao seu valor absoluto e umao sinal, pelo que o seu valor pode variar entre −27+1 = −127 e 27−1 = 127.Logo, o sistema considerado e V F (2, 56,−127, 127). Note-se que, devido acondicao a1 6= 0, no caso do sistema binario obtem-se a1 = 1, qualquer queseja o numero representado. Isto faz com que o primeiro dıgito da mantissaseja superfluo. Na pratica, esse dıgito e usado para representar o sinal damantissa.
Propriedades dos sistemas de vırgula flutuante:
1. Qualquer sistema VF e finito. Determinemos o numero de elementospositivos do sistema V F (β, n, t1, t2). O numero de mantissas diferentes eβn−1(β − 1) (o primeiro dıgito da mantissa nao pode ser 0, por definicao).O numero de expoentes diferentes e t2− t1 +1. Logo, o numero de elementosdo sistema V F (β, n, t1, t2), tendo em conta os numeros negativos e o zero, e2βn−1(β − 1)( t2 − t1 + 1) + 1 (no caso do exemplo 1.1, temos 2× 9× 199×109 + 1 ≈ 3.6 × 1012 elementos, enquanto que no exemplo 1.2 o numero deelementos e 2× 255× 255 + 1 ≈ 1.84× 1019).
Em qualquer sistema de vırgula flutuante existe um elemento maximo,cujo valor e M = (1 − β−n)βt2(no caso do exemplo 1.1, temos M = (1 −10−12)1099 ≈ 1099 , enquanto que para o exemplo 1.2 temos M = (1 −2−155)2127 ≈ 1.70× 1038).
3. Em qualquer sistema de vırgula flutuante existe um elemento com omınimo valor absoluto, igual a m = β−1βt1 = βt1−1(no caso do exemplo1.1, temos m = 10−100 , enquanto que para o exemplo 1.2 temos m =2−128 ≈ 2.9× 10−39).
1.2 Representacao de numeros em sistemas de vırgula flutu-ante
Visto que qualquer sistema VF e finito, ele contem apenas uma pequenaparte dos numeros reais. Quando um numero real nao pertence ao sistemaVF considerado, para o representar nesse sistema e necessario fazer umacerta aproximacao, chamada arredondamento. Denotemos fl(x) a repre-sentacao do real x no sistema VF considerado. Se x ∈ V F (β, n, t1, t2) entaofl(x) = x (diz-se que x tem representacao exacta nesse sistema). Casocontrario, isto e, se x /∈ V F (β, n, t1, t2), ha que escolher fl(x) e essa es-colha pode ser feita de diferentes maneiras. Para melhor compreender este
2
processo, suponhamos que x = σ × 0.a1a2a3...anan+1... × βt. Note-se quequalquer numero real pode ser representado nesta forma, sendo que a man-tissa, regra geral, e uma dızima infinita. Segundo a forma mais simples dearredondamento, o arredondamento por corte, escolhe-se:
fl(x) = σ × 0.a1a2a3...an × βt.
Outra forma mais sofisticada de determinar fl(x) consiste em defini-lo pelaformula
fl(x) =
{σ × 0.a1a2a3...an × βt, se an+1 < β/2;
σ × (0.a1a2a3...an + β−n)× βt, se an+1 ≥ β/2.
Esta forma de aproximacao chama-se arredondamento simetrico.Este tipo de arredondamento envolve um erro igual ao do arredondamentopor corte, no caso de an+1 < β/2, ou menor, no caso em que an+1 ≥ β/2.
Note-se que, ao considerar um certo sistema VF, nem todos os numerosreais podem ser representados nele, mesmo com arredondamento. Os numerosx, tais que |x| > M ou |x| < m, nao tem qualquer representacao no sistema,pelo que ao tentar representa-los ocorrem situacoes de erro. No primeirocaso, essas situacoes designam-se por overflow, enquanto no segundo casosao referidas como underflow .
Exercıcio 1.1 Para cada um dos seguintes numeros reais obtenha (casoseja possıvel) a sua representacao no sistema V F (10, 3,−99, 99),utilizandoarredondamento simetrico:
a) x = 10b) x = 0.001235c) x = 1001d) x = 1/3d) x = 10100
e) x = 10−101
Resposta:
x fl(x)
100 0.100× 103
0.001235 0.124× 10−2
-1001 -0.100× 104
1/3 0.333
10100 nao tem representacao (overflow)
10−101 nao tem representacao (underflow)
3
1.3 Erros de arredondamento
Quando se aproxima um numero real x pela sua representacao em vırgula flu-tuante, fl(x), comete-se um erro geralmente designado por erro de arredonda-mento:
ear = fl(x)− x.
Outras grandezas relacionadas sao o erro de arredondamento absoluto
|ear| = |x− fl(x)|
e o erro de arredondamento relativo:
|δar| =|x− fl(x)||x|
.
Para caracterizar a precisao com que os numeros reais sao aproximadosnum sistema VF utiliza-se o conceito de unidade de arredondamento. Aunidade de arredondamento do sistema VF e um numero real u, tal que
|δar| ≤ u, ∀x ∈ IR,m ≤ |x| ≤M.
O valor de u depende, evidentemente, dos parametros do sistema VF consid-erado, mais precisamente, de n e β. Logicamente, para o mesmo valor de β, aunidade de arredondamento sera tanto mais pequena quanto maior for n, istoe, quanto mais dıgitos utilizarmos para representar os numeros tanto menorsera o erro de arredondamento relativo. Para analisarmos esta questao empormenor, consideremos um numero real x arbitrario e representemo-lo naforma x = σ × 0.a1a2a3...anan+1... × βt. Para simplificar, comecemos porconsiderar o caso do arredondamento por corte. Como ja vimos, neste casofl(x) = σ × 0.a1a2a3...an × βt. Por conseguinte, o erro de arredondamentoabsoluto e
|ear| = |x− fl(x)| = 0.00...0an+1...× βt < βt−n.
No que diz respeito ao erro de arredondamento relativo, temos
|δar| =|x− fl(x)||x|
<βt−n
βt−1= β1−n.
Por conseguinte, qualquer que seja x,tal que m ≤ |x| ≤M,verifica-se
|δar| < β1−n.
4
Logo, podemos afirmar que, no caso do arredondamento por corte, aunidade de arredondamento e
u = β1−n.
Raciocınios semelhantes levam-nos a conclusao que, no caso do arredonda-mento simetrico, a unidade de arredondamento e
u =1
2β1−n.
Por exemplo, no caso do sistema V F (10, 12,−99, 99), e assumindo que oarredondamento e simetrico, temos u = 0.5× 10−11.
Exercıcio 1.2. Para cada um dos numeros referidos no exercıcio 1.3,considerando de novo o sistema V F (10, 3,−99, 99), determine o erro dearredondamento absoluto e o erro de arredondamento relativo (nos casosem que existe representacao). Compare este ultimo com a unidade dearredondamento do sistema.
Resposta:
x fl(x) |ear| |δar|100 0.100× 103 0 0
0.001235 0.124× 10−2 0.5× 10−5 0.004
−1001 -0.100× 104 1 0.001
1/3 0.333 0.33× 10−3 0.001
A unidade de arredondamento, neste caso, e 0.5×10−2 = 0.005, pelo quetodos os numeros considerados tem erro de arredondamento relativo inferiora este valor, como seria de esperar.
1.4 Propagacao dos Erros
Sejam x e y valores aproximados dos numeros reais x e y, respectivamente.Denotaremos por ex e |δx| o erro e o erro relativo de x, repectivamente:
ex = x− x, |δx| =∣∣∣∣ x− xx
∣∣∣∣ .De modo analogo se definem o erro e o erro relativo de y. Suponhamos
que x e y sao dados de um calculo que pretendemos efectuar. O nossoobjectivo e determinar qual o efeito dos erros destes dados no resultado.Para comecar, consideremos o caso das operacoes aritmeticas.
5
Adicao/Subtraccao
Representemos por ex±yo erro de x± y. Note-se que
x± y = (x+ ex)± (y + ey) = (x± y) + (ex ± ey).
Por conseguinte, para o erro de x± y temos
ex±y = ex ± ey
e, para o erro absoluto,|ex±y| ≤ |ex|+ |ey|.
Quanto ao erro relativo, podemos escrever
|δx±y| =|ex ± ey||x± y|
≤ |xδx|+ |yδy||x± y|
. (1.2)
Daqui resulta que, se x ± y for proximo de zero, entao o erro relativo dex±y pode ser muito maior que o dos dados. Voltaremos a este assunto maistarde.
Multiplicacao
No caso da multiplicacao, temos
x.y = (x+ ex).(y + ey) = (x.y) + yex + xey + exey.
Admitindo que ex e ey sao grandezas pequenas, o seu produto pode serdesprezado na expressao anterior, pelo que obtemos
ex.y = x.y − x.y ≈ yex + xey.
Logo, para o erro relativo do produto, resulta
|δx.y| =|ex.y||x.y|
≈ |yex + xey||x.y|
≤ |δx|+ |δy|. (1.3)
6
Divisao
Para deduzir a formula do erro do quociente, suponhamos que os valoresde ex e ey sao desprezaveis em comparacao com x e y, respectivamente.Podemos entao fazer a seguinte aproximacao:
x
y= (x+ ex)
1
y
1
1 +eyy
≈ (x+ ex)1
y
(1− ey
y
)=x
y+yex − xey
y2.
Daqui resulta que
ex/y =x
y− x
y≈ yex − xey
y2.
Quanto ao erro relativo do quociente, obtem-se
|δx/y| = |ex/y||y||x|≈ |yex − xey|
y2
|y||x|≤ |ex||x|
+|ey||y|
= |δx|+ |δy|. (1.4)
Os calculos que acabamos de efectuar mostram que, no caso da multi-plicacao e da divisao, o erro relativo dos resultados e da mesma ordem queo erro relativo dos dados, ou seja, destas operacoes nao resulta uma perdade precisao. Ja no caso da adicao e da subtraccao, como vimos, tal perda deprecisao pode ocorrer. Esse fenomeno designa-se cancelamento subtractivoe e ilustrado pelo exercıcio que se segue.
Exercıcio 1.3. Considere os numeros x = π e y = 2199/700.
1. Determine x e y com 4 dıgitos na mantissa, usando arredondamentosimetrico. Obtenha ainda x− y.
Solucao.
x = 0.3141592 · · · × 10; x = fl(x) = 0.3142 · · · × 10;y = 0.3141428 · · · × 10; y = fl(y) = 0.3141 · · · × 10.
Logo, x− y = 0.1× 10−2.
2. Calcule os erros absolutos e relativos de x e y . Comente.
Solucao.
Dado Erro absoluto Erro relativox |ex| = |x− x| = 0.41× 10−3 |δx| = |ex|/|x| = 0.131× 10−3
y |ey| = |y − y| = 0.43× 10−3 |δy| = |ey|/|y| = 0.137× 10−3.
Como era de esperar, os erros de arredondamento relativos dos dadossao inferiores a unidade de arredondamento do sistema que, neste caso,e u = 0.5× 101−4 = 0.5× 10−3.
7
3. Represente os numeros x e y em ponto flutuante, mas com 6 algarismosna mantissa. Com base nestas novas aproximacoes, calcule de novox− y e comente.
Solucao. Neste caso temos:
x = 0.3141592 · · · × 10; x = fl(x) = 0.314159 · · · × 10;y = 0.3141428 · · · × 10; y = fl(y) = 0.314143 · · · × 10.
Logo, x− y = 0.16× 10−3, o que e muito diferente do valor obtido naalınea 1. Isto sugere que, na alınea 1, houve uma perda de precisao,resultante de cancelamento subtractivo.
4. Tomando como valor exacto da diferenca o resultado da alınea anterior,determine o erro relativo do valor de x−y, obtido na alınea 1. Se usassea estimativa (1.2) para o erro relativo da diferenca, chegaria a mesmaconclusao?
Solucao. Comparando os resultados da alınea 1 e da alınea 3, temos
|δx−y| =|ex−y||x− y|
≈ 0.001− 0.00016
0.00016= 5.25.
Vemos que o erro relativo do resultado da alınea 1 e muito superior aunidade, o que significa uma perda total de precisao. Por outro lado,se usassemos a estimativa (1.2), terıamos
|δx−y| ≤|xδx|+ |yδy||x− y|
≈ 0.00084
0.00016= 5.25,
ou seja, neste caso verifica-se a igualdade na relacao (1.2).
1.5 Estabilidade de algoritmos
Quando se efectua um calculo, geralmente efectua-se por passos. Assim, oerro cometido em cada passo acumula-se com os erros dos passos anteriores.Em resultado, o erro do resultado final pode ser muito maior do que o errode cada passo isolado.
O conjunto dos passos que levam a obtencao de um dado resultadodesigna-se algoritmo. O mesmo resultado pode ser obtido, em princıpio,atraves de algoritmos diferentes. No entanto, os erros propagam-se de formadiferente em cada algoritmo. Por isso, os resultados que se obtem, para omesmo problema, atraves de algoritmos diferentes podem divergir significa-tivamente. Surge assim a definicao de estabilidade numerica.
8
Definicao. Um algoritmo diz-se estavel (ou numericamente estavel)para um certo conjunto de dados se, a pequenos valores dos erros relativosde arredondamento dos dados e da unidade de arredondamento do sistemacorresponderem pequenos valores do erro relativo do resultado.
O exercıcio que se segue ilustra o conceito de estabilidade numerica.Exercıcio 1.3 Considere a funcao real de variavel real
f(x) =1− cosx
x2. (1.5)
1. Supondo que utiliza um sistema de vırgula flutuante com 10 dıgitosna mantissa e arredondamento simetrico, calcule f(10−6) aplicando aformula (1.5).
Solucao. A formula (1.5) pode ser aplicada em 3 passos. Sendox = 10−6, temos
z1 = cosx = 1;z2 = 1− z1 = 0;
z3 = f(x) = z2x2
= 0.
2. Obtenha uma aproximacao de f(10−6), utilizando o desenvolvimentode f em serie de Taylor, em torno de x = 0.
Solucao. Como e sabido, para valores de x proximos de zero, a funcaocos(x) admite o seguinte desenvolvimento em serie de Taylor:
cosx = 1− x2
2+x4
4!+O(x6),
Daqui obtem-se facilmente que
f(x) =1− cos x
x2=
1
2− x2
4!+O(x4). (1.6)
Utilizando a formula (1.6), num sistema VF com 10 dıgitos, obtem-sef(10−6) = 0.5000000000.
3. Sabendo que 1 − cosx = 2 sin2(x/2), calcule f(10−6) utilizando umanova formula para f .
Solucao. Transformando a formula (1.5) obtem-se
f(x) =1− cos x
x2=
2
x2sin2(x/2). (1.7)
9
A semelhanca do que fizemos na alınea 1, apliquemos a formula (1.7)em 5 passos:
w1 = x/2 = 0.5× 10−6;w2 = sen(w1) = 0.5× 10−6
w3 = w22 = 0.25× 10−12;
w4 = w3/x2 = 0.25;
w5 = f(x) = 2w4 = 0.5.
4. Compare os valores obtidos nas alıneas anteriores e classifique os al-goritmos quanto a estabilidade.
Solucao. Verifica-se que o valor obtido na alınea 3 e uma boa aprox-imacao de f(10−6), ja que coincide com o valor dado pela serie de Tay-lor (em todos os 10 dıgitos). Pelo contrario, o valor obtido pelo algo-ritmo da alınea 1 dava uma ma aproximacao (nem um unico dıgito cor-recto). Este facto deve-se apenas aos erros de arredondamento cometi-dos, ja que as formulas usadas sao equivalentes e, se trabalhassemoscom valores exactos, obterıamos em ambos os casos o mesmo resultado.Esta situacao pode interpretar-se do seguinte modo: para valores dex proximos de 0 o algoritmo considerado em 1) e instavel, enquanto oalgoritmo considerado em 3) e estavel.
2 Metodos Numericos paraEquacoes Nao Lineares
2.1 Localizacao de raizes de uma funcao
Seja f uma funcao real, definida num certo intervalo [a, b]. O ponto z ∈[a, b] diz-se um zero ou uma raiz de f se e so se f(z) = 0. No caso def ′(z) 6= 0, z diz-se um zero simples. Se f ′(z) = 0, z diz-se um zero multiplo.Mais precisamente, se f ∈ Ck(z) e se f ′(z) = f ′′(z) = ... = f (k−1)(z) =0, f (k)(z) 6= 0,entao z diz-se um zero de multiplicidade k da funcao f .
Exemplo 2.1. Seja f um polinomio de grau n. Entao, de acordo com oteorema fundamental da algebra, f tem, no total, n raizes em C (somandoas multiplicidades). Em particular, se tivermos f(x) = x2 + 2x + 1, f temuma raiz de multiplicidade dois (raiz dupla) em z = −1. De facto, verifica-se f(−1) = 0; visto que f ′(x) = 2x + 2, temosf ′(−1) = 0; alem disso,como f ′′(x) = 2, f ′′(−1) 6= 0. Se considerarmos o polınomio de terceiro grauf(x) = x3 − x = x(x − 1)(x + 1), verificamos facilmente que ele tem tres
10
raizes simples: z1 = −1, z2 = 0, z3 = 1. Quanto ao polinomio f(x) = x3 + 1,
tem apenas uma raiz real (z1 = −1) e duas raizes complexas (z2,3 = 1±√
3i2 ).
De um modo geral, a determinacao das raizes de um polinomio de graun (ou zeros de uma equacao algebrica) e um problema muito complexoque ocupou os matematicos de varias civilizacoes. Desde o princıpio doseculo XX sabe-se, gracas a Abel, que nao existem formulas resolventes paraequacoes algebricas de qualquer grau. Mais precisamente, para uma equacaoalgebrica de grau superior a 4 nao e possıvel, atraves de uma formula geral,exprimir todas as raizes atraves dos seus coeficientes.
Este exemplo ilustra a importancia dos metodos numericos para a res-olucao de equacoes. Ate no caso de equacoes relativamente simples, como asequacoes algebricas, e impossıvel calcular as suas raizes unicamente atravesde formulas analıticas. Por outro lado, mesmo nos casos em que existemformulas resolventes, estas sao por vezes tao complexas que se torna maiseficiente determinar as raizes a partir de um metodo numerico. Tal e ocaso de algumas equacoes algebricas de terceiro e quarto grau, por exemplo.Naturalmente, isso pressupoe que se escolha um metodo adequado a equacaoconsiderada.
Para tratar o problema do calculo numerico de todas as raizes da funcaof e necessario, em primeiro lugar, localizar as raizes, isto e, determinar, paracada raiz, um intervalo que a contenha e que nao contenha nenhuma outra.Com esse objectivo, recordemos dois teoremas da analise matematica.
Teorema 2.1 (teorema de Bolzano). Se f for contınua em [a, b] e sef(a)f(b) < 0, entao f tem, pelo menos, uma raiz em ]a, b[.
Teorema 2.2 (corolario do teorema de Rolle). Se f for contınuaem [a, b] e continuamente diferenciavel em]a, b[ e se f ′(x) 6= 0 em ]a, b[ entaof tem, no maximo, uma raiz em]a, b[.
Combinando estes dois teoremas e alguns outros resultados da analise, epossıvel, em muitas situacoes, localizar as raizes reais de uma equacao.
Exemplo 2.2. Com base nos teoremas 2.1 e 2.2, determinar o numerode raizes reais da equacao
ex − x2 − 2x = 0.5
e determinar um intervalo que contenha cada uma delas (mas nao contenhaas restantes).
Este problema e equivalente a determinar os zeros da funcao de variavelreal f(x) = ex−x2− 2x− 0.5. Esta funcao e evidentemente contınua em IR,
11
-3 -2 -1 1 2 3x
-2
2
4
f
Figure 1: Grafico relativo ao Exemplo 2.2
assim como todas as suas derivadas de qualquer ordem (ver o seu grafico nafig.1).
Para facilitar a analise do problema, comecemos por calcular uma tabelade valores de f e de f ′.
Tabela 2.1
x −3 −2 −1 0 1 2 3
f(x) −3.45 −0.365 0.868 0.5 −0.782 −1.11 4.59
f ′(x) 4.05 2.14 0.368 −1 −1.28 1.39 12.1
Observando a tabela 2.1, verifica-se imediatamente que o teorema 2.1e aplicavel a f nos intervalos [−2,−1],[0, 1]e [2, 3]. Daqui se conclui que aequacao considerada tem, pelo menos, 3 raizes: z1 ∈ [−2,−1], z2 ∈ [0, 1] ez3 ∈ [2, 3]. Pela aplicacao do teorema 2.2 podemos concluir tambem que, emcada um destes intervalos, a funcao f tem exactamente 1 raiz. Para isso,consideremos a primeira e a segunda derivada de f : f ′(x) = ex − 2x −2, f ′′(x) = ex − 2.
Em relacao a segunda derivada, verifica-se facilmente que ela e positivapara x > ln 2 e negativa para x < ln 2. Temos f ′′(ln 2) = 0, f ′′′(ln 2) = 2, peloque f ′ tem em ln 2 um ponto de mınimo. Assim, no intervalo [−2,−1], f ′
e decrescente. Recorrendo de novo a tabela 2.1, verifica-se que f ′ e semprepositiva neste intervalo, pelo que, pelo teorema 2.2, podemos concluir quef tem uma unica raiz z1 no intervalo [−2,−1]. Do mesmo modo, podemosobservar que a funcao f ′ e crescente em [2, 3] e, de acordo com a tabela,toma sempre valores positivos neste intervalo. Aplicando o teorema 2.2neste intervalo, constata-se que f tem nele uma unica raiz z3. Para aplicar oteorema 2.2 no intervalo [0,1], comecemos por recordar que a funcao f ′ tem
12
um ponto de mınimo em x = ln 2, que pertence a este intervalo. Note-seque f ′(ln 2) = −1.38 < 0; uma vez que, de acordo com a tabela, f ′(0) ef ′(1) tambem sao negativos, podemos concluir que f ′ e negativa em todo ointervalo [0,1]. Logo, o teorema 2.2 e aplicavel neste intervalo e a funcao temnele uma unica raiz z2. Resta esclarecer uma questao: sera que a funcao ftem alguma raiz alem das que acabamos de localizar? Para responder a estapergunta, recordemos que a segunda derivada de f tem uma unica raiz real(ln 2). Daqui, pelo teorema de Rolle, somos levados a concluir que a primeiraderivada de f tem, no maximo, duas raizes reais. Finalmente, aplicando oteorema de Rolle a f ′,conclui-se que f tem, no maximo, 3 raizes reais. Comoja vimos que existem, pelo menos, 3 raizes (z1, z2 e z3), concluımos que estassao as unicas raizes de f.
2.2 Metodo da Bisseccao
Um dos metodos mais simples para o calculo aproximado de raizes e ometodo da bisseccao. Para se poder aplicar este metodo basta conhecerum intervalo que contenha uma unica raiz da funcao considerada (a quale, por condicao, contınua). A ideia do metodo e construir uma sucessao deintervalos encaixados [a, b] ⊃ [a1, b1] ⊃ ... ⊃ [ak, bk] tais que a)cada intervalotem o comprimento igual a metade do intervalo anterior; b)f(ai)f(bi) <0, i = 1, 2, ..., k. Desta ultima condicao, pelo teorema 2.1, resulta que a raize um ponto comum a todos os intervalos da sucessao. Assim, construindoum numero suficientemente grande de intervalos, e possıvel aproximar a raizcom a precisao que se pretender. Vejamos em pormenor o algoritmo destemetodo.
1o Passo. Dado um intervalo [a, b] e uma funcao f tais que f(a)f(b) <0, determina-se o ponto medio desse intervalo: x1 = a+b
2 . Se, por felizcoincidencia, se verificar f(x1) = 0, x1 e a raiz procurada Suponhamos quef(x1) 6= 0,entao verifica-se f(x1)f(a) < 0 ou f(x1)f(a) > 0. No primeirocaso, podemos afirmar que z ∈ [a, x1]; no segundo caso, z ∈ [x1, b]. Entao ointervalo [a1, b1] pode ser definido do seguinte modo:
se f(x1)f(a) < 0 , entao a1 = a; b1 = x1; senao, a1 = x1; b1 = b.Em qualquer dos casos, o novo intervalo [a1, b1] satisfaz f(a1)f(b1) < 0.2o Passo. Repetem-se as accoes do primeiro passo, substituindo o in-
tervalo [a, b] por [a1, b1] e representando por x2 o ponto medio deste inter-valo.Deste passo, resulta o intervalo [a2, b2].
Generalizando, no k-esimo passo (iteracao) realizam-se as seguintesoperacoes.
Determina-se o ponto medio do intervalo anterior: xk =ak−1+bk−1
2 ;
13
Se f(xk)f(ak−1) < 0 , entao ak = ak−1; bk = xk; senao, ak = xk; bk =bk−1. Do k-esimo passo obtem-se o intervalo [ak, bk].
O processo e interrompido quando for satisfeita a condicao de paragem:bk − ak < ε,onde ε e uma tolerancia estabelecida de acordo com a precisaoque se pretende obter.
Note-se que o comprimento do k-esimo intervalo, por construcao, e bk −ak = b−a
2k, pelo que tende para zero, quando k tende para infinito. Logo,
qualquer que seja ε, a condicao de paragem e satisfeita ao fim de um certonumero de passos, que depende do comprimento do intervalo inicial e de ε.Mais precisamente, temos
b− a2k
< ε⇔ b− aε
< 2k ⇐⇒ k > log2(b− aε
).
Assim, o numero de passos do metodo da bisseccao que e necessariorealizar ate satisfazer a condicao de paragem e o mınimo inteiro k, tal quek > log2( b−aε ).
Se tomarmos como k-esima aproximacao da raiz z o valor de xk podemosafirmar que o erro absoluto de xk satisfaz
|z − xk| <bk−1 − ak−1
2=b− a
2k.
E costume nos metodos computacionais representar o erro da k-esimaaproximacao da raiz por ek; usando esta notacao, podemos afirmar que, nometodo da bisseccao e valida a estimativa do erro:
|ek| <b− a
2k.
Exercıcio 2.1. a)Recorrendo ao teorema 2.1, justifique que a raizcubica de 2 pertence ao intervalo [1.2, 1.3].
b) Baseando-se na alınea anterior, efectue tres iteracoes (passos) dometodo da bisseccao com o objectivo de calcular um valor aproximado daraiz cubica de 2.
c) Quantas iteracoes teria que efectuar se pretendesse determinar 3√
2com um erro absoluto inferior a 0.001?
Resposta. Comecemos por observar que determinar a raiz cubica de 2equivale a resolver a equacao f(x) = x3 − 2 = 0.
a) Temos que f(1.2) = 1.23 − 2 = −0.272 < 0 e f(1.3) = 1.33 − 2 =0.197 > 0. Uma vez que a funcao f e contınua, pelo teorema 2.1 concluımosque a raiz procurada esta no intervalo [1.2, 1.3].
14
b) Comecemos com o intervalo [a, b] = [1.2, 1.3]. Entao x1 = a+b2 = 1.25.
Verifica-se que f(1.25) = −0.047 < 0, de onde resulta que f(1.25)f(1.2) > 0.Logo, o intervalo a considerar na iteracao seguinte e [a1, b1] = [1.25, 1.3]. Porconseguinte, x2 = a1+b1
2 = 1.275. Neste caso, f(1.275) = 0.0727 > 0, de onderesulta que f(1.275)f(1.25) < 0. Assim, o intervalo a considerar na terceiraiteracao e [a2, b2] = [1.25, 1.275]. Finalmente, x3 = a2+b2
2 = 1.2625. Nesteponto, temos f(1.2625) = 0.012 > 0, pelo que o intervalo a considerar naiteracao seguinte sera [a3, b3] = [1.25, 1.2625]
c) O comprimento do intervalo inicial e b − a = 0.1. Logo, para atingiruma precisao de ε = 0.001, o numero de iteracoes e dado por log2( b−aε ) =log2( 0.1
0.001) = 6.64. Ou seja, a precisao sera atingida ao fim de 7 iteracoes.
2.3 Metodo do Ponto Fixo
2.3.1 Definicao e exemplos de pontos fixos
Definicao. Seja g uma funcao real, definida num certo intervalo [a, b] ⊂ IR.O numero z ∈ [a, b] diz-se um ponto fixo de g se g(z) = z.
Exemplo 2.3a) Seja g(x) = αx+ β, α 6= 1, α, β ∈ IR.O ponto fixo de g e o que satisfaz αz + β = z, ou seja z = β
1−α . Porexemplo, se for α = 2, β = −3, obtem-se z = 3 (fig.2).
b) Seja g(x) = x2+1. Neste caso, a equacao dos pontos fixos e z2+1 = z.
Por conseguinte, temos z = 12 ±
√12
2 − 1, ou seja, nao existem pontos fixos
reais (fig. 3).c) g(x) = x2. A equacao dos pontos fixos, neste caso, e z2 = z. Logo,
existem dois pontos fixos: z1 = 0, z2 = 1 (fig.4).d) g(x) = cosx. Embora nao seja possıvel determinar analiticamente o
ponto fixo desta funcao, e facil verificar que ela tem um ponto fixo (unico)no intervalo [0, 1]. Com efeito, se definirmos f(x) = cosx − x, verifica-seque f(0) = 1, f(1) = cos 1 − 1 < 0. Logo, sendo a funcao f contınua, peloteorema 2.1, ela tem, pelo menos, um zero z em ]0, 1[. Nesse ponto z verifica-se cos z = z, pelo que z e um ponto fixo de g. Por outro lado, f e uma funcaocontinuamente diferenciavel e a sua derivada, f ′(x) = −senx−1, e negativaem [0, 1]. Logo, pelo Teorema 2.2, f tem uma unica raiz neste intervalo, quee tambem o unico ponto fixo de g (fig.5).
Dada uma funcao g, determinar os seus pontos fixos equivale a calcularas raizes da equacao g(x)− x = 0, ou, por outras palavras, calcular os zerosda funcao f(x) = g(x)−x. Inversamente, se for dada uma equacao f(x) = 0,
15
1 2 3 4 5x
-2
2
4
6
y
z
g
y=x
Figure 2: Exemplo 2.3a
0.5 1 1.5 2x
1
2
3
4
5
y
g
y=x
Figure 3: Exemplo 2.3b
calcular as raizes dessa equacao equivale a determinar os pontos fixos deoutra funcao.
Exemplo 2.4. Consideremos de novo a equacao ex − x2 − 2x = 0.5(exemplo 2.2). Esta equacao pode ser reescrita de varias formas, todas elasequivalentes:
ex − x2 − 0.5
2= x; (2.1)
√ex − 2x− 0.5 = x; (2.2)
ln(x2 + 2x+ 0.5) = x. (2.3)
No caso da equacao (2.1), as raizes da equacao inicial sao vistas como
os pontos fixos da funcao g1(x) = ex−x2−0.52 . Em relacao a equacao (2.2),
ela remete-nos para os pontos fixos de g2(x) =√ex − 2x− 0.5. Note-se
16
0.5 1 1.5 2x
1
2
3
4
y
z2z1
g
y=x
Figure 4: Exemplo 2.3c
0.5 1 1.5 2x
0.5
1
1.5
2
y
zg
y=x
Figure 5: Exemplo 2.3d
que, neste caso, as equacoes so sao equivalentes para valores positivos de x(a funcao g2 nao pode tomar valores negativos). Em particular, a raiz z1,negativa, nao e ponto fixo de g2.Quanto a equacao (2.3), diz-nos que as raizesda equacao inicial sao pontos fixos da funcao g3(x) = ln(x2 +2x+0.5). Nestecaso, a equivalencia tambem nao e valida para qualquer valor de x, ja que odomınio da funcao g3 so inclui os valores de x,para os quais x2+2x+0.5 > 0.Das raizes da equacao inicial apenas z2 e z3 satisfazem esta condicao. Logo,z2 e z3 sao tambem pontos fixos de g3, enquanto z1 nao o e.
O Exemplo 2.4 mostra-nos que as raizes de uma equacao dada podemser tratadas como pontos fixos de diferentes funcoes. Veremos que essefacto pode ser utilizado na escolha de metodos numericos para o calculoaproximado destas raizes.
17
2.3.2 Sucessoes numericas geradas por funcoes
Dada uma funcao real g,com domınio num certo intervalo [a, b], e um numerox0, tal que x0 ∈ [a, b], e possıvel gerar uma sucessao de numeros reais {xn}do seguinte modo:
xk+1 = g(xk), k = 0, 1, ... (2.4)
Se a imagem do intervalo [a,b] estiver contida no proprio intervalo, entaoa relacao (2.4) permite-nos definir uma sucessao infinita de elementos doconjunto considerado. Neste caso, chamaremos a g a funcao iteradora e aostermos xk da sucessao as iteradas. Veremos como as sucessoes geradas dessemodo podem ser utilizadas para aproximar as raizes de uma equacao dada.
Exemplo 2.5. Seja g(x) = x2. O domınio desta funcao e R e a imagemdo intervalo [0,1] por esta funcao e o proprio intervalo.Se tomarmos x0 = 0,a funcao g gera uma sucessao constante:{0, 0, 0, ..}. Se considerarmos 0 <x0 < 1,entao a sucessao gerada e {x0, x
20, x
40, ...} e converge para 0, que e
um dos pontos fixos de g. Se tomarmos x0 = 1, a sucessao das iteradas e denovo constante: {1, 1, 1, ..}( 1 tambem e um ponto fixo de g). Se tomarmosx0 > 1, a sucessao vai ser divergente (tende para infinito).
O exemplo 2.5 sugere-nos que, quando a sucessao gerada por uma funcaog converge, entao o seu limite e um ponto fixo da funcao g. De facto, podeprovar-se que assim e.
Teorema 2.3. Seja {xn} uma sucessao gerada pela funcao g queconverge para um certo limite z. Se g for contınua em z, entao z e pontofixo de g.
Demonstracao . Uma vez que z = limn→∞ xn, temos
z = limn→∞
xn+1 = limn→∞
g(xn).
Da continuidade de g em z resulta que limn→∞ g(xn) = g(limn→∞ xn) =g(z). Obtemos assim que z = g(z),como se pretendia demonstrar.
Exercıcio 2.2. Considere a sucessao gerada pela funcao g(x) = sen(x),com x0 = 1. Prove que esta sucessao converge. Qual e o seu limite?
Resposta. Para provar que a sucessao converge, basta provar que emonotona e limitada. Note-se que, sendo 0 < x < 1, temos 0 < sen(x) < x.Daqui resulta que 1) todos os termos da sucessao considerada pertencemao intervalo [0,1]; 2) a sucessao e mototona decrescente, ja que xk+1 =sen(xk) < xk.Por conseguinte a sucessao e monotona e limitada, logo econvergente. De acordo com o teorema 2.3, a sucessao considerada, sendo
18
convergente, deve convergir para um ponto fixo da funcao iteradora. O unicoponto fixo da funcao g(x) = sen(x) e z = 0, logo e para esse ponto que asucessao converge.
2.3.3 Teorema do Ponto Fixo
O teorema 2.3 diz-nos que uma sucessao gerada por uma funcao iteradorag, a convergir, converge para um ponto fixo daquela funcao. Fica por re-sponder a questao: em que condicoes essa sucessao converge? A resposta aesta questao e dada por um teorema fundamental da Analise Matematica,o teorema do ponto fixo. Embora este teorema possa ser formulado numcontexto mais vasto, por agora, limitar-os-emos a uma versao simplificada,em que g e uma funcao de uma variavel real.
Teorema 2.4 (teorema do ponto fixo). Seja g uma funcao de variavelreal e [a, b] um intervalo tais que:
1. g([a, b]) ⊂ [a, b];
2. g e continuamente diferenciavel em [a, b];
3. maxx∈[a,b] |g′(x)| = L < 1;
Entao sao validas as seguintes afirmacoes:
1. g tem um unico ponto fixo z em [a, b];
2. se x0 ∈ [a, b], a sucessao gerada pela funcao g converge para o pontofixo z.
Demonstracao. 1. Para demonstrar a existencia de, pelo menos, umponto fixo, defina-se a funcao h(x) = g(x) − x. Esta funcao e obviamentecontınua em [a, b]; se tivermos g(a) = a (ou g(b) = b), teremos que a (oub, respectivamente) e ponto fixo de g. Caso contrario, de acordo com acondicao 1), a funcao h satisfaz h(a) = g(a) − a > 0, h(b) = g(b) − b < 0;entao, pelo teorema de Bolzano, existe, pelo menos, um ponto z ∈ [a, b], talque h(z) = 0, ou seja, g(z) = z; logo, z e ponto fixo de g.
Para demonstrar a unicidade, suponhamos que existem em [a, b] doispontos fixos distintos z1,z2. Por definicao, temos g(z1) = z1, g(z2) = z2.Logo, |g(z1) − g(z2)| = |z1 − z2|. Por outro lado, usando o Teorema deLagrange e a condicao 3, temos
|g(z1)− g(z2)| ≤ maxx∈[a,b]
|g′(x)||z1 − z2| = L|z1 − z2|.
19
Obtemos assim que
|z1 − z2| ≤ L|z1 − z2|,
ou seja,|z1 − z2|(1− L) ≤ 0. (2.5)
Mas, de acordo com a condicao 3, L < 1; logo, da desigualdade (2.5)resulta que |z1−z2| = 0, o que contradiz a hipotese de z1 e z2 serem distintos.Desta contradicao conclui-se a unicidade do ponto fixo.
2. Para demonstrar a segunda afirmacao seja x0 um ponto arbitrario de[a, b]. Pela condicao 1 do teorema, temos que x1 = g(x0) tambem pertenceao intervalo [a, b] . Do mesmo modo se conclui que todos os elementos dasucessao, gerada pela funcao g, pertencem aquele intervalo. Vamos agoraprovar que esta sucessao converge para o ponto fixo z. Pela condicao 3 doteorema, temos
|xn − z| = |g(xn−1)− g(z)| ≤ L|xn−1 − z|. (2.6)
Aplicando n vezes a desigualdade (2.6), conclui-se que
|xn − z| ≤ Ln|x0 − z|. (2.7)
Como, por condicao, L < 1, da desigualdade (2.7) resulta que |xn−z| →0, quando n → ∞ (qualquer que seja x0 ∈ [a, b]). Fica assim demonstradaa segunda afirmacao do teorema.
O teorema do ponto fixo nao so nos garante a existencia de um unicoponto fixo z da funcao g num dado intervalo, como nos indica um metodopara obter aproximacoes desse ponto. Na realidade, se tomarmos qualquerponto inicial x0 dentro do intervalo e construirmos a sucessao, gerada pelafuncao g, de acordo com o teorema do ponto fixo essa sucessao convergepara z. O metodo baseado nesta construcao chama-se metodo do pontofixo. Este metodo permite-nos, dada uma funcao iteradora g e um intervalo[a, b] (tal que sejam satisfeitas as condicoes do teorema 2.4), obter umaaproximacao tao precisa quanto quisermos do ponto fixo de g em [a, b]. Oalgoritmo e extremamente simples:
1. Escolher um ponto x0 ∈ [a, b];
2. calcular cada nova iterada pela formula xn = g(xn−1), n = 1, 2, ...;
3. Parar quando se obtiver uma aproximacao aceitavel (sobre os criteriosde paragem falaremos mais abaixo).
20
2.3.4 Estimativa do erro do metodo do ponto fixo
Para efeitos praticos, interessa-nos nao so saber as condicoes em que ummetodo converge, mas tambem estimar o erro das aproximacoes obtidas.No caso do metodo do ponto fixo, a resposta a esta questao e dada peloseguinte teorema.
Teorema 2.5. Nas condicoes do Teorema 2.4, sao validas as seguintesestimativas do erro:
|xn − z| ≤ Ln|x0 − z| (2.8)
(estimativa apriori),
|xn − z| ≤L
1− L|xn − xn−1|, n ≥ 1 (2.9)
(estimativa aposteriori), onde xn−1, xn, sao duas iteradas consecutivas dometodo do ponto fixo e L = maxx∈[a,b] |g′(x)|.
Demonstracao. A formula (2.8) ja foi provada durante a demonstracaodo teorema do ponto fixo (ver (2.7)). Quanto a formula (2.9), pode serprovada do seguinte modo. Comecemos por observar que
|xn−1 − z| = |z − xn−1| ≤ |z − xn|+ |xn − xn−1|. (2.10)
Por outro lado, visto que, de acordo com (2.6),
|xn − z| ≤ L|xn−1 − z|,
da formula (2.10) resulta que
|xn−1 − z|(1− L) ≤ |xn − xn−1|. (2.11)
Observando que 1−L > 0 (pela condicao 3 do teorema 2.4), podem dividir-sepor este valor ambos os membros da desigualdade (2.11), obtendo-se assim
|xn−1 − z| ≤1
1− L|xn − xn−1|. (2.12)
Finalmente, das desigualdades (2.12) e (2.6) obtem-se a estimativa (2.9).
21
Exercıcio 2.3. Considere a equacao cos(x)− 2x = 0.
1. Com base no teorema do ponto fixo, mostre que esta equacao tem umaunica raiz no intervalo [0.4, 0.5] e que o metodo do ponto fixo convergepara essa raiz.
Resolucao. Comecemos por observar que qualquer raiz da equacaodada e um ponto fixo de g(x) = cosx
2 . Mostremos agora que a funcaog satisfaz, no dado intervalo, as condicoes do teorema do ponto fixo.
(a) Para mostrar que g([0.4, 0.5]) ⊂ [0.4, 0.5] comecemos por calcularas imagens dos extremos do intervalo: g[0.4] = cos(0.4)/2 =0.46053 ∈ [0.4, 0.5]; g(0.5) = cos(0.5)/2 = 0.43879 ∈ [0.4, 0.5].Por outro lado, a funcao g e decrescente em [0.4, 0.5] (g′(x) =−senx/2 e negativa naquele intervalo). Daqui se conclui queg([0.4, 0.5]) ⊂ [0.4, 0.5].
(b) g e continuamente diferenciavel em IR e, em particular, no inter-valo considerado.
(c) L = maxx∈[0.4,0.5] |g′(x)| = maxx∈[0.4,0.5]|senx|
2 = sen0.52 = 0.2397 <
1.
Concluimos assim que todas as condicoes do teorema do ponto fixoestao satisfeitas, pelo que o metodo do ponto fixo com a funcao iter-adora g(x) = cosx/2 converge para o ponto fixo.
2. Tomando como aproximacao inicial x0 = 0.4, calcule as duas primeirasiteradas do metodo.
Resolucao. Tomando como aproximacao inicial x0 = 0.4, as duasprimeiras aproximacoes iniciais sao x1 = g(x0) = 0.46053;x2 = g(x1) =0.44791.
3. Obtenha uma estimativa do erro da aproximacao x2, calculada naalınea anterior.
Resolucao.Usando a formula (2.9), obtem-se
|z − x2| ≤L
1− L|x2 − x1| =
0.2397
1− 0.2397|0.44791− 0.46053| = 0.00397.
4. Nas condicoes da alınea anterior, quantas iteracoes e necessario efec-tuar para garantir que o erro absoluto da aproximacao obtida e inferiora 0.001?
22
Resolucao. Para responder a esta questao, e necessario aplicar aestimativa a priori (2.8). De acordo com esta estimativa, temos
|xn − z| ≤ Ln|x0 − z| ≤ 0.2397n|0.5− 0.4| = 0.1× 0.2397n, n ≥ 1.
Logo, para garantir que o erro absoluto da n-esima iterada e inferiora um certo ε, basta escolher n de tal modo que 0.2397n < 10ε. Destainequacao, resulta que n > log0.2397(10ε). No caso de ε = 0.001,obtem-se n > 3.22, de onde se conclui que o erro satisfaz a condicaoexigida ao fim de 4 iteracoes.
2.3.5 Classificacao dos pontos fixos
De acordo com o teorema do ponto fixo, a convergencia das sucessoes geradaspor uma certa funcao g num intervalo [a, b] depende do comportamento dasua derivada g′ nesse intervalo. Isto leva-nos a classificar os pontos fixos zde uma funcao g de acordo com o valor de g′(z). Neste paragrafo iremosassumir que g e continuamente diferenciavel (ou seja, g′ e contınua), pelomenos, numa vizinhanca de cada ponto fixo de g.
Assim, um ponto fixo z de uma funcao g diz-se atractor se |g′(z)| < 1.De facto, se |g′(z)| < 1 e g′ e contınua em z, entao existe uma vizinhancaVε(z) = (z − ε, z + ε) tal que maxx∈Vε(z) |g′(z)| = L < 1. Por outro lado,se x ∈ Vε(z), temos |g(x) − g(z)| ≤ L|x − z| < |x − z| < ε, ou seja, g(x)tambem pertence a Vε(z). Logo, se o intervalo [a, b] estiver contido em Vε(z),nesse intervalo a funcao g satisfaz as condicoes do teorema do ponto fixo.Concluimos portanto que, se z for um ponto fixo atractor, entao existeuma vizinhanca Vε(z) tal que, se x0 ∈ Vε(z) entao a sucessao geradapor g converge para z.
Por outro lado, se z e um ponto fixo de g e |g′(z)| > 1, entao z diz-seum repulsor. Nesse caso, e facil verificar que nenhuma sucessao geradapela funcao g converge para z (excepto a sucessao constante {z, z, . . . } ouqualquer sucessao da forma {. . . , x, z, z, . . . }, onde x e tal que g(x) = z). Defacto, se z e um repulsor, entao existe uma vizinhanca Vε(z) = (z − ε, z + ε)tal que |g′(z)| > 1, ∀x ∈ Vε(z). Assim, seja xk um termo de uma sucessaogerada pela funcao g e suponhamos que xk ∈ Vε(z), xk 6= z. Entao temos|xk+1−z| = |g(xk)−g(z)| ≥ minx∈Vε(z) |g′(x)||xk−z| > |xk−z|. Logo, xk+1
esta mais distante de z do que xk. Se o ponto xk+1 tambem pertencer a Vε(z),entao o mesmo raciocınio se aplica a esse ponto, e vemos que a sucessao seafasta de z. A unica possibilidade de uma sucessao nao constante convergirpara z, sendo z repulsor, e o caso de essa sucessao conter um ponto x (naopertencencente a vizinhanca referida), tal que g(x) = z.
23
0.2 0.25 0.3 0.35 0.4
0.2
0.25
0.3
0.35
0.4
Figure 6: Iteracoes da funcao g(x) = 1.5x(1− x), com x0 = 0.2.
Finalmente, se tivermos |g′(z)| = 1, entao o ponto fixo z diz-se neutro.Neste caso, existem sucessoes geradas pela funcao g que convergem para z eoutras que nao convergem (mesmo que x0 esteja proximo do ponto fixo z).
Exemplo 2.6. Consideremos a funcao g(x) = kx(1 − x), onde k > 0.Esta funcao e conhecida como ”funcao logıstica” e tem diversas aplicacoesem ecologia matematica. Para determinarmos os pontos fixos desta funcao,para um certo valor de k, resolva-se a equacao
kz(1− z) = z. (2.13)
E facil verificar que esta equacao tem 2 raizes: z1 = 0 e z2 = 1− 1/k. Con-sideremos o caso de k = 1.5. Nesse caso, os dois pontos fixos de g sao z1 = 0e z2 = 1/3. Para os classificarmos, observemos que g′(x) = 1.5 − 3x. Logog′(0) = 1.5 e g′(1/3) = 1.5 − 1 = 0.5, ou seja, z1 e repulsor e z2 e atrac-tor. Isto significa que a) nenhuma sucessao gerada pela funcao g poderaconvergir para 0 (excepto a sucessao constante, igual a 0, ou a sucessao{1, 0, 0, . . . }); b) se x0 for suficientemente proximo de 1/3, a sucessao ger-ada por g converge para z2 = 1/3. Mais precisamente pode provar-se que,se 0 < x0 < 1, entao a sucessao {xk} converge para z2. As figuras 6 e 7ilustram esta afirmacao.
Exemplo 2.7. Seja g(x) = x2 + x. Esta funcao tem um ponto fixo(unico) em z = 0. Visto que g′(z) = 2z + 1 = 1, este ponto fixo e neutro.Vejamos qual e o comportamento das sucessoes geradas por esta funcao.
Seja x0 = 0.12. Entao, temos que x1 = x20 +x0 = 0.1344; x2 = x2
1 +x1 =0.152463, etc. E facil verificar que, neste caso, a sucessao e crescente e tendepara +∞.
24
0.3 0.35 0.4 0.45 0.50.3
0.35
0.4
0.45
0.5
Figure 7: Iteracoes da funcao g(x) = 1.5x(1− x), com x0 = 0.5
Mas se escolhermos como iterada inicial x0 = −0.12, ja teremos x1 =x2
0+x0 = −0.1056;x2 = x21+x1 = −0.0945. Neste caso, a sucessao e crescente
e converge para o ponto fixo z = 0. As figuras 8 e 9 ilustram este exemplo.
2.3.6 Monotonia das iteradas do metodo do ponto fixo
Suponhamos que z e um ponto fixo atractor da funcao g. Como vimosno paragrafo anterior, isto verifica-se sempre que |g′(z)| < 1, isto e, −1 <g′(z) < 1. Como vimos, neste caso, qualquer sucessao gerada pela funcaog, com x0 suficientemente proximo de z, converge para z. Neste paragrafo,vamos investigar em que condicoes essa sucessao e monotona (crescente oudecrescente). Tal como antes, admitimos que g e continuamente diferenciavelnuma vizinhanca de z.
Caso 1. Suponhamos que 0 < g′(z) < 1. Entao, da continuidade daderivada de g, resulta que existe uma vizinhanca Vε(z) = (z − ε, z + ε),tal que, se x ∈ Vε(z), entao 0 < g′(x) < 1. Suponhamos que xk e umtermo de uma sucessao, gerada pela funcao g, tal que xk ∈ Vε(z). Parasermos mais especıficos admitamos que z < xk < z + ε. Nesse caso, umavez que xk+1 = g(xk), aplicando o teorema de Lagrange, existe um pontoξk, z ≤ ξk ≤ xk, tal que
xk+1 − z = g(xk)− g(z) = g′(ξk)(xk − z). (2.14)
Por construcao, temos xk− z > 0 e g′(ξk) > 0. Logo, xk+1 > z. Concluimosque, se xk > z, entao tambem xk+1 > z. Por outro lado, uma vez que ze um ponto atractor (g′(ξk) < 1), xk+1 deve estar mais proximo de z doque xk, de onde se conclui que xk+1 < xk. Como o mesmo raciocınio se
25
0.05 0.1 0.15 0.2
0.05
0.1
0.15
0.2
Figure 8: Iteracoes da funcao g(x) = x2 + x, com x0 = 0.12
-0.14-0.12 -0.1 -0.08-0.06-0.04-0.02
-0.14
-0.12
-0.1
-0.08
-0.06
-0.04
-0.02
Figure 9: Iteracoes da funcao g(x) = x2 + x, com x0 = −0.12
26
0.4 0.42 0.44 0.46
0.4
0.42
0.44
0.46
Figure 10: Iteracoes da funcao g(x) = cosx2 , com x0 = 0.39
.
aplica a todas as iteradas subsequentes, podemos concluir que, neste caso,a sucessao {xn} e decrescente (pelo menos, a partir da ordem k). Estasituacao e ilustrada pelo grafico da fig. 7. Se tivermos xk < z, o mesmoraciocınio leva-nos a conlusao que xk+1 > xk. Nesse caso, a sucessao dasiteradas sera crescente. Esta situacao e a que esta representada na fig. 6e na fig. 9. Em qualquer das duas situacoes, as sucessoes sao monotonas.
Caso 2. Suponhamos agora que −1 < g′(z) < 0. Entao, da continuidadeda derivada de g, resulta que existe uma vizinhanca Vε(z) = (z − ε, z + ε),tal que, se x ∈ Vε(z), entao −1 < g′(x) < 0. Admitindo que xk pertence aessa vizinhanca, a igualdade (2.14) e aplicavel. Neste caso, admitindo quexk > z, dessa igualdade resulta que xk+1 < z (uma vez que g′(ξk) < 0). Seaplicarmos o mesmo raciocınio as iteradas seguintes, concluımos que xk+2 >z, xk+3 < z, etc. Se, pelo contrario, tivermos xk < z, entao xk+1 > z,xk+2 < z, etc. Ou seja, neste caso, as iteradas vao ser alternadamentemaiores ou menores que z (uma sucessao deste tipo diz-se alternada). Umapropriedade interessante das sucessoes alternadas e que o limite da sucessaoesta sempre entre dois termos consecutivos, isto e xk < z < xk+1 ou xk+1 <z < xk. Isto permite-nos obter um majorante do erro absoluto de xk+1 alemdaqueles que foram obtidos em 2.3.4:
|xk+1 − z| < |xk+1 − xk|. (2.15)
A sucessao das iteradas do exercıcio 2.3, em que g′(z) < 0, e um exemplo deuma sucessao alternada. Na fig. 10 estao representados graficamente algunstermos desta sucessao.
27
2.3.7 Divergencia do metodo do ponto fixo
O estudo dos pontos repulsores no paragrafo 2.3.5 permite-nos formular oseguinte criterio de divergencia do metodo do ponto fixo.
Teorema 2.6 Se g for continuamente diferenciavel em [a, b] e se |g′(x)| >1,∀x ∈ [a, b], entao nenhuma sucessao gerada pela funcao g pode convergirpara qualquer ponto do intervalo [a, b](excepto a sucessao constante {z, z, ...}ou qualquer sucessao da forma {. . . , x, z, z, ...}, onde z e ponto fixo de g e xe tal que g(x) = z).
Demonstracao. De acordo com as condicoes do teorema e com a clas-sificacao dos pontos fixos, se a funcao g tiver algum ponto fixo em [a, b]esse ponto fixo e repulsor. Por outro lado, se uma sucessao gerada pelafuncao g convergir, ela converge para um ponto fixo de g ( Teorema 2.3).Da conjugacao destes dois factos resulta a afirmacao do teorema.
2.3.8 Ordem de convergencia
Um dos conceitos fundamentais da teoria dos metodos iterativos e o conceitode ordem de convergencia. Este conceito permite-nos comparar a rapidezcom que diferentes metodos convergem e escolher, em cada caso, o metodomais eficiente. Nas definicoes que se seguem representaremos por {xn},n ∈N , um sucessao convergente de numeros reais e por z o seu limite.
Definicao. Diz-se que uma sucessao {xn} tem convergencia, pelomenos, linear (pelo menos,de ordem 1) se existir uma constante k∞ < 1tal que
k∞ = limn→∞
|z − xn+1||z − xn|
;
este limite designa-se coeficiente assimptotico de convergencia.Sempre que k∞ 6= 0, a ordem de convergencia e exactamente 1. Se
k∞ = 0, a ordem de convergencia e superior a 1, caso que estudaremos maisadiante.
Note-se que o coeficiente assimptotico de convergencia nos permite com-parar, quanto a rapidez de convergencia, metodos diferentes que tenhamconvergencia linear. Com efeito, quanto mais pequeno (mais proximo de 0)for k∞ mais rapida sera a convergencia.
Exemplo 2.8. Consideremos a sucessao {xn}, tal que xn+1 = xn/a,a > 1, com x0 ∈ IR. E facil verificar que esta sucessao converge para z = 0,qualquer que seja x0 ∈ IR, ja que este e o unico ponto fixo da funcao iteradorag(x) = x/a. Alem disso, este ponto e atractor, ja que g′(x) = 1/a < 1,∀x ∈ IR. Verifiquemos que esta sucessao tem convergencia linear. Para isso,
28
calculemos o quociente:
k∞ = limn→∞
|z − xn+1||z − xn|
= limn→∞
|xn+1||xn|
=1
a< 1. (2.16)
Concluımos assim que temos convergencia linear. Alem disso, o coefi-ciente assimptotico de convergencia e k∞ = 1
a . Assim, a convergencia seratanto mais rapida quanto maior for a.
Analisemos agora os casos em que a ordem de convergencia e superior a1 (convergencia supralinear).
Definicao. Diz-se que uma sucessao {xn} tem convergencia, pelomenos, de ordem p > 1 se existir uma constante k∞ tal que
k∞ = limn→∞
|z − xn+1||z − xn|p
;
esta constante designa-se coeficiente assimptotico de convergencia. Sek∞ 6= 0, entao a ordem de convergencia e exactamente p.
Exemplo 2.9. Consideremos a sucessao {xn}, tal que xn+1 = bxαn, onde
b 6= 0 , α > 1, com |x0| < |b|−1α−1 . E facil verificar que esta sucessao converge
para z = 0, se x0 satisfizer a condicao indicada. De facto, o ponto z = 0e um ponto fixo da funcao iteradora g(x) = bxα. Alem disso, este ponto e
atractor, ja que g′(0) = 0. Por outro lado, se |x0| < |b|−1α−1 , entao |x1| < |x0|
e, de um modo geral, teremos que |xn+1| < |xn|, ∀n, pelo que a sucessao edecrescente em modulo. Verifiquemos qual e a ordem de convergencia destasucessao. Para isso calculemos o limite:
limn→∞
|z − xn+1||z − xn|p
= limn→∞
|xn+1||xn|p
= limn→∞
|bxαn||xn|p
. (2.17)
Para que este limite seja finito, devemos ter p = α. Neste caso, verifica-sek∞ = |b|. Logo, a ordem de convergencia e α e o coeficiente assimptotico deconvergencia e |b|.
A ordem de convergencia do metodo do ponto fixo depende das pro-priedades da funcao iteradora g. O teorema que se segue diz-nos quais ascondicoes que a funcao g deve satisfazer para garantir que o metodo doponto fixo tenha ordem de convergencia p ≥ 1.
Teorema 2.6 (ordem de convergencia do metodo do ponto fixo). Supon-hamos que g satisfaz as condicoes do teorema do ponto fixo em [a, b] e que,alem disso, g ∈ Cp[a, b], com p ≥ 1. Seja ainda g′(z) = g′′(z) = ... =g(p−1)(z) = 0, g(p)(z) 6= 0, onde z ∈ [a, b] e um ponto fixo de g. Entao
1. g tem um unico ponto fixo z em [a, b];
29
2. se x0 ∈ [a, b], entao a sucessao gerada pela funcao g converge para zcom ordem p;
3. o coeficente assimptotico de convergencia e k∞ = |g(p)(z)|p! .
Demonstracao. A primeira afirmacao resulta do teorema 2.4 (teorema doponto fixo). Resta-nos provar os pontos 2 e 3. Com esse fim, escreva-se odesenvolvimento de g em serie de Taylor, em torno de z:
g(x) = g(z) + g′(z)(x− z) +g′′(z)
2(x− z)2 + · · ·+ g(p)(ξ)
p!(x− z)p, (2.18)
onde ξ ∈ int(z, x). Em particular, se escrevermos a formula (2.18) comx = xm, atendendo as condicoes do teorema, obtem-se
g(xm) = g(z) +g(p)(ξm)
p!(xm − z)p, (2.19)
onde ξm ∈ int(z, xm). Da formula (2.19) resulta imediatamente que
xm+1 − z =g(p)(ξm)
p!(xm − z)p. (2.20)
Dividindo ambos os membros de (2.20) por (xm − z)p e tomando o modulo,obtem-se
|xm+1 − z||xm − z|p
=|g(p)(ξm)|
p!. (2.21)
Calculando o limite quando m→∞ em ambos os membros de (2.21), obtem-se
limm→∞
|xm+1 − z||xm − z|p
=|g(p)(z)|
p!. (2.22)
Da igualdade (2.22) resulta imediatamente que a sucessao {xm} tem ordem
de convergencia p e que k∞ = |g(p)(z)|p! .
Observe-se que, como caso particular do Teorema 2.6, com p = 1, seobtem a seguinte afirmacao: se g satisfizer as condicoes do teorema do pontofixo em [a, b] e se g′(z) 6= 0, entao, qualquer que seja x0 ∈ [a, b], a sucessaogerada pela funcao g converge linearmente para z e o coeficiente assimptoticode convergencia e k∞ = |g′(z)|.
Exercıcio 2.10 Considere a funcao iteradora g(x) = 12
(x+ 1
x
).
30
1. Mostre que os pontos fixos de g sao z1 = 1,z2 = −1.
Resolucao. Da equacao g(z) = 12
(z + 1
z
)= z, multiplicando ambos
os membros por 2z, obtem-se z2 + 1 = 2z2. daqui resulta que z2 = 1,pelo que os pontos fixos de g sao z1 = 1,z2 = −1.
2. Mostre que estes pontos sao atractores.
Resolucao. Visto que g′(x) = 12 −
12x2
, obtem-se g′(1) = g′(−1) = 0.Logo, estes pontos sao atractores.
3. Seja x0 ∈ [1, 2]. Mostre que a sucessao gerada pela funcao g con-verge para z1 = 1, determine a ordem e o coeficiente assimptotico deconvergencia.
Resolucao. Em primeiro lugar, mostremos que a funcao g satisfazas condicoes do teorema do ponto fixo em [1, 2]. Ja sabemos queg′(x) = 1
2 −1
2x2, logo g e continuamente diferenciavel em [1, 2]. Alem
disso, e facil verificar que g′(x) ≥ 0,∀x ∈ [1, 2]. Logo, g e crescenteem [1, 2]. Entao, para mostrar que g([1, 2]) ⊂ [1, 2], basta-nos verificarque g(1) = 1 ∈ [1, 2] e g(2) = 5/4 ∈ [1, 2]. Por outro lado, temosmaxx∈[1,2] |g′(x)| = |g(2)| = 3
8 < 1.
Para determinar a ordem de convergencia e o coeficiente assimptoticode convergencia da sucessao considerada, vamos aplicar o teorema 2.6.Para isso, analisemos as derivadas de g. Ja sabemos que g′(1) = 0.Quanto a segunda derivada, temos g′′(x) = 1
x3. Logo,g′′ e contınua
em [1, 2] e g′′(1) = 1 6= 0. Daqui resulta que o teorema 2.6 e aplicavelcom p = 2 e a ordem de convergencia e 2. Quanto ao coeficiente
assimptotico de convergencia, temos k∞ = |g′′(1)|2 = 1
2 .
2.4 Metodo de Newton
No paragrafo anterior vimos que o metodo do ponto fixo tem um vastodomınio de aplicacao e permite, na maioria dos casos, obter boas aprox-imacoes de raizes de equacoes. No entanto, vimos tambem que, em geral,aquele metodo garante apenas primeira ordem de convergencia (ordens su-periores so se obtem, de acordo com o teorema 2.6, se algumas derivadas dafuncao iteradora se anularem no ponto fixo).
O metodo de Newton, que vamos estudar neste praragrafo, tem a impor-tante vantagem de proporcionar, em geral, convergencia de segunda ordem(quadratica). E um dos metodos mais frequentemente utilizados, ja quecombina a rapidez de convergencia com a simplicidade do processo iterativo.
31
Mais adiante, veremos que o metodo de Newton pode ser encarado comoum caso particular do metodo do ponto fixo. Por agora, vamos introduzi-loatraves da sua interpretacao geometrica.
2.4.1 Interpretacao geometrica do metodo de Newton
Seja f uma funcao continuamente diferenciavel num certo intervalo [a, b].Suponhamos que nesse intervalo a funcao f tem uma unica raiz e que a suaderivada nao se anula ( f ′(x) 6= 0,∀x ∈ [a, b]).
Sendo x0 um ponto arbitrario de [a, b], podemos tracar a tangente aografico de f que passa pelo ponto (x0, f(x0)) (ver fig. 11). Sendo f ′(x0) 6= 0,por condicao, essa recta intersecta o eixo das abcissas num certo ponto(x1, 0). Para determinar x1, comecemos por escrever a equacao da tangenteao grafico de f :
y − f(x0) = f ′(x0)(x− x0). (2.23)
Igualando y a 0 na equacao (2.23), obtem-se a abcissa x1 procurada:
x1 = x0 −f(x0)
f ′(x0). (2.24)
As iteradas seguintes obtem-se do mesmo modo. Mais precisamente, paradeterminar x2, traca-se a tangente ao grafico de f que passa pelo ponto(x1, f(x1)) e procura-se o ponto onde essa recta intersecta o eixo das abcis-sas; e assim sucessivamente. Obtem-se deste modo uma sucessao de pontos{xk}, k = 0, 1, 2, . . . , que podem ser calculados pela formula de recorrencia
xk+1 = xk −f(xk)
f ′(xk). (2.25)
A interpretacao geometrica, representada na Fig.11, sugere-nos que estasucessao converge para a raiz z da equacao considerada. Nos paragrafosseguintes, vamos demonstrar analiticamente que de facto assim e.
2.4.2 Estimativa do erro do metodo de Newton
Em primeiro lugar, vamos deduzir uma formula que nos permite majoraro erro de cada iterada do metodo de Newton, admitindo que e conhecidoum majorante do erro da iterada anterior. Vamos supor que a funcao fsatisfaz no intervalo [a, b] as condicoes formuladas no paragrafo anterior (fe continuamente diferenciavel em [a, b] e a sua derivada nao se anula nesteintervalo). Alem disso, vamos supor que a segunda derivada de f tambem
32
3.5 4 4.5 5
2.5
5
7.5
10
12.5
15
x0x1x2x3
z
Figure 11: Interpretacao geometrica do metodo de Newton.
e contınua neste intervalo. Seja {xn}, n = 0, 1, 2, . . . a sucessao das iteradasdo metodo de Newton (que se consideram pertencentes ao intervalo [a, b]).Se desenvolvermos f em serie de Taylor, em torno de xk, obtem-se
f(x) = f(xk) + (x− xk)f ′(xk) +(x− xk)2
2f ′′(ξk), (2.26)
onde ξk ∈ int(xk, x). Escrevendo o desenvolvimento (2.26) com x = z,obtem-se
f(z) = f(xk) + (z − xk)f ′(xk) +(z − xk)2
2f ′′(ξk) = 0, (2.27)
onde ξk ∈ int(xk, z). Uma vez que, por condicao, f ′(xk) 6= 0, podemosdividir ambos os membros de (2.27) por f ′(xk), obtendo assim
f(xk)
f ′(xk)+ (z − xk) +
(z − xk)2
2f ′(xk)f ′′(ξk) = 0. (2.28)
Atendendo a formula iteradora do metodo de Newton (2.25), da equacao(2.28) resulta que
z − xk+1 = −(z − xk)2
2f ′(xk)f ′′(ξk). (2.29)
A igualdade (2.29) da-nos a relacao que procuravamos entre o erro de xk+1(ek+1)e o erro de xk(ek). No segundo membro desta desigualdade aparece o valorf ′′(ξk), o qual nao podemos calcular exactamente, ja que sabemos apenasque ξk e um ponto situado entre xk e z. Por isso, para podermos majoraro erro absoluto de xk(|ek|), precisamos de majorar o modulo da segundaderivada de f (que se supoe contınua). Introduzindo a notacao
M = maxx∈[a,b]
|f ′′(x)| (2.30)
33
da igualdade (2.29) obtem-se a seguinte relacao:
|ek+1| ≤ |ek|2M
2|f ′(xk)|. (2.31)
Saliente-se que na desigualdade (2.31) |ek+1| e comparado com o quadradode |ek|, o que indica um rapido decrescimento do erro. Seja
µ = minx∈[a,b]
|f ′(x)|. (2.32)
Entao, a desigualdade (2.31) pode ser reforcada, substituindo |f ′(xk)| porµ:
|ek+1| ≤ |ek|2M
2µ. (2.33)
Nesta ultima desigualdade o segundo membro nao depende de k. Na pratica,usam-se frequentemente as formulas (2.31) e (2.33) para obter uma estima-tiva de |ek+1|.
Exemplo 2.10. Consideremos a equacao f(x) = cosx−2x = 0, que ja foianalisada no Exercıcio 2.3. Vamos obter aproximacoes da raiz desta equacao,situada no intervalo [0.4, 0.5], aplicando o metodo de Newton, e majorar oerro dessas aproximacoes. Sendo x0 = 0.4, da formula (2.25) obtem-se quex1 = 0.45066547 e x2 = 0.45018365. Vamos agora obter majorantes de |e1|e |e2|. Em primeiro lugar, note-se que |e0| ≤ 0.5−0.4 = 0.1. Para podermosaplicar a desigualdade (2.31) e necessario majorar |f ′′(x)| e minorar |f ′(x)|.Temos f ′(x) = − sinx− 2 e f ′′(x) = − cosx, logo
µ = minx∈[0.4,0.5]
|f ′(x)| = minx∈[0.4,0.5]
|2 + sinx| = 2 + sin 0.4 = 2.389;
M = maxx∈[0.4,0.5]
|f ′′(x)| = maxx∈[0.4,0.5]
| cosx| = cos 0.4 = 0.921.
Por conseguinte, da desigualdade (2.33) resulta a seguinte majoracao parao erro absoluto de x1:
|e1| ≤M
2µ|e0|2 ≤
0.921
2× 2.3890.01 = 0.001927.
Em relacao ao erro de x2, obtem-se, do mesmo modo:
|e2| ≤M
2µ|e1|2 ≤
0.921
2× 2.3890.001927 = 0.696× 10−7.
Vemos assim que, com duas iteradas apenas, conseguiu-se obter um resul-tado de alta precisao. Para terminar este exemplo apresentamos na Tabela
34
2.2 uma comparacao dos resultados obtidos para esta equacao pelos metodosde Newton e do ponto fixo.
k xk(Ponto Fixo) |ek| xk(Newton) |ek|0 0.4 0.0501 0.4 0.05011 0.46053 0.0105 0.45066547 0.48× 10−3
2 0.44791 0.0022 0.45018365 0.4× 10−7
Tabela 2.2.Resultados numericos do metodo de Newton e do metodo doponto fixo (Exemplo 2.10).
2.4.3 Condicoes suficientes de convergencia do metodo de New-ton
Ate agora, analisamos o erro do metodo de Newton, partindo do princıpioque a aproximacao inicial e tal que as iteradas convergem para a raiz procu-rada. No entanto, nem sempre e facil prever, para uma dada aproximacaoinicial, se o metodo de Newton vai convergir, e para que raiz ele vai convergir(se a equacao tiver varias raizes). Neste paragrafo, vamos enunciar um con-junto de condicoes que, uma vez satisfeitas, garantem que, se a aproximacaoinicial x0 do metodo de Newton, pertencer a um certo intervalo, entao ometodo converge para a raiz da equacao que se encontra nesse intervalo.
Teorema 2.7 Seja f uma funcao de variavel real que satisfaz as seguintescondicoes:
1. f f e contınua em [a, b];f(a)f(b) < 0.
2. f ∈ C1([a, b]),f ′(x) 6= 0 em [a,b].
3. f ∈ C2([a, b]),f ′′(x) ≥ 0 ou f ′′(x) ≤ 0 em [a, b].
4. |f(a)||f ′(a)| < b− a, |f(b)|
|f ′(b)| < b− a.
Entao, qualquer que seja a aproximacao inicial x0 ∈ [a, b], o metodo deNewton converge para a unica raiz z de f em [a, b].
Nalgumas situacoes tem interesse tambem a seguinte variante deste teo-rema:
Teorema 2.8 Suponhamos que f satisfaz as primeiras 3 condicoes doteorema 2.7 . Entao, se a aproximacao inicial x0 for tal que f(x0)f ′′(x) ≥0,∀x ∈ [a, b] , o metodo de Newton converge para a unica raiz z de f em[a, b] e a sucessao das iteradas e monotona.
Nao iremos fazer a demonstracao completa destes teoremas, mas apenasinvestigar o significado e a razao de ser de cada uma das condicoes. As
35
primeiras condicoes, como sabemos pelos Teoremas 2.1 e 2.2, garantem quea funcao considerada tem um unico zero em [a, b]. Alem disso, a segundacondicao e essencial para o metodo de Newton, pois se ela nao se verificar(isto e, se a derivada de f se anular nalgum ponto de [a, b]), o metodo deNewton pode nao ser aplicavel. Quanto a terceira condicao, ela significaque no domınio considerado a segunda derivada de f nao muda de sinalou, por outras palavras, a funcao nao tem pontos de inflexao no intervaloconsiderado. Para entendermos a razao de ser desta condicao, analisemos oseguinte exemplo.
Exemplo 2.11. Consideremos f(x) = x3 − x, no intervalo [−0.5, 05].Neste intervalo, a funcao e continuamente diferenciavel, com f ′(x) = 3x2−1.Alem disso, f tem sinais opostos nos extremos do intervalo (f(−0.5) =3/8, f(0.5) = −3/8) e f ′ nao se anula (e sempre negativa). Por con-seguinte, as duas primeiras condicoes do teorema 2.8 estao satisfeitas nointervalo [−0.5, 0.5]. Em relacao a terceira condicao, temos f ′′(x) = 6x,logo, f ′′(x) muda de sinal em x = 0, pelo que esta condicao nao e satisfeita.Vejamos agora atraves de um exemplo que a convergencia do metodo deNewton nao esta garantida para qualquer aproximacao inicial no intervalo[−0.5, 0.5]. Seja x0 = 1/
√5. Embora este ponto pertenca ao intervalo con-
siderado, verifica-se imediatamente que as iteradas do metodo de Newton,neste caso, formam uma sucessao divergente: x1 = −1/
√5, x2 = 1/
√5, x3 =
−1/√
5, . . . .Vejamos agora um exemplo que ilustra a importancia da condicao 4.Exemplo 2.12. Seja f(x) = lnx, que tem uma unica raiz em z = 1. Se
considerarmos, por exemplo, o intervalo [0.5, 3] vemos que neste intervaloestao satisfeitas as primeiras 3 condicoes dos teoremas 2.7 e 2.8:
1. f(0.5)f(3) < 0;
2. f ′(x) = 1/x 6= 0, ∀x ∈ [0.5, 3];
3. f ′′(x) = −1/x2 < 0, ∀x ∈ [0.5, 3].
No entanto, a convergencia do metodo de Newton nao esta asseguradapara qualquer aproximacao inicial neste intervalo. Se tomarmos, por exem-plo, x0 = 3, temos x1 = 3− 3 ln(3) < 0, pelo que o metodo de Newton naopode ser aplicado (f(x) nao esta definida para x < 0).
Neste caso, e facil ver que falha a condicao 4 do teorema 2.7. Com efeito,temos
|f(3)||f ′(3)|
= 3ln(3) > 3− 0.5 = 2.5.
36
Para que o metodo de Newton convergisse para a raiz procurada, nestecaso, deverıamos escolher, por exemplo, x0 = 0.5. Nesse caso, a convergenciaresultaria do Teorema 2.8, ja que se verifica f(0.5)f ′′(x) > 0, ∀x ∈ [0.5, 3].
Sobre o significado geometrico da condicao 4 podemos dizer o seguinte:se ela se verificar, entao, tomando x0 = a, a iterada x1 satisfaz |x1 − a| =|f(a)||f ′(a)| < |b − a| (ou seja, a distancia de x1 a a e menor que o comprimento
do intervalo [a, b]). Logo, x1 pertence a esse intervalo. Continuando esteraciocınio, pode mostrar-se que todas as iteradas seguintes continuam apertencer ao intervalo [a, b]. Se comecarmos o processo iterativo a partir de
x0 = b e utilizarmos a condicao |f(b)||f ′(b)| < |b − a|, um raciocınio semelhante
leva-nos a mesma conlusao, isto e, a condicao 4 garante que se x0 ∈ [a, b],todas as iteradas do metodo de Newton se mantem dentro desseintervalo.
2.4.4 Ordem de convergencia do metodo de Newton
O metodo de Newton tambem pode ser encarado como um caso particulardo metodo do ponto fixo. Essa abordagem tem a vantagem de nos permi-tir analisar a convergencia do metodo de Newton com base nos resultadosteoricos sobre o metodo do ponto fixo.
Consideremos a equacao f(x) = 0 e suponhamos que ela tem uma unicaraiz num certo intervalo [a, b]. Admitamos ainda que f ∈ C1([a, b]) e quef ′(x) 6= 0, ∀x ∈ [a, b]. Entao a equacao considerada e equivalente a
x− f(x)
f ′(x)= x. (2.34)
Se definirmos a funcao iteradora
g(x) = x− f(x)
f ′(x),
podemos dizer que a equacao (2.34) e a equacao dos pontos fixos de g. Logo,as raizes de f , que tambem sao pontos fixos de g, podem ser aproximadaspelo processo iterativo
xk+1 = g(xk) = xk −f(xk)
f ′(xk). (2.35)
Verificamos portanto que este metodo e identico ao metodo Newton, apli-cado a funcao f(x). Logo, para determinar a ordem de convergencia dometodo de Newton, basta determinar, com base no teorema 2.6, a ordem
37
de convergencia da sucessao gerada por esta funcao iteradora. Para isso,comecemos por calcular as suas derivadas. Temos
g′(x) =f(x)f ′′(x)
f ′(x)2.
Daqui, tomando em consideracao que f(z) = 0 e f ′(z) 6= 0, resulta queg′(z) = 0. Quanto a segunda derivada de g, temos
g′′(x) =(f ′(x)f ′′(x) + f(x)f ′′′(x)) f ′(x)2 − f(x)f ′′(x)(f ′(x)2)′
f ′(x)4.
Logo,
g′′(z) =f ′′(z)
f ′(z).
Podemos assim concluir o seguinte:a) Se f ′′(z) 6= 0, entao g′′(z) 6= 0 (uma vez que, por condicao, f ′(z) 6= 0).
Nesse caso, de acordo com o teorema 2.6, o metodo de Newton (ou seja,
o metodo do ponto fixo com a funcao iteradora g(x) = x − f(x)f ′(x)) tem or-
dem de convergencia 2 (convergencia quadratica). Alem disso, o coeficienteassimptotico de convergencia e dado por
k∞ =|f ′′(z)|2|f ′(z)|
.b) Se f ′′(z) = 0, entao g′′(z) = 0 e o metodo de Newton tem ordem de
convergencia, pelo menos, 3 (para saber qual a ordem concreta e necessarioanalisar as derivadas de ordem superior de g.
Exemplo 2.13. Consideremos a equacao f(x) = x3 − x = 0. Umadas suas raizes e z = 0. Se aplicarmos o metodo de Newton ao calculoaproximado desta raiz, isso equivale a utilizar o metodo do ponto fixo coma funcao iteradora
g(x) = x− f(x)
f ′(x)=
2x3
3x2 − 1.
Analisemos a ordem do metodo neste caso. Para isso comecemos por verificarque f ′(0) = −1 6= 0 e f ′′(0) = 0. Entao, de acordo com a analise queacabamos de realizar, o metodo deve ter ordem, pelo menos, 3. Vejamos
qual e, de facto, a ordem de convergencia. Sabemos que g′′(0) = f ′′(0)f ′(0) = 0.
Para determinar g′′′(0), observemos que a funcao g admite, em torno dez = 0, um desenvolvimento em serie de Taylor da forma
g(x) =−2x3
1− 3x2= −2x3 +O(x5),
38
de onde se conclui que g′′′(x) = −12 + O(x2)(quando x → 0), pelo queg′′′(0) = −12. Temos, portanto, convergencia de ordem 3 . O coeficenteassimptotico de convergencia, de acordo com o Teorema 2.6, e
k∞ =|g′′′(0)|
3!= 2.
2.5 Metodo da Secante
Terminaremos este capıtulo com a analise do metodo da secante. Tal comono caso do metodo de Newton, a formula iteradora deste metodo pode serdeduzida a partir da sua interpretacao geometrica.
2.5.1 Interpretacao geometrica do metodo da secante
Seja f uma funcao de variavel real, contınua num certo intervalo [a, b], esuponhamos que f tem nesse intervalo um unico zero z. Para aplicar ometodo da secante, escolhem-se dois numeros, x0 e x1, no intervalo [a, b], econsidera-se a recta que passa pelos pontos (x0, f(x0)) e (x1, f(x1)) (secanteao grafico de f). A equacao dessa recta e
y − f(x1) =f(x1)− f(x0)
x1 − x0(x− x0). (2.36)
Depois, determina-se o ponto onde esta recta intersecta o eixo das abcissas.Esse ponto existe, desde que f(x0) 6= f(x1), condicao que consideramossatisfeita. Desigando por x2 a abcissa desse ponto, obtem-se a seguinteequacao para x2:
x2 = x1 −x1 − x0
f(x1)− f(x0)f(x1). (2.37)
Deste modo, a nova aproximacao da raiz x2 fica definida a partir de x0 e x1.A formula que nos permite determinar cada aproximacao xk+1 a partir dasduas anteriores (xk e xk−1) e analoga a (2.37):
xk+1 = xk −xk − xk−1
f(xk)− f(xk−1)f(xk), k = 1, 2, . . . (2.38)
Um exemplo de aplicacao do metodo da secante esta representado na fig.12.
39
2.5 2.75 3.25 3.5 3.75 4
-0.5
0.5
1
1.5
2
2.5
x0 x1x2 x3 x4z
Figure 12: Interpretacao geometrica do metodo da secante.
2.5.2 Estimativa do erro do metodo da secante
No caso do metodo de Newton, vimos que o erro de cada iterada pode serestimado a partir do erro da iterada anterior e das propriedades da funcaof . No caso do metodo da secante, e de realcar uma diferenca fundamental:e que cada iterada depende das duas iteradas anteriores, e nao apenas daultima. Neste caso, diz-se que temos um metodo iterativo a dois passos.Sendo assim, e natural que o erro de cada iterada do metodo da secante possaser determinado a partir dos erros das duas ultimas iteradas. Suponhamosentao que xm−1 e xm sao duas iteradas consecutivas do metodo da secante.A iterada seguinte, xm+1, pode ser determinada atraves da formula (2.38).Representemos os erros de xm−1 e xm por em−1 e em, respectivamente; istoe, em−1 = z−xm−1 e em = z−xm. Alem disso, suponhamos que a funcao f eduas vezes continuamente diferenciavel num intervalo I ⊂ [a, b] que contemxm−1, xm, xm+1 e z e que f ′ nao se anula em I. Entao, usando resultadosda teoria da interpolacao (ver, por exemplo, [2]), pode mostrar-se que em+1
(erro de xm+1) satisfaz a desigualdade:
em+1 = − f′′(ξm)
2f ′(ηm)emem−1, (2.39)
onde ξm e ηm representam pontos que pertencem ao intervalo I acimareferido. Note-se que a formula (2.39) e muito parecida com a formula(2.29) para o erro do metodo de Newton. A diferenca e que, como seria deesperar, neste caso, o erro da nova iterada do metodo da secante e avaliadoa partir do produto dos erros das duas ultimas iteradas, enquanto que nometodo de Newton o erro da nova iterada e avaliado a partir do quadradodo erro da iterada anterior.
Na pratica, para usar a formula (2.39), a semelhanca do que fizemosno caso do metodo de Newton, convem majorar em I o modulo da segunda
40
derivada de f e minorar o modulo da sua primeira derivada. Para simplificar,suponhamos que I = [a, b] e seja
M = maxx∈[a,b]
|f ′′(x)|,
µ = minx∈[a,b]
|f ′(x)|.
Entao, da formula (2.39) resulta imediatamente a seguinte majoracao parao erro absoluto do metodo da secante:
|em+1| ≤M
2µ|em||em−1|. (2.40)
Normalmente, os erros absolutos das suas iteradas iniciais, |e0| e |e1|, saomajorados pelo comprimento do intervalo [a, b], isto e, sao evidentes as de-sigualdades |e0| < |b−a| e |e1| < |b−a|. A partir daı, os erros das sucessivasiteradas sao majorados por recorrencia, isto e: o erro |e2| majora-se a partirdos erros |e0| e |e1|; o erro |e3| majora-se a partir dos erros |e1| e |e2|; e assimsucessivamente.
Exemplo 2.14. Consideremos mais uma vez a equacao f(x) = cos(x)−2x = 0, que tem uma raiz no intervalo [0.4, 0.5]. Para aproximar essa raizpretende-se usar o metodo da secante.
1. Tomando como aproximacoes iniciais os pontos x0 = 0.5 e x1 = 0.4,calcule as iteradas x2 e x3 pelo metodo da secante.
Aplicando a formula (2.38), temos
x2 = x1 −x1 − x0
f(x1)− f(x0)f(x1) = 0.449721; (2.41)
x3 = x2 −x2 − x1
f(x2)− f(x1)f(x2) = 0.450188. (2.42)
2. Obtenha majorantes do erro absoluto de x0, x1, x2 e x3.O caminho mais facil seria majorar |e0| e |e1| pelo comprimento do in-
tervalo: |b − a| = |0.5 − 0.4| = 0.1. O majorante pode, no entanto, ser umpouco melhorado, se tivermos em conta o sinal de f em cada um dos pontosxi calculados. Para tal, observemos a tabela:
i xi f(xi)
0 0.5 −0.1221 0.4 0.1212 0.449721 0.00113 0.450188 −0.00001
41
Da observacao da tabela conclui-se que os pontos x1 e x2 se encontram aesquerda da raiz z (onde f e positiva), enquanto x0 e x3 se encontram adireita (onde f e negativa). Sendo assim, para os erros de x0 e x1 obtem-seos seguintes majorantes:
|e0| = |z − x0| ≤ |x2 − x0| = |0.449721− 0.5| = 0.050258,
|e1| = |z − x1| ≤ |x3 − x1| = |0.450188− 0.4| = 0.050188.
Recordemos do exemplo 2.10 que neste caso se tem M = 0.921, µ = 2.389.Daqui, pela estimativa (2.40), obtem-se
|e2| ≤M
2µ|e1||e0| ≤ 0.193× 0.050188× 0.050258 = 0.4868× 10−3, (2.43)
|e3| ≤M
2µ|e2||e1| ≤ 0.193×0.4868×10−3×0.050188 = 0.4715×10−5. (2.44)
Vemos assim que, ao fim de duas iteracoes, o metodo da secante nosproporciona uma aproximacao com um erro da ordem de 10−5. No caso demetodo de Newton, com o mesmo numero de iteracoes, obtem-se um erro daordem de 10−7 (ver exemplo 2.10). Assim, este exemplo sugere que o metodode Newton converge mais rapidamente que o da secante. Por outro lado,no mesmo exemplo, a precisao que se consegue obter com duas iteradas dometodo do ponto fixo e da ordem de 10−2. Logo, a julgar por este exemplo,e de esperar que a ordem de convergencia do metodo da secante esteja entrea ordem do metodo do ponto fixo no caso geral(1) e a do metodo de Newton(2). Esta conjectura e confirmada pelos resultados que vamos expor naproxima seccao.
2.5.3 Convergencia do metodo da secante
Com base na estimativa do erro que foi deduzida no paragrafo anterior podeprovar-se o seguinte teorema sobre a convergencia do metodo da secante (verdemonstracao em [2] ).
Teorema 2.9 Seja f uma funcao duas vezes continuamente diferenciavelnuma vizinhanca de z e tal que f ′(z) 6= 0. Entao, se os valores iniciais x0 e x1
forem suficientemente proximos de z, a sucessao {xm} gerada pelo metododa secante converge para z.
O exemplo numerico que foi analisado na seccao anterior sugere que ometodo da secante e mais rapido que o metodo do ponto fixo (que , nocaso geral, tem ordem 1), mas menos rapido que o de Newton, que temconvergencia quadratica. Assim, se {xm} for uma sucessao gerada pelo
42
metodo da secante, existe um numero real p, tal que 1 < p < 2, para o qualse verifica
limm→∞
|z − xm+1||z − xm|p
= K∞, (2.45)
onde K∞ e uma constante positiva. (De acordo com a designacao intro-duzida na seccao 2.3.8, esta constante designa-se coeficiente assimptotico deconvergencia.) Mais precisamente, pode provar-se (ver detalhes em [1]) que,no caso do metodo da secante, temos
p =1 +√
5
2≈ 1.6,
isto e, a ordem de convergencia e dada pelo chamado numero de ouro (sobre aimportancia dessse numero e as suas numerosas aplicacoes ver, por exemplo,[3] ).
O Teorema 2.9, acima enunciado, tem a desvantagem de nao ser facil-mente aplicavel. Na realidade, o que significa a condicao ”se x0 e x1 foremsuficientemente proximos de z?” Para a pratica, sao bastante mais uteisproposicoes, do tipo dos Teoremas 2.7 e 2.8, que nos dem condicoes su-ficientes para a convergencia do metodo da secante, desde que as aprox-imacoes iniciais pertencam a um dado intervalo. Passamos a enunciar taisteoremas.
Teorema 2.10. Se forem satisfeitas as condicoes do teorema 2.7, entaoo metodo da secante converge para a raiz z de f em [a, b], quaisquer quesejam as aproximacoes iniciais x0,x1, pertencentes a [a, b].
Teorema 2.11. Se forem satisfeitas as primeiras 3 condicoes do teorema2.7, e se as aproximacoes iniciais satisfizerem f(x0)f ′′(x) ≥ 0,f(x1)f ′′(x) ≥0, ∀x ∈ [a, b], entao o metodo da secante converge para a raiz z de f em[a, b].
3 Metodos Numericos para Sistemas de Equacoes
3.1 Normas Matriciais
Neste capıtulo trataremos de metodos computacionais para a resolucao desistemas de equacoes (lineares e nao lineares). Para a analise do erro destesmetodos, necessitaremos frequentemente de recorrer a normas vectoriais ematriciais, pelo que comecaremos por fazer uma breve introducao sobre estetema.
Seja E um espaco linear. Tipicamente, teremos E = IRn (vectores reaisde dimensao n) ou E = IRn×n (matrizes reais de n linhas e n colunas).
43
Definicao 3.1 . Uma aplicacao φ de E em IR+0 diz-se uma norma se
satisfizer as seguintes condicoes:
1. φ(x) ≥ 0, ∀x ∈ E, sendo φ(x) = 0, se e so se x = 0.
2. φ(λx) = |λ|φ(x), ∀x ∈ E, λ ∈ IR.
3. φ(x+ y) ≤ φ(x) + φ(y),∀x, y ∈ E.
Comecaremos por rever alguns exemplos de normas em IRn. Como habit-ualmente, representaremos qualquer elemento de IRn por x = (x1, x2, . . . , xn),onde xi ∈ IR.
1. Norma do maximo. φ(x) = ‖x‖∞ = maxi=1,...,n |xi|.2. Norma 1. φ(x) = ‖x‖1 =
∑ni=1 |xi|.
3. Norma euclidiana. φ(x) = ‖x‖2 =√∑n
i |xi|2.
4. Norma p. φ(x) = ‖x‖p = (∑n
i |xi|p)1p , p ≥ 1.
Note-se que a norma 1 e a norma euclidiana sao casos particulares dasnormas p, respectivamente, com p = 1 e p = 2. A norma do maximo e umcaso limite, quando p→∞. Pode provar-se que todos os exemplos indicadosdefinem normas, isto e, satisfazem as 3 condicoes acima formuladas.
Passamos agora a considerar o caso de E = IRn×n. Neste caso, os ele-mentos de E sao matrizes de n linhas e n colunas, por exemplo:
A =
a11 a12 . . . a1n
a21 a22 . . . a2n
... ... ... ...an1 an2 . . . ann
.Represente-se por ‖.‖v uma norma qualquer em IRn. Entao podemos
definir uma norma ‖.‖M em E pela seguinte igualdade:
‖A‖M = supx∈IRn,x 6=0‖Ax‖v‖x‖v
. (3.1)
Neste caso, diz-se que a norma ‖.‖M e a norma matricial,induzida pelanorma vectorial ‖.‖v.
Esta forma de construir normas matriciais permite-nos associar umanorma matricial a cada uma das nomas vectoriais introduzidas. As normasmatriciais introduzidas deste modo gozam de duas propriedades muito uteis:
1. A norma ‖.‖M e compatıvel com a norma ‖.‖v, isto e,
‖Ax‖v ≤ ‖A‖M ‖x‖v, ∀x ∈ IRn, ∀A ∈ IRn×n. (3.2)
44
Esta propriedade e uma consequencia imediata da formula (3.1).2. A norma ‖.‖M e regular, isto e,
‖AB‖M ≤ ‖A‖M ‖B‖M , ∀A,B ∈ IRn×n. (3.3)
Sao as seguintes as normas matriciais, induzidas pelas normas vectoriaismais correntes.
1. A norma matricial induzida pela norma do maximo chama-se normapor linha e e definida pela seguinte formula:
‖A‖∞ = maxi=1,...,n
n∑j=1
|aij |. (3.4)
2. A norma matricial induzida pela norma 1 chama-se norma por col-una e e definida pela seguinte formula:
‖A‖1 = maxj=1,...,n
n∑i=1
|aij |. (3.5)
3. A norma matricial induzida pela norma euclidiana e definida pelaseguinte formula:
‖A‖2 =√ρ(ATA). (3.6)
Na formula (3.6) o sımbolo ρ(A) representa o raio espectral de A, quese define como o maximo dos modulos dos valores proprios de A:
ρ(A) = maxi=1,...,n
|λi|. (3.7)
Note-se que, se A for uma matriz simetrica, isto e, se AT = A, entaoverifica-se
‖A‖2 =√ρ(ATA) =
√ρ(A2) = ρ(A), (3.8)
isto e, o raio espectral de A coincide com a sua norma euclidiana.Exemplo 3.1. Consideremos a matriz
A =
2 1 −31 3 42 −1 3
.
45
Neste caso, as normas de A sao
‖A‖∞ = max(6, 8, 6) = 8;‖A‖1 = max(5, 5, 10) = 10;
‖A‖2 =√
36.3 ≈ 6.02;
(3.9)
Para calcular ‖A‖2, foram necessarios os seguintes calculos auxiliares:
B = AT A =
9 3 43 11 64 6 34
. (3.10)
Os valores proprios de B sao λ1 = 6.8,λ2 = 10.9, λ3 = 36.3.Para terminar este exemplo, calculemos o raio espectral de A. Os valores
proprios de A sao λ1 = 3.69,λ2,3 = 2.15 ± i3.07. Logo,|λ2| = |λ3| = 3.75.Por conseguinte, ρ(A) = max(3.69, 3.75) = 3.75.
Note-se que, no exemplo 3.1, qualquer das normas de A e maior que oraio espectral da matriz. Tal nao acontece por acaso. Para terminar esteparagrafo, vamos verificar uma relacao existente entre as normas e o raioespectral das matrizes.
Teorema 3.1. Seja A ∈ IRn×n. Entao, qualquer que seja a normamatricial p, associada a uma certa norma vectorial em IRn, verifica-se
rσ(A) ≤ ‖A‖p. (3.11)
Demonstracao. Seja x um vector proprio de A, associado ao valor proprioλ, tal que |λ| = rσ(A). Entao
‖Ax‖p = ‖λx‖p = |λ| ‖x‖p. (3.12)
Entao podemos afirmar que
‖A‖p = supx∈IRn,x 6=0
‖Ax‖p‖x‖p
≥ |λ|, (3.13)
de onde resulta a afirmacao do teorema.
46
3.2 Condicionamento de Sistemas Lineares
Como vimos no capıtulo 1, um dos aspectos mais importantes a ter em con-sideracao quando se quando se analisam metodos numericos para a solucaode um determinado problema e a sensibilidade desses metodos em relacao apequenos erros nos dados. Se for dado um certo sistema linear
Ax = b
os dados sao o segundo membro do sistema ( b ∈ IRn) e a matriz (A ∈ IRn×n),que se supoe nao singular. Vamos, portanto, analisar ate que ponto umpequeno erro, em termos relativos, do vector b ou da matriz A, pode afectara solucao do sistema. Para isso, representemos por A uma perturbacao damatriz A.
A ≈ A. (3.14)
Analogamente, representemos por b uma perturbacao do segundo membrodo sistema:
b ≈ b. (3.15)
Se substituirmos A por A e b por b no sistema inicial, obteremos um novosistema, cuja solucao representaremos por x.
Como vimos no capıtulo 1, para medir a precisao dos resultados numericose essencial o conceito de erro relativo. Vamos, portanto, designar por errorelativo de um vector x (numa certa norma vectorial p) o quociente
‖δx‖p =‖x− x‖p‖x‖p
; (3.16)
analogamente, designaremos por erro relativo de uma matriz A (na normamatricial p)o quociente
‖δA‖p =‖A−A‖p‖A‖p
. (3.17)
O nosso objectivo e, escolhendo uma certa norma vectorial e a norma ma-tricial associada, estimar o erro relativo ‖δx‖p em funcao dos erros relativos‖δb‖p e ‖δA‖p.
Definicao 3.2. Um sistema linear diz-se bem condicionado se e so sea pequenos erros relativos do segundo membro e da matriz correspondempequenos erros relativos da solucao.
47
Para analisarmos o problema do condicionamento, comecemos por con-siderar o caso mais simples em que ‖δA‖p = 0. Nesse caso, temos
A x = b (3.18)
De (3.18) obtem-sex− x = A−1 (b− b) (3.19)
Por conseguinte, qualquer que seja a norma vectorial escolhida, e valida aseguinte estimativa para o erro absoluto de x:
‖x− x‖p ≤ ‖A−1‖p ‖b− b‖p. (3.20)
Por outro lado, da igualdade Ax = b obtem-se
‖b‖p ≤ ‖A‖p ‖x‖p, (3.21)
logo1
‖x‖p≤ ‖A‖p‖b‖p
. (3.22)
Multiplicando cada um dos membros de 3.20 pelo membro correspondentede 3.22, obtem-se
‖x− x‖p‖x‖p
≤ ‖A‖p ‖A−1‖p‖b− b‖p‖b‖p
. (3.23)
Obtivemos assim a estimativa que procuravamos para o erro relativo dasolucao, em funcao do erro relativo do segundo membro. A desigualdade(3.23) leva-nos a definicao de numero de condicao de uma matriz.
Definicao 3.3. Seja A uma matriz invertıvel. Chama-se numero decondicao de A (na norma p) o seguinte valor
condp(A) = ‖A‖p ‖A−1‖p. (3.24)
Como resulta da desigualdade 3.23, o numero de condicao da matriz Ada -nos a relacao entre o erro relativo da solucao de um sistema linear e oerro relativo do seu segundo membro. Assim, se o numero de condicaode A for alto, daı resulta que pequenos erros relativos do segundomembro podem provocar erros muito significativos na solucao, oque significa, por definicao, que o sistema e mal condicionado.
48
Note-se, entretanto, que o numero de condicao de uma matriz e sempremaior ou igual a 1, desde que consideremos normas matriciais induzidas(ver problema 3-b no fim desta seccao). Por conseguinte, um sistema bemcondicionado e aquele que tem o numero de condicao nao muito maior que 1.
Tambem se usa a definicao do numero de condicao atraves do raio es-pectral, sendo nesse caso dado pela expressao
cond∗(A) = rσ(A) rσ(A−1). (3.25)
De acordo com o teorema 3.1, podemos escrever
cond∗(A) ≤ condp(A), (3.26)
qualquer que seja a norma matricial p considerada (p ≥ 1). Daqui re-sulta que, se o numero de condicao cond∗(A) for alto, todos os numeros decondicao da matriz sao altos, pelo que o sistema e mal condicionado. Noentanto, pode acontecer que o sistema seja mal condicionado, mesmo que onumero de condicao cond∗(A) seja pequeno (ver problema 1 desta seccao).
Atendendo a que os valores pro prios de A−1 sao os inversos dos valorespro prios de A, o numero de condicao cond∗(A) pode escrever-se sob a forma
cond∗(A) =maxλi∈σ(A) |λi|minλi∈σ(A) |λi|
. (3.27)
No caso de a matriz A ser simetrica, entao, como foi observado, a sua normaeuclidiana coincide com o raio espectral, pelo que podemos escrever:
cond2(A) = cond∗(A). (3.28)
Consideremos agora o caso mais geral em que o sistema linear podeestar afectado de erros, nao so no segundo membro, mas tambem na propriamatriz. Sobre este caso, e conhecido o seguinte teorema.
Teorema 3.2. Consideremos o sistema linear Ax = b, onde A e umamatriz invertıvel. Sejam δA e δb definidos pelas igualdades (3.17) e (3.16),respectivamente, e suponhamos que
‖A− A‖p ≤1
‖A−1‖p. (3.29)
Entao verifica-se a seguinte desigualdade:
‖δx‖p ≤cond(A)
1 − condp(A) ‖δA‖p{‖δA‖p + ‖δb‖p} . (3.30)
49
Observacao. Note-se que a desigualdade 3.23, obtida anteriormente, eum caso particular de 3.30, que se obtem fazendo ‖δA‖p = 0. A demon-stracao do teorema 3.2 encontra-se, por exemplo, em [1].
A desigualdade (3.30) confirma a nossa conclusao de que os sistemas line-ares com numeros de condicao altos sao mal condicionados. O exemplo quese segue mostra como os problemas de mau condicionamento podem surgirmesmo em sistemas de pequenas dimensoes, com matrizes aparentemente”bem comportadas”.
Exemplo 3.2 Consideremos o sistema linear Ax = b, onde
A =
10 7 8 77 5 6 58 6 10 97 5 9 10
, b =
32233331
. (3.31)
Verifica-se imediatamente, por substituicao, que a solucao deste sistema ex = (1, 1, 1, 1)T . Sabe-se que a matriz A e nao singular. Alem disso, ela e,evidentemente, simetrica e a sua norma, por linhas ou por colunas, e
‖A‖∞ = ‖A‖1 = max(32, 23, 33, 31) = 33.
Entretanto, se substituirmos o vector b por um vector b, tal que
b = (32.1, 22.9, 33.1, 30.9)T ,
a solucao do sistema passa a ser
x = (9.2,−12.6, 4.5,−1.1)T ,
o que e totalmente diferente da solucao do sistema inicial. Por outraspalavras, uma perturbacao do segundo membro, tal que
‖δb‖∞ =0.1
33≈ 0, 3%
leva-nos a uma nova solucao, cuja norma e cerca de 13 vezes maior que a dasolucao original.
50
Observemos ainda o que acontece se a matriz A sofrer uma ligeira per-turbacao das suas componentes, sendo substituıda por
A =
10 7 8.1 7.2
7.08 5.04 6 58 5.98 9.89 9
6.99 5 9 9.98
, (3.32)
mantendo-se o segundo membro inalterado. Neste caso, a solucao do sistemapassa a ser
x = (−81, 137,−34, 22)T .
Neste caso, a diferenca em relacao a solucao inicial e ainda mais acentuada.Entretanto, a norma da perturbacao e relativamente pequena:
‖A− A‖∞ = max(0.3, 0.12, 0.13, 0.03) = 0.3,
de onde resulta
‖δA‖∞ =0.3
33≈ 0, 9%.
Vejamos como interpretar estes factos, com base na teoria do condi-cionamento que expusemos nesta seccao. Para isso, precisamos de conhecera inversa de A :
A−1 =
25 −41 10 −6−41 68 −17 1010 −17 5 −3−6 10 −3 2
. (3.33)
Podemos imediatamente constatar que
‖A−1‖∞ = max(82, 136, 35, 21) = 136.
Assim, o numero de condicao de A, segundo a norma ∞ (que coincide coma norma 1 neste caso), tem o valor
cond∞(A) = cond1(A) = 33× 136 = 4488.
Conhecendo o valor do numero de condicao, ja nao nos surpreende o facto deas pequenas perturbacoes que introduzimos no segundo membro e na matrizterem alterado completamente a solucao. Com efeito, a estimativa (3.23),aplicada a este caso, diz-nos que
‖δx‖ ≤ 4488× 0, 3% = 1346%,
51
o que explica inteiramente os resultados obtidos, no que diz respeito a per-turbacao do segundo membro do sistema. No caso das alteracoes produzidasna matriz A, nao se pode aplicar a estimativa 3.30, uma vez que, para a per-turbacao considerada, nao se verifica a condicao
‖A− A‖∞ ≤1
‖A−1‖∞.
No entanto, dado o alto valor do numero de condicao, seria de esperar, dequalquer modo, que a solucao do sistema ficasse muito alterada, tal como severificou.
3.2.1 Problemas sobre Condicionamento
1. Seja A uma matriz quadrada, de dimensao n, com a forma
A =
1 −1 . . . . . . −10 1 −1 . . . −1. . . . . . . . . . . . . . .0 . . . 0 1 −10 . . . . . . 0 1
.
(a) Calcule A−1.
(b) Determine os numeros de condicao cond1(A) e cond∞(A).
(c) Sejam b1 e b2 dois vectores de IRn tais que
‖δb‖∞ =‖b1 − b2‖∞‖b1‖∞
≤ 10−5.
Sejam x1 e x2, respectivamente, as solucoes dos sitemas Ax = b1e Ax = b2. Determine um majorante de
‖δx‖∞ =‖x1 − x2‖∞‖x1‖∞
no caso de n = 20. Comente.
2. Seja A a matriz real
A =
1 0 11 −1 0a 0 3
,
52
onde a ∈ IR. Suponhamos que, ao resolver o sistema Ax = b, com umcerto valor de a, se obteve a solucao x = (1, 1, 1). Supondo que o valorde a esta afectado de um certo erro, de valor absoluto nao superiora ε, determine um majorante de ‖∆x‖∞, onde ∆x e a diferenca entrea solucao obtida e a que se obteria se fosse conhecido o valor exactode a.
3.3 Metodos Directos
Para resolver um sistema linear, podemos optar por dois caminhos diferentes:
1. ou procuramos a solucao por um metodo de aproximacoes sucessivas,utilizando um metodo iterativo;
2. ou optamos por reduzir o sistema a uma forma mais simples, de modoa obter a solucao exacta atraves de simples substituicoes. Nesse caso,dizemos que estamos a aplicar um metodo directo.
Neste paragrafo, falaremos de metodos directos. Quando se utiliza ummetodo directo, sabe-se que o erro do metodo e nulo, visto que o metodoconduz teoricamente a solucao exacta do problema. Isto, porem nao querdizer que a solucao que se obtem, de facto, seja exacta, visto que, comosabemos, quando se efectuam calculos numericos sao inevitaveis os erros dearredondamento.
3.3.1 Metodo da Eliminacao de Gauss
Um dos metodos mais simples e mais conhecidos para a solucao de sistemaslineares e o metodo da eliminacao de Gauss. A ideia basica deste metodoconsiste em reduzir o sistema dado, Ax = b, a um sistema equivalente,U x = b′, onde U e uma matriz triangular superior. Este ultimo sistemapode ser resolvido imediatamente, por substituicao, comecando pela ultimaequacao. Assim, podemos dizer que a resolucao de um sistema pelo metodode Gauss se divide em tres etapas:
1. Reducao da matriz A a forma triangular superior;
2. Transformacao do segundo membro do sistema;
3. Resolucao do sistema com a matriz triangular superior.
Vejamos com mais pormenor em que consiste cada uma destas etapas eavaliemos, em termo de numero de operacoes aritmeticas, o volume doscalculos correspondentes.
53
1.Reducao da matriz A a forma triangular superior. Suponhamosque a matriz dada tem a forma
A =
a11 a12 . . . a1n
a21 a22 . . . a2n
. . . . . . . . . . . .an1 an2 . . . ann
. (3.34)
Admitindo que a11 6= 0, esta matriz pode ser transformada, somando acada linha (a comecar pela 2a), um multiplo da 1a. Assim, obtem-se umanova matriz A(1), com a forma:
A(1) =
a11 a12 . . . a1n
0 a(1)22 . . . a
(1)2n
. . . . . . . . . . . .
0 a(1)n2 . . . a
(1)nn .
. (3.35)
As componentes de A(1) obtem-se atraves da formula:
a(1)ij = aij − mi1 a1j , (3.36)
ondemi1 =
ai1a11
. (3.37)
Continuando este processo, podemos a seguir transformar a matriz A(1)
numa nova matriz A(2), que sera triangular superior, ate a 2a coluna. Ge-neralizando, vamos efectuar uma sucessao de transformacoes, em que cadatransformada A(k) tem a forma
A(k) =
a11 a12 . . . . . . . . . a1n
0 a(1)22 . . . . . . . . . a
(1)2n
. . . . . . . . . . . . . . . . . .
0 . . . 0 a(k−1)kk . . . a
(k−1)kn
. . . . . . . . . . . . . . . . . .
0 . . . 0 a(k)nk . . . a
(k)nn
. (3.38)
As componentes de A(k) obtem-se a partir das de A(k−1) atraves da formula:
a(k)ij = a
(k−1)ij − mik a
(k−1)kj , i = k + 1, . . . , n; j = k + 1, . . . , n; (3.39)
onde
mik =a
(k−1)ik
a(k−1)kk
(3.40)
54
(neste caso, pressupoe-se que a(k−1)kk 6= 0). Ao fim de n− 1 transformacoes
obtem-se
A(n−1) =
a11 a12 . . . . . . . . . a1n
0 a(1)22 . . . . . . . . . a
(1)2n
. . . . . . . . . . . . . . . . . .
0 . . . 0 a(k−1)kk . . . a
(k−1)kn
. . . . . . . . . . . . . . . . . .
0 . . . . . . . . . 0 a(n−1)nn
. (3.41)
No caso de algum dos coeficientes a(k−1)kk ser igual a 0, para se poder aplicar
a eliminacao de Gauss torna-se necessario alterar a ordem das linhas. Essecaso sera visto em detalhe mais adiante, durante a resolucao do Exemplo 3.3.Note-se que, se a matriz A for nao singular, existe sempre uma permutacaodas suas linhas, depois da qual ela pode ser reduzida a forma (3.41), comtodos os elementos da diagonal principal diferentes de 0.
2. Transformacao do segundo membro. O segundo membro dosistema e sujeito a transformacoes analogas as que se operam sobre a matriz,de modo a garantir a equivalencia do sistema resultante ao inicial. Assim,a transformacao do vector b tambem se realiza em n − 1 passos, sendo aprimeira transformada, b(1), obtida segundo a formula:
b(1)i = bi − mi1 b1, i = 2, 3, . . . , n. (3.42)
Analogamente, a k-esima transformada do segundo membro e obtida pelaformula:
b(k)i = bi − mik b
(k−1)k , i = k + 1, . . . , n. (3.43)
Os coeficientes mik sao dados pelas formulas (3.37) e (3.40).
3. Resolucao do sistema triangular superior. Depois de reduzido osistema inicial a forma triangular superior, com a matriz (3.41), a sua solucao
55
obtem-se facilmente, mediante um processo de substituicoes sucessivas:
xn = b(n−1)n
a(n−1)nn
;
xn−1 =b(n−2)n−1 − a
(n−2)n−1,n xn
a(n−2)n−1,n−1
;
. . . . . . . . .
x1 =b1 −
∑ni=2 a1,i xia1,1
.
(3.44)
Vejamos agora qual o numero de operacoes aritmeticas necessarias paraefectuar cada uma das etapas que acabamos de descrever.
1. Reducao da matriz a forma triangular superior. O numerode operacoes necessarias para a transformacao da matriz A esta relacionadocom o numero de vezes que sao aplicadas as formulas (3.36) e (3.39).
No 1o passo, a formula (3.36) e aplicada (n− 1)2 vezes. Isto implica quese realizem (n− 1)2 multiplicacoes e outras tantas adicoes (ou subtraccoes).Entretanto, para calcular os quocientes da formula (3.37), efectuam-se n−1divisoes. Todas estas operacoes continuam a efectuar-se nos passos seguintesda transformacao, mas em menor numero, de acordo com a quantidade decomponentes que sao alteradas em cada passo. Em geral, no k-esimo passoefectuam-se (n−k)2 multiplicacoes e outras tantas adicoes (ou subtraccoes),assim como n − k divisoes. Assim, o numero total de multiplicacoes efec-tuadas durante a transformacao da matriz, igual ao numero de adicoes (ousubtraccoes), e
M(n) = AS(n) =n−1∑k=1
(n− k)2 =n(n− 1)(2n− 1)
6, (3.45)
Quanto as divisoes, o numero total e
D(n) =n−1∑k=1
(n− k) =n(n− 1)
2. (3.46)
56
Assim, o numero total de operacoes efectuadas durante a transformacao damatriz e, em termos assimptoticos, para valores altos de n,
TO(n) = M(n) +AS(n) +D(n) ≈ 2n3
3+O(n2). (3.47)
2. Transformacao do segundo membro. Quando transformamoso vector b, aplicamos a formula (3.43) as suas componentes. Esta formulae aplicada n − k vezes durante o k-esimo passo da transformacao, o queimplica n−k multiplicacoes e outras tantas adicoes (ou subtraccoes). Assim,o numero total de multiplicacoes efectuadas durante a transformacao dosegundo membro, igual ao numero de adicoes (ou subtraccoes), e
M(n) = AS(n) =n−1∑k=1
(n− k) =n(n− 1)
2, (3.48)
Por conseguinte, o numero total de operacoes exigidas por esta etapa doscalculos e, em termos assimptoticos,
TO(n) = M(n) + AS(n) ≈ n2. (3.49)
3. Resolucao do sistema triangular. Para resolver o sistema triangu-lar, efectuamos a serie de substituicoes (3.44). Como resulta destas formulas,o numero total de multiplicacoes para resolver o sistema e n(n−1)/2, igualao numero total de adicoes (ou subtraccoes). Quanto ao numero de divisoes,e igual a n.
Assim, o numero total de operacoes efectuadas para resolver o sistematriangular e, em termos assimptoticos,
TO(n) = M(n) + AS(n) + D(n) ≈ n2. (3.50)
Exemplo 3.3. Consideremos o sistema linear Ax = b, onde
A =
2 1 3−2 −1 12 4 2
, b =
5−14
. (3.51)
Pretende-se reslolver este sistema pelo metodo da eliminacao de Gauss.Comecemos por reduzir A a forma triangular superior. O primeiro passo
57
para isso, como vimos, consiste em transformar a matriz A na matriz A(1).Neste caso, usando as formulas (3.36) e (3.37), obtem-se:
a(1)22 = a22 − m21 a12 = 0,
a(1)23 = a23 − m21 a13 = 4,
(3.52)
ondem21 =
a21
a11= −1. (3.53)
Verifica-se que o novo elemento da diagonal principal, a(1)22 , e nulo. Como
sabemos, neste caso nao e possıvel aplicar o metodo da eliminacao de Gaussdirectamente a matriz A. Vamos entao proceder a uma troca de linhas; maisprecisamente, troquemos a segunda linha com a terceira. Obtem-se assim onovo sistema A′ x = b′, onde
A′ =
2 1 32 4 2−2 −1 1
, b′ = 5
4−1
. (3.54)
Aplicando o metodo da eliminacao de Gauss a este sistema, usemos de novoas formulas (3.36) e (3.37), que nos dao
a(1)′
22 = a′22 − m′21 a12 = 4− 1 = 3
a(1)′
23 = a′23 − m′21 a13 = 2− 3 = −1,
a(1)′
32 = a′32 − m′31 a12 = −1 + 1 = 0,
a(1)′
33 = a′33 − m′31 a13 = 1 + 3 = 4,
(3.55)
onde
m′21 =a′21
a11= 1, (3.56)
m′31 =a′31
a11= −1. (3.57)
Obtemos assim a matriz transformada da forma
A′ =
2 1 30 3 −10 0 4
. (3.58)
Como esta matriz e triangular superior, nao e necessario transforma-la mais.
58
A segunda etapa da aplicacao do metodo da eliminacao de Gauss consisteem transformar o segundo membro do sistema, isto e, o vector b′. Para isso,utilizamos a formula (3.42), que neste caso nos da
b(1)′
2 = b′2 − m′21 b′1 = 4− 5 = −1;
b(1)′
3 = b′3 − m′31 b′1 = −1 + 5 = 4.
(3.59)
Obtemos assim o vector transformado b(1)′ = (5,−1, 4)T .Por ultimo, resta-nos resolver o sistema triangular superior A(1)′ x =
b(1)′ . Para isso, usamos a substituicao ascendente, isto e, comecamos pordeterminar x3 a partir da ultima equacao, para depois determinar x2 dasegunda, e x1 da primeira. Usando as formulas (3.44), obtem-se
x3 =b(1)′
3
a(1)33
= 1;
x2 =b(2)2 − a
(1)23 x3
a(1)22
= −1+12 = 0;
x1 = b1 − a13 x3 − − a12 x2a1,1
= 5−32 = 1.
(3.60)
Assim, obtemos a solucao do sistema: x = (1, 0, 1)T .
3.3.2 Influencia dos erros de arredondamento
Ate agora, ao falarmos do metodo de Gauss, nao entramos em consideracaocom com os erros cometidos durante os calculos. Na seccao 3.2 ja vimosque pequenos erros nos dados iniciais do sistema podem afectar muito asolucao, se a matriz for mal condicionada. Alem dos erros dos dados iniciais,que sao inevitaveis, ha que ter em conta tambem o erro computacional,resultante dos arredondamentos efectuados durante os calculos. Um dosinconvenientes do metodo de Gauss, assim como de outros metodos directos,de que falaremos adiante, consiste em que esses erros tem frequentementetendencia para se propagar durante os calculos, de tal modo que podemadquirir um peso muito grande na solucao, mesmo que o sistema seja bemcondicionado. No entanto, o efeito destes erros pode ser muito atenuado, sedurante os calculos for usada uma estrategia adequada, a chamada estrategiade pivot. Vamos ver em que consiste esta estrategia.
Ao falarmos da transformacao da matriz A, vimos que, para que ela sepudesse efectuar, e necessario que todos os elementos da diagonal principalda matriz triangular superior U sejam diferentes de 0. Estes elementos sao
59
representados por a(k−1)kk e sao chamados geralmente os pivots, dada a sua
importancia para a aplicacao do metodo de Gauss 1. Vimos tambem que,no caso de um dos pivots ser nulo, se podia mesmo assim aplicar o metodo,desde que se efectuasse uma troca de linhas na matriz. Se o pivot nao fornulo, mas proximo de 0, o metodo continua a ser teoricamente aplicavel,mesmo sem trocas de linhas. So que, ao ficarmos com um denominadormuito pequeno no segundo membro de (3.40), cria-se uma situacao em queos erros de arredondamento podem propagar-se de uma forma desastrosa.A estrategia de pivot tem por objectivo evitar que isto aconteca. Assim,em cada passo da transformacao da matriz, verifica-se o pivot e, caso seconsidere conveniente, efectua-se uma troca de linhas que substitui o pivotpor outra componente maior. A estrategia de pivot tem diversas variantes,das quais falaremos aqui de duas: a pesquisa parcial e a pesquisa total depivot.
Pesquisa parcial de pivot. Neste caso, em cada passo da trans-formacao da matriz, e inspeccionada a coluna k da matriz A(k−1), maisprecisamente, as componentes dessa coluna que se situam da diagonal prin-cipal para baixo. Seja
ck = maxk≤i≤n
∣∣∣a(k−1)ik
∣∣∣ . (3.61)
Se o maximo no segundo membro de (3.61) for atingido para i = k, issosignifica que o actual pivot e, em modulo, a maior componente daquelacoluna. Nesse caso, continuam-se os calculos normalmente. Se o maximofor atingido para um certo i 6= k, entao troca-se de lugar a linha k com alinha i e so depois se prosseguem os calculos. Evidentemente, ao fazer esssatroca, tambem se efectua uma permutacao correspondente nas componentesdo vector b.
Pesquisa total de pivot. De acordo com esta estrategia, mais radical,e inspeccionada nao so a coluna k da matriz A(k−1), mas tambem todas ascolunas subsequentes. Seja
ck = maxk≤i,j≤n
∣∣∣a(k−1)ij
∣∣∣ . (3.62)
Sejam i′ e j′, respectivamente, os valores dos ındices i e j para os quais eatingido o maximo no segundo membro de (3.62). Se i′ nao coincidir com
1em frances, pivot tem o significado de base, apoio
60
k, a linha i′ troca de lugar com a linha k. Se, alem disso, j′ nao coincidircom k, entao a coluna j′ tambem vai trocar de lugar com a coluna k (o quecorresponde a uma troca de lugar entre as incognitas xj′ e xk).
Comparando as duas variantes acima mencionadas da pesquisa de pivot,ve-se que a pesquisa total e bastante mais dispendiosa do que a parcial,uma vez que exige um numero de comparacoes muito maior. Entretanto, apratica do calculo numerico tem demonstrado que, na grande maioria doscasos reais, a pesquisa parcial conduz a resultados praticamente tao bonscomo os da total. Isto explica que, hoje em dia, a pesquisa parcial sejamais frequentemente escolhida quando se elaboram algoritmos baseados nometodo de Gauss.
O exemplo que se segue mostra ate que ponto os erros de arredondamentopodem influir na solucao de um sistema linear, quando e aplicado o metododa eliminacao de Gauss. Vamos ver tambem como a pesquisa parcial depivot pode contribuir para melhorar esta situacao.
Exemplo 3.4. Consideremos o sistema linear Ax = b, onde
A =
10−6 0 11 10−6 21 2 −1
, b =
132
. (3.63)
Se procurarmos resolve -lo, utilizando o metodo da eliminacao de Gauss,chegamos ao sistema equivalente U x = b′,onde
U =
10−6 0 10 10−6 2− 106
0 0 2× 1012 − 5× 106 − 1
, b′ =
13− 106
2× 1012 − 7× 106 + 2
.(3.64)
Na realidade, a matriz U e o vector b′ que vamos obter nao sao exactamenteos da formula (4.71), devido aos erros de arredondamento. Suponhamosque os calculos sao efectuados num computador em que os numeros saorepresentados no sistema decimal, com seis dıgitos na mantissa. Entao , emvez de U e b′, teremos a seguinte matriz e o seguinte vector:
U =
1.00000× 10−6 0 1.000000 1.00000× 10−6 −0.999998× 106
0 0 1.99999× 1012
, b =
1−0.999997× 106
1.99999× 1012
.(3.65)
61
Assim, ao resolvermos o sistema (3.65) por substituicao ascendente, obtemossucessivamente:
x3 = 1.99999× 1012
1.99999× 1012 = 1.00000;
x2 = −0.999997× 106 + 0.999998× 106 x3
1.00000× 10−6 = 1.00000× 106;
∗[0.5cm]x1 = 1.00000− 1.00000 x3
1.00000× 10−6 = 0.
(3.66)Substituindo os valores assim obtidos no sistema dado, verifica-se que elesestao longe de o satisfazer, o que indica que este resultado apresenta umerro relativo muito grande. Este erro, no entanto, nao tem a ver com ocondicionamento do sistema, visto que o numero de condicao de A tem ovalor
cond∞(A) = ‖A‖∞‖A−1‖∞ ≈ 3× 4 = 12;
logo, o sistema nao se pode considerar mal condicionado. Nestas condicoes,temos razoes para pensar que o mau resultado obtido resulta da instabilidadenumerica do metodo, a qual, como vimos, pode ser contrariada atraves dapesquisa de pivot. Vejamos entao que resultado obterıamos, se aplicassemosa pesquisa parcial de pivot. Neste caso, comecarıamos por trocar a primeiralinha de A com a segunda, visto que a21 > a11. Depois da primeira trans-formacao, obterıamos a matriz A(1), com a seguinte forma:
A(1) =
1 10−6 20 −10−12 1− 2× 10−6
0 2− 10−6 −3
. (3.67)
Neste caso, a pesquisa de pivot leva-nos trocar a segunda linha com a ter-ceira, visto que a32 > a22. Depois de efectuar esta troca, realiza-se a se-gunda transformacao da matriz, que nos leva ao sistema A(2)x = b(2). Se oscalculos forem realizados com a precisao acima referida, obtemos a seguintematriz e o seguinte vector:
A(2) =
1.00000 1.00000× 10−6 2.000000 2.00000 −3.000000 0 9.99998× 10−1
, b(2) =
3.00000−1.00000
9.99997× 10−1
.(3.68)
62
Resolvendo o sistema (3.68), obtem-se a seguinte solucao:
x3 = 9.99997× 10−1
1.00000 = 9.99999× 10−1;
∗[0.5cm]x2 = −1.00000 + 3.00000 x32.00000 = 1.00000;
x1 = 3.00000 − 2.00000 x3 − 1.00000× 10−6 x2 = 9.99999× 10−1.(3.69)
Esta solucao e bastante diferente da que obtivemos, quando nao foi uti-lizada a pesquisa de pivot. Se substituirmos estes valores no sistema (3.63),veremos que a nova solucao esta correcta, dentro dos limites da precisaoutilizada. Este exemplo mostra-nos como a pesquisa de pivot pode desem-penhar um papel essencial na eliminacao dos problemas da instabilidadenumerica quando se resolvem sistemas lineares pelo metodo da eliminacaode Gauss.
3.3.3 Metodos de Factorizacao
Neste paragrafo, vamos tratar de metodos directos que se baseiam na facto-rizacao da matriz do sistema linear.
Definicao 3.4.Chama-se factorizacao LU de uma matriz A a sua repre-sentacao sob a forma do produto de de duas matrizes:
A = LU,
onde L e triangular inferior e U e triangular superior.Se for conhecida a factorizacao LU de uma matriz A, entao o sistema
linear Ax = b reduz-se a dois sistemas lineares com matrizes triangulares:
Lg = b,U x = g,
onde g e um vector auxiliar. Alem da solucao de sistemas lineares, a facto-rizacao LU tem outras aplicacoes , como por exemplo o calculo de determi-nantes. Com efeito, o determinante de A e igual ao produto dos determi-nantes de L e U , os quais se calculam imediatamente, dado estas matrizesserem triangulares:
det L = l11 l22 . . . lnn,det U = u11 u22 . . . unn.
63
Note-se que, para calcularmos o determinante a partir da definicao , terıamosde somar, para uma matriz de ordem n, n! parcelas, cada uma das quais eum produto de n componentes da matriz A. Tal calculo significaria, so parauma matriz de ordem 10, efectuar mais de 30 milhoes de multiplicacoes !Compreende-se portanto que tal maneira de calcular um determinante naoe aplicavel na pratica. Em compensacao, se utilizarmos a factorizacao, omesmo determinante pode ser calculado apenas com algumas centenas deoperacoes aritmeticas.
Uma vantagem dos metodos de factorizacao consiste em que, uma vezfactorizada uma matriz, se pretendermos resolver varios sistemas diferentescom essa matriz, basta resolver os sistemas triangulares correspondentes (asmatrizes L e U so precisam de ser determinadas uma vez). Isto e muitoimportante, dado que, como vamos ver, nos metodos de factorizacao a de-terminacao das matrizes L e U e precisamente a etapa mais dispendiosa,em termos de numero de operacoes .
O problema da factorizacao LU de uma matriz A ∈ Ln, nao singular,admite sempre multiplas solucoes. Com efeito, podemos abordar o problemacomo um sistema de n2 equacoes :
aij =
n∑k=1
lik ukj ; i = 1, . . . , n; j = 1, . . . , n; (3.70)
onde lik e ukj sao as incognitas e representam as componentes das matrizes
L e U , respectivamente. Uma vez que cada uma destas matrizes tem n(n+1)2
componentes nao nulas, o numero total de incognitas do sistema e n(n+ 1),superior, portanto, ao numero de equacoes. Isto explica que o sistema sejaindeterminado, isto e, que ele admita muitas solucoes diferentes. A cadauma dessas solucoes corresponde uma certa factorizacao , que se caracterizapor um conjunto de condicoes suplementares. Vamos analisar alguns tiposparticulares de factorizacao, com maior interesse pratico.
3.3.4 Factorizacao de Doolittle
Este tipo de factorizacao caracteriza-se pelo conjunto de condicoes :
lii = 1; i = 1, . . . , n. (3.71)
Vamos mostrar como, a partir destas condicoes, se podem deduzir as formulaspara as componentes das matrizes L e U , que ficam assim completamentedeterminadas.
64
Seja aij uma componente da matriz A, com i ≤ j. Entao , atendendo aforma das matrizes L e U , bem como a condicao 3.71, podemos escrever:
aij =i∑
k=1
lik ukj =i−1∑k=1
lik ukj + uij ; i = 1, . . . , n; j = i, . . . , n. (3.72)
Da formula 3.72 resulta imediatamente que
uij = aij −i−1∑k=1
lik ukj . (3.73)
A fim de deduzir uma formula analoga para a matriz L, consideremos o casode uma componente aij , com i > j. Neste caso, em vez da formula 3.72,temos
aij =
j∑k=1
lik ukj =
j−1∑k=1
lik ukj + lij ujj , i = 1, . . . , n; j = i, . . . , n. (3.74)
Daqui resulta
lij =aij −
∑j−1k=1 lik ukjujj
. (3.75)
Utilizando as formulas 3.73 e 3.75, podem calcular-se todas as componentesdas matrizes L e U . Para isso, basta que todas as componentes da diagonalprincipal de U sejam diferentes de 0. Se, durante processo de calculo, seobtiver alguma dessas componentes igual a 0, entao , tal como acontece nometodo da eliminacao de Gauss, deve-se proceder a alteracoes na matrizU . Neste caso, podemos, por exemplo alterar a ordem das colunas de U ,mantendo a matriz L sem alteracao . Isto corresponde a trocar a ordemdas colunas de A, ou seja, a trocar a ordem das incognitas. Neste caso, aocalcular o determinante de A com base na factorizacao , deve-se entrar emconta com as permutacoes efectuadas. Assim,
detA = (−1)Nt det L detU, (3.76)
onde Nt e o numero de trocas de colunas efectuadas. A troca de colunasde L tambem pode ser aplicada para atenuar os problemas de instabilidadenumerica que podem ocorrer durante a factorizacao . Para isso, usa-se amesma estrategia da pesquisa parcial de pivot que acima descrevemos, parao metodo de Gauss.
E interessante notar que o metodo da eliminacao de Gauss e identico aometodo de Doolittle, podendo, neste sentido, ser considerado tambem um
65
metodo de factorizacao . Para verificar isso, recordemos que no metodo daeliminacao de Gauss se obtem uma matriz triangular superior U , dada pelaformula 3.41. Alem disso, durante o calculo da matriz U sao utilizados oscoeficientes mik, k = 1, . . . , n; i = k + 1, . . . , n, definidos pela formula 3.40.Se dispusermos estes coeficientes numa matriz de ordem n, preenchermosa sua diagonal principal com 1 e as restantes posicoes com 0, obtemos aseguinte matriz triangular inferior:
L =
1 0 . . . 0m21 1 . . . 0. . . . . . . . . . . .. . . mn−1,n−2 1 0. . . mn,n−2 mn,n−1 1
. (3.77)
Teorema 3.3. As matrizes L e U , dadas, respectivamente, pelas formulas3.77 e 3.41,formam uma factorizacao LU da matriz A, identica a factorizacaode Doolittle.
Demonstracao . Vamos demonstrar as igualdades
a(i−1)ij = uij , i = 1, . . . , n; j = i, . . . , n; (3.78)
mij = lij , j = 1, . . . , n; i = j, . . . , n. (3.79)
Para isso, basta comparar as formulas do metodo de Gauss com as da fac-torizacao de Doolittle. Em primeiro lugar, sabemos que
a1j = u1j , j = 1, . . . , n; mi1 =ai1a11
, i = 2, . . . , n. (3.80)
Vamos agora supor, por inducao , que a igualdade 3.78 se verifica para aslinhas da matriz U , com ındice k = 1, . . . , i − 1, e que a igualdade 3.79 severifica para todas colunas de L, com ındice k = 1, . . . , j − 1. Verifiquemosque, nesse caso, as mesmas identidades se mantem validas para a i-esimalinha de U e para a j-esima coluna de L. De facto, de acordo com a formula3.39 do metodo de Gauss, temos
a(k)ij = a
(k−1)ij − mik a
(k−1)kj , (k = 1, . . . , n− 1), (3.81)
66
onde se subentende que a(0)ij = aij , i = 1, . . . , n; j = 1, . . . , n. Aplicando a
formula 3.81, sucessivamente, com k = 1, . . . , i− 1 obtem-se:
a(1)ij = aij − mi1 a1j ;
a(2)ij = a
(1)ij − mi2 a
(1)2j = aij − mi1 a1j − mi2 a
(1)2j ;
. . . . . . . . . . . . . . .
a(i−1)ij = a
(i−2)ij − mi,i−1 a
(i−2)i−1,j = aij −
∑i−1k=1mik a
(k−1)kj .
(3.82)Se, de acordo com a hipotese da inducao , substituirmos os coeficientes mi,k
e a(k−1)k,j , no segundo membro de 3.82, por lik e ukj , obtemos a formula 3.73,
de onde se conclui que a(i−1)ij = uij , com j = i, . . . , n.
Considerando agora a j-esima coluna de L, as suas componentes, deacordo com 3.40, tem a forma
mij =a
(j−1)ij
a(j−1)jj
; i = j, . . . , n. (3.83)
Do mesmo modo como deduzimos a formula 3.82, podemos mostrar que
a(j−1)ij = aij −
j−1∑k=1
mik a(k−1)kj . (3.84)
Se, no segundo membro da formula 3.83, substituirmos o numerador deacordo com 3.84, obtemos
mij =aij −
∑j−1k=1mik a
(k−1)kj
a(j−1)jj
; i = j, . . . , n. (3.85)
Mas, de acordo com a hipotese da inducao , podemos substituir, no segundo
membro de 3.85, a(k−1)kj por ukj ,k = 1, . . . , j e mik por lik,k = 1 . . . , i. Entao
o segundo membro de 3.85 fica igual ao segundo membro de 3.75, de ondese conclui que mij = lij , para todas as componentes da j-esima coluna damatriz L. Fica assim provada, por inducao , a afirmacao do teorema.
Do teorema que acabamos de demonstrar resulta que os metodos deGauss e de Doolittle sao identicos, no sentido em que, dado um certo sis-tema linear, para o resolver segundo cada um desses metodos, efectuam-se
67
exactamente as mesmas operacoes aritmeticas. Em particular, as tres etapasque distinguimos no metodo de Gauss coincidem com as etapas do metodode Doolittle (ou de qualquer outro metodo de factorizacao ):
1. factorizacao da matriz A (reducao da matriz a forma triangular);
2. resolucao do sistema Lg = b (transformacao do segundo membro);
3. resolucao do sistema U x = g (resolucao do sistema triangular supe-rior).
Entao , de acordo com o que dissemos em relacao ao metodo de Gauss,podemos concluir que a etapa mais dispendiosa dos calculos, quando seaplica o metodo de Doolittle, e a primeira, exigindo cerca de 2n3/3 operacoesaritmeticas. As outras duas etapas exigem cerca de n2 operacoes cada uma.Estes numeros tambem se aplicam a factorizacao de Crout, de que falaremosa seguir.
3.3.5 Factorizacao de Crout
Outro tipo, bastante comum, de factorizacao (a factorizacao de Crout)baseia-se nas condicoes
uii = 1, i = 1, . . . , n. (3.86)
As formulas para as componentes das matrizes L e U da factorizacao deCrout deduzem-se da mesma maneira que no caso da factorizacao de Doolit-tle. Assim no caso de i ≥ j, e valida a igualdade
aij =
j∑k=1
lik ukj =
j−1∑k=1
lik ukj + lij , j = 1, . . . , n; i = j, . . . , n. (3.87)
Daqui obtem-se imediatamente
lij = aij −j−1∑k=1
lik ukj . (3.88)
No que diz respeito a matriz L, partimos da igualdade
aij =i∑
k=1
lik ukj =i−1∑k=1
lik ukj + lii uij ; i = 1, . . . , n; j = 1, . . . , i. (3.89)
68
Da igualdade 3.89 resulta
uij =aij −
∑i−1k=1 lik ukjlii
. (3.90)
As formulas 3.88 e 3.90, aplicadas alternadamente (comecando com a primeiracoluna de L e acabando com a ultima linha de U) permitem-nos determinarcompletamente as matrizes L e U da factorizacao de Crout, desde que severifique lii 6= 0, i = 1, . . . , n. Se, durante o processo de factorizacao , severificar que lii = 0, para um certo i, procede-se a uma troca de linhasna matriz L, mantendo U sem alteracao . Esta troca e acompanhada poruma permutacao analoga das linhas da matriz A e das entradas do segundomembro do sistema. Tal como no caso da factorizacao de Doolittle, estaspermutacoes implicam uma troca de sinal no determinante, de acordo coma formula 3.76.
Tambem no caso da factorizacao de Crout e conveniente aplicar a pesquisaparcial de pivot, conduzindo a trocas de linhas quando os elementos diago-nais lii forem pequenos em modulo.
Exemplo 3.5. Consideremos o sistema Ax = b, onde
A =
2 1 3−2 −1 12 4 2
, b =
5−14
. (3.91)
Comecemos por factorizar A segundo o metodo de Doolittle. Como resultada formula 3.73, neste caso a primeira linha de U e igual a primeira linhade A, ou seja,
u11 = 2, u12 = 1, u13 = 3. (3.92)
Calculando os elementos da primeira coluna de L, de acordo com a formula3.75, obtemos
l11 = 1, l21 =a21
u11= −1, l31 =
a31
u11= 1. (3.93)
Passemos ao calculo da segunda linha de U . Temos
u22 = a22 − l21 u12 = 0;u23 = a23 − l21 u13 = 4.
(3.94)
Como sabemos, sendo u22 = 0, nao e possıvel prosseguir os calculos semalterar a matriz A. Assim, e uma vez que u23 6= 0, vamos trocar de lugar
69
a segunda com a terceira coluna de U , fazendo simultaneamente a mesmatroca em A. Sejam U ′ e A′, respectivamente, as matrizes resultantes destaspermutacoes. Entao podemos escrever
u′22 = u23,u′23 = u22.
(3.95)
Continuando o processo de factorizacao com as matrizes U ′ e A′, obtem-se
l32 =a′32 − l31 u
′12
u′22= a33 − l31 u13
u23= −1
4;
u′33 = a′33 − l31 u′13 − l32 u
′23 = a32 − l31 u12 − l32 u22 = 3.
(3.96)Recapitulando, obtivemos a seguinte factorizacao de A:
L =
1 0 0−1 1 01 −1
4 1
, U ′ =
2 3 10 4 00 0 3
. (3.97)
Para obter a solucao do sistema dado, comecemos por resolver o sistemacom a matriz triangular inferior Lg = b, de acordo com o metodo habitual:
g1 = b1 = 5;−g1 + g2 = b2 ⇔ g2 = 4;
g1 − g2/4 + g3 = b3 ⇔ g3 = 0.(3.98)
Ao resolver depois o sistema U ′ x = g, temos de ter em conta que a segundacoluna de U trocou de lugar com a terceira. Isto equivale a uma troca delugares entre x2 e x3. Assim, temos
3x2 = g3 ⇔ x2 = 0;4x3 = g2 ⇔ x3 = 1.
2x1 + 3x3 + x2 = g1 ⇔ x1 = 1.(3.99)
Se, em vez do metodo de Doolittle, quisermos aplicar a factorizacao deCrout, teremos de basear os calculos nas formulas 3.88 e 3.90. Nesse caso,a primeira coluna de L fica igual a primeira coluna de A. Para a primeiralinha de U , obtem-se
u11 = 1; u12 =a12
l11=
1
2; u13 =
a13
l11=
3
2. (3.100)
Na segunda coluna de L, obtemos:
l22 = a22 − l21 u12 = 0;l32 = a32 − l31 u12 = 3.
(3.101)
70
Uma vez que l22 = 0, torna-se necessario trocar a segunda com a terceiralinha de L (e, consequentemente, de A). Feita esta permutacao , obtemos
l′22 = l32 = 3;l′32 = l22 = 0.
(3.102)
Resta calcular as componentes da segunda linha de U e terceira coluna de L :
u23 =a′23 − l′21u13
l′22= −1
3;
l′33 = a′33 − l′31u13 − l′32u23 = 4.
(3.103)
Consequentemente, a factorizacao de Crout da matriz dada tem a forma:
L′ =
2 0 02 3 0−2 0 4
, U =
1 12
32
0 1 −13
0 0 1
. (3.104)
A partir de qualquer uma das factorizacoes de A obtidas, utilizando aformula 3.76, calcula-se facilmente o determinante de A
detA = det L′ (−1)1 = det U ′ (−1)1 = −24.
Se quisermos resolver o sistema dado, com base na factorizacao de Crout,basta considerar o segundo membro b′ = (5, 4,−1)T (uma vez que foi trocadaa segunda com a terceira linha de U), apos o que se resolvem os sistemasL′ g = b′ e U x = g, utilizando, como sempre, a substuicao descendente(para o primeiro sistema) e a substituicao ascendente (para o segundo).
3.3.6 Factorizacao de Cholesky
Os dois tipos de factorizacao de que falamos acima existem para qualquermatriz nao singular (ainda que seja necessario efectuar uma troca de linhasou colunas). Quanto a factorizacao de Cholesky, de que vamos falar a seguir,so existe para matrizes simetricas e definidas positivas. Embora se trate deuma restricao muito forte, este tipo de factorizacao nao deixa de ter inte-resse pratico, visto que estas propriedades sao satisfeitas pelas matrizes quesurgem em muitos problemas de calculo numerico, por exemplo, no metododos mınimos quadrados e em certos problemas de valores de fronteira paraequacoes diferenciais. A grande vantagem desta factorizacao consiste em
71
que so temos de calcular uma matriz, L, visto que uma matriz simetrica edefinida positiva pode ser representada sob a forma:
A = LLT . (3.105)
Isto significa que o numero de operacoes para resolver um sistema linearfica reduzido a cerca de metade, quando se compara o metodo de Choleskycom outros metodos de factorizacao ou com o metodo de Gauss. Para es-tudar melhor as propriedades da factorizacao de Cholesky, comecemos pordemonstrar o seguinte teorema.
Teorema 3.4. Seja A uma matriz simetrica e definida positiva. Entaoexiste uma matriz real, triangular inferior, L, tal que se verifica a factor-izacao 3.105.
Demonstracao . Provemos esta afirmacao por inducao em n (ordemda matriz). Para o caso de n = 1, a afirmacao e trivial, visto que, nessecaso, A e um numero positivo e L, a sua raiz quadrada. Suponhamos agoraque a afirmacao do teorema e verdadeira para qualquer matriz de ordemnao superior a n − 1. Seja A o menor principal de A, de ordem n − 1.Verifiquemos que A tambem e uma matriz definida positiva. Na realidade,se A e definida positiva, entao verifica-se
n∑i=1,j=1
aij xi xj > 0, ∀x ∈ IRn, x 6= 0. (3.106)
Em particular, se considerarmos x = (x1, x2, . . . , xn−1, 0), onde xi, i =1, . . . , n− 1 sao numeros reais arbitrarios, de 3.106 resulta
n−1∑i=1,j=1
aij xi xj > 0; (3.107)
o que significa que A tambem e definida positiva.Vamos procurar uma matriz triangular inferior L, de ordem n, com a
forma
L =
[L 0γT z
], (3.108)
onde L e uma matriz de ordem n− 1, γ ∈ IRn−1, z ∈ IR. Queremos provarque, para essa matriz L, se verifica
LLT =
[L 0γT z
] [LT γ0 z
]= A =
[A ccT d
], (3.109)
72
onde c ∈ IRn−1 e d = ann ∈ IR. Por hipotese da inducao , existe umamatriz real, triangular inferior L, tal que
A = L LT . (3.110)
Resta provar que existe um vector γ ∈ IRn−1 e um numero real z, para osquais se verifica a igualdade 3.109. Quanto a existencia de γ, ela resulta deo sistema L γ = c ter solucao .De facto, a matriz L e nao singular, uma vezque
(det L)2 = det A > 0
(visto que A e definida positiva).Mostremos agora que existe um numeroreal z, para o qual a igualdade 3.109 e verdadeira. Para que a igualdade3.109 seja verdadeira, deve verificar-se
detA = ( det L)2 z2 > 0. (3.111)
Uma vez que detA > 0,a equacao 3.111 tem duas raizes e o valor de zpode ser determinado, escolhendo a raiz positiva. Assim, o teorema ficademonstrado, por inducao .
Observacao . Note-se que o problema da factorizacao de Choleskyadmite mais do que uma solucao . Isso verifica-se pela equacao 3.111, a qualadmite duas raizes reais:
z1,2 = ±√detA
det L(3.112)
Por convencao , escolhe-se sempre a raiz positiva desta equacao , o quesignifica que todos os elementos da diagonal principal de L sao positivos.
Vejamos agora, em termos praticos, como se pode calcular a matriz Lda factorizacao de Cholesky. Seja aij uma componente de A, com i ≥ j.Entao, da igualdade 3.105 resulta
aij =
j∑k=1
lik ljk =
j−1∑k=1
lik ljk + lij ljj , j = 1, . . . , n; i = j, . . . , n. (3.113)
No caso de i = j, da igualdade 3.113 obtem-se a formula para os elementosda diagonal principal de L:
lii =
√√√√aii −i−1∑k=1
l2ik; i = 1, . . . , n. (3.114)
73
De acordo com o teorema 3.2, todos os elementos da diagonal principal deL sao reais, pelo que o segundo membro de 3.114 e sempre real. Uma vezcalculado ljj , podemos obter as restantes componentes da j-esima coluna deL. Da formula 3.113 resulta:
lij =aij −
∑j−1k=1 lik ljkljj
; i = j + 1, . . . , n. (3.115)
Assim, usando as formulas 3.114 e 3.115, alternadamente, pode ser obtida afactorizacao de Cholesky da matriz A.
Exemplo 3.6. Consideremos a matriz de ordem n com a forma geral
A =
4 2 0 . . . 02 5 2 . . . 00 2 5 . . . 0. . . . . . . . . . . . . . .0 . . . 2 5 20 . . . 0 2 5
. (3.116)
Trata-se de uma matriz simetrica tridiagonal, isto e
aij 6= 0 ⇒ |i − j| ≤ 1.
Matrizes com estas caracterısticas aparecem frequentemente nas aplicacoes.Vamos tentar obter a sua factorizacao de Cholesky. Como nao e simplesdeterminar, a partida, se a matriz dada e definida positiva, vamos tentarutilizar as formulas 3.114 e 3.115 e verificar se elas sao aplicaveis. Comece-mos, como sempre, pela componente l11. De acordo com 3.114, o seu valore o seguinte:
l11 =√a11 = 2. (3.117)
As restantes componentes desta coluna sao dadas pela formula 3.115:
l21 = a21l11
= 1;
lk1 = ak1l11
= 0; k = 3, . . . , n.
(3.118)
Vamos provar agora, por inducao , que as restantes colunas da matrix L tema mesma estrutura, isto e para a coluna j, verifica-se:
ljj = 2;lj+1,j = 1;li,j = 0; i = j + 2, . . . , n.
(3.119)
74
Para a primeira coluna, as formulas 3.119 ja estao provadas. Suponhamosagora que estas formulas sao validas para todas as colunas, ate a j − 1.Vejamos o que acontece com a coluna j. De acordo com a formula 3.114,podemos escrever
ljj =
√√√√ajj −j−1∑k=1
l2jk =√ajj − l2j,j−1 = 2. (3.120)
Depois, aplicando a formula 3.115, obtemos
lj+1,j =aj+1,j
ljj= 1;
li,j = 0; i = j + 2, . . . , n.
(3.121)
Fica assim provado que a factorizacao de Cholesky da matriz dada e definidapor uma matriz triangular inferior com a forma
L =
2 0 0 . . . 01 2 0 . . . 00 1 2 . . . 0. . . . . . . . . . . . . . .0 . . . 1 2 00 . . . 0 1 2
. (3.122)
O determinante de A pode ser calculado com base nesta factorizacao :
detA = (det L)2 = (l11 l22 . . . lnn)2 = (2n)2 = 4n. (3.123)
Uma vez que a formula 3.123 e valida para qualquer n, ela pode servir paracalcularmos os menores principais da matriz A dada. Assim, temos
A1 = 4, A2 = 42 . . . , An = detA = 4n.
Fica assim provado que todos os menores principais de A sao positivos, deonde resulta que A e definida positiva.
3.4 Netodos Iterativos para Sistemas Lineares
Neste paragrafo, vamos estudar alguns metodos iterativos, utilizados nocalculo aproximado de solucoes de sistemas lineares. Antes disso, porem,vamos apresentar alguns conceitos gerais, relativos aos metodos iterativos,e que nos voltarao a ser uteis nos capıtulos posteriores.
75
3.4.1 Nocoes Basicas sobre Metodos Iterativos
Em muitos problemas matematicos, quando se revela impossıvel ou muitodifıcil calcular a solucao exacta de uma certa equacao , opta-se por obterum valor aproximado dessa solucao . Esse valor aproximado e geralmenteobtido por um metodo de aproximacoes sucessivas, ou metodo iterativo, emque cada nova aproximacao e obtida a partir da anterior (ou das anteriores),atraves de um algoritmo conhecido, pretendendo-se deste modo tornar o errode aproximacao cada vez mais pequeno.
Definicao 3.5. Seja E um espaco normado e X um subconjunto de E.Chama-se metodo iterativo a p passos em E (definido sobre X) a qualqueraplicacao Ψ, que a cada sucessao de p elementos (ξ0, . . . , ξp−1) ∈ Xp, fazcorresponder uma sucessao infinita {x(k)}∞k=0,x(k) ∈ E , com as seguintes propriedades:
1. Os primeiros p termos da sucessao {x(k)}∞k=0 sao os dados:
x(i) = ξi, i = 0, . . . , p− 1.
2. Os restantes elementos elementos de {x(k)}∞k=0 sao obtidos a partirdestes, de acordo com a formula
x(k+p) = φ(xk, xk+1, . . . , xk+p−1),
onde φ e uma funcao conhecida (chamada a funcao iteradora) comdomınio em Xp e valores em E.
Naturalmente, na pratica so se calcula um numero finito de termos dasucessao {x(k)}∞k=0 (tambem chamados iteradas), tantos quantos necessariospara alcancar a precisao pretendida. Por isso, a cada metodo iterativo estaogeralmente associados criterios de paragem, isto e, regras que nos permitemverificar se uma dada iterada tem ou nao a precisao exigida.
Definicao 3.6. Dizemos que um metodo iterativo a p passos, definidosobre X, e convergente para um certo x ∈ IRn, se, para quaisquer valoresiniciais (ξ0, . . . , ξp−1) ∈ A se verificar x(k) → x, quando k → ∞ (segundo anorma considerada em IRn).
76
Note-se que a convergencia em espacos de dimensao finita nao depende danorma considerada. Daı que, no caso dos metodos iterativos para sistemaslineares, que vamos estudar nos proximos paragrafos, a convergencia numacerta norma e equivalente a convergencia noutra norma qualquer.
Alem da convergencia, outra propriedade importante dos metodos iter-ativos e a sua estabilidade.
Definicao 3.7. Um metodo iterativo Ψ a p passos, definido no conjuntoX, diz-se estavel em A ⊂ X, se existir uma constante C > 0, tal que
maxn∈N‖x(n) − y(n)‖ ≤ C max
i=1,...,p‖ξi − ηi‖ ∀ξ, η ∈ Ap, (3.124)
onde {x(n)} e {y(n)} sao , respectivamente, as sucessoes geradas a partir deξ = (ξ0, ξ1, . . . , ξp−1) e η = (η0, η1, . . . , ηp−1).
Para representar o erro da k-esima iterada usaremos a notacao e(k),ou seja, e(k) = x(k) − x. Nos proximos paragrafos, vamos analisar algunsmetodos iterativos para o calculo aproximado da solucao do sistema linear
Ax = b, (3.125)
onde A ∈ IRn×n e b ∈ IRn. Suponhamos, como habitualmente, que amatriz A e nao singular, pelo que o sistema (3.125) tem uma unica solucao.
Com o objectivo de construir um metodo iterativo, comecamos por re-duzir o sistema (3.125) a forma equivalente
x = G(x) = Cx+ g, (3.126)
onde C e uma certa matriz (a que chamaremos matriz de iteracao), e g e umvector auxiliar (g ∈ IRn). Uma vez escrito o sistema na forma (3.126), pode-mos dizer que a sua solucao e um ponto fixo da funcao G (funcao definidaem IRn e com valores no mesmo espaco). A ideia e determinar o pontofixo de G por um metodo analogo ao metodo do ponto fixo, utilizado nocapıtulo anterior para aproximar os pontos fixos de funcoes de uma variavel.Assim, dada uma certa aproximacao inicial x(0) ∈ IRn,vamos construir umasucessao de vectores atraves da seguinte formula de recorrencia:
x(k+1) = G(x(k)) = Cx(k) + g, k = 0, 1, . . . (3.127)
Naturalmente, esta transformacao do sistema pode ser feita de muitas maneirasdiferentes, dando origem a diferentes metodos iterativos. Vamos descreverdois deles.
77
3.4.2 Metodo de Jacobi
Para deduzir a formula do metodo de Jacobi, comecamos por reescrever osistema (3.125) do seguinte modo:
x1 = b1−a12x2−a13x3−···−a1nxna11
x2 = b2−a21x1−a23x3−···−a2nxna22
· · ·xn =
bn−an1x1−an2x2−···−an,n−1xn−1
ann
. (3.128)
Ficamos asim com o sistema escrito na forma (3.126), de tal modo queo sistema (3.128) e equivalente ao inicial. Para poder escrever o sistemanesta forma, basta que todos os elementos da diagonal principal da matrizA sejam diferentes de 0 (aii 6= 0,i = 1, . . . , n). Se iterarmos agora a funcaoG correspondente ao sistema (3.128), obtem-se a seguinte formula iteradora:
x(k+1)1 =
b1−a12x(k)2 −a13x(k)3 −···−a1nx
(k)n
a11
x(k+1)2 =
b2−a21x(k)1 −a23x(k)3 −···−a2nx
(k)n
a22
· · ·
x(k+1)n =
bn−an1x(k)1 −an2x(k)2 −···−an,n−1x
(k)n−1
ann
, k = 0, 1, 2, . . . (3.129)
A formula (3.129) pode escrever-se na seguinte forma compacta:
x(k+1)i =
biaii−∑n
j=1,...,n,j 6=i aijx(k)j
aii, (3.130)
i = 1, . . . , n, k = 0, 1, 2, . . .Exemplo 3.7. Consideremos o sistema Ax = b, onde
A =
2 1 0−1 2 10 −1 2
, b =
221
. (3.131)
1. Efectuar uma iteracao do metodo de Jacobi, tomando como aproxi-macao inicial x(0) = (0.5, 0.8, 1).
78
Resolucao:
x(1)1 =
b1−a12x(0)2 −a13x(0)3
a11= 1
2 (2− 0.8− 0) = 0.6;
x(1)2 =
b2−a21x(0)1 −a23x(0)3
a22= 1
2(2 + 0.5− 1) = 0.75;
x(1)3 =
b3−a31x(0)1 −a32x(0)2
a33= 1
2(1− 0 + 0.8) = 0.9;
. (3.132)
2. Sabendo que a solucao exacta do sistema e x = (0.583, 0.833, 0.917),calcule ‖e(0)‖1 e ‖e(1)‖1.
Resolucao:
e(0) = x− x(0) = (0.083, 0.033,−0.083) ‖e(0)‖1 = 0.199;
e(1) = x− x(1) = (−0.017, 0.083, 0.017) ‖e(1)‖1 = 0.117.(3.133)
A norma do erro da primeira iterada e menor, o que significa que elaesta mais proxima da solucao exacta do que a aproximacao inicial. Noentanto, daqui nada podemos concluir sobre a convergencia do metodo.Esta sera analisada mais tarde, segundo os metodos que vamos estudarnos proximos paragrafos.
3.4.3 Metodo de Gauss-Seidel
O metodo de Gauss-Seidel e um dos metodos iterativos mais comuns pararesolucao de sistemas lineares. Para deduzirmos a sua forma iteradora, par-timos de novo do sistema na forma (3.128).
Neste caso, porem, a formula iteradora tem o seguinte aspecto:
x(k+1)1 =
b1−a12x(k)2 −a13x(k)3 −···−a1nx
(k)n
a11
x(k+1)2 =
b2−a21x(k+1)1 −a23x(k)3 −···−a2nx
(k)n
a22
· · ·
x(k+1)n =
bn−an1x(k+1)1 −an2x(k+1)
2 −···−an,n−1x(k+1)n−1
ann
, k = 0, 1, 2, . . . (3.134)
A diferenca em relacao ao metodo de Jacobi esta em que, para determinar a
componente x(k+1)i da iterada (k+1) (com i > 1) utilizamos as componentes
x(k+1)1 , ..., x
(k+1)i−1 dessa mesma iterada, enquanto que no metodo de Jacobi
79
as componentes de x(k+1) eram calculadas apenas a partir das componentesde x(k) (iterada anterior). A formula (3.134) pode ser escrita na seguinteforma compacta:
x(k+1)i =
biaii−∑i−1
j=1 ai,jx(k+1)j +
∑nj=i+1 ai,jx
(k)j
aii, (3.135)
i = 1, . . . , n, k = 0, 1, 2, . . . . Vejamos um exemplo de aplicacao.Exemplo 3.8. Consideremos de novo o sistema (3.131).
1. Efectuar uma iteracao do metodo de Gauss-Seidel, tomando comoaproximacao inicial x(0) = (0.5, 0.8, 1).
Resolucao:
x(1)1 =
b1−a12x(0)2 −a13x(0)3
a11= 1
2(2− 0.8− 0) = 0.6;
x(1)2 =
b2−a21x(1)1 −a23x(0)3
a22= 1
2(2 + 0.6− 1) = 0.8;
x(1)3 =
b3−a31x(1)1 −a32x(1)2
a33= 1
2(1− 0 + 0.8) = 0.9;
. (3.136)
2. Sabendo que a solucao exacta do sistema e x = (0.583, 0.833, 0.917),calcule ‖e(0)‖1 e ‖e(1)‖1.
Resolucao:
e(0) = x− x(0) = (0.083, 0.033,−0.083) ‖e(0)‖1 = 0.199;
e(1) = x− x(1) = (−0.017, 0.033, 0.017) ‖e(1)‖1 = 0.067.(3.137)
Tal como acontecia no caso do metodo de Jacobi, tambem aqui anorma do erro diminui da aproximacao inicial para a primeira iterada,o que significa que esta esta mais proxima da solucao exacta do que aaproximacao inicial. Na caso do metodo de Gauss-Seidel, a diferencaentre os dois erros e mais acentuada. Embora, como ja foi dito, nadase possa concluir destes factos, eles estao de acordo com a analise quevamos fazer mais adiante sobre a convergencia dos metodos iterativos.
3.4.4 Forma Matricial dos Metodos Iterativos
Para podermos escrever as formulas iteradoras dos metodos iterativos de ummodo mais compacto e estudar a sua convergencia, convem representa-los
80
na forma matricial. Para isso, dada uma certa matriz A, comecamos pordefinir as matrizes L,D,U , tais que
L =
0 0 . . . 0a21 0 . . . 0. . . . . . . . . . . .an1 an2 . . . 0
, D =
a11 0 . . . 00 a22 . . . 0. . . . . . . . . . . .0 0 . . . ann
, U =
0 a12 . . . a1n
0 0 . . . a2n
. . . . . . . . . . . .0 0 . . . 0
.(3.138)
Obviamente, verifica-se A = L + D + U . Nesta seccao, vamos supor quetodas as entradas diagonais da matriz A sao diferentes de zero, ou seja,aii 6= 0, i = 1, . . . , n. Daqui resulta que a matriz D e invertıvel.
Metodo de Jacobi na Forma MatricialUtilizando estas notacoes, vejamos como se pode escrever a formula it-
eradora do metodo de Jacobi na forma (3.127), identificando o vector g e amatriz de iteracao correspondentes.
Comecamos por escrever a formula (3.130) com o auxılio das matrizesL,D e U :
x(k+1) = D−1(b− Lx(k) − Ux(k)
). (3.139)
Esta equacao tambem pode ser escrita na seguinte forma:
x(k+1) = D−1b−D−1(L+ U)x(k). (3.140)
Comparando esta equacao com a formula geral para os metodos iterativos(3.127), concluımos que, no caso do metodo de Jacobi, o vector auxiliar g ea matriz de iteracao tem a forma
gj = D−1b, Cj = −D−1(L+ U). (3.141)
Uma vez que a matriz D e diagonal e que todas as entradas da sua diagonalsao diferentes de zero, a sua inversa pode ser determinada imediatamente:
D−1 =
1a11
0 . . . 0
0 1a22
. . . 0
. . . . . . . . . . . .0 0 . . . 1
ann
. (3.142)
Por conseguinte, a matriz de iteracao Cj tem a seguinte forma:
Cj = −D−1(L+ U) =
0 −a12
a11. . . −a1n
a11
−a21a22
0 . . . −a2na22
. . . . . . . . . . . .− an1ann
− an2ann
. . . 0
. (3.143)
81
Metodo de Gauss-Seidel na Forma MatricialVejamos agora como se pode escrever a formula iteradora do metodo de
Gauss-Seidel na forma (3.127).Com o auxılio das matrizes L,D e U , a formula (3.135) tem o seguinte
aspecto:
x(k+1) = D−1(b− Lx(k+1) − Ux(k)
). (3.144)
Multiplicando ambos os membros de (3.144) a esquerda por D, obtem-se:
Dx(k+1) = b− Lx(k+1) − Ux(k). (3.145)
Juntando no primeiro membro os termos que contem x(k+1), temos a seguinteequacao:
(L+D)x(k+1) = b− Ux(k). (3.146)
Uma vez que a matriz D, por condicao, e invertıvel, L + D tambem o e (oseu determinante e igual ao determinante de D). Assim, podemos escrever
x(k+1) = (L+D)−1b− (L+D)−1Ux(k). (3.147)
Finalmente, comparando a equacao (3.147) com a formula geral para osmetodos iterativos, concluımos que, no caso do metodo de Gauss-Seidel, ovector auxiliar g e a matriz de iteracao tem a forma
ggs = (L+D)−1b, Cgs = −(L+D)−1U. (3.148)
Neste caso, nao e possıvel encontrar uma forma explıcita para a inversa de(L + D). Tudo o que se pode dizer e que e uma matriz triangular inferiore que os seus elementos diagonais sao os inversos dos elementos diagonaisde A. Logo, tambem nao e possıvel encontrar uma forma explıcita para amatriz de iteracao Cgs.
Exemplo 3.9. Determinemos os vectores g e as matrizes de iteracao dosmetodos de Jacobi e do metodo de Gauss-Seidel para o sistema do exemplo3.7. No caso do metodo de Jacobi, o vector g tem a forma
ggs = D−1b = (b1a11
,b2a22
,b3a33
) = (1, 1,1
2). (3.149)
A matriz CJ obtem-se da aplicacao directa da formula (3.143):
CJ = −D−1(L+ U) =
0 −12 0
12 0 −1
20 1
2 0
. (3.150)
82
No caso do metodo de Gauss-Seidel, para poder determinar o vector g e amatriz de iteracao comecamos por calcular a inversa de L+D:
(L+D)−1 =
12 0 014
12 0
18
14
12
. (3.151)
Daqui, pelas formulas (3.148), resulta que
ggs = (1,3
2,5
4), Cgs =
0 −12 0
0 −14 −1
20 −1
8 −14
. (3.152)
3.4.5 Convergencia dos Metodos Iterativos para Sistemas Lin-eares
Uma vez definido um metodo iterativo para a aproximacao da solucao deum sistema, e fundamental saber em que condicoes esse metodo gera umasucessao que converge para a solucao exacta. Nos teoremas que se seguem,estabelecem-se condicoes sobre a matriz A que garantem a convergencia dosmetodos iterativos considerados.
Conforme resulta das formulas (3.126) e (3.127), os erros das iteradas,que representaremos por e(k), k = 0, 1, . . ., satisfazem a igualdade
e(k+1) = x(k+1) − x = C(x(k) − x); k = 0, 1, 2, . . . (3.153)
A igualdade (3.153) tambem pode ser escrita na forma
e(k+1) = C e(k); k = 0, 1, 2, . . . (3.154)
onde C e a matriz de iteracao do metodo considerado. No paragrafo anteriorja foi analisada a forma das matrizes de iteracao dos metodos de Jacobi eGauss-Seidel.
Vejamos agora que propriedades deve apresentar a matriz C para garan-tir a convergencia de um metodo iterativo. Em primeiro lugar, notemos queda igualdade (3.154) resulta imediatamente uma formula que exprime o errode qualquer iterada atraves do erro da aproximacao inicial:
e(k) = Ck e(0); k = 0, 1, 2, . . . (3.155)
Definicao 3.7 Uma matriz C ∈ IRn×n, que representa uma aplicacaolinear no espaco vectorial IRn, diz-se convergente se e so se
limk→∞
Ck x = 0, ∀x ∈ IRn. (3.156)
83
Estamos agora em condicoes de formular o teorema que nos da umacondicao necessaria e suficiente para a convergencia dos metodos iterativosbaseados na formula (3.127).
Teorema 3.5. Seja {x(k)}∞k=0 uma sucessao em IRn, gerada pela formula (3.127),com uma certa matriz C. Entao esta sucessao converge para a solucao dosistema Ax = b , qualquer que seja a aproximacao inicial x0, se e so se amatriz C for convergente.
Demonstracao .Condicao suficiente. Seja C uma matriz convergente eseja e(k) o erro da k-esima iterada. Entao , de acordo com as formulas (3.155)e (3.156), temos
limk→∞
e(k) = limk→∞
Ck e(0) = 0, (3.157)
qualquer que seja o vector e(0) ∈ IRn, independentemente da norma consid-erada. Isto significa que o metodo iterativo converge, qualquer que seja aaproximacao inicial x(0) ∈ IRn .
Condicao necessaria. Suponhamos que C nao e convergente. Entaoexiste um vector v ∈ IRn tal que a sucessao {Ck v}∞k=0 nao converge para ovector nulo. Seja x(0) = x+v, onde x e a solucao exacta do sistema. Entao, de acordo com (3.155), temos e(k) = Ck v e, de acordo com a definicaode v, a sucessao {e(k)}∞k=0 nao tende para o vector nulo. Isto, por sua vez,significa que o metodo iterativo nao e convergente, no caso da aproximacaoinicial x(0) = x+ v .
Na pratica, pode nao ser facil averiguar se a matriz C e ou nao conver-gente. Vamos a seguir apresentar dois teoremas que nos ajudam a verificaresta propriedade.
Teorema 3.6. Seja C ∈ IRn×n. Se existir uma norma matricial M ,associada a uma certa norma vectorial V , tal que
‖C‖M < 1 (3.158)
entao a matriz C e convergente.Demonstracao . Seja x um vector arbitrario de IRn. Entao , de acordo
com as propriedades das normas matriciais, temos
‖Ck x‖V ≤ ‖Ck‖M ‖x‖V ≤ (‖C‖M )k ‖x‖V . (3.159)
84
Da desigualdade (3.159) resulta imediatamente que, se ‖C‖M < 1, entao
limk→∞
‖Ck x‖V = 0, (3.160)
o que significa, por definicao , que a matriz C e convergente.
Teorema 3.7. Para que a matriz C ∈ IRn×n seja convergente e necessarioe suficiente que
rσ(C) < 1. (3.161)
Demonstracao .Condicao suficiente. Se tivermos rσ(C) = ρ < 1, deacordo com o teorema 2.1.31 de [4], para qualquer ε > 0, existe uma normamatricial N(ε) tal que
‖C‖N(ε) ≤ ρ + ε. (3.162)
Se considerarmos ε = 1−ρ2 , obtemos
‖C‖N(ε) ≤ρ + 1
2< 1. (3.163)
Da desigualdade (3.163) resulta, pelo teorema 3.5, que a matriz C e conver-gente.
Condicao necessaria. Suponhamos que a condicao (3.161) nao se veri-fica, isto e, que rσ(C) ≥ 1. Entao existe, pelo menos, um valor proprio λde C, tal que |λ| = ρ ≥ 1. Seja v um vector proprio de C, associado ao valorproprio λ. Entao , para qualquer norma vectorial, verifica-se
‖Ck v‖ = ‖λk v‖ = |λ|k ‖v‖. (3.164)
Da igualdade (3.164), visto que |λ| = ρ ≥ 1, resulta que a sucessao {Ck v}∞k=0
nao converge para o vector nulo, pelo que a matriz C nao e convergente.
Se aplicarmos os teoremas 3.5 e 3.6 aos metodos de Jacobi e Gauss-Seidel,podemos obter outros criterios de convergencia para estes metodos. A fimde formularmos estes criterios numa forma mais concisa, vamos introduzirmais algumas definicoes .
Definicao 3.8. Diz-se que a matriz A ∈ IRn×n tem a diagonal estrita-mente dominante por linhas, se forem satisfeitas as condicoes :
n∑j=1,j 6=i
|aij | < |aii|, i = 1, . . . , n. (3.165)
85
Definicao 3.9 . Diz-se que a matriz A ∈ IRn×n tem a diagonal estrita-mente dominante por colunas, se forem satisfeitas as condicoes :
n∑i=1,i 6=j
|aij | < |ajj |, j = 1, . . . , n. (3.166)
Do teorema 3.6 resulta o seguinte criterio de convergencia para o metodode Jacobi.
Teorema 3.8. Se a matriz A tiver a diagonal estritamente dominantepor linhas, entao o metodo de Jacobi converge para a solucao do sistemaAx = b, qualquer que seja a aproximacao inicial x(0) ∈ IRn.
Demonstracao . Se a matriz A tiver a diagonal estritamente dominantepor linhas, das desigualdades (3.165) resulta
n∑j=1,j 6=i
|aij ||aii|
< 1 i = 1, . . . , n. (3.167)
De acordo com a forma da matriz CJ , dada por (3.141), as desigualdades(3.167) implicam
‖CJ‖∞ = maxi=1,...,n
n∑j=1,j 6=i
|aij ||aii|
< 1, (3.168)
De acordo com o teorema 3.6, a condicao (3.168) garante que a matriz CJ econvergente. Isto, por sua vez, implica, de acordo com o teorema 3.5, queo metodo de Jacobi e convergente, qualquer que seja a aproximacao inicial.
Um criterio analogo e valido no caso de a matriz A ter a diagonal estri-tamente dominante por colunas.
Teorema 3.9. Se a matriz A tiver a diagonal estritamente domi nantepor colunas, entao o metodo de Jacobi converge para a solucao do sistemaAx = b, qualquer que seja a aproximacao inicial x(0) ∈ IRn.
Demonstracao . Suponhamos que a matriz A satisfaz as condicoes(3.166). Seja D uma matriz diagonal com a forma D = diag(a11, . . . , ann).
86
Entao podemos definir uma norma matricial M pela igualdade
‖C‖M = ‖DCD−1‖1, ∀C ∈ Ln. (3.169)
Das condicoes (3.166) obtem-se facilmente que
‖CJ‖M = ‖DCJ D−1‖1 < 1. (3.170)
De acordo com o teoremas 3.5 e 3.6, da desigualdade (3.170) resulta que ometodo de Jacobi converge para a solucao do sistema Ax = b, qualquerque seja a aproximacao inicial x(0) ∈ IRn. .
Exemplo 3.10 Voltemos ao sistema do exemplo 3.7 . Recordemos quea matriz A daquele sistema tem a forma
A =
2 1 0−1 2 10 −1 2
. (3.171)
Verifica-se facilmente que esta matriz nao tem a diagonal estritamente dom-inante por linhas, uma vez que, neste caso,
|a21| + |a23| = 2 = |a22|.
Do mesmo modo se pode verificar que aquela matriz tambem nao tem adiagonal estritamente dominante por colunas. Por conseguinte, os teoremas3.8 e 3.9 nao sao , neste caso, aplicaveis. Vejamos se e possıvel aplicardirectamente o teorema 3.7.
A matriz CJ , neste caso, tem a forma:
CJ =
0 −1/2 01/2 0 −1/20 1/2 0
. (3.172)
Os valores proprios desta matriz sao as raizes da equacao
λ3 +λ
2= 0.
Determinando estas raizes, obtem-se
λ1 = 0, λ2 =i√2, λ3 = − i√
2.
87
Por conseguinte, o raio espectral de CJ e
rσ(CJ) = |λ2| =1√2< 1.
Logo, pelo teorema 3.7, podemos concluir que o metodo de Jacobi convergepara a solucao do sistema considerado, qualquer que seja a aproximacaoinicial.
Vamos agora averiguar as condicoes que garantem a convergencia dometodo de Gauss-Seidel. Representemos por Cgs a matriz
Cgs = −(L+D)−1U (3.173)
De acordo com o teorema 3.5, o metodo de Gauss-Seidel converge, qualquerque seja a aproximacao inicial, se e so se a matriz CS for convergente.Para isso, de acordo com o teorema 3.7, e necessario e suficiente que o raioespectral daquela matriz seja menor que 1.
Por outro lado, pode mostrar-se que o metodo de Gauss-Seidel, tal comoo de Jacobi, converge sempre que a matriz do sistema tem a diagonal es-tritamente dominante por linhas. Para isso, comecemos por introduzir asseguintes notacoes :
αi =
i−1∑j=1
∣∣∣∣aijaii∣∣∣∣ , βi =
n∑j=i+1
∣∣∣∣aijaii∣∣∣∣ . (3.174)
Sendo conhecidos αi e βi,i = 1, . . . , n, define-se a grandeza η atraves daformula
η = maxi=1,...,n
(βi
1− αi
). (3.175)
Note-se que, se a matriz A tiver a diagonal estritamente dominante, porlinhas, entao αi + βi < 1, i = 1, . . . , n, de onde resulta que η < 1.
Teorema 3.10. Seja A uma matriz com a diagonal estritamente dom-inante por linhas. Entao o metodo de Gauss-Seidel converge, qualquer queseja a aproximacao inicial e e valida a estimativa do erro
‖e(k)‖∞ ≤ ηk ‖e(0)‖∞. (3.176)
88
Demonstracao . Da formula (3.135) deduz-se facilmente que o erro dak-esima iterada do metodo de Gauss-Seidel satisfaz a igualdade
e(k+1)i =
1
aii
− i−1∑j=1
aij e(k+1)j −
n∑j=i+1
aij e(k)j
, i = 1, 2, . . . , n; k = 0, 1, . . . .
(3.177)Tomando o modulo de ambos os membros da igualdade (3.177) e entrandoem conta com as definicoes das grandezas αi e βi, obtem-se
|e(k+1)i | ≤ αi ‖e(k+1)‖∞ + βi‖e(k)‖∞, i = 1, 2, . . . , n; k = 0, 1, . . . . (3.178)
Seja m o ındice para o qual se verifica |e(k+1)m | = ‖e(k+1)‖∞. Entao , es-
crevendo a desigualdade (3.178), com i = m, obtem-se
‖e(k+1)‖∞ ≤ αm ‖e(k+1)‖∞ + βm‖e(k)‖∞, k = 0, 1, . . . . (3.179)
Daqui resulta
‖e(k+1)‖∞(1 − αm) ≤ βm‖e(k)‖∞, k = 0, 1, . . . . (3.180)
Visto que αm < 1, podemos dividir ambos os membros de (3.180) por1 − αm e obter
‖e(k+1)‖∞ ≤βm
1 − αm‖e(k)‖∞ ≤ η‖e(k)‖∞, k = 0, 1, . . . . (3.181)
Da desigualdade (3.181) resulta a estimativa do erro (3.176). Por outro lado,uma vez que a matriz tem a diagonal estritamente dominante, por linhas,η < 1. Logo, a desigualdade (3.176) implica que
limk→∞
‖e(k)‖∞ = 0,
o que garante a convergencia do metodo de Gauss-Seidel, qualquer que sejaa aproximacao inicial .
Exemplo 3.11 Consideremos o mesmo sistema linear dos exemplos an-teriores. Neste caso, iremos debrucar-nos sobre o convergencia do metodo deGauss-Seidel. Ja vimos que a matriz deste sistema nao tem a diagonal estri-tamente dominante por linhas. Por conseguinte, o teorema 3.10 nao e, nestecaso, aplicavel. Vejamos se e possıvel aplicar directamente o teorema 3.7.
89
A matriz Cgs, de acordo com (3.151), tem a forma:
Cgs =
0 −1/2 00 −1/4 −1/20 −1/8 −1/4
. (3.182)
Os valores proprios desta matriz sao as raizes da equacao
λ3 +λ2
2= 0.
Determinando estas raizes, obtem-se
λ1 = λ2 = 0, λ3 = −1
2.
Por conseguinte, o raio espectral de Cgs e
rσ(Cgs) = |λ3| =1
2.
Logo, pelo teorema 3.7, podemos concluir que o metodo de Gauss-Seidelconverge para a solucao do sistema considerado, qualquer que seja a aprox-imacao inicial.
3.4.6 Rapidez de Convergencia dos Metodos Iterativos. Analisedo Erro
Nos paragrafos anteriores, estudamos as condicoes que garantem a con-vergencia dos metodos iterativos de Jacobi e de Gauss-Seidel. Atendendo aosresultados ja obtidos, vamos agora classificar estes metodos e compara-losquanto a rapidez de convergencia. Considerando qualquer norma vectorial Ve a norma matricial M , a ela associada, podemos afirmar que, para qualquermetodo iterativo que verifique a igualdade (3.154), se verifica
‖e(k+1)‖V ≤ ‖C‖M ‖e(k)‖V . (3.183)
A rapidez de convergencia depende, naturalmente, das propriedades da ma-triz C e da aproximacao inicial escolhida. Nalguns casos especiais, podeacontecer que a solucao exacta seja obtida com um numero finito de it-eracoes . No entanto, na maioria dos casos com interesse pratico, verifica-seque a ordem de convergencia dos metodos aqui analisados e precisamente 1,ou seja, a convergencia e linear. Como sabemos, a rapidez de convergencia
90
de metodos da mesma ordem e caracterizada pelo factor assimpotico deconvergencia. Para avaliar esse factor, recorre-se frequentemente ao limite
c1 = limk→∞
‖e(k+1)‖V‖e(k)‖V
. (3.184)
A existencia do limite c1 depende das propriedades da matriz C e da normaV considerada. Alem disso, para a mesma matriz C, o limite pode ter difer-entes valores, conforme a aproximacao inicial escolhida. Pode-se mostrarque, se a matriz C tiver um unico valor proprio λ ∈ IR, tal que |λ| = rσ(C) ,entao, para a maioria das aproximacoes iniciais, o limite c1 existe e verifica-se c1 = rσ(C) . Logo, se o limite c1 existir e o metodo iterativo convergir,entao 0 < c1 < 1 e este valor pode ser tomado como o factor assimptotico deconvergencia. Isto significa que, para valores de c1 proximos de 0, teremosconvergencia rapida, enquanto que para valores de c1 proximos de 1 teremosconvergencia lenta (isto e, sao necessarias muitas iteracoes para atingir umadada precisao ). Na pratica, o valor de c1 nao pode ser obtido directamenteda formula (3.184), uma vez que os valores ‖e(k+1)‖V e ‖e(k)‖V nao sao, emgeral, conhecidos para nenhuma iterada. Por isso, recorre-se frequentementea igualdade
x(k+1) − x(k) = −e(k+1) + e(k) = −C e(k) + C e(k−1) = C(x(k) − x(k−1)).(3.185)
Desta igualdade depreende-se que a diferenca entre iteradas sucessivas variacom k do mesmo modo que o erro e(k) (ambas estas grandezas satisfazemuma relacao do tipo (3.154). Logo, se o limite (3.184) existir, tambem existeo limite
c′1 = limk→∞
‖x(k+1) − x(k)‖V‖x(k) − x(k−1)‖V
. (3.186)
e os dois limites (c1 e c′1) tem o mesmo valor, para a maioria das aprox-imacoes iniciais. A formula (3.186) pode ser utilizada na pratica para avaliarc′1 e, logo, c1. Com esse fim, calcula-se, para sucessivos valores de k, a razao
r(k) =‖x(k+1) − x(k)‖V‖x(k) − x(k−1)‖V
,
ate que o seu valor estabilize. O numero assim obtido e tomado comouma estimativa de c1. Por outro lado, os valores de r(k) tambem podem serutilizados para obter estimativas do erro e(k). Na realidade, se considerarmosum valor c2 tal que r(k) ≤ c2, ∀k > k0 (aqui k0 representa a ordem a partir
91
da qual o valor de r(k) estabiliza), podemos esperar que, para k > k0, severifique
‖e(k+1)‖V = ‖x(k+1) − x‖V ≤ c2 ‖x(k) − x‖V . (3.187)
Por outro lado, da desigualdade triangular, temos
‖x(k) − x‖V ≤ ‖x(k) − x(k+1)‖V + ‖x(k+1) − x‖V (3.188)
Das desigualdades (3.188) e (3.187) resulta
‖x(k) − x‖V ≤ ‖x(k) − x(k+1)‖V + c2‖x(k) − x‖V , (3.189)
de onde‖x(k) − x‖V (1 − c2) ≤ ‖x(k) − x(k+1)‖V . (3.190)
Uma vez que c2 < 1, por construcao , da desigualdade (3.190) obtem-se
‖e(k)‖V = ‖x(k) − x‖V ≤‖x(k) − x(k+1)‖V
1 − c2. (3.191)
De (3.191), utilizando (3.187), obtem-se tambem
‖e(k+1)‖V = ‖x(k+1) − x‖V ≤c2
1 − c2‖x(k) − x(k+1)‖V . (3.192)
A desigualdade (3.192) permite-nos majorar o erro de uma dada iterada,bastando para isso conhecer a diferenca entre as duas ultimas iteradas e ovalor de c2.
Exemplo 3.12 . Regressando, mais uma vez, ao sistema linear quefoi analisado no exemplo 3.10, vamos efectuar uma analise do erro, para osmetodos de Jacobi e de Gauss-Seidel. Com este fim, foi aplicado cada umdestes metodos ao sistema considerado. Partindo da aproximacao inicialx(0) = (0.5, 0.8, 1.0), foram efectuadas iteracoes ate que fosse satisfeita acondicao
‖x(k) − x(k+1)‖2 ≤ 0.01.
Em cada iteracao foi avaliada a norma ‖x(k) − x(k+1)‖2 e, a partir da 2a
iteracao , a razao r(k) correspondente. Os resultados obtidos no caso dometodo de Jacobi sao dados na tabela 3.1, enquanto os resultados obtidospara o metodo de Gauss-Seidel se encontram na tabela 3.2. Nestas tabelasverifica-se nitidamente que os valores de r(k) tendem para c1 = 0.7071, no
92
k x(k)1 x
(k)2 x
(k)3 ‖x(k+1) − x(k)‖2 r(k)
1 0.6 0.75 0.9 0.152 0.625 0.85 0.875 0.106066 0.70710643 0.575 0.875 0.925 0.07500 0.70710664 0.5625 0.825 0.9375 0.05303 0.70710695 0.5875 0.8125 0.9125 0.03750 0.70710686 0.59375 0.8375 0.90625 0.02652 0.70710837 0.58125 0.84375 0.91875 0.01875 0.70710758 0.578125 0.83125 0.921875 0.01326 0.70710619 0.584375 0.828125 0.915625 0.00938 0.7071068
Table 1: Resultados do metodo de Jacobi para o exemplo 3.12
k x(k)1 x
(k)2 x
(k)3 ‖x(k+1) − x(k)‖2 r(k)
1 0.6 0.8 0.9 0.1414212 0.6 0.85 0.925 0.055902 0.39528463 0.575 0.825 0.9125 0.037500 0.67081874 0.5875 0.8375 0.91875 0.018750 0.55 0.58125 0.83125 0.915625 0.009375 0.5
Table 2: Resultados do metodo de Gauss-Seidel para o exemplo 3.12
93
caso do metodo de Jacobi, e para c1 = 0.5, no caso do metodo de Gauss-Seidel. Estes valores coincidem com os raios espectrais das matrizes CJ eCgs, respectivamente (ver exemplos 3.10 e 3.11). Com base nestes valores,podemos obter estimativas do erro para cada um dos metodos.
Para o metodo de Jacobi, de acordo om a formula (3.191), considerandoc2 = 0.70711, temos
‖e(9)‖2 ≤c2
1 − c2‖x(9) − x(8)‖2 ≤ 0.0242.
No caso do metodo de Gauss-Seidel, tomando c2 = 0.5, temos
‖e(5)‖2 ≤c2
1 − c2‖x(5) − x(4)‖2 ≤ 0.01.
No exemplo que acabamos de ver, constatamos que o metodo de Gauss-Seidel convergia mais rapidamente que o de Jacobi, o que resulta de o raioespectral da matriz Cgs ser inferior ao da matriz CJ . Vamos ver que este casonao e uma excepcao , no que diz respeito a relacao entre os dois metodos.
A fim de compararmos o metodo de Gauss-Seidel com o de Jacobi, quantoa rapidez de convergencia, consideremos o caso em que a matriz A do sistematem a diagonal estritamente dominante por linhas. Neste caso, de acordocom os teoremas 3.8 e 3.10, ambos os metodos convergem para a solucaoexacta, qualquer que seja a aproximacao inicial escolhida. Alem disso, parao metodo de Jacobi e valida a estimativa do erro
‖e(k)‖∞ ≤ µk ‖e(0)‖∞, k = 1, 2, . . . (3.193)
onde µ = ‖CJ‖∞ . Recordando a forma da matriz CJ , dada por (3.143), eas definicoes das grandezas αi e βi, dadas por (3.174), podemos concluir que
µ = maxi=1,...,n
(αi + βi). (3.194)
Por outro lado, para o metodo de Gauss-Seidel, segundo o teorema 3.10, evalida a estimativa do erro
‖e(k)‖∞ ≤ ηk ‖e(0)‖∞, k = 1, 2, . . . (3.195)
Para estabelecer uma relacao entre a rapidez de convergencia dos dois metodos,basta-nos portanto comparar o parametro µ da formula (3.193) com o parametroη da formula (3.195).
94
Exemplo 3.13. Consideremos o sistema Ax = b, onde A e uma matrizde ordem n com a forma geral
A =
5 2 0 . . . 02 5 2 . . . 0. . . . . . . . . . . . . . .0 . . . 2 5 20 . . . 0 2 5
.Neste caso, verifica-se imediatamente que a matriz A tem a diagonal estri-tamente dominante, por linhas, pelo que tanto o metodo de Gauss-Seidelcomo o de Jacobi convergem, qualquer que seja a aproximacao inicial. Alemdisso, atendendo as formulas (3.174), temos
αi = βi = 2/5, i = 1, 2, . . . , n.
Daqui, pelas formulas (3.194) e (3.175), obtem-se imediatamente
µ = 4/5, η = 2/3.
Assim, neste exemplo verifica-se η < µ. Por conseguinte, e de esperar que,no caso deste sistema, o metodo de Gauss-Seidel convirja mais rapidamenteque o de Jacobi.
Note-se, porem, que esta comparacao entre os dois metodos so e validapara matrizes com a diagoal estritamente dominante por linhas. No casogeral, nem sempre o metodo de Gauss-Seidel e mais rapido que o de Jacobi,havendo mesmo casos particulares em que o segundo e convergente e oprimeiro nao .
Um aspecto importante a salientar e que os metodos iterativos parasistemas lineares, quando convergem para qualquer aproximacao inicial,sao estaveis (ver Definicao 3.7). Ou seja, partindo de dois vectores ini-ciais proximos, ξ0 e η0 , obtem-se sempre duas sucessoes {x(n)} e {y(n)}igualmente proximas, convergindo para o mesmo vector x (solucao exacta).Esta propriedade e de grande importancia pratica, uma vez que no calculonumerico sao inevitaveis os erros se arredondamento, os quais, como ja vi-mos, podem propagar-se ao longo de uma sucessao de operacoes , conduzindoa erros muito grandes no resultado final. Esta situacao verifica-se, por ex-emplo, na resolucao de sistemas lineares por metodos directos, mesmo queeles sejam bem condicionados. Os metodos iterativos, desde que sejam apli-cados a sistemas bem condicionados, sao sempre estaveis, ou seja, quando seusam estes metodos nao ha perigo de os erros de arredondamento cometi-dos nos calculos poderem resultar em erros significativos no resultado final.
95
Isto representa, portanto, uma importante vantagem dos metodos iterativossobre os directos, sobretudo quando se trata de resolver sistemas de grandesdimensoes .
3.4.7 Metodo das relaxacoes sucessivas (SOR)
Neste paragrafo estudaremos uma generalizacao do metodo de Gauss-Seidel,muito utilizada no calculo numerico, que e conhecida como metodo das re-laxacoes sucessivas ou metodo SOR (sucessive over-relaxations). Nestecaso, estamos mais propriamente a falar de uma famılia de metodos, quedepende de um parametro (ω), cuja matriz de iteracao e dada pela seguinteformula:
C(ω) = −M−1(ω)N(ω),
onde
M(ω) = L+1
ωD, N(ω) = U + (1− 1
ω)D, (3.196)
sendo as matrizes L,D e U definidas como no caso dos metodos de Jacobi ede Gauss-Seidel.
E facil verificar que, no caso de ω = 1, se obtem M(1) = L + D eN(1) = U , pelo que C(1) = −(L + D)−1U , ou seja, neste caso o metodoSOR reduz-se ao metodo de Gauss-Seidel.
Na pratica, cada iterada x(k+1) do metodo SOR calcula-se atraves daseguinte formula:
x(k+1) = ωz(k+1) + (1− ω)x(k), (3.197)
ondez(k+1) = D−1
(b− Lx(k+1) − Ux(k)
). (3.198)
Comparando as formulas, verifica-se que z(k+1) coincide com a (k+1)-esimaiterada do metodo de Gauss-Seidel. Assim, podemos dizer que cada iteradado metodo SOR e uma media ponderada entre a nova iterada (pelo metodode Gauss-Seidel) e a iterada anterior, sendo ω o peso da nova iterada. Maisuma vez se confirma que, no caso de ω = 1, o metodo SOR se reduz aometodo de Gauss-Seidel. Ao introduzir o parametro ω, isso permite-nosacelerar a convergencia do metodo iterativo, pois e possıvel, para cada sis-tema em particular, escolher esse parametro de tal modo que a convergenciaseja a mais rapida possıvel. Isso consegue-se teoricamente (estudando o raioespectral da matiz C(ω) como funcao de ω) ou experimentalmente, testandodiferentes valores de ω.
Pode ser demonstrado o seguinte teorema sobre a convergencia do metodoSOR, no caso de um sistema linear Ax = b:
96
Teorema 3.11. Sendo A uma matriz simetrica e definida positiva, ometodo SOR converge para a solucao do sistema se e so se ω ∈]0, 2[.
Deste teorema resulta, em particular, que o metodo de Gauss-Seidel(ω = 1) e sempre convergente para matrizes simetricas e definidas positivas.
3.5 Metodos Iterativos para Sistemas Nao- lineares
Consideremos agora um sistema de n equacoes nao lineares, da formaF1(x1, x2, . . . , xn) = 0,F2(x1, x2, . . . , xn) = 0,
. . .Fn(x1, x2, . . . , xn) = 0,
(3.199)
onde cada uma das funcoes Fi e uma funcao real de n variaveis reais. Estesistema pode ser escrito na forma vectorial:
F(x) = 0,
onde F = (F1, F2, . . . , Fn), x = (x1, x2, . . . , xn). O ponto z ∈ IRn diz-sesolucao (ou raiz) do sistema (3.199) se F(z) = 0. O problema da deter-minacao das razes de um sistema nao linear e complexo (muito mais com-plexo do que no caso de sistemas lineares). Em primeiro lugar, nao hanenhum criterio simples que nos permita verificar se o sistema (3.199) temou nao solucao. No caso de varias solucoes, nao e facil isolar cada uma,isto e, definir um conjunto em IRn que contenha essa raiz, e nao contenhaoutras.
Uma das abordagens para conseguir localizar raizes de um sistema naolinear e baseada no teorema do ponto fixo (que pode ser reformulado parafuncoes de IRn em IRn). O mesmo teorema permite-nos definir um metodoiterativo (metodo do ponto fixo em IRn) para aproximar as raizes do sistema.Finalmente, veremos que o metodo de Newton (estudado em 2.4, no caso deuma equacao) tambem pode ser generalizado para sistemas de n equacoes.
3.5.1 Metodo do Ponto Fixo em IRn
A fim de investigar as condicoes que garantem a convergencia do metododo ponto fixo em IRn, vamos formular uma generalizacao do teorema doponto fixo, estudado em 2.3. Com esse objectivo, necessitamos de introduziralgumas definicoes.
97
Definicao 3.9. Seja E um espaco normado, X ⊂ E e G uma funcaode X em E. A funcao G diz-se lipschitziana em X, se e so se existir umaconstante q, tal que
‖G(x1) − G(x2)‖ ≤ q ‖x1 − x2‖, ∀x1, x2 ∈ X. (3.200)
Ao ınfimo de todas as constantes q, para as quais a desigualdade (3.200)se verifica, chama-se constante de Lipschitz de G em X e representa-se porLG,X .
Definicao 3.10. Diz-se que G e uma contraccao (ou uma funcao con-tractiva) em X se e so se G for lipschitziana em X e LG,X < 1.
Exemplo 3.14. Seja E = IR e G(x) = x2. Verifiquemos para quevalores de r a funcao G e contractiva em X = [−r, r]. Temos
|G(x1) − G(x2)| = |x21 − x2
2| = |x1 − x2|.|x1 + x2|. (3.201)
Se x1 e x2 pertencerem a X, podemos escrever
|x1 + x2| ≤ r + r = 2r. (3.202)
Substituindo (3.202) em (3.201), obtem-se
|G(x1) − G(x2)| ≤ 2r |x1 − x2|, (3.203)
donde se conclui que G e lipschitziana em X, com a constante LG,X = 2r.Por conseguinte, se r < 1/2, podemos afirmar que G e contractiva em X.
No caso de funcoes de uma variavel real, a condicao de contractividadepode ser expressa noutros termos, tornando-se mais facil a sua verificacao.
Teorema 3.12. Seja G uma funcao com domınio em X = [a, b] e valoresreais e suponhamos que G ∈ C1([a, b]). Entao G e contractiva em X se eso se
maxx∈[a,b]
|G′(x)| < 1. (3.204)
98
Demonstracao. Pelo teorema de Lagrange, quaisquer que sejam x1 ex2, pertencentes a [a, b], existe ξ ∈ [x1, x2], tal que
|G(x1) − G(x2)| = |G′(ξ)||x1 − x2|. (3.205)
Entao podemos afirmar que a constante de Lipschitz de G e
LG = maxx∈[a,b]
|G′(x)| < 1, (3.206)
donde se conclui que G e contractiva em [a, b].Demonstremos agora a inversa. Suponhamos que existe em [a, b] um
certo y, tal que |G′(y)| ≥ 1. Entretanto, pelo teorema de Lagrange, paraqualquer h > 0, existe θ ∈ [0, 1] tal que
|G(y + h) − G(y)| = |G′(y + θh)|h. (3.207)
Visto que G′ e contınua em [a, b],para qualquer ρ < 1, existe h0 tal que|G′(y+ θh0)| > ρ. Escrevendo a desigualdade (3.207) com h = h0, obtem-se
|G(y + h0) − G(y)| = |G′(y + θh0)|h0 > ρh0. (3.208)
A desigualdade (3.208) implica que G nao e contractiva em [a, b], ficandoassim demonstrado o teorema.
O teorema 3.12 permite-nos substituir a condicao de contractividade pelacondicao (3.204), quando se consideram funcoes de uma variavel. Foi issoprecisamente o que fizemos na sec. 2.3. No caso em que G e uma funcao deIRn em IRn, com derivadas parciais contınuas, a contractividade pode serverificada atraves da sua matriz jacobiana.
Definicao 3.11. Seja G uma funcao vectorial, tal que
G(x) = (G1(x), G2(x), . . . , Gn(x)) ,
onde Gi e uma funcao escalar com domınio em X ⊂ IRn. Se em X exis-tirem as derivadas parciais ∂Gi
∂xj, i, j ∈ {1, ..., n} chama-se jacobiana de G e
representa-se por JG a matriz
JG(x) =
∂G1∂x1
. . . ∂G1∂xn
∂G2∂x1
. . . ∂G2∂xn
. . . . . . . . .∂Gn∂x1
. . . ∂Gn∂xn
. (3.209)
99
Teorema 3.13. Seja X um conjunto convexo 2 em IRn e G uma funcaovector com domınio em X e valores em IRn. Entao, se as derivadas parciais∂Gi∂xj
, i, j ∈ {1, ..., n} forem contınuas em X e se
supx∈X‖JG(x)‖∞ < 1, (3.210)
G e contractiva em X (segundo a norma do maximo).
Demonstracao . Sejam x1 e x2 dois elementos de X. Segundo oteorema de Lagrange para funcoes de n variaveis, para cada funcao Gi existeum ponto ξi, pertencente ao segmento [x1, x2], tal que
Gi(x1)−Gi(x2) = (∇Gi(ξi), x1 − x2) , (3.211)
onde
∇Gi(x) =
(∂Gi∂x1
, . . . ,∂Gi∂xn
), i ∈ {1, . . . , n}, (3.212)
e os parentesis curvos no segundo membro de (3.211) denotam o produtointerno de vectores. Note-se que todos os pontos ξi pertencem a X, uma vezque este conjunto e convexo. De (3.211) e (3.212) obtem-se
|Gi(x1)−Gi(x2)| ≤ maxj=1,...,n |x1,j − x2,j |∑n
j=1
∣∣∣∂Gi∂xj(ξi)∣∣∣ = ‖∇Gi(ξi)‖1 ‖x1 − x2‖∞,
i ∈ {1, . . . , n}.(3.213)
Seja i′ o ındice para o qual se verifica
|Gi′(x1)−Gi′(x2)| = ‖G(x1)−G(x2)‖∞.
No caso de i = i′, a desigualdade (3.213) toma o aspecto
‖G(x1)−G(x2)‖∞ ≤ ‖∇Gi′(ξi′)‖1 ‖x1 − x2‖∞. (3.214)
Atendendo a que
‖∇Gi′(ξi′)‖1 ≤ max1≤i≤n
n∑j=1
∣∣∣∣∂Gi∂xj(ξi′)
∣∣∣∣ = ‖JG(ξi′)‖∞ < 1, (3.215)
2Um conjunto X diz-se convexo se, para quaisquer x1, x2 pertencentes a X, todos ospontos do segmento [x1, x2] tambem pertencerem a X.
100
de (3.214) resulta que G e contractiva em X, segundo a norma do maximo.
Nalguns casos, pode ser mais comodo considerar em IRn outras normasque nao a do maximo, por exemplo, a norma ‖.‖1. Por isso enunciamos aquio seguinte teorema, analogo a 3.13.
Teorema 3.14. Seja X um conjunto convexo em IRn e G uma funcaovectorial com domınio em X e valores em IRn. Entao , se as derivadasparciais ∂Gi
∂xj, i, j ∈ {1, ..., n} forem contınuas em X e se
supx∈X‖JG(x)‖1 < 1, (3.216)
G e contractiva em X (segundo a norma ‖‖1).
Para nao nos tornarmos repetitivos, nao apresentamos a demonstracaodeste teorema, que pode ser obtida por um raciocınio semelhante ao doteorema anterior.
Estamos agora em condicoes de formular o teorema do ponto fixo, paraespacos normados de dimensao finita, por exemplo, os espacos IRn.
Teorema 3.15 (Teorema do ponto fixo em Rn). Seja E um espaconormado de dimensao finita e X um sub-conjunto fechado de E. Seja Guma funcao contractiva em X, tal que G(X) ⊂ X. Entao sao verdadeirasas afirmacoes
1. G tem um unico ponto fixo z em X.
2. Se {x(k)}k∈IN for uma sucessao de termos em E tal que x(0) ∈ X ex(k+1) = G(x(k)),∀k ∈ IN, entao x(k) → z.
3. Se G satisfaz em X a desigualdade (3.200) com a constante q < 1,entao sao validas as desigualdades
‖x(n+1) − z‖ ≤ q‖x(n) − z‖, ∀n ∈ IN. (3.217)
e
‖x(m) − z‖ ≤ qm
1− q‖x(1) − x(0)‖,∀m ∈ IN. (3.218)
101
Demonstracao . Em primeiro lugar, note-se que se x(0) ∈ X, entaox(k) ∈ X, ∀k, visto que G(X) ⊂ X. Comecemos por provar que a sucessaode que se fala no ponto 2 e convergente. Para tanto, basta provar que ela euma sucessao de Cauchy. Uma vez que G e contractiva em X, existe umaconstante q < 1, tal que
‖G(x1)−G(x2)‖ ≤ q ‖x1 − x2‖, ∀x1, x2 ∈ X. (3.219)
Em particular, se tivermos dois termos consecutivos da sucessao consideradaverifica-se
‖x(k+1) − x(k)‖ ≤ q ‖x(k) − x(k−1)‖, ∀k ∈ IN. (3.220)
Sejam agora x(m) e x(n) dois termos quaisquer da sucessao , com n > m.Entao podemos escrever
‖x(n) − x(m)‖ = ‖x(n) − x(n−1) + x(n−1) − x(n−2) + . . . + x(m+1) − x(m)‖ ≤‖x(n) − x(n−1)‖ + ‖x(n−1) − x(n−2)‖ + · · ·+ ‖x(m+1) − x(m)‖.
(3.221)Das desigualdades (3.220) e (3.221) obtem-se
‖x(n) − x(m)‖ ≤(qn−m−1 + · · ·+ q + 1
)‖x(m+1) − x(m)‖ ≤
qm(qn−m−1 + · · ·+ q + 1)
)‖x(1) − x(0)‖. (3.222)
A soma que figura no segundo membro de (3.222) e a soma de uma progressaogeometrica de razao q. Como q < 1, podemos escrever
qmn−m−1∑k=0
qk <qm
1− q, ∀n ∈ IN. (3.223)
Substituindo (3.223) em (3.222), obtem-se
‖x(m) − x(n)‖ < qm
1− q‖x(1) − x(0)‖, ∀n > m. (3.224)
Da desigualdade (3.224) resulta que ∀ε > 0, existe n0 ∈ IN tal que
‖x(m) − x(n)‖ < ε, ∀m,n > n0. (3.225)
Daqui resulta, por definicao que a sucessao considerada e uma sucessao deCauchy e, logo, converge. Representemos por z o seu limite. Uma vez que Xe fechado, z ∈ X. Provemos agora que z e um ponto fixo de G. Utilizandoo facto de G ser contractiva, podemos escrever
‖x(m+1) −G(z)‖ = ‖G(x(m))−G(z)‖ ≤ q‖x(m) − z‖, ∀m. (3.226)
102
Logo ‖x(m+1) − G(z)‖ → 0, de onde resulta que x(m) → G(z), quandom → ∞. Por conseguinte, G(z) = z, tal como se pretendia demonstrar.Fica assim demonstrado o ponto 2 do teorema. A desigualdade (3.217), porsua vez, resulta de (3.226). Quanto a desigualdade (3.218), ela obtem-se de(3.224), se fizermos n tender para infinito. Resta-nos provar que z e o unicoponto fixo de G em X. Suponhamos que existem dois pontos fixos de G emX e representemo-los por z e z′. Uma vez que G e contractiva, temos
‖G(z′)−G(z)‖ = ‖z′ − z‖ ≤ q‖z′ − z‖, (3.227)
donde‖z′ − z‖ (1 − q) ≤ 0. (3.228)
Uma vez que 1− q > 0, de (3.228) resulta que z′ = z.
Exemplo 3.15. Consideremos o sistema de duas equacoes{3x1 + x2
2 = 0x2
1 + 3x2 = 1(3.229)
Vamos utilizar o teorema do ponto fixo para provar que este sistema temuma unica raiz no conjunto
X = {(x1, x2) ∈ IR2 : −1/3 ≤ x1 ≤ 0; 0 ≤ x2 ≤ 1/3} (3.230)
Com esse fim, o sistema (3.229) pode ser reescrito na forma x = G(x), onde
G1(x1, x2) = −x223
G2(x1, x2) =1−x21
3 .
(3.231)
Verifiquemos se a funcao G, definida por (3.231), satisfaz as condicoes doteorema do ponto fixo em X. Em primeiro lugar, constata-se que o conjuntoX e um quadrado, contendo a sua fronteira, pelo que e convexo e fechado.Alem disso, vemos que as derivadas parciais de G1 e G2 sao contınuas emX, de tal modo que a jacobiana de G tem a forma
JG(x1, x2) =
0 −2x23
−2x13 0
(3.232)
103
Assim, podemos escrever
‖JG(x1, x2)‖∞ = max
(2|x2|
3,2|x1|
3
)(3.233)
de onde resulta que
‖JG(x1, x2)‖∞ ≤2
9< 1 ∀(x1, x2) ∈ X. (3.234)
Com base no teorema 3.13 podemos, portanto, afirmar que G e contractivaem X (segundo a norma do maximo) com a constante q = 2
9 . Para podermosaplicar o teorema do ponto fixo, precisamos tambem de verificar queG(X) ⊂ X.
Seja x = (x1, x2) ∈ X. Entao
G1(x1, x2) = −x223 ∈ [−1/3, 0]
G2(x1, x2) =1−x21
3 ∈ [0, 1/3].
(3.235)
Por conseguinte, temos (G1(x1, x2), G2(x1, x2)) ∈ X, de onde se conclui queG(X) ⊂ X.
Logo, a funcao G satisfaz as condicoes do teorema do ponto fixo em X.Assim, podemos garantir que aquela funcao tem um unico ponto fixo emX, que, por construcao , sera a unica raiz do sistema (3.229) no referidoconjunto.
Podemos tambem aplicar o metodo do ponto fixo para aproximar a raizconsiderada. Tomemos como aproximacao inicial qualquer ponto do con-junto X, por exemplo, a origem das coordenadas: x(0) = (0, 0). Entaoobtem-se as seguintes aproximacoes sucessivas:
x(1)1 = G1(0, 0) = 0
x(1)2 = G2(0, 0) = 1
3
, (3.236)
x(2)1 = G1(0, 1/3) = − 1
27
x(2)2 = G2(0, 1/3) = 1
3
. (3.237)
Tentemos agora obter uma estimativa do erro da iterada x(2). De acordocom a desigualdade (3.218), podemos escrever
‖x(2) − z‖∞ ≤q2
1− q‖x(1) − x(0)‖∞, (3.238)
104
onde q = 29 . Neste caso, temos ‖x(1) − x(0)‖∞ = 1/3; logo
‖x(2) − z‖∞ ≤4
63
1
3=
4
189. (3.239)
Esta ultima estimativa poderia ser refinada se, em vez da desigualdade(3.218), aplicassemos a desigualdade
‖x(m+1) − z‖∞ ≤q
1− q‖x(m+1) − x(m)‖∞, (3.240)
que tambem se pode deduzir do teorema do ponto fixo. Nesse caso obterıamos
‖x(2) − z‖∞ ≤q
1− q‖x(2) − x(1)‖∞ =
2
189. (3.241)
3.5.2 Metodo de Newton para Sistemas de Equacoes
Como vimos em 2.3, o metodo de Newton pode ser considerado como umcaso particular do metodo do ponto fixo. No caso de uma funcao f (de umavariavel real), a funcao iteradora do metodo d Newton, como sabemos, tema forma:
g(x) = x− f(x)
f ′(x). (3.242)
Ao considerar o sistema (3.199), em vez da funcao f temos uma funcaovectorial F (de n variaveis). Admitindo que todas as derivadas parciais deF existem e sao contınuas num certo conjunto D, onde se procura a raiz dosistema, podemos definir do seguinte modo a matriz jacobiana de F :
JF (x) =
∂F1∂x1
∂F1∂x2
. . . ∂F1∂xn
∂F2∂x1
∂F2∂x2
. . . ∂F2∂xn
. . . . . . . . . . . .∂Fn∂x1
∂Fn∂x2
. . . ∂Fn∂xn
. (3.243)
Admitamos ainda que JF e invertıvel no domınio considerado. Quandose generaliza o metodo de Newton para sistemas de equacoes, e naturalsubstituir na formula (3.242) a expressao 1/f ′(x) pela inversa da matrizjacobiana. Obtem-se assim formalmente a seguinte funcao iteradora para ometodo de Newton, aplicado a sistemas nao-lineares:
G(x) = x− J−1F (x).F (x). (3.244)
Daqui resulta a formula iteradora do metodo de Newton para sistemas nao-lineares:
x(k+1) = G(x(k)) = x(k) − J−1F (x(k)).F (x(k)), (3.245)
105
onde x(k) ∈ IRn representa, como habitualmente, a k-esima iterada dometodo (k = 0, 1, 2 . . . ). A formula (3.245), embora definindo o metodode Newton para sistemas, nao e a que geralmente se aplica nos calculos.Do ponto de vista computacional, nao e vantajoso aplicar directamente estaformula, ja que isso nos obrigaria, em cada iteracao do metodo, a inverteruma matriz de ordem n (a jacobiana de F ), o que seria muito dispendiosoem termos de numero de operacoes. Na realidade, isso pode ser evitado,como vamos ver. Comecemos por rescrever a formula iteradora na forma
x(k+1) − x(k) = −J−1F (x(k)).F (x(k)), (3.246)
Introduzindo a notacao ∆x(k) = x(k+1) − x(k) e multiplicando ambos osmembros de (3.246), a esquerda, por JF (x(k)), obtem-se
JF (x(k))∆x(k) = −F (x(k)). (3.247)
Esta formula, juntamente com
x(k+1) = ∆x(k) + x(k), (3.248)
define um processo iterativo, equivalente ao da formula (3.246), mas ondenao aparece a inversa da jacobiana. Na pratica, em vez de inverter a Jaco-biana de F, basta-nos em cada iteracao resolver o sistema linear (3.247), cujamatriz e essa jacobiana. Esse sistema linear pode ser resolvido por qualquerdos metodos directos ou iterativos que estudamos nas seccoes anteriores e,como sabemos, a sua resolucao implica muito menos operacoes do que a in-versao da matriz correspondente. Uma vez resolvido o sistema, a sua solucao∆x(k) da-nos a ”correccao” que, somada a iterada anterior, permite obter anova iterada x(k+1) (formula 3.248). Este processo e repetido ate que se veri-fique uma das seguintes condicoes (ou as duas): ‖∆x(k)‖ < ε, ‖F (x(k))‖ < ε,sendo ε uma margem de erro previamente estabelecida. Nestas condicoes deparagem pode usar-se qualquer das normas vectoriais estudadas.
Exemplo 3.16. Consideremos de novo o sistema de duas equacoes{3x1 + x2
2 = 0x2
1 + 3x2 = 1(3.249)
Partindo da aproximacao inicial x(0) = (0, 0), vamos efectuar duas iteracoesdo metodo de Newton para aproximar a sua solucao. Neste caso, temos
F1(x1, x2) = 3x1 + x22,
F2(x1, x2) = x21 + 3x2 − 1.
(3.250)
106
A matriz jacobiana de F tem a forma
JF (x1, x2) =
[3 2x2
2x1 3
]. (3.251)
Logo, para a primeira iteracao, temos
JF (x(0))∆x(0) = −F (x(0)), (3.252)
onde
JF (x(0)) =
[3 00 3
], (3.253)
F (x(0)) = (F1(0, 0), F2(0, 0)) = (0,−1). (3.254)
Resolvendo o sistema (3.252), obtem-se
∆x(0) = (0, 1/3). (3.255)
Logo,x(1) = ∆x(0) + x(0) = (0, 1/3). (3.256)
Passemos a segunda iteracao, com o sistema linear
JF (x(1))∆x(1) = −F (x(1)), (3.257)
onde
JF
(x(1)
)=
[3 2
30 3
], (3.258)
F (x(1)) = (F1(0, 1/3), F2(0, 1/3)) = (1/9, 0). (3.259)
Resolvendo agora o sistema (3.257), obtem-se
∆x(1) = (−1/27, 0). (3.260)
Finalmente, a segunda iterada obtem-se somando:
x(2) = ∆x(1) + x(1) = (−1/27, 1/3). (3.261)
Note-se que, embora neste caso as duas primeiras iteracoes do metodo deNewton coincidam com as do metodo do ponto fixo, isto nao e o que acon-tece em geral. Em regra, tal como acontece no caso de n = 1, o metodo deNewton, quando converge, define uma sucessao de aproximacoes com con-vergencia quadratica, enquanto o metodo do ponto fixo apresenta apenasconvergencia linear. Assim, de uma maneira geral, o metodo de New-ton, com o mesmo numero de iteracoes, permite atingir um resultado maispreciso.
107
4 Aproximacao de funcoes
4.1 Interpolacao Polinomial
Neste capıtulo dedicar-nos-emos a aproximacao de funcoes. Um dos metodosmais classicos e simples de aproximacao de funcoes e a interpolacao. Oobjectivo da interpolacao e reconstruir num certo intervalo [a, b] uma funcaof , cujos valores sao conhecidos apenas num nunmero finito de pontos desseintervalo. Esses pontos sao os chamados nos de interpolacao e vamosrepresenta-los genericamente por xi (i = 0, 1, 2, ..., n). Assim, pode dizer-seque os dados para a interpolacao sao sempre constituıdos por uma tabelade valores de f :
x0 x1 . . . xnf0 f1 . . . fn
onde fi = f(xi) representa o valor de f no no de interpolacao xi (i =0, 1, . . . , n).
Para exemplificar com uma aplicacao bem simples, consideremos a seguintetabela, que representa os valores da populacao de uma especie (Ni, em mil-hares) determinados em diferentes instantes ti:
Tabela 4.1ti 10 12 14 16
Ni 10 15 22 18
Na fig. 13 esta representada graficamente esta tabela. Suponhamos que o
nosso objectivo e reconstruir a funcao N no intervalo [10, 16]. Claro que,de um modo geral, existem muitas diferentes maneiras de ”interpolar” estesdados. Na mesma figura, estao representadas duas das possıveis funcoesinterpoladoras. Se usarmos a funcao interpoladora N1, por exemplo, a es-timativa para a populacao em t = 15 e N1(15) = 20, enquanto que seusarmos a funcao N2, o valor correspondente e N2(15) = 22.188. Como foimencionado no inıcio, neste capıtulo iremos estudar um tipo especıfico deinterpolacao, a interpolacao polinomial, o que significa que iremos consid-erar sempre como funcoes interpoladoras polinomios. Na proxima seccao,definiremos polinomio interpolador e demonstraremos a existencia e unici-dade de tal polinomio.
4.1.1 Existencia e Unicidade de Polinomio Interpolador
Comecamos por formular a seguinte definicao:
108
11 12 13 14 15 16t
12
14
16
18
20
22
T
Figure 13: Exemplo de interpolacao : funcao N1 (cinzento), funcao N2
(preto)
Definicao 4.1. Chama-se polinomio interpolador da funcao f nos nosx0, ..., xn o polinomio P , de grau nao superior a n, que satisfaz as relacoesP (xi) = fi, i = 0, 1, . . . , n.
Dada esta definicao, a primeira questao que se poe e saber se, paraqualquer tabela, existe sempre um polinomio interpolador, e se este e unico.
No caso de dois pontos (x0, x1) e simples responder a esta questao. Comefeito, segundo a definicao 4.1, neste caso o polinomio interpolador tem graumenor ou igual a um, ou seja, e uma funcao linear. Como o grafico de umafuncao linear e uma recta, e obvio que o polinomio interpolador neste casoexiste e e unico: e a funcao P1 que tem como grafico a recta que passapelos pontos (x0, f0) e (x1, f1). Quando se considera um grau n arbitrario,o problema ja nao e tao simples, mas a resposta continua a ser positiva.
Para analisarmos o caso geral, recordemos que um polinomio de graunao superior a n pode ser escrito na forma
Pn(x) = a0 + a1x+ a2x2 + · · ·+ anx
n, (4.1)
onde os coeficientes ai (i = 0, 1, . . . , n) sao certos numeros reais. Paraconstruir o polinomio interpolador precisamos de calcular esses coeficientes.Recorrendo de novo a definicao, podemos escrever que o polinomio Pn sa-tisfaz as igualdades
Pn(x0) = a0 + a1x0 + a2x20 + · · ·+ anx
n0 = f0,
Pn(x1) = a0 + a1x1 + a2x21 + · · ·+ anx
n1 = f1,
. . .Pn(xn) = a0 + a1xn + a2x
2n + · · ·+ anx
nn = fn.
(4.2)
109
Observando o conjunto das igualdades (4.2), verificamos que elas formamum sistema de n+ 1 equacoes lineares em ordem as incognitas a0,a1,. . . ,an.Escrevendo esse sistema na forma matricial, obtem-se
1 x0 x20 . . . xn0
1 x1 x21 . . . xn1
. . . . . . . . . . . . . . .1 xn x2
n . . . xnn
a0
a1
. . .an
=
f0
f1
. . .fn
. (4.3)
A matriz deste sistema e conhecida como a matriz de Vandermonde. Va-mos representa-la por V (x0, x1, . . . , xn). Para nos certificarmos de que osistema (4.3) tem sempre solucao unica, basta verificar que esta matriz einvertıvel, ou seja, que o seu determinante e diferente de zero. Mais umavez, comecemos pelo caso de n = 1. Neste caso, e evidente que
Det V (x0, x1) = x1 − x0 6= 0,
ja que, por construcao, os nos de interpolacao sao distintos. Passando aocaso geral, pretendemos mostrar que Det V (x0, x1, . . . , xn) 6= 0, para n =1, 2, . . . . Para isso, recorremos a um resultado auxiliar da algebra (ver, porexemplo, [?], p.45) que nos diz o seguinte:
Det V (x0, x1, . . . , xn) =
n∏i,j=0
(xi − xj) (4.4)
(consideram-se todos os pares xi, xj tais que i > j). Mais uma vez se verificaque o determinante da matriz de Vandermonde e nao nulo, para qualquer n,e por conseguinte o sistema (4.3) tem sempre uma unica solucao (desde queos nos de interpolacao sejam distintos entre si). Daqui se conclui tambemque, para qualquer tabela de valores de uma funcao f , existe um unicopolinomio interpolador.
4.1.2 Formula Interpolaora de Lagrange
Uma vez esclarecida a questao teorica da existencia de polinomio interpo-lador, poe-se a questao pratica de o calcular. Para isso, sao conhecidasdiferentes formulas, que iremos estudar e comparar, do ponto de vista daeficiencia computacional. Nunca e demais lembrar que o polinomio inter-polador, para cada tabela, e unico (como acabamos de provar). Por isso,independentemente da formula que usarmos para o construir, o resultado fi-nal sera sempre o mesmo (desprezando eventuais erros de arredondamento).
110
Uma das formulas mais simples para a construcao do polinomio interpo-lador e a formula interpoladora de Lagrange. Esta formula baseia-seno facto de que os polinomios de grau nao superior a n formam um espacolinear de dimensao n + 1 (o espaco Pn). Como sabemos da albegra linear,isto significa que, sendo dados n + 1 polinomios de grau nao superior a n,linearmente independentes, qualquer outro polinomio do mesmo espaco seexprime como uma combinacao linear dos primeiros. Sendo assim, o metodode Lagrange para construir o polinomio interpolador comeca com a definicaode n+ 1 polinomios que formam uma base em Pn. Vamos representar essespolinomios por Li(x), i = 0, 1, . . . , n e designa-los polinomios de Lagrange.
Os polinomios de Lagrange sao construıdos para um certo conjunto denos de interpolacao xi, i = 0, 1 . . . , n (distintos entre si), de tal modo que enatural estabelecer uma correspondencia entre xi e Li, para cada i. Estacorrespondencia estabelece-se do seguinte modo: designamos Li o polinomiode grau n, tal que Li(xi) = 1, e Li(xj) = 0, se j ∈ {0, 1, ..., n},j 6= i.Como como construir tal polinomio? Uma vez que ele se anula nos pontosx0, . . . , xi−1, xi+1, . . . , xn, e facil concluir que ele tem a forma:
Li(x) = Ai(x−x0) . . . (x−xi−1)(x−xi+1) . . . (x−xn), i = 0, ..., n. (4.5)
onde Ai e uma certa constante real (nao depende de x). Para definir o valordesta constante, basta ter em conta a condicao Li(xi) = 1. De acordo com(4.5), temos
Li(xi) = Ai(xi−x0) . . . (xi−xi−1)(xi−xi+1) . . . (xi−xn) = 1, i = 0, ..., n.(4.6)
Daqui imediatamente se conclui
Ai =1
(xi − x0) . . . (xi − xi−1)(xi − xi+1) . . . (xi − xn), i = 0, ..., n.
(4.7)Introduzindo (4.7) em (4.6) e escrevendo tudo numa forma mais compacta,obtem-se finalmente
Li(x) =(x− x0) . . . (x− xi−1)(x− xi+1) . . . (x− xn)
(xi − x0) . . . (xi − xi−1)(xi − xi+1) . . . (xi − xn)=
∏nj=0,j 6=i(x− xj)∏nj=0,j 6=i(xi − xj)
,
(4.8)i = 0, ..., n. E obvio que cada uma das funcoes Li e um polinomio de graun. Para provar que eles formam uma base no espaco dos polinomios degrau menor ou igual a n, basta verificar que sao linearmente independentes.
111
Considere-se uma combinacao linear com a forma
n∑j=0
CjLj(x) (4.9)
onde Cj sao coeficientes reais arbitrarios. Queremos provar que, quaisquerque sejam estes coeficientes (nao simultaneamente nulos), a soma (4.9) naose transforma numa funcao nula. Admitamos o contrario. Isto e, suponha-mos que, para um certo conjunto de valores de Cj , se tem
n∑j=0
CjLj(x) ≡ 0, ∀x ∈ IR. (4.10)
Nesse caso teremos, em particular
n∑j=0
CjLj(xi) = 0, i = 0, 1, ..., n. (4.11)
Por outro lado, temos, por definicao de Li,
n∑j=0
CjLj(xi) = Ci, i = 0, 1, ..., n. (4.12)
Comparando (4.11) com (4.12), conclui-se que Ci = 0,i = 0, 1, .., n, istoe, a identidade (4.10) so se verifica se todos os coeficentes Ci se anularemsimultaneamente. Logo, os polinomios de Lagrange sao linearmente inde-pendentes, e formam, portanto uma base em Pn. Sendo assim, dada umatabela de valores de uma certa funcao f nos pontos xi, o polinomio interpo-lador de f nesses pontos, Pn, tambem pode ser representado (de uma unicamaneira) na forma
Pn(x) =
n∑j=0
djLj(x). (4.13)
Resta-nos determinar o valor dos coeficientes dj nesta formula, o que e bas-tante facil, dada a definicao dos polinomios de Lagrange. Com efeito e facilverificar que, para que o polinomio Pn, dado por (4.13), seja o polinomiointerpolador de f , basta escolher dj = f(xj),j = 0, ..., n, isto e, tem-se
Pn(x) =
n∑j=0
f(xj)Lj(x). (4.14)
112
Para provar a formula (4.14), basta recordar a definicao dos polinomios deLagrange. De facto, calculando Pn em xi, pela formula (4.14), obtem-se
Pn(xi) =
n∑j=0
f(xj)Lj(xi) = f(xi)Li(xi) = f(xi), i = 0, 1, ..., n. (4.15)
A igualdade (4.15), verificada em todos os nos xi, comprova que o polinomioPn, definido por (4.14), e o polinomio interpolador de f nestes nos.
A formula (4.14) e conhecida como a formula interpoladora de Lagrange,sendo os polinomios de Lagrange definidos por (4.8).
4.1.3 Escolha dos nos de interpolacao
Por vezes, ao resolver um problema, dispomos de mais informacao do que aestritamente necessaria para aplicar um determinado metodo. Por exemplo,se quisermos aproximar uma funcao por um polinomio de grau 2 e conhecer-mos os seus valores em quatro pontos, e obvio que teremos de descartar umdos pontos.
Sao possıveis varios criterios para a seleccao de pontos. Em primeirolugar, se a funcao considerada apresentar uma ou mais descontinuidades, fazsentido aproxima-la por trocos. Isto e, se por exemplo ela for descontınuaem x = 0, a interpolacao deve ser feita separadamente para valores de xpositivos e negativos. Nao faria sentido aproximar uma tal funcao usandodois nos de sinais opostos.
Se nao for este o caso, isto e, se a funcao a interpolar for contınua em todoo domınio considerado, entao o criterio mais comum para a escolha dos nosde interpolacao e a proximidade. Isto e, se quisermos aproximar a funcaonum certo ponto x, devem escolher-se primeiro os dois pontos mais proximosde x, sendo os pontos seguintes escolhidos pelo mesmo criterio. Embora oerro de interpolacao, como veremos mais adiante, dependa de varios factores,na ausencia de outras informacoes sobre a funcao, esta e a melhor escolhapossıvel para o minimizar. No final do capıtulo da interpolacao, refrir-nos-emos a outros possıveis criterios de escolha dos pontos, relacionados com aminimizacao do erro de interpolacao.
Exemplo 4.1. Consideremos a funcao f , dada pela tabela 4.1. O nossoobjectivo e obter valores aproximados de N(15) (valor da populacao noinstante t = 15), por interpolacao polinomial, aplicando a formula interpo-radora de Lagrange.
1. Obter um valor aproximado de N(15), por interpolacao linear.
113
Resolucao. Para aplicar interpolacao linear (isto e, por um polinomiode grau nao superior a 1) devemos considerar os valores de N em doispontos. De acordo com o que ficou dito, neste caso os pontos deveraoser os nos mais proximos de x = 15: x0 = 14 e x1 = 16. Note-se quea ordem dos pontos escolhidos e arbitraria, nao influindo no resultadoda interpolacao. Seja entao P1 o polinomio que interpola N em x0
e x1. Para o calcularmos, comecamos por construir os respectivospolinomios de Lagrange. De acordo com a formula (4.8), temos
L0(x) =x− x1
x0 − x1=x− 16
−2, L1(x) =
x− x0
x1 − x0=x− 14
2. (4.16)
Entao, de acordo com a formula interpoladora de Lagrange, temos
P1(x) = f(x0)L0(x) + f(x1)L1(x) = 22x− 16
−2+ 18
x− 14
2. (4.17)
Logo, a aproximacao procurada e P1(15) = 11 + 9 = 20.
2. Obter uma nova aproximacao de N(15), por interpolacao quadratica.
Resolucao. Neste caso, sao necessarios 3 nos de interpolacao. Usandode novo o criterio de proximidade, o terceiro ponto a considerar ex2 = 12. Os polinomios de Lagange correspondentes sao
L0(x) = (x−x1)(x−x2)(x0−x1)(x0−x2) = (x−16)(x−12)
(−2)(2) ,
L1(x) = (x−x0)(x−x2)(x1−x0)(x1−x2) = (x−14)(x−12)
(2)(4) ,
L2(x) = (x−x0)(x−x1)(x2−x0)(x2−x1) = (x−14)(x−16)
(−2)(−4) .
(4.18)
Aplicando de novo a formula interpoladora de Lagrange, temos
P2(x) = f(x0)L0(x) + f(x1)L1(x) + f(x2)L2(x) =
22 (x−16)(x−12)−4 + 18 (x−14)(x−12)
8 + 15 (x−14)(x−16)8 .
(4.19)
Logo, a aproximacao procurada e P2(15) = 22 ∗ 3/4 + 18 ∗ 3/8 + 15 ∗(−1/8) = 21.375.
3. Obter uma nova aproximacao de N(15), usando todos os pontos databela.
Resolucao. Se usarmos todos os pontos da tabela, estaremos a fazerinterpolacao cubica (de grau 3). Uma vez que, como ja foi dito, aordem dos pontos e irrelevante para o resultado, designemos x3 = 10,
114
mantendo as designacoes dos restantes pontos de interpolacao. Ospolinomios de Lagange correspondentes, neste caso, sao
L0(x) = (x−x1)(x−x2)(x−x3)(x0−x1)(x0−x2)(x0−x3) = (x−16)(x−12)(x−10)
(−2)(2)(4) ,
L1(x) = (x−x0)(x−x2)(x−x3)(x1−x0)(x1−x2)(x1−x3) = (x−14)(x−12)(x−10)
(2)(4)(6) ,
L2(x) = (x−x0)(x−x1)(x−x3)(x2−x0)(x2−x1)(x2−x3) = (x−14)(x−16)(x−10)
(−2)(−4)2 ,
L3(x) = (x−x0)(x−x1)(x−x2)(x3−x0)(x3−x1)(x3−x2) = (x−14)(x−16)(x−12)
(−2)(−6)(−4) .
(4.20)
Neste caso, aplicando a formula interpoladora de Lagrange, temos
P3(x) = f(x0)L0(x) + f(x1)L1(x) + f(x2)L2(x) + f(x3)L3(x) =
22 (x−16)(x−12)(x−10)−16 + 18 (x−14)(x−12)(x−10)
48 + 15 (x−14)(x−16)(x−10)16 + 10 (x−14)(x−16)(x−12)
−48 .
(4.21)Logo, a aproximacao por interpolacao cubica e P3(15) = 22.1875.
4.1.4 Formula Interpoladora de Newton
No exemplo anterior a formula interpoladora de Lagrange foi utilizada paracalcular uma sucessao de aproximacoes do mesmo valor. Deste calculoressalta uma desvantagem da formula de Lagrange: cada vez que se de-termina uma nova aproximacao, e necessario calcula-la do princıpio, semaproveitar nada dos calculos efectuados anteriormente. Vamos agora es-tudar uma formula de calculo alternativa, em que o polinomio interpo-lador de um certo grau e obtido como uma correccao do polinomio do grauanterior. Numa sucessao de calculos como a que acabamos de efectuar,este metodo, conhecido como formula interpoladora de Newton, e ni-tidamente mais eficiente, diminuindo substancialmente o numero total deoperacoes aritmeticas necessarias.
Para estudar a formula interpoladora de Newton, comecemos por formu-lar o seguinte problema. Seja Pn o polinomio de grau menor ou a igual an que interpola uma certa funcao f nos nos x0,x1,..,xn. Se a tabela consi-derada de f acrescentarmos mais um no, xn+1, o resultado da interpolacaopassara a ser o polinomio Pn+1, que interpola f tambem neste ponto. Vamosconstruir Pn+1 a partir de Pn. Comecemos por escrever
Pn+1(x) = Pn(x) + Cn+1(x). (4.22)
Nesta formula, assumindo que Pn+1 e diferente de Pn, Cn+1 e um polinomiode grau n+ 1 (o mesmo grau de Pn+1). Facilmente se verifica que as raizes
115
deste polinomio coincidem com os nos de interpolacao x0,x1,...,xn. Comefeito, da igualdade (4.22) resulta imediatamente
Cn+1(xi) = Pn+1(xi)−Pn(xi) = f(xi)−f(xi) = 0, i = 0, 1, . . . n. (4.23)
Por conseguinte, Cn+1 pode ser escrito na forma
Cn+1(x) = An(x− x0)(x− x1) . . . (x− xn), (4.24)
onde An nao depende de x. Podemos entao rescrever a formula (4.22) como
Pn+1(x) = Pn(x) +An+1(x− x0)(x− x1) . . . (x− xn). (4.25)
O problema de calcular Pn+1 ficou pois reduzido a determinar An+1, umaconstante que depende dos valores de f em x0, x1,...,xn. Note-se que, nocaso de Pn+1(x) coincidir com Pn(x) (o que acontece se tivermos Pn(xn+1) =f(xn+1)), resulta que An+1 = 0. Se excluirmos este caso, Pn+1 e umpolinomio de grau n+ 1, que se pode escrever na forma
Pn+1(x) = An+1xn+1 + . . . , (4.26)
ou seja, An+1 e o coeficiente do termo em xn+1 (termo principal) de Pn+1.Podemos assim introduzir a seguinte definicao:
Definicao 4.2 Chama-se diferenca dividida de ordem k da funcao f nosnos x0, x1, . . . , xk o coeficiente Ak = f [x0, x1, . . . , xk] do termo em xk dopolinomio Pk que interpola f nos nos considerados.
Vamos deduzir agora uma formula para calcular as diferencas divididas.Para facilitar, comecemos por considerar as diferencas divididas de primeiraordem. Seja P0 o polinomio que interpola f em x0, P0(x) ≡ f(x0). Sendo x1
um novo ponto de interpolacao, de acordo com a formula (4.25), o polinomioP1, que interpola f em x0 e x1, e dado por
P1(x) = P0(x) +A1(x− x0) = f(x0) +A1(x− x0). (4.27)
Neste caso, o valor de A1 = f [x0, x1] (diferenca dividida de f em x0 e x1)deduz-se facilmente a partir da condicao P1(x1) = f(x1). De acordo com(4.27), obtem-se
P1(x1) = f(x0) +A1(x1 − x0) = f(x1). (4.28)
Assim, resulta que
A1 = f [x0, x1] =f(x1)− f(x0)
x1 − x0. (4.29)
116
No caso de f(x0) = f(x1), temos A1 = 0. Este e o unico caso em que opolinomio P1 coincide com P0, ou seja, tem grau 0. Generalizando a formula(4.29) para quaisquer dois nos de interpolacao xi, xj podemos escrever aformula da diferenca dividida de primeira ordem:
f [xi, xj ] =f(xj)− f(xi)
xj − xi. (4.30)
Esta formula justifica a razao da designacao diferenca dividida. A diferencadividida de primeira ordem tem um significado geometrico simples: e odeclive da recta que passa pelos pontos (x0, f(x0)) e (x1, f(x0)) (e que etambem o grafico do polinomio P1).
Para construir as diferencas divididas de ordem superior (correspon-dentes aos polinomios interpoladores de grau maior que um) torna-se necessariodeduzir uma formula de recorrencia. Suponhamos que e conhecido o polinomioPn, que interpola f em x0, x1, ..., xn. Isto significa que e conhecida a diferencadividida f [x0, x1, ..., xn], ja que ela e o coeficiente do termo principal dePn. Para construir o polinomio Pn+1 precisamos de introduzir mais umno de interpolacao xn+1. Em geral, temos Pn+1(xn+1) = f(xn+1), masPn(xn+1) 6= f(xn+1), ja que Pn+1 interpola f neste ponto (ao contrariode Pn). 3 Vamos definir o polinomio auxiliar Qn do seguinte modo: Qninterpola f em x1, x2, . . . , xn+1. Logo, Qn e um polinomio de grau nao su-perior a n, tal como Pn, mas cujo termo termo principal tem o coeficientef [x1, x2, ..., xn+1]. Vamos mostrar que Pn+1 pode ser obtido a partir de Pne Qn, atraves da formula
Pn+1(x) =Pn(x)(xn+1 − x) +Qn(x)(x− x0)
xn+1 − x0. (4.31)
Para tanto, basta provar que Pn+1(xi) = f(xi), i = 0, 1, ..., n+ 1. Se i = 0 ,temos
Pn+1(x0) =f(x0)(xn+1 − x0)
xn+1 − x0= f(x0). (4.32)
Por outro lado, se i ∈ {1, 2, ..., n}, verifica-se
Pn+1(xi) =f(xi)(xn+1 − xi) + f(xi)(xi − x0)
xn+1 − x0=f(xi)(xn+1 − x0)
xn+1 − x0= f(xi).
(4.33)
3Pode dar-se o caso de Pn(xn+1) = f(xn+1). Nesse caso, Pn+1 coincide com Pn e adiferenca dividida f [x0, x1, ..., xn+1] e nula.
117
Finalmente, para i = n+ 1, obtem-se
Pn+1(xn+1) =f(xn+1)(xn+1 − x0)
xn+1 − x0= f(xn+1). (4.34)
Com isto, acabamos de provar que Pn+1, definido pela formula (4.31), eo polinomio que interpola f nos pontos x0, x1, ..., xn+1. Por definicao, adiferenca dividida f [x0, x1, ..., xn+1] e o coeficiente do termo principal destepolinomio. Assim, ela pode ser calculada atraves da formula
f [x0, x1, ..., xn+1] =−f [x0, x1, ..., xn] + f [x1, x2, ..., xn+1]
xn+1 − x0, (4.35)
onde, como ja sabemos, f [x0, x1, ..., xn] e o coeficiente do termo principal dePn e f [x1, ..., xn+1] e o coeficiente do termo pricipal de Qn . Assim, a formula(4.35) permite-nos calcular uma diferenca dividida de ordem n+ 1 a partirde duas diferencas divididas de ordem n. Aplicando sucessivamente estaformula de recorrencia, podemos calcular diferencas divididas de qualquerordem (desde que, evidentemente, se diponha de valores suficentes da funcaof).
Recapitulando, para construir o polinomio interpolador Pn pela formulade Newton, num certo conjunto de nos, devemos proceder do seguinte modo:
1. Calcular as diferencas divididas de f nos pontos considerados, ate aordem n, com base na formula (4.35);
2. Determinar P0(x) ≡ f(x0);
3. Obter os polinomios P1, P2,..., Pn, atraves da aplicacao sucessiva daformula (4.25), onde An+1 representa uma diferenca dividida da ordemcorrespondente.
Este processo pode ser facilmente programado e realizado num computador.Quando os calculos sao efectuados manualmente, para facilitar, e costumerepresentar as diferencas divididas numa tabela. Apresentamos a seguir umexemplo.
Exemplo 4.2 Construir a tabela de diferencas divididas correspondentea tabela 4.1, ordenando os nos de interpolacao segundo a sua proximidadeao ponto x = 4 (a semelhanca do exemplo 4.1).
Resolucao. Uma tabela de diferencas divididas tem a forma de umtriangulo e a sua construcao comeca sempre pela base, que neste caso, con-siste em duas colunas que contem os dados do problema: uma com os valoresde xi, e outra, com os de f(xi) = fi. Na coluna seguinte, sao calculadas
118
as diferencas divididas de primeira ordem. No caso concreto da tabela 4.1,temos 4 pontos, logo podemos calcular tres destas diferencas:
f [x0, x1] = f1−f0x1−x0 = 18−22
16−14 = −2.
f [x1, x2] = f2−f1x2−x1 = 15−18
12−16 = 34 .
f [x2, x3] = f3−f2x3−x2 = 10−15
10−12 = 52 .
Segue-se a coluna com as diferencas de segunda ordem:
f [x0, x1, x2] = f [x1,x2]−f [x0,x1]x2−x0 = −11
8 .
f [x1, x2, x3] = f [x2,x3]−f [x1,x2]x3−x1 = − 7
24 .
Finalmente, o vertice do triangulo e constituıdo pela diferenca dividida deterceira ordem:
f [x0, x1, x2, x3] =f [x1, x2, x3]− f [x0, x1, x2]
x3 − x0= −13
48. (4.36)
O aspecto grafico da tabela e o seguinte
xi fi ordem1 ordem2 ordem3
14 22−2
16 18 −118
34 −13
4812 15 − 7
2452
10 10
A localizacao de cada entrada da tabela sugere-nos a maneira como ela secalcula: o numerador da fraccao da formula (4.35) e a diferenca entre asduas entradas adjacentes da coluna anterior; o denominador dessa fraccao ea diferenca entre os extremos da base do triangulo cujo vertice se encontrana entrada a calcular.
Exemplo 4.3 Obter aproximacaoes de N(15) por interpolacao linear,quadratica e cubica, recorrendo a formula interpoladora de Newton.
Resolucao. A tabela de diferencas divididas para este problema ja foicalculada no exemplo anterior. Logo, para calcular as aproximacoes pedidas,basta utilizar sucessivamente a formula (4.25).
Para comecar P0(x) ≡ f(x0) = 22. Logo, aplicando a formula (4.25)com n = 0 obtem-se o polinomio interpolador de primeiro grau:
P1(x) = P0(x) + f [x0, x1](x− x0) = 22− 2(x− 14). (4.37)
119
ano 1991 1992 1993 1994 2004 2010
P 6.5 220 320 415 792.5 996.85S 200 222.5 237 246.5 374 475
Table 3: P - valor medio das propinas das licenciaturas (em euros), S- salariomınimo nacional (em euros)
Utilizando este polinomio, obtem-se a aproximacao por interpolacao linear:
N(15) ≈ P1(15) = 22− 2(15− 14) = 20. (4.38)
Aplicando agora a formula (4.25) com n = 1 obtem-se o polinomio interpo-lador de segundo grau:
P2(x) = P1(x)+f [x0, x1, x2](x−x0)(x−x1) = 22−2(x−14)−11/8(x−14)(x−16).(4.39)
Usando este polinomio, obtem-se agora a aproximacao por interpolacaoquadratica:
N(15) ≈ P2(15) = 20− 11/8(15− 14)(15− 16) = 171/8 = 21.375. (4.40)
Finalmente, para obter o polinomio interpolador de grau 3, aplica-se amesma formula com n = 2:
P3(x) = P2(x) + f [x0, x1, x2, x3](x− x0)(x− x1)(x− x2) =22− 2(x− 14)− 11/8(x− 14)(x− 16)− 13/48(x− x0)(x− x1)(x− x2).
(4.41)Obtem-se assim a aproximacao por interpolacao cubica:
N(15) ≈ P3(15) = 21.375− 13/48(15− 14)(15− 16)(15− 12) = 22.1875.(4.42)
Exemplo 4.4. Na tabela 3, e dada a evolucao do valor medio daspropinas das Licenciaturas em Portugal, no perıodo entre 1991 e 2010, bemcomo da evolucao do salario mınimo nacional no mesmo perıodo. Parasimplificar, quando nos referimos as propinas no ano lectivo N/N + 1, con-sideramos que as mesmas correspondem ao ano N .
Vamos ilustrar a aplicacao da interpolacao polinomial com alguns calculosbaseados nesta tabela.
1. Em primeiro lugar, vamos responder a seguinte questao: sera que aevolucao do valor das propinas entre 1991 e 1994 permitia prever ovalor que elas iriam atingir em 2004 e em 2010? Responderemos aesta pergunta por etapas.
120
(a) Utilizando interpolacao quadratica e aplicando a formula de New-ton, obter uma estimativa do valor das propinas em 2004, ou sejaP (2004) (baseado apenas nos valores das propinas em anos ante-riores).
Resolucao. Para realizar interpolacao quadratica, devemos basear-nos nos valores da funcao P em 3 anos anteriores a 2004; comodevemos escolher os 3 anos mais proximos de 2004, temos 1992,1993 e 1994. Calculemos as diferencas divididas:
P [1992, 1993] = (P (1993)− P (1992))/(1993− 1992) = 100,
P [1993, 1994] = (P (1994)− P (1993))/(1994− 1993) = 95,
P [1992, 1993, 1994] = (P [1993, 1994]−P [1992, 1993])/(1994−1992) = −2.5.
Aplicando a formula interpoladora de Newton, tem-se:
P2(t) = P (1992)+P [1992, 1993](t−1992)+P [1992, 1993, 1994](t−1992)(t−1993).
Finalmente, substituindo t pelo ano em causa, 2004, obtem-se
P2(2004) = 1090.
(b) Aplicando de novo a formula de Newton,obtenha uma estimativade P (2010), tambem por interpolacao quadratica.
Resolucao. Para realizar interpolacao quadratica, neste caso,devemos basear-nos nos valores da funcao P nos 3 anos maisrecentes, anteriores a 2010: 1993, 1994 e 2004. Calculemos asdiferencas divididas:
P [1994, 2004] = (P (2004)− P (1994))/(2004− 1994) = 37.76;
P [1993, 1994, 2004] = (P [1994, 2004]−P [1993, 1994])/(2004−1993) = −5.20364.
Aplicando a formula interpoladora de Newton, tem-se:
Q2(t) = P (1993)+P [1993, 1994](t−1993)+P [1993, 1994, 2004](t−1993)(t−1994).
Finalmente, substituindo t pelo ano em causa, 2010, obtem-se
Q2(2010) = 519.611.
121
(c) Compare as estimativas obtidas com os valores reais dados natabela.
Resolucao. A primeira estimativa, comparada com o valor deP (2004) = 792.5, tem um erro por excesso de 297.5, enquantoque a segunda, comparada o valor de P (2010) = 996.85, temum erro por defeito de 273.5. Erros tao significativos nao sao deestranhar, ja que neste caso estamos a fazer uma extrapolacao,isto e, estamos a basear-nos em valores da funcao P em certosintervalos para obter estimativas do valor dessa funcao em pontosexteriores a esses intervalos.
2. Utilizando interpolacao cubica, e com base na mesma tabela, que pre-visao se pode fazer para o valor das propinas em 2012?
Resolucao. Para realizar interpolacao cubica, neste caso, devemosbasear-nos nos valores da funcao P nos 4 anos tabelados anterioresa 2012: 1993, 1994 , 2004 e 2010. Vamos aproveitar as diferencasdivididas que ja calculamos: as que nao envolvem 2010. Alem disso,preciso de calcular mais tres diferencas:
P [2004, 2010] = (P (2010)− P (2004))/(2010− 2004) = 34.0417;
P [1994, 2004, 2010] = (P [2004, 2010]−P [1994, 2004])/(2010−1994) = −0.232396.
P[1993,1994,2004,2010]=(P[1994,2004,2010]-P[1993,1994,2004])/(2010-1993)=0.292426. Aplicando a formula interpoladora de Newton, tem-se:
P3(t) = Q2(t)+P [1993, 1994, 2004, 2010](t−1993)(t−1994)(t−2004).
Finalmente, substituindo t pelo ano em causa, 2012, obtem-se
P3(2012) = 1145.43.
Nas restantes perguntas desta serie sugere-se que use a formula inter-poladora de Lagrange.
3. Atendendo a que a razao P/S e um importante ındice sobre o grau deacesso ao ensino superior, analise a evolucao desta razao a partir databela dada e, por interpolacao quadratica, obtenha uma previsao dovalor que esse ındice atingira em 2012.
122
Resolucao. Comecemos por calcular as razoes P/S nos anos consid-erados:
P (1991)/S(1991) = 0.0325, P (1992)/S(1992) = 0.988764,
P (1993)/S(1993) = 1.35021, P (1994)/S(1994) = 1.68357,
P (2004)/S(2004) = 2.11925, P (2010)/S(2010) = 2.09863.
Ou seja, durante o perıodo considerada o valor das propinas passou decerca de 3 por cento para mais do dobro do valor do salario mınimo.Para fazer a interpolacao pedida, tenhamos em conta os tres ultimosvalores da tabela, correspondentes a t0 = 1994, t1 = 2004 e t2 = 2010.Os polinomios de Lagrange para estes pontos sao:
l0(t) =(t− t1)(t− t2)
(t0 − t1)(t0 − t2))=
(t− 2004)(t− 2010)
(1994− 2004)(1994− 2010),
l1(t) =(t− t0)(t− t2)
(t1 − t0)(t1 − t2))=
(t− 1994)(t− 2010)
(2004− 1994)(2004− 2010),
l2(t) =(t− t0)(t− t1)
(t2 − t0)(t2 − t1)=
(t− 1994)(t− 2004)
(2010− 1994)(2010− 20004).
Para construir o polinomio interpolador, basta entao usar a formula:
Π2(t) = P (t0)/S(t0)l0(t) + P (t1)/S(t1)l1(t) + P (t2)/S(t2)l2(t) =
= 1.68357l0(t) + 2.11925l1(t) + 2.09863l2(t).
Finalmente, para responder a pergunta, basta calcular
Π2(2012) = 2.04475.
4. Com base na evolucao do salario mınimo nacional, obtenha, por in-terpolacao quadratica, uma previsao do salario mınimo nacional em2012.
Resolucao. Aproveitando os polinomios de Lagrange ja calculados,podemos escrever a seguinte formula para o polinomio S2, que inter-pola a funcao S(t) em 1994, 2004 e 2010:
S2(t) = S(1994)l0(t) + S(2004)l1(t) + S(2010)l2(t).
Substituındo entao t por 2012, obtem-se a previsao pedida:
S2(2012) = 512.75.
123
5. Com base nas respostas as duas ultimas alıneas, obtenha uma novaprevisao do valor das propinas em 2012 e compare com a previsaoobtida na alınea 2.
Resolucao. Neste caso, baseamo-nos na previsao do valor do salarioS2(2012) e na da relacao propinas/salario para o mesmo ano Π2(2012) =2.04475. Nesse caso, obtem-se
P (2012) = S2(2012)Π2(2012) = 1048.45.
6. Admitindo agora que o valor das propinas e uma funcao do salariomınimo nacional, obtenha, por interpolacao quadratica e cubica, esti-mativas do valor das propinas quando o salario mınimo nacional atingir500 euros.
Resolucao. Neste caso, devemos considerar a seguinte tabela:P (237) = 320,P (246.5) = 420,P (374) = 792.5 e P (475) = 976.85.Para obter as estimativas pedidas , e necessario interpolar essa tabela.No caso da interpolacao quadratica, utilizam-se os 3 ultimos valoresde P (aqueles cuja abcissa e mais proxima de 500). Representandopor P2 o polinomio interpolador correspondente (que se pode obterpor qualquer uma das formulas ja utilizadas) obtem-se
P2(2012) = 1034.46.
No caso da interpolacao cubica, utilizam-se todos os pontos da tabelaacima representada. Neste caso obtem-se
P3(2012) = 1193.04.
4.1.5 Erro de Interpolacao
Nesta seccao vamos debrucar-nos sobre o erro de interpolacao. Chama-seerro de interpolacao a diferenca
en(x) = f(x)− Pn(x), (4.43)
onde Pn e o polinomio que interpola f em n+1 pontos (x0, x1, . . . , xn.) Esteestudo permite-nos, entre outras coisas, decidir qual o grau do polinomiointerpolador que melhor aproxima a funcao considerada num certo ponto.
Para comecar, assumindo que se pretende aproximar f num certo in-tervalo [a, b] (ao qual pertencem os nos de interpolacao), seja x um pontogenerico deste intervalo. Naturalmente, se x coincidir com um dos nos de
124
interpolacao xi teremos en(x) = en(xi) = f(xi) − Pn(xi) = 0. Suponhamosagora que x nao e nenhum dos nos. Para avaliar en(x), neste caso, vamosconstruir o polinomio Pn+1 que interpola f em x0, x1, . . . , xn, x. De acordocom a formula inerpoladora de Newton, temos
Pn+1(x) = Pn(x) + f [x0, x1, ..., xn, x]n∏i=0
(x− xi). (4.44)
Em particular,
Pn+1(x) = Pn(x) + f [x0, x1, ..., xn, x]n∏i=0
(x− xi). (4.45)
Visto que, por construcao, Pn+1(x) = f(x), temos en(x) = Pn+1(x)−Pn(x).Logo, de (4.45) resulta que
en(x) = Pn+1(x)− Pn(x) = f [x0, x1, ..., xn, x]
n∏i=0
(x− xi). (4.46)
Visto que x e um ponto generico de [a, b], a formula (4.46) pode ser utilizadapara estimar o erro de interpolacao em qualquer ponto deste intervalo.
Contudo, esta formula nao e facilmente aplicavel, ja que a estimativa doerro que ela proporciona depende depende de f [x0, x1, ..., xn, x], grandezaque geralmente nao e conhecida (alias, em geral, nem a funcao f e suposta-mente conhecida no ponto x). Assim, para que a formula (4.46) possa teralguma utilidade pratica, e necessario relacionar as diferencas divididas deuma funcao f com as suas derivadas (assumindo que estas sao conhecidas).
No caso de n = 1 existe uma relacao simples entre as diferencas divididasde uma funcao e a sua primeira derivada. De facto, se f for uma funcao con-tinuamente diferenciavel em [x0, x1], de acordo com o teorema de Lagrange,verifica-se
f [x0, x1] =f(x1)− f(x0)
x1 − x0= f ′(ξ), (4.47)
onde ξ ∈ [x0, x1]. Surge naturalmente a pergunta: sera que existe umarelacao semelhante entre as diferncas divididas de uma certa ordem k e aderivada de f da mesma ordem? A resposta a esta pergunta e positiva e edada pelo seguinte teorema, que constitui uma generalizacao do teorema deLagrange.
Teorema. Seja f ∈ Ck([a, b]) uma funcao dada e sejam x0, x1, ..., xk umconjunto de k + 1 pontos distintos de [a, b](k ≥ 1). Entao existe ξ ∈ [a, b]
125
tal que
f [x0, x1, . . . , xk] =f (k)(ξ)
k!. (4.48)
Demonstracao. Sejaek(x) = f(x)− Pk(x) (4.49)
o erro de interpolacao de f por Pk, onde Pk representa o polinomio interpo-lador de f em x0, x1, . . . , xk. Por definicao, temos
ek(xi) = 0, i = 0, 1, ..., k, (4.50)
ou seja, ek tem, pelo menos, k + 1 zeros distintos em [a, b]. Alem disso,ek tem k derivadas contınuas em [a, b], segundo resulta das condicoes do
teorema. Aplicando k vezes o teorema de Rolle, conclui-se que e(k)k se anula,
pelo menos, uma vez em [a, b]. Logo, existe ξ ∈ [a, b], tal que e(k)k (ξ) = 0.
Mostremos que para este ponto ξ e valida a igualdade (4.48). Com efeito,pela definicao de diferenca dividida de ordem k, temos
e(k)k (ξ) = f (k)(ξ)− P (k)
k (ξ) = f (k)(ξ)− k!f [x0, . . . , xk] = 0. (4.51)
De (4.51) resulta que
f [x0, x1, . . . , xk] =f (k)(ξ)
k!, (4.52)
como se pretendia demonstrar.Da igualdade (4.52) resulta imediatamente que
|f [x0, x1, . . . , xk]| ≤1
k!maxx∈[a,b]
|f (k)(x)|. (4.53)
Combinando esta formula com (4.46) obtem-se que
|en(x)| = |f [x0, x1, ..., xn, x]|n∏i=0
|x−xi| ≤1
(n+ 1)!maxx∈[a,b]
|f (n+1)(x)|n∏i=0
|x−xi|.
(4.54)Exemplo 4.5. Suponhamos que se aproxima a funcao f(x) = cos(x) no
intervalo [−1, 1] por um polinomio interpolador nos nos x0 = −1,x1 = 0, x2 = 1.
1. Determine o polinomio interpolador P2(x).
126
A fim de aplicar a formula interpoladora de Newton, comecemos porcalcular as diferencas divididas de f :
f [x0, x1] =cos(x1)− cos(x0)
x1 − x0= 1− cos(−1);
f [x1, x2] =cos(x2)− cos(x1)
x2 − x1= cos(1)− 1.
f [x0, x1, x2] =f [x1, x2]− f [x1, x2]
x2 − x0=
2 cos(1)− 2
2= cos(1)− 1.
De acordo com a formula (4.25), o poli´nomio interpolador e dado por:
P2(x) = f(−1) + f [−1, 0](x+ 1) + f [−1, 0, 1](x+ 1)x =cos(−1) + (1− cos(−1))(x+ 1) + (cos(1)− 1)(x+ 1)x.
(4.55)
2. Determine um majorante de e2(x) sendo x ∈ [−1, 1].
De acordo com (4.54), temos
|e2(x)| ≤ 1
3!max
x∈[−1,1]|f (3)(x)||x+ 1||x||x− 1|.
Alem disso,
maxx∈[−1,1]
|f (3)(x)| = maxx∈[−1,1]
| sin(x)| = sin(1).
Por conseguinte,
|e2(x)| ≤ sin(1)
3!|x+ 1||x||x− 1|. (4.56)
3. Determine um majorante do maximo de erro de interpolacao no inter-valo [−1, 1].
Pretende-se majorar maxx∈[−1,1] |e2(x)|. Para isso, baseando-nos naresposta anterior, basta obter
maxx∈[−1,1]
|ω(x)|,
ondeω(x) = x(x− 1)(x+ 1) = x3 − x.
127
-1.0 -0.5 0.5 1.0
0.01
0.02
0.03
0.04
0.05
Figure 14: Grafico relativo ao Exemplo 4.5. O erro de interpolacao absolutode facto cometido esta representado a preto. A cinzento esta representadoo majorante do erro dado pela formula (4.56).
Para determinar os pontos de extremo de ω(x) resolve-se a equacao
ω′(x) = 3x2 − 1 = 0,
a qual tem duas raizes: y1 = − 1√3, y2 = 1√
3. E facil verificar que a
primeira destas raizes corresponde a um maximo local de ω, enquantoa segunda, a um mınimo local. Por outro lado, sendo ω uma funcaoımpar, facilmente se deduz que o mınimo local e o simetrico do maximolocal. Concluındo, temos
maxx∈[−1,1]
|ω(x)| = |ω(y1)| = |ω(y2)| = |y2(y2−1)(y2+1)| = 2
3√
3. (4.57)
Finalmente, combinando (4.56) com (4.57), obtem-se
maxx∈[−1,1]
|e2(x)| ≤ sin(1)
3!max
x∈[−1,1]|ω(x)| = sin(1)
3!
2
3√
3≈ 0.054 (4.58)
Este exemplo e ilustrado pela Fig. 14.
4.2 Metodo dos Mınimos Quadrados
A interpolacao polinomial, estudada na seccao anterior, tem o inconvenientede ser extremamente sensıvel a erros nos dados iniciais. Na realidade, amatriz de Vandermonde, de que falamos na Sec. 4.1.1 torna-se mal condi-cionada com o aumento de n, podendo o seu numero de condicao atingirvalores muito elevados, mesmo para valores de n nao muito altos. Isto
128
significa que um pequeno desvio num dos valores de f pode resultar numpolinomio interpolador totalmente diferente. Esta propriedade e altamenteindesejavel para as aplicacoes, ja que na vida real e frequente trabalharmoscom tabelas de valores que sao afectados de erros (resultantes de medicoesou de arredondamentos). Por outro lado, como foi observado quando es-tudamos o erro de interpolacao, este erro pode, nalgun casos, aumentar como aumento do grau do polinomio (recorde-se o exemplo de Runge). Tudoisto nos alerta para o facto de que a interpolacao polinomial pode nao ser,nalguns casos, uma boa forma de aproximar funcoes, sobretudo quando onumero de dados e elevado.
Nesta seccao, vamos estudar um metodo alternativo para aproximacaode funcoes. o metodo dos mınimos quadrados. Entre as vantagens destemetodo contam-se: 1) permmitir uma grande variedade de funcoes ajusta-doras, sem que a forma da funcao dependa do numero de dados; 2) ser menossensıvel aos erros dos dados (em comparacao com a interpolacao); 3) com oaumento do numero de dados, a qualidade da aproximacao tende a aumen-tar. Por outro lado, a soma dos quadrados dos desvios, sendo o criterio paraa escolha da funcao ajustadora, constitui um ındice para avaliar a qualidadeda aproximacao.
Na proxima seccao iremos descrever o metodo dos mınimos quadrados,com ajustamentos lineares, no caso discreto. Nas seccoes seguintes abor-daremos ajustamentos nao lineares e o caso contınuo.
4.2.1 Ajustamentos lineares no caso discreto
O caso discreto caracteriza-se pela forma como e dada a funcao a aproxi-mar: atraves de uma tabela de pontos (tal como no caso da interpolacaopolinomial). Seja pois fi = f(xi) o valor de f no ponto xi (i = 0, 1, . . . , n).O nosso objectivo e construir uma funcao g de uma certa foma (dita afuncao ajustadora), definida num intervalo que contem esses pontos e queconstiui, num certo sentido, a melhor aproximacao de f entre as funcoes daforma considerada. A funcao ajustadora depende de um certo numero deparametros, que representaremos genericamente por a0,a1,. . . ,am. No casodos ajustamentos lineares, de que trataremos em primeiro lugar, a funcaoajustadora pertence a um espaco linear de dimensao m + 1, podendo serescrita na forma
g(x) =
m∑i=0
aiφi(x), (4.59)
onde φ sao funcoes dadas, chamadas as funcoes de base . Essas funcoesdevem estar definidas em todos os pontos xi e devem, alem disso, ser linear-
129
mente independentes, no seguinte sentido. Se fizermos corresponder a cadafuncao φi um vector φi tal que φi = (φi(x0), . . . , φi(xn)), os vectores φi saolinearmente independentes em IRn+1,i = 0, . . . ,m.
Nos problemas da vida real as funcoes de base sao escolhidas de acordocom as propriedades conhecidas da funcao a aproximar.
Exemplo 4.6. Se a funcao a aproximar for uma funcao linear, entaoas funcoes de base deverao ser φ0(x) = 1,φ1(x) = x, de tal modo que oespaco linear onde se procura a funcao ajustadora e o espaco das funcoesda forma g(x) = a0 + a1x, ou seja, o dos polinomios de grau nao superior a1. Num contexto mais geral, se quisermos usar como funcao ajustadora umpolinomio de grau m, entao as funcoes de base a utilizar serao φi(x) = xi,i = 0, . . . ,m. Note-se que estas funcoes sao linearmente independentes nosentido acima mencionado, quaisquer que sejam os pontos xi,i = 0, 1, . . . , n,ja que os vectores φi neste caso tem a forma φi = (xi0, x
i1, ..., x
in), i = 0, ...,m
que formam sempre um sistema linearmente independente.Exemplo 4.7 No caso da aproximacao de funcoes periodicas e comum
tambem usarem-se bases de funcoes trigonometricas, como por exemploφ0(x) = 1,φi(x) = cos(ix),i = 0, 1, . . . ,m. Com funcoes deste tipo, o sis-tema podera ser ou nao linearmente independente, consoante a escolha dospontos xj e o numero de funcoes de base. Se tivermos, por exemplo,xj =jπ/4, j = 0, .., 4, os vectores φi neste caso tem a forma φ0 = (1, 1, 1, 1, 1),¯phii = (1, cos(i ∗ π/4), cos(2i ∗ π/4), cos(3i ∗ π/4) cos(i ∗ π)), o que garante
que sao linearmente independentes, para m ≤ 4.Uma vez escolhidas as funcoes de base, determinar a funcao ajustadora
corresponde a determinar os coeficientes ai da formula (4.59). Estes coefi-cientes sao determinados com base no criterio dos mınimos quadrados, ouseja, de forma a minimizar a soma
Q(a0, a1, ..., am) =
n∑i=0
(f(xi)− g(xi))2. (4.60)
Visto que Q e uma funcao de m+ 1 variaveis, a solucao deste problemade minimizacao obtem-se resolvendo o sistema
∂Q∂a0
= 0,∂Q∂a1
= 0,
. . .∂Q∂am
= 0.
, (4.61)
dito o sistema normal.
130
Para construir o sistema normal e discutir as suas propriedades, e con-veneiente, antes de mais, introduzir a notacao de produto interno. Para osnossos fins, o produto interno de duas funcoes u, v, definidas nos pontos xida tabela e dado por
(u, v) =n∑i=0
u(xi)v(xi). (4.62)
Nesta notacao a funcao Q pode ser dada pelo produto
Q(a0, a1, ..., am) = (f − g, f − g). (4.63)
Por outro lado, usando as propriedades do produto interno, as derivadasparciais de Q podem ser representadas do seguinte modo:
∂Q
∂ai=∂(f − g, f − g)
∂ai= 2
(∂(f − g)
∂ai, f − g
). (4.64)
Alem disso, utilizando a expressao (4.59) para g, de (4.64) obtem-se
∂Q
∂ai= −2
∂(∑m
j=0 ajφj
)∂ai
, f − g
= −2(φi, f − g). (4.65)
Sendo assim, cada uma das equacoes do sistema (4.61) pode ser escritasimplesmente na forma:
(φi, f − g) = 0, i = 0, . . . ,m, (4.66)
ou seja,(φi, g) = (φi, f), i = 0, . . . ,m. (4.67)
Usando mais uma vez a represententacao (4.59) e a propriedade distributivado produto interno, obtem-se finalmente
m∑j=0
aj(φi, φj) = (φi, f), i = 0, . . . ,m, (4.68)
que constitui a forma compacta do sistema normal. Conluımos assim que osistema normal e um sistema linear de m+ 1 equacoes lineares que pode serescrito na forma
(φ0, φ0) (φ0, φ1) . . . (φ0, φm)(φ1, φ0) (φ1, φ1) . . . (φ1, φm). . . . . . . . . . . .
(φm, φ0) (φm, φ1) . . . (φm, φm)
a0
a1
. . .am
=
(φ0, f)(φ1, f). . .
(φm, f)
. (4.69)
131
Este sistema tem solucao unica se e so se a sua matriz, que representare-mos por S, for invertıvel. Pode provar-se (ver, por exemplo, [2]), que ainversibilidade da matriz S resulta do facto de as funcoes de base formaremum sistema linearmente independente, no sentido acima referido. Sendoassim, nessas condicoes,o sistema normal tem sempre solucao unica.
Alem disso, verifica-se facilmente que S e uma matriz simetrica, dadoque
Sij = (φi, φj) = (φj , φi) = Sji, ∀i, j ∈ {0, ...,m}.
Isto facilita o calculo de S, pois basta calcular as entradas da diagonalprincipal e as que se encontram acima delas.
Vale a pena tambem referir que, nas condicoes dadas, S e uma matrizdefinida positiva.
Em geral, o sistema (4.69) e resolvido numericamente usando, por ex-emplo, um dos metodos estudados no capıtulo 3. O conhecimento das pro-priedades de S e pois util para a escolha do metodo a aplicar. Uma dasescolhas mais frequentes e o metodo de Cholesky, ja que e aplicavel a ma-trizes simetricas e definidas positivas.
O exemplo seguinte ilustra a aplicacao do metodo acima descrito.Exemplo 4.8 Consideremos a funcao f , dada pela seguinte tabela:
xi 1 2 3 4
fi 7 4.5 3 2
Pretende-se aproximar f atraves de uma funcao ajustadora da forma
g(x) = a0 +a1
x.
Da forma da funcao ajustadora resulta que se trata de um ajustamentolinear com duas funcoes de base:
φ0(x) = 1, φ1(x) = x.
Para resolver o problema, os valores de a0 e a1 devem ser obtidos pelaresolucao do sistema normal, que neste caso e um sistema de duas equacoes:[
(φ0, φ0) (φ0, φ1)(φ1, φ0) (φ1, φ1)
] [a0
a1
]=
[(φ0, f)(φ1, f)
]. (4.70)
Calculemos os produtos internos necessarios:
(φ0, φ0) =3∑i=0
φ0(xi)2 = 1 + 1 + 1 + 1 = 4;
132
(φ0, φ1) =3∑i=0
φ0(xi)φ1(xi) =1
1+
1
2+
1
3+
1
4=
25
12;
(φ1, φ1) =3∑i=0
φ1(xi)2 =
1
1+
1
22+
1
32+
1
42=
205
144;
(φ0, f) =3∑i=0
φ0(xi)f(xi) = f(1) + f(2) + f(3) + f(4) = 16.5;
(φ1, f) =
3∑i=0
φ1(xi)f(xi) = f(1) +f(2)
2+f(3)
3+f(4)
4= 10.75.
Substituindo estes valores no sistema (4.70) obtem-se[4 25/12
25/12 205/144
] [a0
a1
]=
[16.510.75
]. (4.71)
A solucao do sistema (4.71) e a0 = 0.8077,a1 = 6.369, pelo que a funcaoajustadora procurada e
g(x) = 0.8077 +6.369
x.
Para concluir este exemplo, vamos calcular o mınimo
min(a0,a1)∈IR2
Q(a0, a1) = min(a0,a1)∈IR2
3∑i=0
(f(xi)−g(xi))2 =
3∑i=0
(f(xi)−a0−a1/xi)2.
De acordo com os calculos ja efectuados, este mınimo e atingido quandoa0 = 0.8077,a1 = 6.369, pelo que basta calcular Q(0.8077, 6.369). O resumodos calculos e apresentado na tabela que se segue.
xi fi g(xi) (fi − g(xi))2
1 7 7.177 0.0312 4.5 3.992 0.2583 3 2.931 0.0054 2 2.400 0.160
O valor procurado e a soma dos valores da ultima coluna da tabela:
Q(0.8077, 6.369) = 0.454.
Conforme resulta da definicao do metodo, quaisquer que sejam a0, a1 ∈ IRverifica-se
Q(a0, a1) ≥ 0.454.
Este Exemplo e ilustrado pela Fig. 15.
133
1.5 2.0 2.5 3.0 3.5 4.0
1
2
3
4
5
6
7
Figure 15: Grafico relativo ao Exemplo 4.8. Os pontos representam osvalores de f e a linha contınua e o grafico da funcao ajustadora.
References
[1] K.E. Atkinson, An Introduction to Numerical Analysis, John Wiley &sons, New York, 1978.
[2] M. Carpentier,Analise Numerica, Seccao de Folhas, IST, 1993.
[3] M. Graca, P. Lima, Matematica Experimental, ISTPress, 2007.
[4] P.Lima, Metodos Numericos da Algebra, disponıvel emwww.math.ist.utl.pt/∼plima/LMAC/mna.pdf .
[5] H. Pina, Metodos Numericos, Mc Graw Hill, 1995.
134
