eBook Calculo Numerico

Embed Size (px)

Citation preview

Leonardo F. Guidi

Notas da disciplina Clculo Numrico A

1 Instituto de Matemtica Universidade Federal do Rio Grande do Sul Av. Bento Gonalves, 9500 Porto Alegre - RS

SumrioCaptulo 1. Representao de nmeros em mquinas . . . . . . . . . . . . . . . . . . . . . . . . . 1.1. 1.2. Sistema de numerao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1. 1.2.1. 1.2.2. 1.2.3. Mudana de base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Representao de nmeros inteiros . . . . . . . . . . . . . . . . . . . . . . . . . . . . Representao de nmeros com parte fracionria . . . . . . . . . . . . . . . . . . . . Representao de ponto utuante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Propriedades do conjunto F (b, n, e1 , e2 ) . . . . . . . . . . . . . . . . . . . . . . . . . Representao em ponto utuante IEEE754 . . . . . . . . . . . . . . . . . . . . . . . 1.3. Erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1. 1.3.2. 1.3.3. 1.3.4. 1.3.5. 2.1. Origem dos erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conceitos iniciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Arredondamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Propagao dos erros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Instabilidade numrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Aritmtica de mquina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 6 8 8 9 10 11 12 13 13 14 14 16 16 20 21 21 23 24 25 27 28 29 29 31 32 32 33 33 35

Captulo 2. Equaes no lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodos de quebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1. 2.1.2. 2.2. 2.2.1. 2.2.2. 2.3. 2.3.1. Mtodo da bisseco . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Limitaes do mtodo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo da falsa posio ou regula falsi . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo da iterao linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo da secante . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodos de ponto xo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Mtodos de mltiplos pontos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Captulo 3. Sistemas de equaes lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Captulo 4. Interpolao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Algoritmo de Horner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Interpolao polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1. 4.1.2. Interpolao pelos polinmios de Lagrange . . . . . . . . . . . . . . . . . . . . . . . Interpolao de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Sumrio 4.1.3. 4.2. Erros de truncamento na interpolao por polinmios . . . . . . . . . . . . . . . . . . Fenmeno de Runge . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interpolao spline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1. Interpolao spline cbica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spline natural . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Spline com as mesmas condies de f nos extremos . . . . . . . . . . . . . . . . . . Captulo 5. Ajuste de mnimos quadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1. 5.2. Ajuste de mnimos quadrados polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Caso particular: regresso linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Ajuste de mnimos quadrados por uma combinao linear de funes . . . . . . . . . . . . . . Problemas de condicionamento no mtodo de ajuste de funes pelo mtodo dos mnimos quadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3. 5.4. Ajuste de mnimos quadrados por uma combinao linear de funes ortogonais . . . . . . . . Ajustes no lineares redutveis ao caso linear . . . . . . . . . . . . . . . . . . . . . . . . . .

338 38 40 41 43 44 46 46 49 50 51 51 53 55 55 56 57 58 62 62 65 65 65 67 68 68 68 68 69 69 70 71 73 76 77 81 82 84 89

Captulo 6. Derivao numrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Aproximao da derivada por diferenas nitas . . . . . . . . . . . . . . . . . . . . . Erros de truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Erros de arredondamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1. 7.1. 7.2. Extrapolao de Richardson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quadratura por interpolao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quadraturas newtonianas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.1. 7.2.2. 7.2.3. Regra do trapzio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Erro de truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Regra de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Regras de ordem superior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Regra 3/8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Regra de Bode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.4. Regras compostas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Erro de truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Regra de Simpson composta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.5. 7.3. 8.1. 8.2. 8.3. Mtodo de Romberg . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Quadratura gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo da srie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Erro de truncamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mtodo Runge-Kutta Clssico (4 estgios) . . . . . . . . . . . . . . . . . . . . . . . . Regra do trapzio composta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Captulo 7. Integrao numrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Captulo 8. Equaes Diferenciais Ordinrias . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Sumrio 8.4. Mtodos de mltiplos passos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.1. Mtodo de Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8.4.2. Mtodo de Adasm-Moulton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

491 91 91

Captulo 1

Representao de nmeros em mquinas1.1. Sistema de numeraoUm sistema de numerao formado por uma coleo de smbolos e regras para representar de maneira consistente conjuntos de nmeros. Para desempenhar favoravelmente o seu propsito, um sistema de numerao deve possuir as seguintes propriedades: 1. Deve ser capaz de representar um conjunto de nmeros distintos de maneira padronizada. 2. Deve reetir as estruturas algbricas e aritmticas dos nmeros O sistema mais utilizado o sistema de numerao posicional de base10 ou sistema de numerao base-10 o sistema decimal. Em um sistema de numerao posicional a posio relativa dos smbolos guarda informao sobre o nmero que se quer representar. Mais especicamente, posies adjacentes esto relacionadas entre si por uma constante multiplicativa denominada base do sistema de numerao. Assim, cada smbolo em uma determinada posio contribui aditivamente com uma quantidade dada pela multiplicao do valor numrico do smbolo pela base elevada ao expoente dado pelo posio. Dessa forma, o numeral, formado por uma seqncia de smbolos, representa o nmero calculado a partir da soma da contribuio conjunta de cada smbolo e de sua posio. Denio (Sistemas de numerao base-b). Sistemas de numerao base-b. Dado um natural b > 0, a coleo de smbolos , , e os algarismos 1 {0, 1, . . . , b 1}, o numeral 2 dn dn1 d1 d0 , d1 representa o nmero real positivo dj bj .jn b

(1.1)

(1.2)

Aplicando o sinal frente do numeral representamos os reais negativos. Observao. Quando no h possibilidade de equvoco, dispensamos o subscrito 10 na representao decimal. Dessa forma, por exemplo, o numeral 1, 23 10 , ou simplesmente 1, 23, representa oQuando a base maior do que 10 utilizamos as letras maisculas A, B, C, . . . para representar os algarismos essa a notao utilizada na base hexadecimal (base-16). 2 Por se tratar de nmeros reais, a seqncia sempre limitada esquerda e portanto n Z.1

Captulo 1. Representao de nmeros em mquinas

6

dos se no houver mais algum outro algarismo esquerda. O mesmo ocorre para os algarismos zero preciso de alguma medida. Assim, evitamos notaes da forma 00 . . . 02, 3 para indicar o racional23 10 ;

nmero 1 100 +2 101 +3 102 . Os algarismos 0 antes da vrgula tambm no so representa-

direita aps a vrgula com exceo do uso nas cincias exatas quando faz-se necessrio registrar a o mesmo para notaes da forma 2, 30 . . . 00, com exceo dos casos em que se quer indicar no23 10

o racional

mas uma medida com determinada preciso. Quando a representao de um nmero

possuir dzimas peridicas utiliza-se como notao uma barra sobre a seqncia que se repete. Por 67 exemplo, 1 representado como 0, 3 e 495 representado como 0, 135. 3 Os algarismos esquerda da vrgula formam a parte inteira do numeral, os demais formam a sua parte fracionria. Os sistemas de numerao base-b. O sistemas de numerao base-b so capazes de representar unicamente os nmeros inteiros, no entanto, o mesmo no verdade para nmeros racionais. Por exemplo, o racional11 10

possui as representaes 1, 1 e 1, 09.

1.1.1. Mudana de base A partir das expresses (1.1) e (1.2) possvel construir o procedimento de mudana de bases. Seja, ento, X um nmero representado na base b como d n dn1 d1 d0 , d1 . . . b . Encontrar a sua representao em uma outra base g signica encontrar a seqncia de algarismos dm dm1 d1 d0 , d1 . . . tal que X=jm

dj g j .

Inicialmente devemos separar as partes inteira e fracionria de um numeral. Denominamos a parte inteira de X como X i e, de acordo com (1.2), na nova base gm

X =j=0

i

dj g j ,

(1.3)

de maneira semelhante, denominamos X f a parte fracionria de x:

Xf =j=1

dj g j .

(1.4)

Ao dividirmos X i por g temos Xi d0 = + g g

m

dj g j1 .j=1

O que nos fornece o primeiro dgito, dividindo o segundo membro da igualdade anterior temosm j1 j=1 dj g

g

=

d1 + g

m

dj g j2 .j=2

Captulo 1. Representao de nmeros em mquinas

7

E assim sucessivamente at encontrarmos o dgito dm o que ocorre quando o quociente da ltima diviso for nulo. Convm lembrar que ao iniciar o processo, conhecemos a representao de X i na base, portanto as operaes de diviso devem ser feitas na base b. Em resumo, iniciamos o procedimento com a representao de X i na base b, em seguida devemos encontrar a representao de g na base b. Uma vez encontrada essa representao, dividimos X i por g; o resto da diviso o primeiro algarismo da representao de X i em base g o algarismo em si est representado na base b e o quociente dever ser dividido novamente no passo seguinte. Ao contrrio do ocorre no procedimento para a mudana de base na parte inteira, a mudana de base para a parte fracionria envolve operaes de multiplicao. A partir da expresso (1.4) possvel observar que g X f = d1 +j=2

dj g j+1 ,

o que nos fornece o primeiro dgito da parte fracionria note que no lado direito da expresso anterior o primeiro termo um inteiro e o segundo um fracionrio. Repetindo a operao com a parte fracionria de g X f temos

gj=2

dj g j+1 = d2 +j=3

dj g j+2 .

E assim por diante recuperamos os termos da parte fracionria de X na base g. Aqui vale a mesma observao do pargrafo anterior, as operaes de multiplicao devem ser realizadas na base b e os dgitos obtidos esto representados na base b. Exemplo: Vamos representar o numeral 53, 205 6 no sistema base-8. Nesse caso X i = 536 e X f = 0, 136 . Para encontrar os dgitos da parte inteira na base 8 devemos realizar sucessivas operaes de diviso por 8, como X i est representado na base 6, devemos realizar todas as operaes nessa base3 , ou seja vamos dividir por 126 (= 810 ): 16 536 = 46 + , 126 126 o resto da diviso, 16 o primeiro dgito (na base 6), na base 8 temos o mesmo dgito, ou seja, 18 . Em seguida, vamos dividir o quociente 4 6 : 46 46 = 06 + , 126 1263

Para facilitar as operaes conveniente utilizar uma tabela de multiplicaes em base 6 (tabuada em base 6): 16 26 36 46 56 16 16 26 36 46 56 26 26 46 106 126 146 36 36 106 136 206 236 46 46 126 206 246 326 56 56 146 236 326 416

Captulo 1. Representao de nmeros em mquinas

8

aqui o resto da diviso 46 = 48 e o procedimento termina pois o quociente da diviso nulo. Portanto a parte inteira representada pelo numeral 41 8 . Encontramos a representao da parte fracionria atravs de operaes de multiplicao por 8 na base 6, ou seja, 126 : 0, 205 6 126 = 2, 506 = 26 + 0, 506 , a parte inteira da multiplicao 26 o primeiro dgito aps a vrgula, em base 8, 2 8 . Em seguida multiplicamos a parte fracionria 0, 50 6 por 8, ou melhor, 126 : 0, 506 126 = 10, 50 6 = 106 + 0, 506 , a parte inteira que resulta da multiplicao 10 6 = 68 e a parte fracionria novamente assume a mesma forma, 0, 506 . Isto nos permite concluir que a partir deste ponto sempre obteremos o mesmo dgito. Portanto 0, 205 6 representa o mesmo nmero que 0, 26 8 . Combinando a parte inteira e fracionria temos nalmente que 53, 205 6 = 41, 26 8 .

1.2. Aritmtica de mquinaNo ocidente a utilizao de um sistema posicional base-10 com algarismos indo-arbicos corrente desde pelo menos o sculo XIV, no entanto, em outras culturas comum encontrarmos sistemas ou pelo menos o seu reexo na linguagem base-5, base-8, base-9, base-12 e mesmo a utilizao matemtica de sistemas posicionais base-20 (civilizao maia) e base-60 utilizada pelos sumrios e babilnios com reexos at hoje na notao para medir ngulos em graus minutos e segundos e nas unidades de tempo minuto e segundo. Os computadores digitais atuais, em quase sua totalidade, representam internamente os nmeros em base-2 (base binria) e realizam operaes aritmticas nessa base. No incio da computao digital chegaram a ser construdas mquinas que representavam nmeros em base ternria (base-3) e em base decimal mas as idias foram abandonadas devido ao fato de que os componentes de hardware so de fabricao mais complexa do que os necessrios para uma mquina que trabalhe em base binria. No entanto muitas das calculadoras atuais realizam operaes utilizando uma aritmtica de base no binria apesar de armazenar os nmeros dessa forma.

1.2.1. Representao de nmeros inteiros Tipicamente um nmero inteiro armazenado em um computador como uma seqncia de dgitos binrios de comprimento xo, dessa forma possvel representar os nmeros de zero a 2 d 1 se forem utilizados d dgitos ou bits exclusivamente para representar um nmero positivo. Se houver

Captulo 1. Representao de nmeros em mquinas

9

necessidade de representar os inteiros negativos ento um dos bits deve ser reservado para o sinal. Por exemplo, um registro de mquina formado por uma seqncia de 32 bits : s d30 d29 ... d2 d1 d0

mos representar os nmeros uma seqncia binria de 32 dgitos (ou bits) s d 30 d29 . . . d2 d1 d0 pode representar todos os inteiros entre 2 31 + 1 e 231 1 (2147483647 e 2147483647) na forma (1)s d0 20 + d1 21 + . . . + d30 230 . Nesse caso o maior inteiro representvel dado numeral I max = 11 . . . 1112 . Portanto Imax = 1 + 21 + 22 + . . . + 230 = 231 1 = 2.147.483.647. Levando em conta os nmeros negativos e que o zero possui duas representaes possveis (como 0 e 231 1. 0), essa disposio de dados na registro permite representar os 2 32 1 inteiros entre 231 1 e De maneira geral essas tcnica de registro com n bits permite armazenar exatamente em uma Se em uma operao de soma ou subtrao o resultado for um nmero que no pode ser armaze-

pode ser utilizado para representar o algarismo binrio ((1) s d30 d29 . . . d2 d1 d0 )2 . Ou seja, podere-

mquina os 2n 1 nmeros inteiros entre 2n1 1 e 2n1 1.

nado nos registro ocorre um erro conhecido como overow, nesse caso a mquina deve ser capaz de reconhecer o evento e enviar uma mensagem de erro se no o zesse, poderia retornar um nmero truncado que no corresponde ao resultado correto da operao programada. Observao. Uma forma alternativa de registrar inteiros positivos consiste em utilizar todos os bits do registro para representar os inteiros, neste caso, os n bits seriam capazes de registrar os 2 n inteiros entre 0 e 2n 1. Se houver necessidade de um nmero de inteiros maior do que o suportado sempre possvel combinar dois ou mais registros para armazenar o nmero em questo. O nico impedimento a quantidade de espao existente na mquina para os registros de nmeros. 1.2.2. Representao de nmeros com parte fracionria possvel estender a representao utilizada para os inteiros para representar com preciso nita, nmeros que possuam parte fracionria. Nesse caso, alm do sinal, uma parte do registro utilizada para a parte inteira do nmero e o restante para a parte fracionria. Assim a seqncia de 32 bits s d15 d14 . . . d13 d14 d15 pode ser interpretada, por exemplo, como um registro onde os 16 bits d15 , . . . , d0 representam a parte inteira e os demais 15 bits, d 1 , . . . , d15 , a parte fracionria: (1)s d15 215 + d14 214 + . . . + d0 20 + d1 21 + . . . + d15 215 .

Captulo 1. Representao de nmeros em mquinas

10

Utilizando essa representao o maior nmero representado 2 16 215 e o menor nmero 216 tipo de representao de nmeros conhecido com representao de ponto-xo.

215 . Nesse intervalo s podemos representar exatamente os nmeros em intervalos de 2 15 . Esse De forma geral, utilizando p bits para representar a parte inteira, q bits para a parte fracionria e

1 bit para o sinal podemos representar os nmeros fracionrios no intervalo [2 p + 2q , 2p 2q ] em intervalos igualmente espaados de 2 q . O menor numero representvel tambm possui esse valor. Se em alguma operao o resultado for um nmero menor em mdulo que 2 q dizemos que ocorreu um erro de underow, em particular a regio compreendida pelo intervalo (2 q , 2q ) denominada regio de underow. Da mesma maneira, se o resultado de alguma operao for maior em mdulo que 2p 2q , dizemos que ocorreu um erro de overow e a regio (, 2 p + 2q ) (2p 2q , +) denominada regio de overow. De maneira geral, a cardinalidade do conjunto de nmeros representados por uma implementao de ponto xo com sinal, p dgitos para a parte inteira e o zero negativo 0). e q dgitos para a parte fracionria em base natural b dada por 2 b p+q (contado o zero positivo +0 Alm da limitao a uma quantidade nita de nmeros imposta representao em mquina, a

representao em ponto xo possui a desvantagem de representar nmeros distintos com preciso diferente. Por exemplo os nmeros 9999, 1234 e 0, 0012113 so representados em base 10 com 4 dgitos para a parte inteira e 4 dgitos para a parte fracionria da seguinte forma: 9999, 1234 e 0, 0012. Enquanto que no primeiro caso o nmero representado com oito dgitos o segundo dispe de apenas dois dgitos para represent-lo. A representao em ponto xo possui essa assimetria, representa com maior preciso os nmeros maiores em valor absoluto e com baixa preciso os nmeros mais prximos de zero. Para contornar essa assimetria utilizamos a representao de ponto utuante. 1.2.3. Representao de ponto utuante Denio (ponto utuante). A representao de um nmero real x denominada ponto utuante normalizado na base b, b N, se forem satisfeitas as propriedades 1. x = m be , onde 2. m = 0, d1 d2 . . . dn 3. 1 d1 b 1 e 0 di b 1 n N, para i = 2, 3, . . . n,

4. e1 e e2 , onde e, e1 , e2 Z.

m denominada mantissa, e expoente e n o nmero de dgitos de preciso. Exemplo: O nmero 9999, 1234 em representao de ponto utuante em base 10 com 8 dgitos de preciso 0, 99991234 104 . Utilizando essa mesma prescrio, o nmero 0, 0012113 Em uma representao de ponto utuante normalizado, o primeiro algarismo aps a vrgula necessariamente maior ou igual a 1. Portanto o nmero zero est fora dos casos cobertos pela denio de ponto utuante. Usualmente o inclumos em um conjunto denominado sistema de ponto utuante representado como 0, 12113131 102 .

Captulo 1. Representao de nmeros em mquinas

11

onde possui a representao 0, 00 . . . 0 b e . O sistema de ponto utuante F (b, n, e 1 , e2 ) denido como o conjunto de nmeros que inclui o zero e os pontos utuantes em base b com n dgitos de preciso e expoente que pode variar entre e1 e e2 inclusive. Propriedades do conjunto F (b, n, e1 , e2 ) Ao contrrio do que ocorre com os nmeros representados por um esquema de ponto xo, os elementos do conjunto F (b, n, e1 , e2 ) no so igualmente espaados. Para exemplicar essa proprieelemento positivo mais prximo de zero o numeral 0, 10 10 10 , o numeral seguinte 0, 11 1010 dade vamos considerar o sistema dado pelo conjunto F (10, 2, 10, 10). De acordo com a denio, o

e assim por diante at o numeral 0, 99 10 10 . O espaamento entre eles4 de 0, 01 1010 . Aps o numeral 0, 99 1010 , vm os numerais 0, 1 109 , 0, 11 109 , . . . , 0, 99 109 . Agora o espaa-

zero.

espaamento 0, 01 1010 . Portanto os elementos esto mais densamente acumulados em torno do A cardinalidade de um sistema de ponto utuante F (b, n, e 1 , e2 ) dada por |F (b, n, e1 , e2 )|: |F (b, n, e1 , e2 )| = 1 + 2(b 1)bn1 (e2 e1 + 1).

mento j 0, 01 109 e assim por diante at os maiores numerais 0, 10 10 10 , . . . , 0, 99 1010 cujo

(1.5)

A cardinalidade do conjunto calculada a partir de todas as combinaes possveis para a representao de um numeral como elemento de F (b, n, e 1 , e2 ), somadas ao elemento zero que no pode ser representado segundo a denio de ponto utuante usual. O primeiro termo do lado direito de (1.5) deve-se ao zero. O fator 2 deve-se ao sinal. O fator (b 1) deve-se aos possveis valores que o dgito assumir. d1 pode assumir. O fator bn1 deve-se combinao dos b possveis valores que os dgitos d 2 , . . ., dn

podem assumir. E nalmente, o fator (e 2 e1 + 1) deve-se aos possveis valores que o expoente pode Um outra propriedade importante dos pontos utuantes diz respeito s relaes algbricas que ao

contrrio dos reais, racionais e inteiros, em geral no so vlidas. Vamos representar as operaes de adio, subtrao, multiplicao e diviso em ponto utuante, respectivamente, pelos smbolos , , , . Dados trs numeros com respresentao em ponto utuante x, y e z, em geral x y = x + y, x y = x y, (x y) z = x (y z), x (y z) = (x y) (x z).

Note que a diferena entre esses primeiros numerais menor do que o menor numeral representvel pelo sistema, ou seja, se for realizada uma operao de subtrao em ponto utuante entre quaisquer dois elementos consecutivos o resultado ser nulo.

4

Captulo 1. Representao de nmeros em mquinas

12

Ento

Exemplo: Sejam x = 0, 1 103 e y = z = 0, 1 1010 , elemntos do conjunto F (10, 3, 10, 10). z) = 0, 1 103 (0, 1 1010 = 0, 1 103 0, 0 100 = 0, 1 103 , 0, 1 1010 )

x (y

por outro lado (x y) z = (0, 1 103 0, 1 1010 ) = 0, 1 1010 = 0, 0 100 . 0, 1 1010 0, 1 1010

Representao em ponto utuante IEEE754 O padro IEEE754 (a sigla se refere ao Institute of Electrical and Electronics Engineers) foi desenvolvido com o objetivo de unicar as diversas implementaes em mquina de registros e operaes em ponto utuante. A maioria dos processadores e compiladores atuais o suportam ou pelo menos suportam um subconjunto obrigatrio das denies. interessante observar que alm dos numerais subnormais (menores do que o usualmente suportado em uma notao F (2, n, e 2 , e1 )) e os N aN (nota a number reservado para operaes ilegais como razes de nmeros negativos). O padro prev quatro tipos de registros: registros de 32 bits denomindos pontos utuantes de preciso simples, registros de 43 ou mais bits para preciso simples estendida, registros de 64 bits para preciso dupla e registros de 79 ou mais bits para preciso dupla estendida. A implementao de preciso simples obrigatria, as demais so opcionais. A ttulo de ilustrao, vale a pena estudar os registros de preciso simples. Os 32 bits do registro de um numeral em preciso simples so divididos de acordo com o diagrama a abaixo, s e07 ... e00 m22 ... m00 em ponto utuante, um mesmo registro pode conter informao sobre +, , +0, 0, numerais

o bit s responsvel pelo sinal, os bits e 07 e06 . . .e00 so responsveis pelo expoente e nalmente os bits m22 m21 . . .m00 so responsveis pela mantissa. Com os 8 bits do expoente podemos representar inteiros entre 0 e 255, no entanto os registros relativos ao 0 (00000000) e 255 (11111111) so reservados para uso especial, sobram portanto os inteiros entre 1 e 254. Segundo o padro, o inteiro relativo ao expoente est deslocado de 127, ento os 8 bits permitem representar os valores inteiros entre 126 e 127. Os 23 bits restantes so utilizados para representar a mantissa com 24 dgitos binrios (j que o primeiro dgito sempre igual a 1 em uma base binria, no h necessidade explcita de

armazen-lo no registro e com isso ganhamos um bit extra) com uma diferena: no padro IEEE754

Captulo 1. Representao de nmeros em mquinas

13

os pontos utuantes normalizados comeam com o primeiro dgito esquerda posicionado antes da vrgula. Portanto o registro de 32 bits capaz de armazenar os elementos no nulos do sistema F (2, 24, 126, 127) e mais os casos especiais: 1. zeros: bits do expoente e da mantissa todos nulos. O bit de sinal pode ser igual a 0 ou 1, ou seja, 2. subnormais: bits do expoente todos nulos e os da mantissa guardam informao sobre o subnormal. 3. innitos: bits do expoente todos iguais a 1 e os da mantissa iguais a 0. O bit de sinal pode ser 4. N aN : bits do expoente todos iguais a 1 e os demais bits contm informao de diagnstico. Exemplo: Vamos encontrar o registro equivalente ao numeral 1345, 875. O primeiro passo representar o nmero na base binria: 1345, 875 = 10101000001, 1112 . Em seguida vamos reescrev-lo como um ponto utuante normalizado: 1, 0101000001111 2 10 . igual a 0 ou 1, ou seja, h uma representao para + e . h uma representao para +0 e 0.

Ignoramos o 1 antes da vrgula e adicionamos tantos 0 direita quantos forem necessrios para preencher os 23 bits, dessa forma encontramos os bits da mantissa: 01010000011110000000000.

O expoente vale 10, com o deslocamento de 127, o inteiro a ser representado pelos bits do expoente o 137 = 100010012 . O bit de sinal igual a 0 pois o nmero positivo. O registro de 32 bits completo dado por 01000100101010000011110000000000.

1.3. ErrosNa anlise numrica de fundamental importncia saber estimar a exatido do resultado de um clculo computacional. Esta seo contempla a origem dos erros que podem alterar os resultados de uma computao e os mtodos utilizados para estimar os erros presentes nos resultados. 1.3.1. Origem dos erros De maneira geral, podemos classicar os erros que afetam o resultado de um procedimento computacional como erros nos dados de entrada, erros de arredondamento e erros de truncamento. Os erros nos dados de entrada so aqueles relacionados alguma medida fsica. Os aparelhos utilizados para medio sempre possuem uma preciso nita em geral no muito grande quando

Captulo 1. Representao de nmeros em mquinas

14

comparada preciso que possvel ser obtida na representao de nmeros em mquinas e nem sempre possvel melhor-la consideravelmente. Os erros de truncamento so os mais comuns em algoritmos numricos. Ocorrem quando, de alguma maneira, necessrio aproximar um procedimento formado por uma seqncia innita de passos atravs de um outro procedimento nito. Os erros de arredondamento so aqueles relacionados s limitaes que existem na forma de representar nmeros em mquinas. Qualquer que seja a natureza do erro, o mtodo numrico utilizado em um dado problema deve ser capaz de estimar as suas conseqncias sobre o resultado: os dados de sida. Vamos discutir nesta seo como propagar e controlar essa incerteza presente nos valores nas diversas operaes prprias a um procedimento computacional. Devido ao seu carter, os erros de truncamento sero estudados em conjunto com os algoritmos que os geram. 1.3.2. Conceitos iniciais Seja x um nmero que se conhece exatamente e x uma representao nita, ou aproximao, de 1 x, por exemplo x = e x = 3, 14159, ou ainda x = e x = 0, 333. Ento denimos erro absoluto e 3 erro relativo como: Denio (Erro absoluto e erro relativo). O erro absoluto na representao x denido por |x|. x |x x| O erro relativo denido como . |x|

Em uma mquina, denominamos preciso o nmero de dgitos na mantissa. Por outro lado, como o prprio nome diz, a exatido de uma aproximao uma medida de quanto ela est prxima do valor exato. Uma maneira de estimar a exatido utiliza o conceito de dgitos exatos de uma representao. Denio (Dgitos exatos de x). Seja a aproximao x de um nmero x em base b. Dizemos que x possui k dgitos exatos se k for o maior inteiro tal que |x x| 1 bk . |x| 2 Atravs dessas denies podemos apreciar a diferena entre os conceitos de preciso e exatido de uma representao x. Exemplo: Seja x = e x = 3, 141675 = 0, 3141675 10 1 . Nesse caso a preciso da represen |x x| tao x de 7 dgitos, no entanto < 0, 3 104 < 0, 5 104 , ou seja, x possui apenas 4 |x| dgitos exatos. 1.3.3. Arredondamento O operao de arredondamento consiste em encontrar um representao x para um nmero x com uma determinada preciso. Essa operao usualmente realizada em mquinas para representar

Captulo 1. Representao de nmeros em mquinas

15

internamente os nmeros e o resultado de operaes aritmticas com eles. No h uma nica forma de realizar o arredondamento e, de acordo com a aplicao, existem algoritmos mais convenientes. Vamos discutir os 5 mais comuns. Por comodidade e para melhor ilustrar, vamos discuti-los a partir dos mesmos exemplos, sempre utilizando a mantissa 0, d 1 d2 . . . dn . Entre os possveis algoritmos para o arredondamento vamos nos xar nos seguintes: Seja a sequncia {0, 45 ; 0, 44 ; 0, 36 ; 0, 35 ; 0, 34 ; 0, 34 ; 0, 35 ; 0, 36 ; 0, 44 ; 0, 45 ; 0, 46}.

1. No sentido do zero: o ltimo dgito descartado. O algoritmo possui esse nome pois o efeito dele aproximar os nmeros do zero. {0, 4 ; 0, 4 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 4 ; 0, 4 ; 0, 4}. 2. Se afastando do zero: o penltimo dgito arredondado para o algarismo seguinte se o ltimo dgito for no nulo. O algoritmo possui esse nome pois o seu efeito afastar os nmeros do zero. {0, 5 ; 0, 5 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 5; 0, 5; 0, 5}. 3. No sentido de +: o penltimo dgito arredondado para o algarismo seguinte se o numeral for positivo e o ltimo dgito for no nulo. Se o numeral for negativo e o ltimo dgito no nulo, o penltimo algarismo arredondado para o algarismo anterior. O algoritmo possui esse nome pois o seu efeito deslocar os nmeros de um valor maior ou igual a zero (ou seja no sentido de +). {0, 4 ; 0, 4 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 5 ; 0, 5 ; 0, 5}. 4. No sentido de : o penltimo dgito arredondado para o algarismo anterior se o numeral for penltimo algarismo arredondado para o algarismo seguinte. O algoritmo possui esse nome pois o seu efeito deslocar os nmeros de um valor menor ou igual a zero (ou seja, no sentido de ). {0, 5 ; 0, 5 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 3 ; 0, 3 ; 0, 3 ; 0, 4 ; 0, 4 ; 0, 4}. 5. Arredondamento par (base decimal): o penltimo digito arredondado para o algarismo seguinte se o ltimo dgito for maior ou igual a 6. Se o ltimo dgito for menor ou igual a 4 o penltimo dgito arredondado para o algarismo anterior. Se o ltimo dgito for igual a 5 e o penltimo dgito for mpar, o penltimo digito arredondado para o algarismo seguinte. Caso o penltimo seja par, ele mantido. O algoritmo possui esse nome pois como resultado os numerais sempre terminam com um algarismo par. {0, 4 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 3 ; 0, 3 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 4 ; 0, 5} positivo e o ltimo dgito for no nulo. Se o numeral for negativo e o ltimo dgito no nulo, o

Captulo 1. Representao de nmeros em mquinas

16

1.3.4. Propagao dos erros Vamos assumir que x uma aproximao de um nmero x, sem perda de generalidade, vamos supor que x > x. Se quisermos encontrar o valor da funo f calculada em x mas s dispormos de x devemos aproximar f (x) por f (). Se a funo f for diferencivel, pelo teorema do valor mdio x existe um [, x]: x f (x) f () = f ()(x x). x

Porm o teorema no diz nada sobre alm de sua existncia. Para estimar o erro cometido na aproximao f () devemos utilizar mais informaes sobre a funo f . Por exemplo, se soubermos x que a derivada de f limitada no intervalo [, x], max y[,x] |f (y)| = M , ento temos que x x |f (x) f ()| = |f ()| |x x| M |x x|. x Mesmo que no sejamos capazes de determinar o valor mximo da derivada da funo f no intervalo . [, x], ainda podemos considerar o caso em que | x| = pequeno. Nesse caso, se f () = 0 e f x x x x

no varia muito rapidamente na vizinhana de x ento como muito pequeno podemos considerar x M f () na expresso acima, dessa forma x . f () = |f (x) f ()| f ()|x x| = f (). x x x x x

O mesmo procedimento pode ser utilizado para o caso de funes que dependam de duas ou mais variveis. Nesse caso seja o erro na i-sima varivel x i de um funo f que dependa de n variveis, ento o erro propagado na funo f resultado dos erros nas n variveis dado por f (x 1 , x2 , . . . , xn ): n

f (x1 , x2 , . . . , xn ) = i=1

f (x1 , x2 , . . . , xn ) xi . xi

1.3.5. Instabilidade numrica Alguns problemas matemticos e algoritmos numricos possuem a propriedade de ampliar drasticamente os erros presentes nos dados de entrada e assim invalidar a sada ou resposta. No contexto do clculo numrico, esse fenmeno denominado instabilidade numrica. A instabilidade numrica pode estar relacionada s propriedades matemticas do problema ou ento estrutura do algoritmo utilizado para resolv-lo. De qualquer maneira, ao estudar um problema que pretendemos reescrever numericamente imprescindvel a anlise de estabilidade do algoritmo ou o condicionamento do problema matemtico. Dizemos que um problema mal condicionado quando pequenas variaes nos dados de entrada resultam em grandes variaes na resposta. O seguinte exemplo ilustrativo. Exemplo: Seja o polinmio P (x) = (x 1)(x 2) . . . (x 10) = x 10 55 x9 + . . . + 10!. Ou

seja, P um polinmio em x com razes inteiras 1, 2, 3, . . . , 10.

Captulo 1. Representao de nmeros em mquinas

17

Vamos considerar agora o polinmio P : P (x) = P (x) + 104 x9 . Ou seja, P igual P a menos de um erro relativo de 1, 8 10 6 no coeciente do termo x9 . As razes do polinmio P so (com quatro casas decimais) x1 = 1, 0000 x2 = 2, 0000 x3 = 3, 0002 x4 = 3, 9940 x5 = 5, 0803 x6 = 5, 7283 x7 = 7, 4078 0, 6753 i x8 = 7, 4078 + 0, 6753 i x9 = 9, 6907 0, 3513 i x10 = 9, 6907 0, 3513 i

Um erro relativo de 1, 8 106 no coeciente x9 foi capaz de alterar drasticamente parte das razes. Como entender esse fenmeno? Vamos avaliar quanto varia cada raiz x j (j = 1, 2, . . . , 10) quando alteramos o coeciente 55 do termo x9 . Para tanto vamos considerar o novo polinmio P (, x) = P (x) x 9 . Esse

novo polinmio possui dez razes xj (), j = 1, 2, . . . , 10, cada uma delas depende de . A dxj derivada contm informao sobre quanto a i-sima raz de P (, x) varia com relao d =0 a quando = 0. Como P (0, x) P (x), essas derivadas possuem a informao que buscamos

sobre as razes de P . dxj no ser necessrio encontrar explicitamente a dependncia das razes Para calcular d =0 em , basta lembrar que por denio P (, x j ()) = 0 (xj () raiz do polinmio P (, x)) e portanto de acordo com a regra da cadeia d P P dxj (P (, xj ())) = (, xj ()) + (, xj ()) =0 d x d e assim,P dxj (, xj ()) = P . d x (, xj ())

O numerador da expresso acima simplesmente

P P 9 (, xj ()) = x . Para calcular x (, xj ()), basta notar que P (, x) = (x1)(x2) . . . (x10)x 9 , portanto P (, x) =(x2) . . . (x x 10)+(x 1)(x 3) . . . (x 10)+. . . + (x 1) . . . (x 9)9x 8 , ou seja

P (, x) = 9x8 + x Assim, dxj = d 9x8 +

10

10

k=1 l=1,l=k

(x l).

10 k=1

x9 () j

10 l=1,l=k (xj ()

l)

.

Captulo 1. Representao de nmeros em mquinas

18

Quando = 0, P (0, x) = P (x) ento xj (0) = j para j = 1, 2, . . . , 10. Ou seja, dxj d = j 9 = j 9 = j .

=0

10 k=1

10 l=1,l=k (j

l)

10 l=1,l=j (j

l)

Na expresso anterior, o somatrio foi retirado pois o nico termo que contribua era o k = j (Por que?). As derivadas que buscamos so representadas pelo smbolo j : 1 = 2, 8 106 3 = 2, 0 6 = 3, 5 103 8 = 1, 3 104 9 = 9, 6 103 10 = 2, 8 103 7 = 9, 3 103

2 = 1, 3 102 4 = 6, 1 101 5 = 6, 8 102

E assim podemos notar que enquanto que a variao no valor das primeiras razes pequena, o mesmo no pode ser dito para as demais. O valor das derivadas explica o mal condicionamento deste problema. O prximo exemplo contempla a instabilidade numrica devido ao algoritmo. . Exemplo: Seja o seguinte algoritmo para calcular a integral I n = 1, 2, . . . Podemos integrar por partes I n : In = xn ex1x=1 x=0 1 1 n x1 , 0 dx x e

onde n =

n

dx xn1 ex1 (1.6)

0

= 1 n In1 .

Como In = 1 nIn1 , se conhecermos I1 , poderemos encontrar I2 , I3 e assim por diante. Antes de dar prosseguimento vamos examinar um pouco mais a integral I n . No intervalo de integrao [0, 1], o integrando sempre positivo, alm disso, o termo e x1 sempre menor ou igual a 1, isto implica a desigualdade 0 < In 0 1

dx xn =

1 . n+1

(1.7)

Portanto devemos esperar que os termos I n decresam com n e sejam sempre positivos. O primeiro termo I1 pode ser calculado explicitamente utilizando integrao por partes,1

I1 =0

dx x ex1 = 1

1 0

dx ex1 = e1 .

Iremos representar I1 aproximadamente por I1 com cinco casas aps a vrgula: I1 = 0, 36788.

Captulo 1. Representao de nmeros em mquinas

19

Nesse caso = |I1 I1 | 5, 6 107 . Se calcularmos os demais termos a partir da iterao (1.6) encontraremos I1 = 0, 36778 I2 = 0, 26424 I3 = 0, 20728 I4 = 0, 17088 I5 = 0, 14560 I6 = 0, 12640 I7 = 0, 11520 I8 = 0, 07840 I9 = 0, 29440 (?) I10 = 1, 944 (???)

O que ocorreu aqui? A estimativa (1.7) prev o decrescimento de I n com n e In > 0. Vamos analisar a iterao levando em conta a aproximao inicial. Partimos portanto de I1 = I1 + , ento segundo a iterao (1.6) I2 = 1 2(I1 ) = 1 2(I1 + ) = I2 2 e ento I3 = 1 3(I2 ) = 1 3(I2 2) = I3 + 3 2 . Por induo temos ento que In = In + (1)n+1 n!. Ou seja, o erro cresce com o fatorial de n. Dessa forma I9 I9 + 9! 5, 6 107 I9 + 0, 2 e I I10 10! 5, 6 107 I10 2. 10 Felizmente uma pequena alterao do algoritmo permite resolver a questo da instabilidade. Baste reescrever a iterao (1.6) como In1 = 1 (1 In ) n (1.8)

isolando In1 . Utilizando os mesmos critrios possvel vericar que o erro cometido no ponto inicial mais e mais diludo a cada passo do algoritmo (verique!) Alm dessa propriedade, a diferena da iterao (1.8) reside no fato de que o ponto inicial deve ser um valor de n sucientemente grande para que os ndices abaixo possam ser calculados pelo algoritmo. Neste ponto a estimativa (1.7) til, pois ela nos diz que I n decresce com n e 1 1 mais In < n+1 . Portanto podemos utilizar um nmero n sucientemente alto, In = 0 e = n+1 .

Captulo 2

Equaes no linearesVamos estudar mtodos numricos para resolver o seguinte problema. Dada uma funo f contnua, real e de uma varivel, queremos encontrar uma soluo x que satisfaa a equao no linear: f (x ) = 0. (2.1)

Em geral a equao (2.1) no pode ser resolvida exatamente, isto , a soluo x no pode ser funes elementares (polinmios, razo entre polinmios, potncias racionais, e as funes transcendentais: exp,log, trigonomtricas, hiperblicas). H casos em que a prpria funo f no conhecida explicitamente: pode ser denida a partir de uma srie innita, ou a partir de uma integral ou ainda ser soluo de uma equao diferencial. nesses casos utilizamos mtodos numricos para resolver a equao. Idealmente, poderamos dividir o o procedimento nas seguintes etapas: inicialmente devemos encontrar uma regio de interesse onde possam existir solues da equao; em seguida, quando possvel, isolar os intervalos que contm apenas 1 soluo; feito isso, determinamos pelo menos 1 aproximao inicial x0 da soluo (de acordo com o mtodo utilizado, pode ser necessrio utilizar mais de uma aproximao inicial) para cada intervalo; nalmente, a partir das aproximaes iniciais, isto , o mtodo numrico consiste na construo de uma seqncia {x n } que converge para a soluo , n=0n+

descrita a partir de uma combinao nita de operaes algbricas simples (+, , /, , exp, log) e

lim xn = x

soluo da equao (2.1). Portanto os mtodos numricos para encontrar a soluo de equaes no lineares so mtodos iterativos. A cada iterao, utilizamos um subconjunto das aproximaes xn1 , xn2 , . . . , x0 , obtidas anteriormente, para determinar a prxima aproximao x n . Estudaremos os mtodos separados em trs classes principais: Mtodos de quebra: o ponto de partida encontrar um intervalo que contenha pelo menos 1 soluo. Segundo o teorema de Bolzano, basta determinar um intervalo em que a funo f muda de sinal. Os mtodos de quebra consistem na descrio de como subdividir o intervalo inicial em intervalo cada vez menores que ainda contenham a mesma soluo. Nesse caso, a seqncia x0 , x1 , x2 , . . . , xn formada pelos extremos dos intervalos. A soluo numrica ser encontrada

Captulo 2. Equaes no lineares

21

quando a largura do intervalo em uma m-sima iterao for pequeno o suciente para satisfazer as exigncias de exatido. Mtodos de ponto xo: dos mtodos. Mtodos de mltiplos passos: Uma generalizao do mtodo anterior onde a funo depende de mais de uma aproximao anterior, i. e., x n+1 = (xn , xn1 , . . . , xnp ) para algum p n. xn+1 = (xn ). A convergncia do mtodo garantida pelo teorema do ponto xo, da o nome A seqncia {x i }i=n construda a partir da sucessiva iterao i=0

2.1. Mtodos de quebraOs mtodos de quebra utilizam como primeira aproximao um intervalo que contenha pelo menos 1 soluo da equao no linear. As iteraes consistem em seguidas subdivises dos intervalos de maneira que o novo intervalo sempre contem a soluo. O uso do teorema comum a todos os mtodos de quebra, ele fornece condies para que os intervalos contenham pelo menos 1 soluo da equao. Apresentamos o teorema sem sua prova. Teorema (Bolzano) intervalo e [f (a), f (b)] f (I) ou [f (b), f (a)] f (I). Seja I = [a, b] R e uma funo f : I R contnua. Ento o conjunto imagem f (I) tambm um Portanto, se encontrarmos um intervalo [a, b] tal que, por exemplo, f (a) < 0 e f (b) > 0, ento pelo teorema de Bolzano existe, um ponto x [a, b] tal que f (x ) = 0. O que difere os mtodos de quebra entre si a maneira com que os intervalos so subdivididos.

2.1.1. Mtodo da bisseco A aproximao inicial consiste em um intervalo [x 0 , x1 ] tal que f (x0 )f (x1 ) < 0. Ento dividimos x0 + x 1 o intervalo ao meio, ou seja, no ponto x 2 = . Entre os intervalos [x0 , x2 ] e [x2 , x1 ] escolhe2 mos aquele que possui pelo menos 1 soluo, ou seja, aquele em que f (x 2 )f (xi ) < 0 para i = 0, 1. Pode ainda ocorrer que x2 seja soluo, ou seja, x2 = x ou ainda que ambos os intervalos sejam tais que f (x0 )f (x2 ) < 0 e f (x2 )f (x1 ) < 0, nesse caso, ambos subintervalos contm pelo menos 1 soluo e podemos continuar o procedimento em cada um deles em separado. Observao. Se a aproximao inicial [x 0 , x1 ] for tal que f (x0 )f (x1 ) > 0 isto no quer dizer que no exista soluo nesse intervalo, apenas o teorema no permite uma concluso sobre a existncia ou no de soluo nesse intervalo. Nesse caso necessrio escolher outro intervalo ou ento realizar um diviso adicional. Por exemplo, se f (x) = x(1 x), a equao no linear f (x ) = 0 possui solues x = 0 e x = 1, porm f (1)f (3) > 0. Devemos adotar um critrio de parada no processo de subdiviso dos intervalos. Em geral adotamos dois valores pequenos 1 2 de maneira que a soluo escolhida uma vez que uma das trs desigualdades seguintes seja satisfeita: |x n xm | < 1 , f (xn ) < 2 ou f (xm ) < 2 .

Captulo 2. Equaes no lineares

22

O seguinte algoritmo descreve com maior detalhe todos os passos. A entrada do programa consiste nos extremos do intervalo inicial [a, b], a funo f , os parmetros de exatido 1 e 2 e o nmero mximo de passos aceitvel N . A sada pode ser mensagens de erro, o intervalo mnimo (de comprimento 1 ) onde a soluo se encontra (se no for possvel obter a exatido pretendida, 2 ) ou a soluo x com exatido dada pelo parmetro 2 . Por comodidade, os comentrios esto colocados aps o algoritmo. 1. entrada {a, b, f, 1 , 2 , N }

2. x0 a; x1 b; f0 f (x0 ); f1 f (x1 ); i 1 3. se f0 f1 > 0, ento: final 4. enquanto sada{Erro

nos dados

de entrada};

v para:

a) se |x0 x1 | < 1 ento: d) i i + 1 5. se i > N , ento: v para: para: final 6. se f0 2 ento: final 7. final: Comentrios: termina c) se f2 f0 < 0, ento:

|f0 | > 2 e |f1 | > 2 e i < N

sada{ x , x1 }; v para: 0

final x0 x2 ; f0 f2 exigida em ,N, passos.};

b) x2 0, 5 (x0 + x1 ); f2 f (x2 )

x1 x2 ; f1 f2 , seno: atingiu a exatido final,

sada{No

sada{ x }; v para: 0 o programa.

seno:

sada{ x}; v 1

linha 2: as variveis x0 e x1 guardam o valor dos extremos dos intervalos, f 0 e f1 guardam o valor da funo nesses extremos. A varivel i um contador. linha 4: lao em que os intervalos so divididos. Se alguma das condies for satisfeita encontramos a soluo com exatido desejada ou excedemos o nmero mximo de passos. linha 4(a): os intervalo atingiu o menor valor admissvel. O intervalo que contem a soluo retornado como sada, porm a exatido no foi atingida. linha 4(c): teste para encontrar o intervalo que contem pelo menos 1 soluo com certeza. linha 4(d): incremento do contador linha 5: se o cdigo est nesse ponto, o lao 4 terminou, ou seja, ou encontramos a soluo ou excedemos o nmero mximo de passos. armao verdadeira. linha 6: no necessrio testar f1 pois se i > N falso e f0 2 falso, s resta f1 2 como

dada em termos da funo especial W de Lambert:

Como exemplo do mtodo, vamos estudar a equao no linear para f (x) = x e x . A soluo x = W (1) 0, 5671432904097839 . . .

Captulo 2. Equaes no linearesf x 0.5 0.25 1 4 1 x 2 3 4 1 x

23

0.25 0.5 0.75 1

Figura 2.1. Grco da funo f (x) = x ex , no intervalo x [0, 1].

com intervalo inicial (0.0 , 1.0):

A tabela seguinte ilustra o comportamento dos extremos do intervalo para a equao x e x = 0 iterao 1 2 3 4 5 6

x0 0, 5 0, 5 0, 5 0, 5625 0, 5625 0, 5625

x1 1,0 0, 75 0, 625 0, 625 0, 59375 0, 578125

Tabela 2.1. Tabela das primeiras iteraes para o mtodo da bisseco.

Aps 20 iteraes chegamos ao intervalo (0.567142, 0.567143). O valor 0, 567143 satisfatrio como soluo com 6 casa decimais exatas. Limitaes do mtodo No caso de razes mltiplas de polinmios, pode no existir um intervalo onde a funo troca de sinal. Ou pelo menos pode ser difcil encontrar tal intervalo. Veja os grcos abaixo.f x f x

x

x

Captulo 2. Equaes no lineares

24

2.1.2. Mtodo da falsa posio ou regula falsi A diferena bsica entre este mtodo e o mtodo da bisseo est na forma de dividir o intervalo. O mtodo da falsa posio utiliza como ponto intermedirio para diviso do intervalo (x 0 , x1 ), o ponto dado pela interseco entre o eixo x e a reta que une os pontos (x 0 , f (x0 )) e (x1 , f (x1 )). A reta que une esses dois pontos possui equao (x): (x) = f (x0 ) f (x1 ) x0 f (x1 ) x1 f (x0 ) x+ . x0 x 1 x0 x 1

Portanto, o ponto intermedirio xm dado por xm = x 0 (x0 x1 ) f (x0 ). f (x0 ) f (x1 )

f x x x0 x xm x1

Figura 2.2. A reta que une os pontos (x0 , f (x0 )) e (x1 , f (x1 )) est pontilhada. Ela cruza o eixo x no ponto que divide o intervalo, xm

com intervalo inicial (0.0 , 1.0):

A tabela seguinte ilustra o comportamento dos extremos do intervalo para a equao x e x = 0 iterao 1 2 3 4 5 6 x0 0, 0 0, 0 0, 0 0, 0 0, 0 0, 0 x1 0,6127 0, 572181 0, 567703 0, 567206 0, 567150 0, 567144

Tabela 2.2. Tabela das primeiras iteraes para o mtodo da falsa posio.

Aps 7 iteraes chegamos ao resultado nas mesmas condies (6 casas decimais de exatido) utilizadas no mtodo anterior.

Captulo 2. Equaes no lineares

25

2.2. Mtodos de ponto xoOs mtodos de ponto xo so caracterizados por reescrever o equao no linear f (x ) = 0 na forma (x ) = x e utilizar o teorema do ponto xo que veremos logo adiante para garantir a convergncia da seqncia xn+1 = (xn ) para o ponto xo x que soluo de (2.2). Seja portanto a funo : [a, b] R (x) = x + (x)f (x), onde (x) = 0 no intervalo [a, b]. Nesse caso, se (x ) = x , ento como (x) = 0 para todo x [a, b], f (x ) = 0. A soluo x ser ento determinada atravs da convergncia da seqncia {x n } , limn xn = n=0 (2.2)

x , onde xn+1 = (xn ). A garantia da convergncia tratada pelo teorema do ponto xo: Teorema (ponto xo)

Seja uma funo contnua, denida em um intervalo I = [a, b] e tal que as seguintes condies sejam satisfeitas:

(I) I, obs: (a notao indica x I, g(x) I) x I, | (x)| L < 1 obs:( uma contrao)

Ento dado qualquer x0 I, a seqncia xn+1 = (xn ) converge para o ponto x = (x ). O ponto x existe e nico no intervalo I.

Demonstrao: Vamos tratar inicialmente a questo da convergncia (existncia de ponto xo). Seja a distncia entre dois pontos consecutivos x n+1 e xn da seqncia : |xn+1 xn | = |(xn ) (xn1 )|. | (c)||xn1 x |. De acordo com as hipteses, tal que (I) I, ento se x 0 , a aproximao inicial, pertence a I ento (xn ) e (xn1 ) tambm lhe pertencem. Como c [xn , xn1 ] I Segundo o teorema do valor mdio, existe um c [x n , xn1 ] tal que |(xn ) (xn1 )| =

Captulo 2. Equaes no lineares

26

ento segundo as hipteses | (c)| L < 1, ento |xn+1 xn | = |(xn ) (xn1 )| = | (c)||xn xn1 | L|xn xn1 |. Utilizando recursivamente a desigualdade (2.3) temos |xn+1 xn | L|xn xn1 | L2 |xn1 xn2 | . . . Ln |x1 x0 |, portanton

(2.3)

lim |xn+1 xn | lim Ln |x1 x0 | = |x1 x0 | lim Ln .n n

Novamente segundo as hipteses, L < 1 e |x 1 x0 | um nmero nito pois x1 e x0 pertencem ao intervalo nito [a, b]. Assim lim n Ln = 0 en

lim |xn+1 xn | = 0.

Ou seja, a seqncia converge para um x = (x ). Dessa forma, existe pelo menos um ponto x no intervalo [a, b] que satisfaz a equao x = (x ). A seguir vamos vericar que esse ponto nico. Sejam x1 e x2 dois pontos distintos no intervalo I = [a, b] que satisfazem a equao x = (x), ou seja, x1 = (x1 ) ex2 = (x2 ). Ento, de acordo com o teorema do valor mdio, existe um c [x1 , x2 ] tal que |x1 x2 | = | (c)||x1 x2 |. Segundo as hipteses, x1 , x2 I, ento (c) L < 1, ou seja |x1 x2 | < |x1 x2 | o que uma contradio. Portanto, no intervalo I = [a, b] h um e somente um ponto x = (x).

Observao. Note que na demonstrao do teorema do ponto xo, fundamental que a derivada de seja estritamente menor do que 1 em alguma vizinhana I que contm a soluo. Caso contrrio, se uma seqncia cclica de pontos sem convergir para a soluo x , ou mesmo a possibilidade de haver | (x)| 1 em um intervalo I, no podemos excluir a possibilidade de que as iteradas transitem por

Captulo 2. Equaes no lineares

27

mais de uma soluo nesse intervalo. Naturalmente isto no quer dizer que esses comportamentos ocorrem sempre que as hipteses do teorema no forem vlidas.

2.2.1. Mtodo da iterao linear Trata-se de encontrar uma funo que satisfaa as hipteses do teorema do ponto xo para alguma vizinhana em torno da soluo x da equao f (x ) = 0. Como a funo construda a partir de uma outra funo (x) = 0 em um intervalo que contem a soluo de f (x ) = 0, encontr-la signica determinar (x) = 0. A condio de convergncia garantida ento pelo teorema do ponto xo se as suas hipteses forem satisfeitas. f (x) = 0 x = ex = (x). Portanto, como por denio, (x) = x+(x)f (x), no nosso exemplo do ponto xo so vlidas apenas nos intervalos 1 I [W (1), 1], onde (I) I e | (x)| < 1. termos por iterao n 1 2 3 4 5 6 xn 0, 606531 0, 545239 0, 579703 0, 560065 0, 571172 0, 565863 Vamos considerar o exemplo que j estudamos anteriormente, f (x) = x e x . Nesse caso

(x) 1. Assim (x) = 0 para qualquer valor de x. Como | (x)| = ex , as hipteses do teorema Vamos escolher ento a aproximao inicial x 0 = 0, 5. a seqncia dada em seus primeiros

Tabela 2.3. Tabela das primeiras iteraes para o mtodo da iterao linear com (x) = e x .

A soluo com 6 dgitos exatos alcanada aps 22 iteraes. Uma outra possibilidade para a funo seria a escolha (x) = ln x que corresponde a (x) = x + ln x 1 que sempre negativa no intervalo (0, +). No entanto, | (x)| = maior do que a x ex |x| unidade no intervalo (0, 1) que contem a soluo e assim, o teorema do ponto xo no d garantias de convergncia. Podemos perceber que logo nas primeiras iteraes, a seqncia toma valores negativos e, dessa forma, como (x) = ln x, a seqncia no estar denida apenas nos nmeros reais. Em particular essa seqncia no converge para nenhuma soluo de f (x) no plano complexo (a equao possui innitas solues l).A funo W funo W de Lambert e esta relacionada soluo da equao x = ex . Na prtica, no procuramos garantir a hiptese (I) I pois muitas vezes, determinar esse intervalo exatamente equivale a resolver a equao no linear.1

Captulo 2. Equaes no lineares

28

2.2.2. Mtodo Newton-Raphson A partir da demonstrao do teorema do ponto xo, podemos notar que quanto menor for o limite superior L < 1 para o valor absoluto da derivada de na vizinhana da soluo x mais rapidamente a seqncia converge para a soluo da equao no linear. O mtodo de Newton-Raphson um mtodo iterativo que utiliza essa propriedade da convergncia das seqncias para garantir uma convergncia rpida para a soluo a partir do instante que x n+1 se aproxima de uma vizinhana sucientemente prxima de x . Portanto, a idia determinar uma funo (x) tal que (x ) = 0 e assim garantir que, em uma vizinhana prxima de x , a funo tal que | | Tomando a derivada de , por denio temos: 1.

(x) = 1 + (x)f (x) + (x)f (x) e em x = x , soluo da equao f (x ) = 0, temos (x ) = 1 + (x )f (x ). Portanto, a escolha (x) = escolha (2.4) para a funo , a funo dada por

1 f (x)

(2.4)

implica (x ) = 0 de maneira que na vizinhana de x , | | assume pequenos valores. A partir da f (x) . f (x) (2.5)

(x) = x

Vamos novamente utilizar o exemplo f (x) = x e x , nesse caso a iterao dada pela funo : (x) = x Partindo da aproximao inicial x0 = 0, 5: iterao n 1 2 xn 0, 566311 0, 567143(x+1) 1+ex .

(x + 1) x ex = . x 1+e 1 + ex

Tabela 2.4. Tabela das primeiras iteraes para o mtodo Newton-Raphson com (x) =

a seqncia converge para a soluo exata at a 6 a casa decimal em duas iteraes. Se utilizarmos x0 = 1, 0 como aproximao inicial obteramos o mesmo resultado aps trs iteraes. Vamos analisar com um pouco mais de detalhe a questo da convergncia. Se a funo sucientemente bem comportada a ponto de admitir uma expanso em srie de Taylor em torno da soluo

Captulo 2. Equaes no lineares

29

x , ento |xn+1 x | = |(xn ) x | = 1 (x ) + (x )|xn x | + (x )|xn x |2 x 2 (2.6)

1 | (x )||xn x | + | (x )||xn x |2 + O(|xn x |3 ). 2 Como a derivada dada por (x) = f (x)f (x) , (f (x))2

(2.7)

podemos concluir que se f (x ) = 0 ento (x ) = 0. E assim, a desigualdade (2.6) assume a forma |xn x | com | (x )| = 1 | (x )||xn1 x |2 + O(|xn1 x |3 ) 2

f (x ) . Ou seja, se f (x ) = 0 ento a convergncia quadrtica pelo menos. No f (x ) entanto se f (x ) = 0 (por exemplo, no caso de razes mltiplas), a derivada de no ponto x no se anula. Se realizarmos uma expanso de Taylor para (2.7) encontraremos, no caso em que f (x ) = 0 e f (x ) = 0, (xn ) = ou seja (x ) = 1 1 f (x ) + (xn x ) + O((xn x )2 ), 2 3! f (x )

1 se f (x ) = 0. E assim a desigualdade (2.6) assume a forma 2 1 |xn x | |xn1 x | + O(|xn1 x |2 ) 2

e a convergncia linear como nos mtodos de iterao linear.

2.3. Mtodos de mltiplos pontos

2.3.1. Mtodo da secante O mtodo da secante similar ao mtodo da falsa posio, diferem entre si pelo fato de que no mtodo da secante no h diviso e escolha de intervalos, a seqncia de aproximaes calculada a partir das duas ltimas aproximaes e portanto, devemos iniciar com duas aproximaes para a soluo. Ao contrrio do mtodo da falsa posio, no h necessidade de que a soluo esteja entre as duas aproximaes iniciais.

Captulo 2. Equaes no lineares

30

A seqncia montada a partir da regra para iterao 2 xn+1 = xn (xn xn1 ) f (xn ). f (xn ) f (xn1 )

De maneira semelhante que ocorre nos mtodos de ponto xo, para que ocorra convergncia, em geral, as duas primeiras aproximaes devem estar em uma vizinhana sucientemente prxima da soluo. possvel demonstrar que existe um constante K tal que |xn+1 x | = K, n |xn x | lim

1+ 5 1, 618. Ou seja, apesar de ser mais lenta que no mtodo Newton-Raphson, a onde = 2 convergncia mais rpida que a convergncia linear do mtodos de ponto xo. iterao n 1 2 3 4 xn 0, 544221 0, 568826 .567150 .567143

Tabela 2.5. Tabela das primeiras iteraes para o mtodo da secante para f (x) = x e x , com aproximaes iniciais x0 = 0, 9 e x1 = 1, 0.

a seqncia converge para a soluo exata at a 6 a casa decimal em quatro iteraes. Se utilizarmos x0 = 0, 5 e x1 = 1, 0 como primeiras aproximaes obteramos o mesmo resultado aps trs iteraes.

2

comum utilizar as seguintes variaes para minimizar os efeitos de arredondamento: xn

f (x ) (xn xn1 ) f (x n ) n1 f (x ) 1 f (x n ) n1

xn+1 =

,

se |f (xn )| < |f (xn1 )| se |f (xn1 )| < |f (xn )|

xn

(xn xn1 )n1 1 f (x ) n f (x )

,

Captulo 3

Sistemas de equaes linearesOs sistemas de equaes lineares fazem parte da descrio matemtica dos mais diversos fenmenos em todas as reas das cincias naturais e tambm so pea fundamental de diversos algoritmos utilizados em computao. Por exemplo, mais adiante na disciplina, veremos que atravs da discretizao dos domnios onde est denida uma equao diferencial possvel reduzi-la a um sistema de equaes lineares. Em outras aplicaes como mecnica dos udos e mecnica estrutural, no incomum trabalhar com sistemas de ordem 10 5 ou mais. Devido a sua quase onipresena , muito importante poder contar com tcnicas que encontre a soluo para os sistemas de modo eciente e acurado. Atualmente podemos contar com software de alta qualidade para resolver sistemas de equaes lineares e ainda hoje este assunto est sendo ativamente pesquisado. Principalmente a resoluo de sistemas muito grandes em clusters de computadores e computadores com processadores vetoriais. Para melhor compreender e utilizar esses softwares essencial conhecer as propriedades dos algoritmos mais simples. Vamos nos concentrar no estudo de sistemas de n equaes e n incgnitas x 1 , x2 , . . . , xn : a11 x1 a21 x1 + a12 x2 + a22 x2 + . . . + a1n xn + . . . + a2n xn . . . = b1 = b2

(3.1)

an1 x1 + an2 x2 + . . . + ann xn = bn

onde os coecientes aij e as constantes bi so nmeros reais. Denominamos o conjunto de todas as possveis solues de um sistema linear de conjunto soluo. Dados dois sistemas lineares, dizemos que os dois so equivalentes se possuirem o mesmo conjunto soluo. Sob um ponto de vista geomtrico, podemos considerar cada varivel independente x i nas equaes como a i-sima componente de um vetor no R n , dessa forma cada equao do sistema representa um subespao de dimenso n 1 contido no R n . O conjunto soluo formado pela interseco de todos esses subespaos. Por exemplo, no caso de um sistema de duas variveis e duas incgnitas, cada apenas uma soluo se as retas se cruzarem ou ainda innitas solues se as retas se sobreporem. equao representa um reta no R2 e o conjunto soluo ser vazio se as retas forem paralelas, possuir

Captulo 4

InterpolaoNeste captulo estudaremos mtodos que permitem encontrar um valor aproximado para uma funo f calculada em um ponto x do intervalo I, atravs do conhecimento de uma coleo de pares ordenados (pontos) {(xi , f (xi ))}N tais que xi I. Seja P a funo que aproxima f no intervalo I. i=1 Ento, para o conjunto de pontos xi , i = 1, . . . , N

P (xi ) = f (xi ), dizemos que P interpola a funo f nos valores x 1 , x2 , . . ., xN . Ento podemos utilizar a funo P interpolao. Se x estiver fora do intervalo I e ainda assim utilizarmos a funo P para encontrar o valor aproximado de f nesse ponto, o procedimento denominado extrapolao. A escolha de polinmios para realizar essa tarefa possui os seguintes motivos: possvel aproximar uma grande variedade de funes, os polinmios so de fcil manipulao matemtica (principalmente derivao e integrao) e o teorema de Weierstrass Teorema (Weierstrass) Seja f uma funo contnua denida no intervalo fechado limitado [a, b] e seja um nmero positivo. Ento existe um polinmio p, tal que para todo x [a, b], |f (x) p(x)| < . No entanto, da mesma forma que o teorema de Weierstrass garante uma representao de f por um polinmio p to prximo quanto queiramos, ele nada diz sobre o grau de p. Em algumas situaes, o problema de encontrar p que desempenhe esse papel pode ser extraordinariamente do ponto de vista numrico. Antes de discutirmos o procedimento de interpolao por polinmios, vale a pena mencionar um algoritmo til no clculo do valor de p em um ponto x. Trata-se do algoritmo de Horner. Algoritmo de Horner Batizado com o nome do matemtico ingls Willian George Horner mas j conhecido por Isaac Newton em 1669 e mesmo pelo matemtico chins Qin Jiunshao no sc. XIII, o algoritmo consiste para encontrar uma aproximao para o valor de f no ponto x I, esse procedimento denominado

Captulo 4. Interpolao

33

em uma maneira otimizada de calcular p(x) = a m xm + am1 xm1 + . . . + a1 x + a0 atravs de m multiplicaes e m adies. Basta reescrever o polinmio na forma concatenada: p(x) = ((. . . ((am x + am1 )x + ama )x + . . . + a2 )x + a1 )x + a0 . Assim, p(x) pode ser calculado iterativamente: se denominarmos b m = am , am x + am1 = bm x + am1 = bm1 , ento obtemos uma recurso para os b i de modo que p(x) = b0 . Por exemplo, o b1 = b2 x 1 = (3x + 8)x 1 e nalmente p(x) = b0 = b1 x + 1. polinmio p(x) = 3x3 +8x2 x+1 = ((3x+8)x1)+1. Nesse caso b3 = 3, b2 = b3 x+8 = 3x+8,

4.1. Interpolao polinomialSeja fi , i = 1, 2, . . . , n, o valor da funo f calculada nos n pontos de interpolao x i . Encontrar o polinmio de grau m que interpola f nesses pontos consiste em resolver o sistema de equaes lineares fi f (xi ) = p(xi ), ou seja o sistema m1 + . . . + a 1 x1 + a 0 = f 1 am xm + am1 x1 1 am xm + am1 xm1 + . . . + a1 x2 + a0 = f2 2 2 . . . . . . . . . . . . . . . . . . m1 + m + a + a 1 xn + a 0 = f n am xn m1 xn

(4.1)

As m + 1 incgnitas so os coecientes do polinmio, a 0 , a1 , . . . , am e o sistema possui n equaes. Portanto, tipicamente, o sistema no possui soluo se m + 1 < n, possui innitas solues se m + 1 > n e ser unicamente determinado se m + 1 = n. Observao. Vale lembrar que esse procedimento no adequado para determinar polinmios interpolantes se os dados de entrada possurem componentes aleatrios (medidas experimentais), ou de alguma forma, o valor da funo que se quer interpolar no conhecida exatamente ou com suciente preciso.

4.1.1. Interpolao pelos polinmios de Lagrange Como veremos adiante, resolver o sistema (4.1) no a maneira mais simples ou menos sujeita a erros de arredondamento para determinar o polinmio interpolante. O seguinte teorema garante a unicidade do polinmio interpolante, o que nos permite buscar maneiras alternativas de constru-lo. Por ser nico, o resultado ser independente da construo.

Captulo 4. Interpolao

34

Teorema (unicidade do polinmio interpolante) Sejam x1 , . . . , xn , pontos distintos. Para um conjunto arbitrrio de valores f 1 , . . . , fn existe um e somente um polinmio p de grau menor ou igual a n 1 tal que p(xi ) = fi , para i = 1, 2, . . . , n. Vamos supor que para cada 1 j n exista um polinmio de grau n 1, l j (x) tal que para cada

1 k n, o valor de lj no ponto de interpolao xk tal que lj (xk ) = j,k ,

onde j,k o delta de Kronecker1 . Nesse caso, os polinmios lj permitem reescrever o polinmio interpolante p(x):n

p(x) = f1 l1 (x) + f2 l2 (x) + . . . + fn ln (x) =j=1

fj lj (x), = fk . Portanto se formos

podemos trivialmente vericar que p(x k ) = partir das seguintes consideraes.

capazes de construir os polinmios l j a interpolao estar determinada. Vamos ento constru-los a Segundo a sua denio lj (xk ) = 0 para todo xk tal que k = j, ento os pontos xk so razes de lj , se j = k e portanto, a menos de uma constante multiplicativa, C j , o polinmio lj determinado pelo produtrio lj (x) = Cj (x x1 )(x x2 ) . . . (x xj1 )(x xj+1 ) . . . (x xn )n

n j=1 fj lj (xk )

=

n j=1 fj j,k

= Cji=1i=j

(x xi ).

Por m, a constante Cj pode ser determinada atravs da propriedade l j (xj ) = 1:n

lj (xj ) = 1

Cji=1i=j

(xj xi ) = 1,

1

O delta de Kronecker denido pela expresso

j,k = onde j e k so dois nmeros inteiros.

0 1

, ,

j=k , j=k

Captulo 4. Interpolao

35n

ou seja Cj =

i=1i=j

1 . (xj xi )

Dessa forma, os polinmios lj (x), denominados polinmios de Lagrange so determinados a partir do seguinte produtrio lj (x) =i=1i=j

n

x xi xj x i

e a interpolao de Lagrange p(x) =

n

fj lj (x).j=1

Exemplo: Seja a funo f (x) = sen(x) a partir da qual construmos a interpolao nos trs pontos x1 = 0, x2 = 1 e x3 = 2. Ser ento um polinmio de segundo grau. Os pontos de interpolao so dados por j 1 2 3 xj 0 1 2 fj = sen(xj ) 0 sen(1) sen(2)

os polinmios de Lagrange so ento dados por l1 (x) = (x 1)(x 2) x2 3x + 2 = , (0 1)(0 2) 2 (x 0)(x 2) = x2 + 2x (1 0)(1 2) x2 x (x 0)(x 1) = . (2 0)(2 1) 2

l2 (x) = e

l3 (x) = A interpolao dada por p(x) =

sen(2) sen(2) sen(1) x2 + 2sen(1) 2 2

x

4.1.2. Interpolao de Newton De acordo com o teorema da unicidade do polinmio interpolante, todo interpolao de n pontos existem outras maneiras de construir o polinmio p(x) que podem ser mais convenientes. Uma dessas por um polinmio de grau n 1 nica e pode ser obtida pelo mtodo de Lagrange. No entanto,

Captulo 4. Interpolao

36

maneiras a interpolao de Newton, que permite a insero de pontos adicionais de maneira mais simples e menos sujeita deteriorao por erros de arredondamento. O mtodo consiste em determinar o polinmio p(x) = a0 + a1 (x x1 ) + a2 (x x1 )(x x2 ) + . . . + an1 (x x1 ) . . . (x xn1 ). Por construo, p(x1 ) = a0 como p(x) o polinmio interpolante, p(x 1 ) = f1 , ou seja, a0 = f 1 . Da mesma forma, p(x2 ) = a0 + a1 (x2 x1 ) = f2 = f1 + a1 (x2 x1 ) = f2 , ou seja, a1 = f2 a 0 x2 x 1

e assim por diante, os coecientes so determinados recursivamente e o k-simo coeciente determinado em termos dos pontos de interpolao e os coecientes anteriores pela expresso ak = fk+1 a0 k1 j=1 aj (xk+1 k j=1 (xk+1

xj )

x1 ) . . . (xk+1 xj )

.

(4.2)

A frmula de recorrncia (4.2) pode ser convenientemente descrita atravs da notao de diferenas divididas. Seja a funo f [xk , xk+1 , . . . , xl+1 ] denida pela relao de recorrncia f [xk , xk+1 , . . . , xl+1 ] = e Assim, podemos vericar que f [xk , xk+1 ] = f [xk+1 ] f [xk ] xk+1 xk f [xk+1 , xk+1 , . . . , xl+1 ] f [xk , xk+1 , . . . , xl ] xl+1 xk . f [xk ] = fk = f (xk ).

Captulo 4. Interpolao

37 f [xk+1 , xk+2 ] f [xk , xk+1 ] . xk+2 xk

e f [xk , xk+1 , xk+2 ] =

Nessa notao, os coecientes do polinmio so dados por a0 = f [x1 ], a1 = f [x1 , x2 ],

a2 = f [x1 , x2 , x3 ], . . . . . . . . . an1 = f [x1 , x2 , . . . , xn ]. Diagramaticamente, os coecientes so calculados a partir da seqncia de diferenas divididas calculadas recursivamente: x1 x2 x3 . . . . . . f [x1 ] f [x2 ] f [x3 ] . . . . . . f [x1 , x2 ] f [x2 , x3 ] f [x3 , x4 ] . . . . . . f [x1 , x2 , x3 ] f [x2 , x3 , x4 ] f [x3 , x4 , x5 ] . . . . . . f [x1 , . . . , xn1 ] f [x1 , . . . , xn ] ... ... f [x2 , . . . , xn ]

xn2 f [xn2 ] f [xn2 , xn1 ] f [xn2 , xn1 , xn ] xn1 f [xn1 ] xn f [xn ] f [xn1 , xn ]

Exemplo: Vamos realizar a interpolao da funo sen(x) no intervalo x [0, 2] atravs de um polinmio de segundo grau nos pontos x 1 = 0, x2 = 0 e x3 = 0. Neste caso, j 1 2 3 xj 0 1 2 fj = sen(xj ) 0 sen(1) sen(2)

e f [x1 ] = 0, f [x2 ] = sen(1) e f [x3 ] = sen(2). As prximas diferenas divididas so dadas por f [x1 , x2 ] = e f [x2 , x3 ] = f [x2 ] f [x1 ] sen(1) 0 = x2 x 1 10

sen(2) sen(1) f [x3 ] f [x2 ] = . x3 x 2 21

Captulo 4. Interpolao

38

Finalmente, f [x1 , x2 , x3 ] = f [x2 , x3 ] f [x1 x2 ] sen(2) sen(1) + sen(1) = . x3 x 1 20

Portanto, o polinmio interpolante p(x) = f [x1 ] + f [x1 , x2 ] (x x1 ) + f [x1 , x2 , x3 ] (x x1 )(x x2 ) p(x) = sen(1) x + sen(2) sen(1) x(x 1) 2

Exerccio. 1) Inclua o ponto x4 = 1/2 na interpolao anterior e encontre o polinmio interpolante de terceiro grau. 2) Encontre o polinmio interpolante de terceiro grau nos mesmos pontos do exemplo anterior (incluindo o ponto x4 = 1/2) para as funes cos(x), xsen(x) e e x 1. 4.1.3. Erros de truncamento na interpolao por polinmios Seja f uma funo contnua e nvezes diferencivel no intervalo (a, b) que contm os pontos x1 , x2 , . . . , xn e seja p o polinmio de grau n 1 que interpola f nesses pontos. Ento possvel mostrar2 que para cada x (a, b), existe um (x) (a, b) tal que f (x) p(x) = 1 (n) f () n!n i=1

(x xi ).

(4.3)

Poderamos supor que para uma f contnua e sucientemente suave, a seqncia de polinmios interno intervalo (a, b). No entanto, como o exemplo a seguir ilustra, isto nem sempre ocorre. Fenmeno de Runge A seguinte funo, proposta por Carle D. T. Runge ao estudar o comportamento dos erros na interpolao polinomial, f (x) = 1 , 1 + 25x2 x [1, 1] (0.727, 1). polantes {pn }n1 convergiria para f conforme aumentssemos o nmero de pontos de interpolao

tal que a seqncia de polinmios interpolantes {p n }n construdos a partir de pontos de interpolao igualmente espaados no converge 3 para f (x) no intervalo de valores x (1, 0.727)A demonstrao pode ser encontrada nas referncias: Eldn, L.; Wittmeyer-Koch, L. Numerical Analysis (1990), Claudio, D. M.; Marins, J. M. Clculo Numrico Computacional - teoria e prtica 3a ed. (2000). 3 A demonstrao pode ser encontrada na referncia : Isaacson, E. ; Keller, H. Analysis of Numerical Methods (1966).2

Captulo 4. Interpolao

39

Na realidade possvel demonstrar que lim max |f (x) pn (x)| = +.

n+ 1x1

Podemos analisar esse comportamento no regular da interpolao a partir do termon i=1

(x xi )

(4.4)

contido na expresso (4.3). Esse produtrio possui uma utuao para os valores do argumento prximos fronteira do intervalo (1, 1) que progressivamente ampliada conforme aumentamos o nmero de pontos se os mesmos orem igualmente espaados. Os grco seguintes ajudam a ilustrar o comportamento do produtrio (4.4).

6 10 0.0003 4 10 0.0002 2 10 0.0001

-7

-7

-7

x x -1 -0.5 -0.0001 0.5 1 -0.6 -0.4 -0.2 -2 10-7

0.2

0.4

0.6

-0.0002

-4 10

-7

-0.0003

-6 10

-7

Figura 4.1. a) comportamento do produtrio (4.4) com 20 pontos igualmente espaados no intervalo [1, 1]. b) recorte do mesmo produtrio no intervalo [0.7, 0.7].

Esse comportamento pode ser minimizado atravs da escolha de pontos no igualmente espaados. Na realidade possvel demonstrar que a variao do termo (4.4) mnima em valor absoluto quando os pontos xi esto espaados em um intervalo (a, b) segundo a seguinte expresso xi = a+b ab + cos 2 2 2i 1 2n

para i = 1, 2, . . . , n. Esses pontos so denominados pontos de Chebyshev. Utilizando os pontos de Chebyshev no intervalo [1, .1] podemos controlar o comportamento dos n +. polinmios interpolantes para a funo de Runge e garantir a convergncia p n1 (x) f (x) quando

Captulo 4. Interpolao2 10-6

40

1 10

-6

x -1 -0.5 0.5 1

-1 10

-6

-2 10

-6

Figura 4.2. O produtrio (4.4) com 20 pontos de Chebyshev

Ainda assim, existem funes contnuas que requerem um nmero impraticvel de pontos para que a interpolao se aproxime da funo original. Por exemplo, a funo |x| no intervalo [1, 1] requer um polinmio de grau maior que 10 6 para que a interpolao seja exata at 10 3 . em geral, quando utilizamos polinmios de grau maior ou igual a 100, a maior diculdade lidar com os erros de arredondamento.

4.2. Interpolao splineSplines so funes formadas por diferentes polinmios de grau menor ou igual a um m, denidos para cada intervalo entre os pontos de interpolao de modo que em cada ponto de interpolao o spline contnuo,assim como todas as derivadas at ordem m 1.f

sn s1 x s2 x sn2

1

x

x

x x1 x2 x3 xn2

xn

1

xn

Figura 4.3. Interpolao spline

Nas situaes em que o nmero de pontos de interpolao grande (por exemplo, em aplicaes CAD computer-aided design), a inexatido na aproximao obtida com um polinmio de grau ele-

Captulo 4. Interpolao

41

vado dominada pelos erros de arredondamento. Ou ento quando a funo que se quer interpolar possui derivadas de valor numrico elevado em alguma regio do intervalo de interpolao, a aproximao prejudicada em todo o intervalo. Nessas situaes, a interpolao por spline pode auxiliar a tarefa de interpolao. O procedimento de construir splines anlogo qualquer que seja o grau dos polinmios utilizados, como o spline de maior interesse (veremos porque) aquele formado por polinmios de grau 3, nos concentraremos nesse caso apenas.

4.2.1. Interpolao spline cbica Sejam x1 < x2 < . . . < xn os pontos e interpolao. um spline cbico uma funo s(x), denida no intervalo [x1 , xn ] com as seguintes propriedades: 1. s(x), s (x) e s (x) so funes contnuas no intervalo (x 1 , xn ). i = 1, 2, . . . , n. Portanto, s composto por n1 polinmios cbicos, cada polinmio determinado por 4 coecientes (ai ,bi , ci e di ) o que d um total de 4n 4 coecientes a determinar, ou seja 4n 4 incgnitas. Cada polinmio deve satisfazer a condio de continuidade nos pontos de interpolao alm, claro, de interpolar o ponto xi , ou seja, si (xi ) = fi (interpolao), si (xi+1 ) = fi+1 (continuidade de s ), dades de s (x) e s (x): para i = 1, 2, . . . , n 1. As condies acima implicam 2(n 1) equaes. Faltam ainda as continuisi (xi+1 ) = si+1 (xi+1 ) (continuidade de s ), si (xi+1 ) = si+1 (xi+1 ) (continuidade de s ), para i = 1, 2, . . . , n 2. Cada condio equivale a n 2 equaes. Portanto temos at agora um total Essas duas ltimas equaes relacionam-se com as condies de fronteira do spline. Com relao ao comportamento de s(x) no extremo do intervalo, h duas possibilidades a se considerar: i) spline natural, s1 (x1 ) = 0 sn1 (xn ) = 0 de 4n 6 equaes. Restam duas equaes para que seu nmero seja igual ao nmero de incgnitas.

. 2. Em cada subintervalo [xi , xi+1 ], s(x) um polinmio cbico tal que s(x i ) = fi = f (xi ) para

Captulo 4. Interpolao

42

possui esse nome por ser a condio equivalente aproximao por rguas elsticas (uso mais tradicional do spline). ii) spline com mesmas condies de f na extremidade, s1 (x1 ) = f (x1 ) sn1 (xn ) = f (xn ) essa escolha pressupe que a informao sobre o valor da derivada de f nos extremos do intervalo seja conhecida. A aproximao obtida com essa escolha possui uma maior exatido do que a obtida com o spline natural. os coecientes ai ,bi , ci e di dos n 1 polinmios que compe o spline: Nos prximos pargrafos montaremos o sistema de equaes lineares para determinarmos 4n 4 si (x) = ai + bi (x xi ) + ci (x xi )2 + di (x xi )3 . vista da equao (4.5) a interpolao implica fi = si (xi ) = ai para i = 1, 2, . . . , n 1. O que determina o valor dos coecientes a i . para i = 1, 2, . . . , n 2, ou seja, ai + bi (xi+1 xi ) + ci (xi+1 xi )2 + di (xi+1 xi )3 = ai+1 , fi + bi (xi+1 xi ) + ci (xi+1 xi )2 + di (xi+1 xi )3 = fi+1 . (4.6) pode ser reescrita como (4.6) (4.5)

Por ser uma interpolao, a cada xi , temos que s(xi ) = fi , ou seja, si (xi ) = fi . Portanto, em

A continuidade do spline s(x) nos pontos de interpolao implica a equao s i (xi+1 ) = si+1 (xi+1 )

Para aliviar a notao, vamos introduzir a notao h i = (xi+1 xi ). Dessa forma, a equao anterior fi + bi hi + ci h2 + di h3 = fi+1 i i A continuidade na primeira e na segunda derivadas implicam bi + 2ci hi + 3di h2 = bi+1 i e ci + 3di hi = ci+1 para i = 1, 2, . . . , n 2. (4.9) (4.8) (4.7)

Captulo 4. Interpolao

43

Isolando di na equao (4.9) e substituindo em (4.7) e (4.8) encontramos respectivamente fi+1 = fi + bi hi + e bi+1 = bi + hi (ci + ci+1 ) para i = 1, 2, . . . , n 2. (4.11) h2 i (2ci + ci+1 ) 3 (4.10)

Isolando bi na equao (4.10) podemos determin-lo em termos dos valores conhecidos f i , hi e

da incgnita ci (o mesmo acontece com os coecientes d i , a partir da equao (4.9)), bi = para i = 1, 2, . . . , n 2. fi+1 fi hi (2ci + ci+1 ), hi 3 (4.12)

A substituio de bi e bi1 dados pela equao (4.12) na equao (4.11) com os ndices deslocados

de uma unidade, ou seja, bi = bi1 + hi1 (ci1 + ci ), permite encontrar uma equao para os coecientes ci em termos dos valores conhecidos fi e hi : hi1 ci1 + 2(hi1 + hi )ci + hi ci+1 = 3 fi+1 fi hi 3 fi fi1 hi1 , (4.13)

para i = 2, 3, . . . , n1. A equao anterior dene um sistema de equaes lineares para as incgnitas ci . Note que alm dos coecientes c1 , c2 , . . . , cn1 , o sistema envolve um coeciente cn que no est diretamente relacionado a algum dos n 1 polinmios s i . Na realidade, cn est relacionado s condies no extremo do intervalo de interpolao e sua determinao depende do tipo de spline que estamos construindo, se um spline natural ou um spline que satisfaz as mesmas condies de f nos extremos do intervalo de interpolao. As n2 equaes (4.13) envolvem n variveis (as incgnitas c i ), para que o sistema (tipicamente) tenha soluo nica devemos incluir as duas ltimas equaes que descrevem o comportamento do spline nos extremos do intervalo de interpolao. Vamos estudar inicialmente o caso do spline natural.

Spline natural O spline natural deve satisfazer as condies s (x1 ) = 0 e s (xn ) = 0, estas duas equaes implicam respectivamente c1 = 0 e 2cn1 + 6dn1 hn1 = 0. A equao (4.14) implica em termos da equao para os coecientes d i (4.9) que cn = 0. (4.14)

Captulo 4. Interpolao

44

Colecionando esses resultados temos ento a seguinte situao: resolvendo o sistema de equaes c1 = 0 h c + 2(hi1 + hi )ci + hi ci+1 = 3 i1 i1 c = 0 n

fi+1 fi hi

3

fi fi1 hi1

,

i = 2, . . . , n 1

encontramos o valor dos coecientes c i . A partir desses coecientes determinamos o valor dos coecientes bi atravs das equaes (4.12) bi = fi+1 fi hi (2ci + ci+1 ), hi 3

para i = 1, 2, . . . , n 1; e o valor dos coecientes d i atravs da equaes di = ci+1 ci , 3hi (4.15)

para i = 1, 2, . . . , n 1, obtida a partir de (4.9). Os coecientes a i = fi como j havamos determinado anteriormente.

Spline com as mesmas condies de f nos extremos Nesse caso o spline deve satisfazer as condies s (x1 ) = f (x1 ) f1 e s (xn ) = f (xn ) (4.16) (4.17)

vamente e

fn .Para determinar o spline, f1 e fn devem ser valores conhecidos. As condies implicam respectib1 = f 1 bn1 + 2cn1 hn1 + 3dn1 h2 = fn . n1 Como os coecientes bi satisfazem a equao (4.12), a equao (4.16) implica f1 = b 1 = f2 f 1 h1 (2c1 + c2 ), h1 3 f2 f 1 h1

ou seja, 2h1 c1 + h1 c2 = 3 3f1 . (4.18)

Da mesma forma, no caso da equao (4.17), as equaes (4.12) e (4.15) implicam hn1 cn1 + 2hn1 cn = 3 fn fn1 hn1 + 3fn . (4.19)

Captulo 4. Interpolao

45

Em resumo, devemos resolver o sistema formado pelas equaes (4.13), (4.18) e (4.19) f 2h1 c1 + h1 c2 = 3 f2h1 1 3f1 f f f fi1 h c + 2(hi1 + hi )ci + hi ci+1 = 3 i+1i i 3 ihi1 , h i1 i1 fn1 hn1 cn1 + 2hn1 cn + 3fn = 3 fnhn1 bi = fi+1 fi hi (2ci + ci+1 ), hi 3 di = ci+1 ci , 3hi

i = 2, . . . , n 1

e ento determinar os coecientes b i e di atravs das equaes (4.12) e (4.15):

para i = 1, 2, . . . , n 1. Naturalmente, os coecientes a i = fi .

Captulo 5

Ajuste de mnimos quadrados5.1. Ajuste de mnimos quadrados polinomialNo captulo anterior estudamos como encontrar um polinmio de grau m que interpola um conentanto, ainda assim podemos encontrar uma aproximao p(x) para a funo desconhecida que dejunto de n pontos {{xi , fi }}n . Tipicamente1 quando m < n 1 esse polinmio no existe. No i=1

ne os n pontos. Se p(x) fosse uma interpolao, seria necessrio que p(x i ) = fi . A tcnica de ajuste de funes substitui essa exigncia pela minimizao de uma funo Q({r i }) sobre os resduos . ri = p(xi ) fi . A escolha mais comum para Q a soma quadrtica dos resduos, ou sejan

Q({ri }) =

2 ri . i=1 i |ri |,

(5.1) a importncia

Existem outras escolhas possveis para Q, por exemplo, Q = max i |ri | ou Q =

de (5.1) deve-se ao teorema de Gauss-Markov para modelos lineares em estatstica.

Inicialmente vamos tratar o caso em que a funo ajuste um polinmio de grau m. Como veremos a seguir, a condio de que (5.1) seja mnimo fornecer as equaes para determinarmos os coecientes do polinmio. Esse procedimento de encontrar os coecientes que minimizam (5.1) denominado ajuste de mnimos quadrados. Seja ento p(x) o polinmio p(x) = a0 + a1 x + . . . + am xm , uma condio necessria sobre o valor dos coecientes a i para quen

Q=i=1

m j=0

fi

aj xj i

2

claro que existem casos em que a interpolao possvel, por exemplo, poderamos ter um conjunto de trs pontos alinhados (n = 3). Nesse caso, h a interpolao dos mesmos por uma reta (m = 1), mas se fossem trs pontos distintos no alinhados, a interpolao por uma reta no mais seria possvel.

1

Captulo 5. Ajuste de mnimos quadrados

47 Q =0 ak

seja mnimo

para k = 0, 1, . . . , m. Isto implica k = 0, 1, . . . , mn

2 pois

i=1

m j=0

fi

aj xj xk = 0, i i

(5.2)

j m a0 + a1 xi + . . . + ak xk + . . . + am xm = xk . A equao (5.2) = i i i j=0 aj xi ak ak pode ser reescrita comon m j+k aj xi i=1 j=0 m j=0 n j+k xi i=1 n

=i=1 n

fi xk i (5.3)

aj =i=1

fi xk i

para cada k = 0, 1, . . . , m. Ou seja, um sistema de m + 1 equaes lineares cujas incgnitas so os m + 1 coecientes ai do polinmio p(x). O sistema de equaes (5.3) pode ser convenientemente descrito em notao matricial como X T Xa = X T f , onde (5.4)

f1 a1 . XT . e f = . . As equaes dadas pelo sistema (5.4) so . . . fn am denominadas equaes normais. Essa nomenclatura deve-se ao seguinte fato: o sistema (5.4) pode sua transposta, a = ser escrito como X T (Xa f ) = 0 o vetor entre parnteses, Xaf o vetor cujas componentes so dadas pelos resduos da aproximao

X = a0

1 x1 1 x2 . . . . . .

x2 1 x2 2 . . .

xm 1 xm 2 . . . xm n

.. .

1 x n x2 n

,

(5.5)

Captulo 5. Ajuste de mnimos quadrados

48

e, segundo a equao anterior esse vetor normal (ortogonal) aos vetores formados pelos elementos l x1 l x2 das linhas da matriz X T que so da forma . para l = 0, 1, 2, . . . , m. . . xl n

vamos determinar o ajuste de mnimos quadrados para um polinmio de segundo grau p(x) = a0 + a 1 x + a 2 x2 . Os coecientes do polinmio so a soluo do seguinte sistema na representao matricial X T Xa = X T f , onde a matriz X dada por 1 x 1 x2 1 1 2 4

Exemplo: Seja o conjunto de pontos {{x i , fi }}5 : {{2, 0}, {1, 1}, {0, 2}, {1, 1}, {2, 0}} i=1

portanto

1 x2 X = 1 x3 1 x4 1 x5

x2 2 x2 = 3 x2 4 x2 5 5 10 0 10 0

1 1 1 1 0 0 , 1 1 1 1 2 4 10

O vetor de constantes f = f2 e o vetor de incgnitas a = a1 compe o sistema que a2 f3 possui matriz completa 5 0 10 4 0 0 10 0 10 0 34 2 e soluo 58/35 0 3/7

XT X = 0 f1

0 . 34

ao

Assim, o polinmio que ajusta os dados p(x) =

a=

58 3 2 x . 35 7

.

Captulo 5. Ajuste de mnimos quadradosf x 2

49

1

x 2 1 1 2

Figura 5.1. Ajuste de mnimos quadrados

Caso particular: regresso linear Quando o ajuste de mnimos quadrados realizado para um polinmio de 1 o grau as expresses so mais simples, em particular, possvel determinar o valor dos coecientes a 0 e a1 exclusivamente em termos de mdias. Para um conjunto de n pontos {{xi , fi }}n , a soma quadrtica dos desvios, Q, dada por i=1n

Q=i=1

(fi a0 a1 x)2 .

Da mesma forma que no caso mais geral, condies necessrias para que Q seja mnimo so dadas pelas derivadas de Q com relao a a0 e a1 que devem se anular: Q =2 a0 e Q =2 a1 A partir da equao (5.6) temosn n n i=1 n i=1

(fi a0 a1 xi )(1) = 0

(5.6)

(fi a0 a1 xi )(xi ) = 0.

(5.7)

na0 + a1i=1

xi =i=1

fi

a0 + a1

1 n

n

xi =i=1

1 n

n

fi .i=1

Da mesma forma, a equao (5.7) implican n n

a0i=1

xi + a 1i=1

x2 i

=i=1

fi xi

1 a0 n

n i=1

1 xi + a 1 n

n

x2 ii=1

1 = n

n

fi xi .i=1

Captulo 5. Ajuste de mnimos quadrados

501 n n i=1

. Atravs da notao x =

reescrever as equaes (5.6) e (5.7) como o sistema

1 n

n i=1 xi ,

x2 =

. x2 , f = i

1 n

n i=1 fi

. e fx =

1 n

n i=1 fi xi

podemos

a0 + x a 1 x a0 + x2

=

f

a1 = f x

cuja soluo determina os coecientes a 0 e a1 : a0 = a1 = fx f x ,

x2 (x)2 x2 (x)2

f x2 x f x

.

5.2. Ajuste de mnimos quadrados por uma combinao linear de funesA maneira atravs da qual determinamos o ajuste de mnimos quadrados por um polinmio pode ser imediatamente generalizada para o caso da combinao linear de um conjunto de funes conhecidas {i (x)}i . Podemos vericar prontamente que o polinmio p(x) = a 0 + a1 x + . . . + am xm um caso particular da combinao linear de m distintas funes j (x)m

(x) =j=1

aj j (x),

no caso em que j (x) = xj . A soma quadrtica dos resduos Q(f, ) dada por . Q(f, ) =n i=1 n

(fi (xi ))2 =

i=1

m j=1

De maneira totalmente anloga ao caso do ajuste para polinmios, o conjunto de coecientes que minimizam Q(f, ) determinado pela soluo do sistema de equaes lineares Q a0 Q a1 Q am

fi

aj j (xi ) .

2

= 0 = 0 . . . = 0

que equivalente ao sistema de equaes normais T a = T f , (5.8)

Captulo 5. Ajuste de mnimos quadrados

51

onde

1 (x2 ) 2 (x2 ) = . . .. . . . . . 1 (xn ) 2 (xn ) a= a1 a2 ef = . . . am f1 f2 . . . fn .

1 (x1 )

2 (x1 )

m (x1 ) m (x2 ) . . . m (xn )

,

nm

Problemas de condicionamento no mtodo de ajuste de funes pelo mtodo dos mnimos quadrados De modo geral, ao aplicarmos o mtodo de ajuste de mnimos quadrados com um polinmio de grau maior ou igual a 8, a tarefa de resolver o sistema de equaes normais (5.4) muito dicultada por erros de arredondamento. A diculdade est relacionada ao fato de que as matrizes da forma X T X presentes no sistema (5.4) so mal condicionadas. Essa propriedade independe dos valores f i no conjunto de n pontos {{xi , fi }}n para ajuste. Note que a matriz (m + 1) (m + 1) X T X i=1 depende apenas dos valores xi e do grau m do p