UNIVERSIDADE FEDERAL DO RIO GRANDE DO SULINSTITUTO DE MATEMATICADEPARTAMENTO DE MATEMATICA PURA E APLICADAIntrodu caoaoCalculoNumerico2a. Edi caoAlvaro Luiz de BortoliCarolina CardosoMaria Paula Gon calves FachinRudnei Dias da CunhaPorto Alegre, abril de 2003.Alvaro Luiz de Bortoli e professor adjunto da UFRGS,desempenhando suas atividades junto ao Departamento deMatematica Pura e Aplicada, do Instituto de Matematica, desde1996.E formado emEngenharia Mec anica pelaUFRGS(1987);Mestre em Engenharia Mecanica pela UFSC (1990); e Doutor emAerodin amicaeAeroelasticidadepelaUFSCeDeutschesLuft-und Raumfahrt Institut, Alemanha (1995).Carolina Cardoso e Bacharel em Matematica (enfaseMatematicaAplicadaeComputacional) pelaUFRGS(1998) eMestre em Matematica Aplicada pela UFRGS (2001).Maria Paula Goncalves Fachin e professora adjunta daUFRGS, desempenhando suas atividades junto ao Departamentode MatematicaPurae Aplicada, doInstitutode Matematica,desde1990.ELicenciadaemMatematicapelaUFRGS(1986);Mestre em Matematica pela UFRGS (1990); e Doctor ofPhilosophy in Computer Science pela University of Kent atCanterbury, Reino Unido (1994).Rudnei Dias da Cunha e professor adjunto da UFRGS,desempenhando suas atividades junto ao Departamento deMatematica Pura e Aplicada, do Instituto de Matematica, desde1994.E Bacharel em Ciencias de Computa cao pela UFRGS (1988)eDoctorof PhilosophyinComputerSciencepelaUniversityofKent at Canterbury, ReinoUnido(1992). Exerceuas fun coesde programador de computadores e analista de sistemas noCentro de Processamento de Dados da UFRGS (1983-1994). Foicoordenador do Programa de Pos-Gradua cao emMatematicaAplicada da UFRGS (1999-2000) e atualmente ocupa o cargo deVice-Diretor do Instituto de Matem atica da UFRGS.Sumario1 Aritmetica noComputador 61.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 Representa cao em bin ario e decimal . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.1 Bits, bytese palavras . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.2 Convers ao entre representa coes . . . . . . . . . . . . . . . . . . . . . . . . . 71.3 Representa cao de n umeros em um computador . . . . . . . . . . . . . . . . . . . . 101.3.1 Representacao de n umeros inteiros . . . . . . . . . . . . . . . . . . . . . . . 111.3.2 Representacao de n umeros reais . . . . . . . . . . . . . . . . . . . . . . . . . 121.3.2.1 Representa cao racional de n umeros reais . . . . . . . . . . . . . . 121.3.2.2 Representa cao de n umeros reais em ponto-xo . . . . . . . . . . . 131.3.2.3 Representa cao de n umeros reais em ponto-utuante . . . . . . . . 131.3.2.4 Tratamento do zero . . . . . . . . . . . . . . . . . . . . . . . . . . 141.3.3 Caracteriza cao de uma representa cao . . . . . . . . . . . . . . . . . . . . . . 141.3.4 Arredondamentos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161.3.5 Opera coes aritmeticas de ponto-utuante . . . . . . . . . . . . . . . . . . . 201.3.5.1 Erros em opera coes aritmeticas de ponto-utuante. . . . . . . . . 201.4 Perda de dgitos signicativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221.4.1 Subtra cao de valores quase identicos . . . . . . . . . . . . . . . . . . . . . . 221.4.2 Teorema sobre a perda de precisao . . . . . . . . . . . . . . . . . . . . . . . 241.5 Condicionamento de um problema . . . . . . . . . . . . . . . . . . . . . . . . . . . 251.6 Computa coes estaveis e instaveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261.7 Desastres causados por erros aritmeticos no computador . . . . . . . . . . . . . . . 271.7.1 Falha do sistema de msseis Patriot . . . . . . . . . . . . . . . . . . . . . 271.7.2 Explos ao do foguete Ariane 5 . . . . . . . . . . . . . . . . . . . . . . . . . . 281.8 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282 CalculodeRazesdeFun coesNao-Lineares 292.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.2 Metodo da Bissec cao. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.3 Metodo da posi cao falsa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332.3.1 Melhorando o metodo da posi cao falsa. . . . . . . . . . . . . . . . . . . . . 362.3.2 An alise do erro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.4 Metodo de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.4.1 An alise do erro . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 402.5 Deriva cao numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 412.5.1 O metodo de Newton-Raphson e as razes complexas def(x) . . . . . . . . 442.6 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4423 CalculodeRazesdePolinomios 453.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.2 Resultados teoricos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.3 Enumera cao e localizacao de razes de polin omios. . . . . . . . . . . . . . . . . . . 463.3.1 Regra de Descartes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.3.2 Regra de Du Gua . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.3.3 Regra da lacuna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.3.4 Cota de Laguerre-Thibault . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.5 Cota de Fujiwara. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.6 Cota de Kojima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.3.7 Cota de Cauchy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 483.4 Metodo de Newton-Viete . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 493.5 Metodo de Horner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.5.1 C alculo do quociente e do resto. . . . . . . . . . . . . . . . . . . . . . . . . 513.5.2 Dea cao de um polin omio. . . . . . . . . . . . . . . . . . . . . . . . . . . . 523.5.3 Calcular a expans ao de Taylor de um polin omio . . . . . . . . . . . . . . . . 523.5.3.1 O metodo de Horner e sua rela cao com a derivada dep(z) . . . . . 533.5.3.2 O metodo de Newton-Raphson usado em conjunto com o algoritmoparcial de Horner . . . . . . . . . . . . . . . . . . . . . . . . . . . 543.6 Razes complexas de equacoes polinomiais . . . . . . . . . . . . . . . . . . . . . . . 563.6.1 Metodo de Bairstow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.7 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 604 Resolu cao deSistemas deEquacoesLineares 614.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 614.2 Resolucao de Sistemas Triangulares de Equa coes Lineares . . . . . . . . . . . . . . 624.3 Resolucao de Sistemas de Equa coes Lineares por Eliminacao Gaussiana . . . . . . 644.3.1 Diculdades . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.3.2 Elimina cao Gaussiana e a Fatora cao LU. . . . . . . . . . . . . . . . . . . . 674.3.3 O Custo Computacional da Fatora cao LU. . . . . . . . . . . . . . . . . . . 694.3.4 Resolucao de sistemas com m ultiplos termos independentes . . . . . . . . . 704.3.4.1 C alculo da inversa de uma matriz . . . . . . . . . . . . . . . . . . 704.4 Resolucao Iterativa de Sistemas de Equa coes Lineares . . . . . . . . . . . . . . . . 734.4.1 Normas de vetores e de matrizes . . . . . . . . . . . . . . . . . . . . . . . . 744.4.2 Normas de matrizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.4.3 N umero de condi cao de uma matriz . . . . . . . . . . . . . . . . . . . . . . 754.4.4 Erros computacionais e condicionamento . . . . . . . . . . . . . . . . . . . . 764.4.5 Metodos iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.4.6 Renamento iterativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.4.7 Metodo iterativo de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.4.8 Metodo iterativo de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . 824.4.9 Extrapola cao de um metodo iterativo . . . . . . . . . . . . . . . . . . . . . 854.5 Metodo do Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 864.5.1 Forma Quadr atica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 874.5.2 Descri cao do metodo do Gradiente . . . . . . . . . . . . . . . . . . . . . . . 904.6 Metodo das Dire coes-Conjugadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 944.7 Metodo dos Gradientes-Conjugados . . . . . . . . . . . . . . . . . . . . . . . . . . . 984.8 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1005 Resolu cao deSistemas deEquacoesNao-Lineares 1025.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.2 Metodo de Newton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1025.3 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10936 AutovaloreseAutovetores 1106.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1106.2 Teoremas de limites sobre autovalores . . . . . . . . . . . . . . . . . . . . . . . . . 1136.3 Calculo de autovalores e autovetores via determinantes. . . . . . . . . . . . . . . . 1156.4 Autovalores de uma matriz tridiagonal simetrica . . . . . . . . . . . . . . . . . . . 1166.5 Metodos para aproxima cao de autovalores e autovetores . . . . . . . . . . . . . . . 1206.5.1 Metodo da potencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1206.5.2 O metodo da potencia com transla cao da origem. . . . . . . . . . . . . . . 1236.5.3 Metodo da itera cao inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . 1246.5.4 O metodo da itera cao inversa e o quociente de Rayleigh . . . . . . . . . . . 1266.6 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1267 Interpolacao 1287.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1287.2 Interpola cao polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1297.3 Forma de Newton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1317.4 Forma de Lagrange. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1327.5 Forma de Newton com diferen cas divididas . . . . . . . . . . . . . . . . . . . . . . 1347.6 Forma de Newton com diferen cas simples . . . . . . . . . . . . . . . . . . . . . . . 1377.7 Interpola cao inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1387.8 Interpola cao por splines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1397.9 Estudo do erro na interpola cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1427.9.1 Estimativa para o erro. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1437.10Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1448 Ajustededadosexperimentais 1468.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1468.2 Mnimos quadrados - domnio discreto . . . . . . . . . . . . . . . . . . . . . . . . . 1488.3 Ajuste linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1488.4 Ajuste polinomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1498.5 Ajustamento por fun coes nao lineares nos par ametros lineariza cao . . . . . . . . 1508.5.1 Ajustamento por uma fun cao exponencial . . . . . . . . . . . . . . . . . . . 1508.5.2 Ajustamento por uma fun cao potencia. . . . . . . . . . . . . . . . . . . . . 1518.5.3 Ajustamento por uma fun cao hiperb olica . . . . . . . . . . . . . . . . . . . 1518.5.4 Ajustamento por uma fun cao do tipoy =xa0+a1 x. . . . . . . . . . . . . . . 1518.5.5 Ajustamento por uma fun cao do tipoy =1a0+a1 x+a2 x2. . . . . . . . . . . . 1518.5.6 Ajustamento por uma fun cao do tipoy = a eb x+c x2. . . . . . . . . . . . . 1518.6 Escolha do melhor ajuste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1518.7 Mnimos quadrados - domnio contnuo . . . . . . . . . . . . . . . . . . . . . . . . . 1548.7.1 Polin omios ortogonais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1578.8 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1599 Integra caoNumerica 1619.1 Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1619.2 Integra cao numerica via interpola cao polinomial . . . . . . . . . . . . . . . . . . . 1619.2.1 Regra do Trapezio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1629.2.2 Metodo dos Coecientes a Determinar . . . . . . . . . . . . . . . . . . . . . 1659.2.3 Regra de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1669.2.4 Regra de Simpson com exatid ao crescente . . . . . . . . . . . . . . . . . . . 1679.2.5 Mudan ca do intervalo de integra cao . . . . . . . . . . . . . . . . . . . . . . 1689.2.6 Quadratura Gaussiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1699.3 Integra cao de fun coes mal comportadas . . . . . . . . . . . . . . . . . . . . . . . . 1739.4 Intervalos de integra cao innitos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1749.5 Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174410 Solu caoNumericadeEqua coesDiferenciaisOrdinarias 17810.1Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17810.2Problema de Valor Inicial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18010.2.1 Existencia da Solucao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18110.2.2 Erros na solucao numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18110.2.3 Metodo da Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18110.2.3.1 Vantagens e desvantagens. . . . . . . . . . . . . . . . . . . . . . . 18310.2.4 Metodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18310.2.5 Metodo de Heum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18410.2.5.1 Erro de truncamento para o metodo de Heum . . . . . . . . . . . 18510.2.6 Metodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18610.2.6.1 Metodo modicado de Euler . . . . . . . . . . . . . . . . . . . . . 18710.2.6.2 Metodo de Runge-Kutta de 4aOrdem. . . . . . . . . . . . . . . . 18710.2.6.3 Erros do metodo de Runge-Kutta . . . . . . . . . . . . . . . . . . 18710.2.6.4 Avalia cao da Fun cao versus Ordem do Metodo Runge-Kutta . . . 18710.2.6.5 Metodo Adaptativo de Runge-Kutta-Fehlberg . . . . . . . . . . . 18810.2.7 Metodos de passo m ultiplo . . . . . . . . . . . . . . . . . . . . . . . . . . . 18910.2.7.1 Convergencia, Estabilidade e Consistencia . . . . . . . . . . . . . . 19210.2.7.2 Erros de truncamento . . . . . . . . . . . . . . . . . . . . . . . . . 19310.2.7.3 Erros de truncamento globais . . . . . . . . . . . . . . . . . . . . . 19310.2.8 Sistemas de Equa coes Diferenciais Ordin arias . . . . . . . . . . . . . . . . . 19510.2.8.1 Metodo da Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . 19610.2.8.2 Metodo de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . 19610.2.9 Solu cao via decomposi cao em autovalores e autovetores . . . . . . . . . . . 19710.2.9.1 O expoente de uma matriz . . . . . . . . . . . . . . . . . . . . . . 19810.2.10Equa coes rgidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19910.3Problemas de Valor de Fronteira . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20010.3.1 Metodo do disparo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20110.3.2 Metodo de Newton. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20310.3.3 Metodo da colocacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20310.3.4 Deriva cao numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20510.3.5 Solu cao por diferen cas-nitas . . . . . . . . . . . . . . . . . . . . . . . . . . 20810.3.5.1 O caso linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20810.4Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20911 Solu caoNumericadeEqua coesDiferenciaisParciais 21211.1Introdu cao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21211.2Equa coes parab olicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21311.2.1 Metodo explcito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21311.2.2 Metodo de Crank-Nicolson . . . . . . . . . . . . . . . . . . . . . . . . . . . 21611.2.2.1 Aproxima cao ponderada . . . . . . . . . . . . . . . . . . . . . . . . 21711.2.3 Condi coes de fronteira . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21811.3Equa coes diferenciais parciais elpticas . . . . . . . . . . . . . . . . . . . . . . . . . 21911.4Exerccios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2225Captulo1AritmeticanoComputador1.1 Introdu caoHojeemdia, oscomputadoresutilizamumsistemadenumera caoembase2, emsuagrandemaioria. Essesistemaechamadodebin arioeutilizaosalgarismos0e1pararepresentarosn umeros (apesar de que quaisquer outros dois smbolos poderiam ser usados).A nossa sociedade, ao contr ario, utiliza um sistema de numeracao decimal, ou base 10; muitoprovavelmente, pelo fato dos seres humanos terem dez dedos nas maos, os quais eram utilizados como uma crianca os utiliza para contar quantidades. A palavra dgito, sin onimo para algarismo,vem do latim digitus, dedo.Ambos os sistemas citados decimal e bin ario s ao sistemas posicionais, i.e., os n umeros s aoformados por somas de potencias, convenientemente multiplicadas pelos algarismos.Por exemplo,o n umero(420, 325)10 = 4 102+ 2 101+ 0 100+ 3 101+ 2 101+ 5 103(1.1)e representado no sistema decimal como a soma das potencias de 10 mostradas acima.Aprincipal caracterstica deumsistemadenumeracaoposicional eanecessidadedarepre-sentacaodozeroporumsmbolo. Aparentemente, ozeroj aerautilizadopelosmaiasepelosbabil onios,esses porvolta de300 A.C.Onosso sistema denumeracao decimalfoiinventado naIndiaporvoltadoano600D.C. e, tendosidousadopormuitosseculospelospovos arabesnoOrienteMedio, foiintroduzidonaEuropaduranteasinvas oesmouras, noperodoentre1200e1600 (da o nome de algarismos ar abicos).Vejacomoozeroeimportantenumsistemaposicional, comparadocomumsistemanao-posicional, como o romano, por exemplo. Nesse ultimo, o n umero 401 e representado como CCCCI(os romanosn aoutilizavamas abrevia coescomoIVpararepresentar4). Porem, nosistemadecimal, o 0 e necessario para distinguir 401 de 41 ele efetivamente serve como um espa cadordos algarismos, em termos das potencias de 10.Einteressantenotarqueosistemadecimal erautilizadopararepresentarapenasn umerosinteiros, e n ao fra coes decimais, ate o seculo XVII. Em pases de lngua inglesa, ate hoje persisteousodefra coes inteirascomo1/4,1/8,1/16,3/4,comoporexemploemplacasdesinaliza caorodovi aria e na especicacao dos di ametros de ferramentas.Como dissemos ao iniciarmos esse captulo, os computadores utilizam normalmente um sistemade numeracao bin ario, ou base 2. Esse sistema nao e, no entanto, t ao recente quanto os computado-res; na verdade, j a era usado como base para um algoritmo de multiplica cao no Papiro Matem aticode Rhind, escrito ha 4,000 anos atr as [12, p ag. 7].6Introdu c ao ao Calculo Numerico Aritmetica no Computador1.2 RepresentacaoembinarioedecimalGenericamente, podemosdizerqueumsistemadenumeracaonumabaseadmiteapenasosdgitos 0, 1, . . ., 1. Assim, o n umero (1001, 11101)2 representa o n umero (9, 90625)10, onde ossubscritos indicam a base do sistema de numera cao utilizado:(1001, 11101)2= 1 23+ 0 22+ 0 21+ 1 20+1 21+ 1 22+ 1 23+ 0 24+ 1 25= 8 + 0 + 0 + 1 + 0, 5 + 0, 25 + 0, 125 + 0 + 0, 03125= (9, 90625)10Assim como existem n umeros reais, em decimal, que tem parte fracion aria com n umero innitodedgitos ditosirracionais tambemexistemn umerosreais embin ariocomamesmacarac-terstica. Mais ainda, existem n umeros reais, em decimal, cuja parte fracion aria tem um n umeronito dedgitos, para os quais a sua representa cao em bin ario apresenta um n umero innitodedgitos. Por exemplo, o n umero 1/10 n ao tem uma representacao bin aria nita:110= (0, 0001100110011 . . .)2 =116 +132 +064 +0128 +1256 +1512 +01024 + . . .mas o conjunto de dgitos 0011 repete-se.1.2.1 Bits, bytesepalavrasAmemoriadeumcomputador podeser descritacomoumconjuntodepalavras. Amaioriadoscomputadorestemsuamemoriaestruturadadetal formaquecadaacessodeleituraouescritaefeitoemtermosdeumaoumaispalavras, asquaissaoacessadasporumendere co unico. Tipicamente,uma palavra e composta por 32 bits;processadores de ultima gera cao paramicrocomputadores pessoais com palavras de 64 bitsj a sao uma realidade hoje.Um bit(contracao em ingles de binary digit) e a menor unidade de informa cao armazenadaem um computador, podendo representar os valores 0 ou 1. Um conjunto de 8 bits e chamado debyte; nele, podemos armazenar 28= 256 diferentes valores inteiros, atraves das combinacoes de 0e 1 entre os diferentes bits.1.2.2 ConversaoentrerepresentacoesEconvenientesabercomoconverterumn umerodecimal parasuarepresentacaoembin arioevice-versa. Emalgumasaplica coes envolvendo ousodecomputadores, enecessario sabercomoconverter para decimal um valor armazenado de forma bin aria.Podemos efetuar a conversao de um n umero real decimal x para bin ario convertendo separa-damenteaspartesinteira efracion aria dexip(x)efp(x)edepoisjustaporarepresenta caobin aria dessas duas partes, separando-as por um ponto. Nos algoritmos apresentados a seguir, aconversao de decimal para bin ario resulta em um string de caracteres 0 e 1, e nao num n umeroformado pelos mesmos algarismos, pois nao sao representa coes equivalentes.Suponhaent ao o n umerox = (401, 640625)10. Arepresenta cao bin aria deip(x)= (401)10eobtida, inicialmente, dividindo-se 401 por 2; essa divis ao devolve um quociente e um resto. O restoe, necessariamente, 0 ou 1. Ap os, divide-se esse quociente por 2, obtendo-se um outro quocientee resto. Esse processo e repetido ate que o quociente seja 1;a representa cao bin aria e formada,ent ao, pelo ultimoquocienteepelosrestos, tomadosnaordeminversa aqueforamobtidos. Atabela 1.1 mostra esse processo. A representacao em bin ario de ip(x) = (401)10 e, portanto,(401)10= (110010001)2 == 1 28+ 1 27+ 0 26+ 0 25+ 1 24+ 023+ 0 22+ 0 21+ 1 20== 256 + 128 + 16 + 1 = 401A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 7Introdu c ao ao Calculo Numerico Aritmetica no Computadordividendo quociente resto401 200 1200 100 0100 50 050 25 025 12 112 6 06 3 03 1 1Tabela1.1: Processodeconversaoparabinariodaparteinteirade(401, 640625)10; osdgitossublinhados comporao a representacao binaria.Para convertermos a parte fracion aria fp(x) = (0, 640625)10, fazemos um processo de multipli-cac oessucessivaspor2. Inicialmente,multiplicamos 0, 640625 por 2, resultando em 1, 28125. Odgito ` a esquerda do ponto decimal ser a um dos dgitos da representa cao bin aria de fp(x);comoesse dgito e igual a 1, subtramos 1 do n umero, resultando em 0, 28125. Esse n umero e, novamente,multiplicado por 2, resultando em 0, 5625; como o dgito ` a esquerda do ponto decimal e o 0, bastamultiplicar novamente esse n umero por 2. O processo continua ate que o n umero multiplicandoseja igual a 1, 0; os dgitos 0 e 1, ` a esquerda do ponto decimal, formam a representa cao bin aria defp(x), agora na mesmaordem em que foram obtidos, conforme mostrado na tabela 1.2multiplicando resultado0, 640625 1, 281250, 28125 0, 56250, 5625 1, 1250, 125 0, 250, 25 0, 50, 5 1, 0Tabela 1.2: Processo de conversao para binario da parte fracionaria de (401, 640625)10; os dgitossublinhados comporao a representacao binaria.Logo, a representa cao em bin ario de fp(x) = 0, 640625 e(0, 640625)10= (0, 101001)2 == 1 21+ 0 22+ 1 23+ 0 24+ 0 25+ 1 26== 0, 5 + 0, 125 + 0, 15625 = 0, 640625e, portanto, podemosjustaporasrepresenta coesbin ariasdeip(401, 640625)=(110010001)2efp(401, 640625) =(0, 101001)2, obtendoarepresenta caobin ariade (401, 640625)10, aqual e(110010001, 101001)2.Osalgoritmosapresentadosaseguirsumarizamosprocessosdeconvers aodedecimal parabin ario e de bin ario para decimal.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 8Introdu c ao ao Calculo Numerico Aritmetica no ComputadorAlgoritmo1.2.1Convers ao decimal para binarioproc conv dec para bin(input:x; output:b)%x eh um numero real decimal em modulo eb eh a sua% representacao binaria, armazenada como um% string de caracteres.x [x[b convdecparabinip(ip(x))[convdecparabinfp(fp(x))endprocproc conv dec para bin ip(input:x; output:b)%x eh um numero inteiro eb eh a sua representacao% binaria, armazenada como um string de caracteres.d x/2|r x 2dx db num2str(r)whilex 1d x/2|r x 2dx db num2str(r)[bendwhileendprocproc conv dec para bin fp(input:x; output:b)%x eh um numero fracionario menor do que 1 e%b eh a sua representacao binaria, armazenada% como um string de caracteres.x [x[b

.

whilex < 1d x 2i d|ifd > 1 thenx d 1elsex dendifb b[num2str(i)endwhileendprocA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 9Introdu c ao ao Calculo Numerico Aritmetica no ComputadorAlgoritmo1.2.2Convers ao binario para decimalproc conv bin para dec(input:b; output:x)%b eh um numero binario, armazenado como% um string de caracteres ex eh a sua% representacao em decimal.l length(b)p ndstr(b,

.

) % Localiza onde estah o ponto decimalembifp ,= 0 then %b eh um numero binario com parte fracionariax convbinparadecip(substr(b, 1, p 1))+convbinparadecfp(substr(b, p, l))else %b eh um numero inteirox convbinparadecip(b)endprocproc conv bin para dec ip(input:b; output:x)%b eh um numero binario inteiro, armazenado como% um string de caracteres ex eh a sua% representacao em decimal.l length(b)x 0k lfori = 1, 2, . . . , lx x + str2num(b[i]) 2kk k 1endforendprocproc conv bin para dec fp(input:b; output:x)%b eh um numero binario menor do que 1, armazenado% como um string de caracteres (com . aa frente) e%x eh a sua representacao em decimal.l length(b)x 0k 1fori = 2, 3, . . . , l % Desconsidera o 1o. caracter (.)x x + str2num(b[i]) 2kk k + 1endforendprocNote,noentanto,que sepor umlado umn umero inteiro decimalpodeser facilmente repre-sentadoporumn umerointeirobin ario, arepresenta caodapartefracion arianessesistemadenumera cao exige, normalmente, um n umero bastante elevado de dgitos. Por exemplo, (0, 249)10 o qual difere de (0, 25)10 = (0, 01)2 por apenas (0, 001)10 n ao tem representa cao bin aria nita(pois essa diferenca nao pode ser representada de forma nita): seus primeiros 20 dgitos sao(0, 0011111110111110011100 . . .)2evejaque, como0, 249 0.Note que o padrao de bits0 00000000 00000000000000000000000n ao representa 0, mas sim 1 (uma vez queb0n ao e armazenado).Temos, entao, duas op coes para representar o zero:1. Representar explicitamenteb0: comisso,reduzimosaprecisao,poisobit b23damantissan ao poder a ser representado;2. Escolherumcertovalorde Eoqual, quandoopadr aodebits de Mfor 00 . . . 00, ser aconsiderado como representando o 0. Essa e a estrategia utilizada no padr ao IEEE-754,Note que, em ambas opcoes, persiste a representacao dupla para o zero, +0 e 0;usualmente osinal e desconsiderado, nessa situa cao.1.3.3 Caracteriza caodeumarepresentacaoAmdecaracterizarmosumarepresenta caoden umeros reais, sejaemponto-xoouponto-utuante, podemos denir algumas quantidades, as quais s ao:1. A precis ao,p, e a quantidade de bitsdisponvel para representar o n umero;2. O menor n umero representavel, em modulo, MINR;3. O maior n umero representavel, em modulo,MAXR;4. O menor n umero representavel, , tal que 1+ ,= 1, tambem chamado de epsilon da m aquinaou unidade de arredondamento da m aquina;5. A menor separa cao possvel entre dois n umeros representaveis, ULP (do ingles units-in-the-last-place).Para um sistema de ponto-xo, teremos ent ao1. p = [ e [ +[ d [;2. MINR e obtido fazendo-ses = 0,e = 0 e colocando 1 no bit menos signicativo ded;A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 14Introdu c ao ao Calculo Numerico Aritmetica no Computador3. MAXR e obtido fazendo-ses = 0, ee ed tendo todos os seus bitsiguais a 1;4. MINR;5. ULPMINR.Usando-se uma palavra de32 bits,dividida em campose com 15 bitsed com 16 bits,podemoscalcular essas quantidades, conforme mostra a tabela 1.5.p 31MINR 216 0, 000015MAXR (2151) +

16i=1 2i= 32767, 9999847412109375 215 216ULP 216Tabela 1.5: Valores caracterizadoresde uma representa cao em ponto-xo.J a para um sistema de ponto-utuante, essas quantidades s ao obtidas de forma diferente:1. p = [ M[ + 1 (poisb0 = 1 n ao e armazenado);2. MINR, e obtido fazendo-ses = 0,E = 128 eM= 1/2.3. MAXR, e obtido fazendo-ses = 0,E = 127 eMtendo todos os seus bitsiguais a 1;4. = 2(p1), para arredondamento por corte, ou =122(p1)= 2p, para arredondamentopor adi cao (ver 1.3.4);5. ULP=(0, 00 . . . 01)22E=2(p1) 2E=2E, paraqualquern umeroxnaformaM 2E.Usandoumapalavrade 32bits divididaconformeexpressoacima, essas quantidades temosseguintes valores, conforme mostra a tabela 1.6, onde foicalculado usando-se arredondamentoporcorte. Comparando comatabela1.5, ef acilnotarquearepresentacao emponto-utuantep 24MINR 212128 0, 146937 1038MAXR (

24i=1 2i) 2127 0, 170141 1039 2(241)= 0, 119209 106ULP 2E= 223+ETabela 1.6: Valores caracterizadoresde uma representa cao em ponto-utuante.oferece um intervalo muito maior de n umeros representaveis; alem disso, a separacao entre essesn umeros e bem menor.Na representacao em ponto-utuante, e importante estabelecer o intervalo de valores possveisparaoexpoente. Os limites desse intervalos aoomenor e omaior expoente, MINEe MAXE,respectivamente, esuadeni caodependedecomoosexpoentessaoarmazenados, conformeatabela 1.7.Cabe, aqui, uma observa cao referente ao . Suponha que se desconhe cam as caractersticas dosistema de ponto-utuante de um computador ou calculadora; nesse caso, e possvel estimar o,usando o algoritmo 1.3.1, o qual baseia-se na deni cao 1 + ,= 1:sinal-e-modulo complemento-de-2MINE |E|11 |E|1MAXE +|E|11 +|E|11Tabela 1.7: Denicao dos valores do menor e maior expoentes num sistema de ponto-utuante.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 15Introdu c ao ao Calculo Numerico Aritmetica no ComputadorAlgoritmo1.3.1Estimacao deproc macheps(output: )s 1, 0t 2, 0while(t > 1, 0)s 0, 5 st s + 1, 0if (t 1, 0) then 2, 0 sendifendwhileendprocPor exemplo, executando-se esse algoritmoemuma calculadoraHP-48SX, teremos comoresultado=7, 27595761411012. OprocessadorSaturndaHP-48SXutiliza38bits pararepresentar a mantissa e, portanto, o valor calculado para e uma boa aproxima cao.1.3.4 ArredondamentosConforme salientado anteriormente, qualquer representa cao de um n umero real, num computador,ser ainexata, salvoalgumaspoucasexce coes. Esseerronarepresentacaoest aassociado` abaseutilizadapararepresenta cao eaofatodeque, necessariamente,existeumn umeronitodebitspara armazenar o n umero.Quanto ` a base, apesar de alguns fabricantes de computadores terem utilizado outras bases (16,no caso doIBM360 e 8, noBurroughsB-6700), em 1960-1970, tipicamente se utiliza a base2, porsermaisf acildeseimplementaroscircuitosdoprocessador,usandoumalogicabin aria.Dessa forma, temos de conviver com o problema de certos n umeros, com representa cao exata embase decimal, nao poderem ser representados de forma exata em base binaria.Ooutroproblemaalimita caonotamanhodapalavrapararepresentarumn umerorealpodesermitigadaaoseaumentaraprecis ao. Qualquerlinguagemdeprograma cao cientcaoferece a possibilidade de se utilizar vari aveis em precis ao dupla e, algumas, em precisao qu adrupla,i.e., utilizamosduasouquatropalavraspararepresentarumn umeroreal. Outraalternativaeutilizar um processador cuja palavra contenha um maior n umero de bits: o recem-lancado IntelItaniumtemumapalavrade64bits. Noteadiferen casutil entreessasduasalternativas:umprogramaqueutilizevari aveisemprecis aosimples permitir asetrabalharcomn umerosdediferentes precis oes, se utilizarmos dois computadores com palavras de tamanhos diferentes. Porexemplo, um programa em Fortran90 com variaveis de precis ao simples tipo REAL tera umaprecis ao dep = 24 num computador que utilize oIntelPentiumII, mas, num IntelPentium4, o mesmo programa tera uma precisao dep = 53.Assim, com as limitacoes impostas pela escolha da base e precisao da representacao, pode-seperceber que existe um n umero nito de n umeros representaveis ou de m aquina; isso em marcantecontrastecomosn umerosreais, cujaquantidadeeinnita. Maisainda, entrequaisquerdoisn umeros reais, existem innitos outros n umeros; j a em qualquer das representacoes (ponto-xo ouponto-utuante), se tomarmos dois n umeros representaveis consecutivos, i.e., h a uma diferen ca de1 no bit menos signicativo, nao h a qualquer outro n umero. Assim, se posicionarmos os n umerosrepresentaveis sobre a reta dos reais, veremos que existem espacos entre cada n umero representavel.Para demonstrar isso, considere um computador hipotetico com uma palavra de 7 bits, e duasrepresentacoes de n umeros reais:1. Ponto-xo: [ s [ = 1, [ e [ = 2, [ d [ = 4;2. Ponto-utuante: [ s [ = 1, [ M[ = 5, [ E[ = 2.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 16Introdu c ao ao Calculo Numerico Aritmetica no ComputadorAs guras 1.1 e 1.2 mostram a distribui cao dos n umeros de maquina ao longo da reta dos reais.Paraa representacao em ponto-xo, note que a separacao entre os n umeros representaveis e constante:isso eexplicadopoisovalordeULP econstante. Janarepresenta caoemponto-utuante, essaseparacao aumenta ` a medida que um n umero torna-se maior, ou seja, o expoente cresce. Isso podeser vericado na expressao paraULP, 2E, a qual depende do valor do expoente.Outra caracterstica que pode ser observada e o menor intervalo de representacao no sistemadeponto-xo, bemcomoumaregi aodeunderowmuitomaiordoquenosistemadeponto-utuante.Observando asguras1.1e1.2, podemosvericaraexistencia deespa cosentreosn umerosrepresentaveis. Suponha, ent ao, umn umerocomo, porexemplo, 2/3=0, 666 . . ., cujarepre-sentacao em bin ario e (0, 101010. . .)2. Para ns de explana cao consideremos apenas um sistemadeponto-utuante, apesardasobserva coesaseguirseremv alidasparaumarepresenta caodeponto-xo tambem.Como 2/3 n ao tem representa cao nita em 24 bits, temos duas opcoes para armazena-lo:1. (0, 101010 . . . 1010)22. (0, 101010 . . . 1011)2o primeiro n umero e obtido descartando-se os bitsb24b25 . . .; j a o segundo obtem-se descartando-seos bitsem excesso e somando 1 ab24.Asitua caoquetemose, portanto, aseguinte: semprequeumn umeroenao-represent avel,devemosescolheroutroentreosdoisn umerosdemaquina, maispr oximosdaquele. Essesdoisn umeros s ao arredondamentosdo n umero original, chamados de arredondamento porcortee poradi c ao, correspondentes aos itens 1 e 2 acima, respectivamente.Como estamos aproximando um n umero n ao-representavel por um outro, o mais pr oximo dele,nossa representacao daquele n umero traz associada a si um erro, o qual pode ser medido de duasformas. Quando um n umero real x e aproximado por um n umero x, o erro e x x. O erro absolutoe denido como[ x x [ (1.5)e o erro relativo e dado porx xx. (1.6)Esse ultimo eomaisutilizadoporpermitirumacompara cao maisjustaentrequantidades comdiferentes relacoes de magnitude2.Vejamos formalmente, agora, quais os errosassociados a esses arredondamentos. Chamemosdexcexaos n umeros de maquina correspondentes aos arredondamentos por corte e por adi cao,respectivamente, de um n umero n ao-representavel x, os quais satisfazem a rela caoxc< x < xapoisxc= 0, 1010 . . . 1010[x = 0, 1010 . . . 1010[1010 . . .xa= 0, 1010 . . . 1011[porem,x pode estar mais pr oximo dexcou dexa, conforme mostrado na gura 1.3.Escrevendo os n umerosxcexacomoxc= (0, b0b1 . . . b22b23)22E(1.7)xa= _(0, b0b1 . . . b22b23)2 + 224_2E(1.8)podemos calcular os erros absoluto e relativo associados aos dois arredondamentos como2Porexemplo, umerrode1mnamedicaodadistanciaentreaTerraeJ upiterepequeno; porem, umerrode5cmnumaincisaonumcorpohumanopodeserconsideradobastantealto.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 17Introdu c ao ao Calculo Numerico Aritmetica no ComputadorFigura1.1: Distribui c aoden umerosdem aquinaemumsistemadeponto-xo; observequeadist ancia entre dois n umeros representaveis e a mesma.Figura1.2: Distribui c aode n umeros de m aquinaemumsistemade ponto-utuante; aqui, adist ancia entre dois n umeros representaveis aumenta ` a medida que se afastam do 0.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 18Introdu c ao ao Calculo Numerico Aritmetica no ComputadorFigura 1.3: Umn umero n ao-represent avel x pode ser melhor representado por xcou xa.Aqui, consideramos umsistemadeponto-utuantecomquatrodgitos namantissa, comxc=0, 6875 =(0, 1011)2e xa=0, 75 =(0, 1100)2; no diagrama ` a esquerda, x =0, 70625 =(0, 10110100110011 . . .)2e, ` a direita,x = 0, 734375 = (0, 101111)2.Errosassociados aoarredondamentopor corteSe o arredondamento por corte foi escolhi-do, ent aox encontra-se `a esquerda do ponto medio do intervalo [xc, xa]. Ent ao:[ x xc[ 12[ xaxc[ =12[ (M + 224) 2EM 2E[ = (1.9)21 224 2E= 2E25(1.10)x xcx2E25M 2E=225M22521= 224... M 12(1.11)Errosassociados aoarredondamentopor adi caoSe o arredondamento por adi cao foi esco-lhido, ent aoxencontra-se`adireitadopontomediodointervalo[xc, xa]. Poranalogia,escrevemos[ x xa[ 2E25(1.12)x xax 224(1.13)Generalizando, podemos dizer que o erro relativo entrex e seu arredondamentox ex xx (1.14)oux = (x) = x(1 + ), [ [ , =x xx(1.15)Uma medida muito utilizada para se determinar a qualidade numerica de um valor e o n umerode dgitossignicativos,DIGSE. Aplicando logaritmos aos dois lados da expressao 1.14, temoslog10x (x)x log10 (1.16)A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 19Introdu c ao ao Calculo Numerico Aritmetica no ComputadoreDIGSE e denido comoDIGSE(x, (x)) = log10x (x)x(1.17)Paraumaprecis aop=24, log10 2247, ouseja, temosnomnimosetecasasdecimaisdeprecis ao.1.3.5 Opera coesaritmeticasdeponto-utuanteNa soma e subtra cao, os expoentes dos dois operandos devem ser iguais. Para tal, seleciona-se omaior dos dois expoentes, e a mantissa e expoente do outro operando s ao ajustados de tal formaa coincidir os expoentes. Por isso, as opera coes aritmeticas s ao sempre efetuadas com o dobro debitsutilizados para armazenar os n umeros. Uma vez feito o ajuste dos expoentes, basta calcular(a rp) (b rp) = (a b) rpA multiplica cao e a divis ao sao calculadas como(arp) (b rq) = ab rp+q(arp) (b rq) = a b rpqComo essas operacoes sao efetuadas em v arias etapas, a cada parcela do processo, deslocam-se osbits `a esquerda, de forma a sobrar bits menos signicativos; a cada deslocamento, o expoente deveser modicado adequadamente.1.3.5.1 Errosemopera coesaritmeticasdeponto-utuanteSempre que dois n umeros de ponto-utuante sofrerem o efeito de uma das quatro opera coes arit-meticas, as seguintes etapas s ao efetuadas:1. Aopera cao efeitadeformacorreta, i.e., comodobrodon umerodebits usadosparaarmazenar cada operando;2. O resultado e normalizado;3.E feito o arredondamento, de forma que o resultado normalizado possa ser armazenado napalavra.O exemplo a seguir mostra por que deve-se efetuar a normaliza cao antesdo arredondamento.Exemplo1.1Considerex = 0, 45230102ey= 0, 25470103,emumsistemadeponto-utuante com cinco casas na mantissa. O resultado dex y, sem normaliza cao, ex y = 0, 0000011520Se efetuarmos o arredondamento, sem normaliza-lo, o resultado ser a 0!Comomostramosanteriormente,arepresentacao emponto-utuantetrazassociada asi umerro. O exemplo a seguir mostra, no entanto, que uma unica opera cao aritmetica simples tem umerro que n ao excede a.Exemplo1.2Considere x = 0, 31426103e y = 0, 92577105e, calculando as quatro operacoesaritmeticas, temos:x + y = 0, 9289100000 105x y = 0, 9226274000 105x y = 0, 2909324802 108x y = 0, 3394579647 102e, arredondando para cinco casas decimais, por adicao, temosA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 20Introdu c ao ao Calculo Numerico Aritmetica no Computadorerro relativo (x + y) 0, 92891 1058, 5 106< 105 (x y) 0, 92263 1052, 3 106< 105 (x y) 0, 29093 1082, 8 106< 105 (x y) 0, 33946 1026, 0 106< 105onde 105e o dessa representa c ao em cinco casas decimais.De forma generica, podemos dizer que, sex eysao n umeros representaveis, ent ao para umaopera cao aritmetica qualquer , (x y) = (x y)(1 + ), [ [ (1.18)e, sex ey n ao sao n umeros representaveis, ent ao ( (x) (y)) = (x(1 + 1) y(1 + 2)) (1 + 3), [ 1,2,3[ (1.19)A equa cao (1.19) nos diz que o erro associado ao encadeamento de operacoes aritmeticas podeser maior do que. Considere o exemplo abaixo: (x(y + z)) = (x( (y + z))(1 + 1), [ 1[ 224= (x(y + z)(1 + 2))(1 + 2), [ 2[ 224= x(y + z)(1 + 1 + 2 + 12) x(y + z)(1 + 1 + 2)... [ 1 + 2[ 223 x(y + z)(1 + 3), [ 3[ 223, 12 3Esse exemplo nos leva a supor que, caso a quantidade de opera coes aritmeticas a serem feitasseja muito grande, entao o erro crescer a proporcionalmente. O teorema a seguir mostra que essahip otese e verdadeira.Teorema1.3.1Sejamx0, x1, . . ., xnn umerosrepresentaveispositivos, eaunidadedearre-dondamentodamaquina. Entao, oerrorelativodearredondamentoaosecalcular ni=0 xinaordem natural, i.e.x0 + x1 + . . . + xn, e de no maximo (1 + )n1 ou, aproximadamente,n.Prova: SejaSk=x0 + x1 + . . . + xke (Sk),asquaispodemserrepresentadaspelasf ormulasde recorrencia_S0= x0Sk+1= Sk + xk+1, k 0(1.20)_ (S0) = x0 (Sk+1) = ( (Sk) + xk+1), k 0(1.21)e chamemos dekekaos erros relativos associadosa (Sk) e (Sk+1),k= (Sk) SkSk(1.22)k= (Sk+1) ( (Sk) + xk+1) (Sk) + xk+1(1.23)de ondek+1= (Sk+1) Sk+1Sk+1=( (Sk) + xk+1)(1 + k) (Sk + xk+1)Sk+1... (Sk) = Skk + Sk(por (1.22))...=(Sk(1 + k) + xk+1)(1 + k) (Sk + xk+1)Sk+1A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 21Introdu c ao ao Calculo Numerico Aritmetica no Computadore, rearranjando os termos, obtemosk+1 = k + kSkSk+1(1.24)Como, por deni cao,Sk< Sk+1e [ k[ ,[ k+1[ +[ k[(1 + ) = + [ k[, = 1 + podemos entao escrever:[ 0[ = 0[ 1[ [ 2[ + [ 3[ + ( + ) = + + 2...ou[ n[ + ( + ) = + + 2 + . . . + n1 == (1 + + . . . + n1) == n1 1== (1 + )n1== (1 + )n1e, pelo bin omio de Newton, tem-se(1 + )n1 = 1 +_n1_ =_n2_2+ . . . 1 n.1.4 PerdadedgitossignicativosApesar de erros de arredondamento serem inevitaveis e difceis de controlar, existem alguns tiposde erros, em computa coes numericas, que podem ser evitados.Por exemplo, suponha a subtra cao de dois n umerosx ey pr oximos entre si:x = .3721478693y = .3720230572x y = .0001248121Se o computador utilizado oferecer apenas cinco dgitos decimais na mantissa, teramos: (x) = .37215 (y) = .37202 (x) (y) = .00013Nesse caso, o erro relativo e bastante grande, da ordem de 4%.1.4.1 Subtra caodevaloresquaseidenticosComo regra, devemos evitar situacoes que levem `a subtra cao de valores quase identicos - normal-mente causados por express oes inadequadas do ponto de vista numerico.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 22Introdu c ao ao Calculo Numerico Aritmetica no ComputadorExemplo1.3Considere a express aoy _x2+ 1 1 (1.25)Ora, parax 0 entaox2e calculado por (1.28) ex1 =2cb b24ac=cax2(1.29)2. Seb24ac eb < 0 entaox1e calculado por (1.27) ex2 =2cb +b24ac=cax1(1.30)Porexemplo, sejaaequa caox2 106x + 1 = 0. Calculandox1ex2atravesdasexpress oesusuais, utilizandoosoftwareMaplecomprecis aode10dgitosdecimais, emumcomputadorPentiumII, temos:x1=106+101242=106+ 1062= 106x2=106101242=1061062= 0A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 23Introdu c ao ao Calculo Numerico Aritmetica no ComputadorAraizx2foi calculadadeformaerrada, poissofreucancelamentocatastr oco. Noentanto, serecalculamo-la usando (1.30), temosx2 =2106+10124=2106+ 106= 106Ao substituirmos esse valor na equa cao, teremos 10121 +1, o qual pode ser considerado comoaproximando 0 para a precis ao utilizada. Note que, para esse exemplo particular,x1nao e raiz daequacao, pois 1012106106+ 1 ,= 0.1.4.2 TeoremasobreaperdadeprecisaoUmaquestaoquesurgeapartirdosexemplosanteriores easeguinte: Quantosdgitosbinariossignicativos sao perdidos na subtra c aoxy quandox y? O teorema a seguir nos da limitantesextremos para o n umero de dgitos bin arios perdidos nessa situa cao, baseado na rela cao [ 1y/x[,o que nos d a uma medida de quao pr oximox e dey.Teorema1.4.1Sexeys aon umerosemponto-utuantebin ariospositivos, normalizados, talquex > ye2q 1 yx 2pentao no maximoq e no mnimop dgitos binarios signicativos s ao perdidos na subtra caox y.Prova. Considerando apenas o extremo superior da desigualdade,temos quex eys ao da formax = r 2n,_12 r < 1_y = s 2m,_12 s < 1_Comox>y, porhip otese, oexpoentedeydever aserigualadoaodexantesdeserealizarasubtra cao (note que, como nao pode haver dgitos `aesquerda do ponto binario, sempre se faz comque o menor n umero iguale seu expoente ao maior, introduzindo zeros imediatamente ` a direita doponto binario). Logo,ydeve ser escrito comoy = (s 2mn) 2ne, da,x y = (r s 2mn) 2nA mantissa desse n umero satisfaz a seguinte rela cao:r s 2mn= r_1 s 2mr 2n_= r_1 yx_< 2pPara normalizara representacaodex y, um deslocamento deaomenosp bitspara a esquerdae necessario. Entao, ao menosp zeros s ao inseridos ao nal da mantissa, efetivamente perdendop bits de precisao.O exemplo a seguir ilustra a utiliza cao desse resultado.Exemplo1.5Suponhaumamantissade5dgitos decimais eque x=0, 31457105e y =0, 31453 104e que os c alculos sejam efetuados com o dobro de dgitos. Ora, para calcularx y,temos:x = 0, 31457 00000105y = 0, 03145 30000105x y = 0, 28311 70000105 (x y) = 0, 28311 105A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 24Introdu c ao ao Calculo Numerico Aritmetica no Computadorcomo foi necess ario inserir um dgito 0 apos o ponto decimal dey, espera-se que se perder a umdgito ao fazer a normaliza cao, como pode-se vericar pelo resultado. O erro relativo, no caso, e0, 283117 1050, 28311 1050, 283117 105= 2, 4724760435105, 105o que demonstra que, efetivamente, houve perda de dgitos signicativos na subtracao.Exemplo1.6Considerea express aoy = x senx. Comosenx x parax 1, ocorrer a perdadedgitos signicativosemy. Proponhaumaformaalternativaparacalcular yeestipuleumintervalo parax no qual pode-se utilizar a express ao original.Solu cao:Usando a serie de Taylor para senx, temos:y = x senx= x _x x33!+x55! x77!+ . . ._=_x33! x55!+x77! . . ._e, parax 0, podemos truncar a serie como segue, utilizando apenas quatro termos: y = (x3/6)(1 (x2/20)(1 (x2/42)(1 x2/72)))UsandooTeoremadaPerdadePrecis ao, podemosexigirqueapenas umbitsejaperdidoseescrevermos1 senxx12, x > 0,Essa desigualdade e satisfeita se [ x[ 1, 9 e, nessa situacao, podemos usar a express ao original.Para 0 < x < 1, 9, devemos usar a expressao baseada na serie truncada de Taylor, pois ela eliminao problema.Suponhax = 106em uma calculadora HP-48SX. Ent ao, teremos:y = 106sen(106) = 0, 000001 0, 000001 = 0 y = 1, 66666666667 10191.5 CondicionamentodeumproblemaO condicionamento de um problema diz respeito a quao exato podemos resolve-lo em uma dadaprecis ao de ponto-utuante, independentementedo algoritmo utilizado para resolve-lo.Seguindoaderiva caoem[12], esupondoquedesejamosavaliarumafun caoy =f(x), j asabemosquequalqueropera cao deponto-utuanteacarretar a aexistencia deumerro. Logo,oque efetivamente calcula-se e uma aproximacao (y) = (f( (x)))mas, por simplicidade, assumimos que (f) =f. Podemos calcular, ent ao, o erro relativo emycomo (y) yy=f( (x)) f(x) (x) x

xf(x) (x) xx.O termof((x))f(x)(x)xe uma aproximacao paraf

(x). Assim, (y) yy f(x) (x) xx(1.31)A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 25Introdu c ao ao Calculo Numerico Aritmetica no Computadorondef(x) = [ x[ [ f

(x) [[ f(x) [(1.32)a qual e chamada de n umero de condi cao defemx. Esse fator e que mede o quanto os erros dearredondamento emx sao amplicados ao se avaliarf(x).Ent ao,para seavaliar on umero dedgitos corretos em (y), aplicamos logaritmos aosdoismembros da equacao (1.31),log10_ (y) yy_ log10_ (x) xx_log10 f(x)Note que log10_(x)xx_ e aproximadamente igual a 7 parap = 24 e a 16 parap = 53.1.6 Computa coesestaveiseinstaveisUm processo numerico e dito inst avel se pequenos erros ocorridos num passo sao ampliados nospassos seguintes, degradando a exatid ao do processo.Considere, por exemplo, a seq uencia de n umeros dada por___x0 = 1x1 =13xn+1 =133xn43xn1, n 1a qual gera os n umerosxn =_13_npoisx0 =130= 1,x1 =131=13. Paran = m + 1, temos:xm+1=133xm 43xm1 =133_13_m 43_13_m1=_13_m1_133 43_=_13_m+1Utilizando a forma de recorrencia acima para gerar os n umeros, teremos:n xn (1/3)nerrorelativo0 1,000000000000000000000 1,000000000000000000000 0,0000000000000000000001 0,333333333333333310000 0,333333333333333310000 0,0000000000000000000002 0,111111111111110940000 0,111111111111111100000 0,0000000000000014988013 0,037037037037036258000 0,037037037037037028000 0,0000000000000207958654 0,012345679012342514000 0,012345679012345677000 0,0000000000002561544735 0,004115226337435884400 0,004115226337448558300 0,0000000000030797552016 0,001371742112432145600 0,001371742112482852900 0,0000000000369655985517 0,000457247370624785240 0,000457247370827617560 0,0000000004435942960618 0,000152415789464541850 0,000152415790275872500 0,0000000053231404445499 0,000050805260179967644 0,000050805263425290837 0,00000006387769640489010 0,000016935074827137338 0,000016935087808430279 0,00000076653236686197211 0,000005644977344304949 0,000005645029269476759 0,00000919838841059637612 0,000001881468722471661 0,000001881676423158920 0,00011038066093717230013 0,000000626394671637267 0,000000627225474386307 0,00132456793125653350014 0,000000205751947132609 0,000000209075158128769 0,01589481517508941800015 0,000000056398875391618 0,000000069691719376256 0,19073778210108555000016 -0,000000029940802813135 0,000000023230573125419 2,28885338521303970000017 -0,000000204941979379076 0,000000007743524375140 27,466240622556484000000Osvaloresacimademonstramqueoalgoritmo einstavel, equalquererropresenteemxnemultiplicado por 13/3 emxn+1. Portanto, h a a possibilidade de queo erro existente emx3(daordemde1015)sejapropagado parax17porumfator(13/3)14109; ouseja, oerroemx17devido unicamente ax3pode ser de 104, que n ao e desprezvel. Alem disso, os erros devido aosdemais n umerosx4, x5, . . . , xksao propagados parax17por fatores da forma (13/3)k.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 26Introdu c ao ao Calculo Numerico Aritmetica no Computador1.7 Desastres causados por erros aritmeticos no computadorEsse captulo apresentou os conceitos ligados ` a computa cao numerica; dentre estes, certamente oconceito de erroassociado a qualquercalculo numerico e o mais importante. Como nao h a comoevit a-los, enecessario queoprogramador e/ou analista numerico saiba como trata-los deformaque a ocorrencia deles nao leve a falhas catastr ocas. Infelizmente, isso nem sempre e levado emconta, no dia-a-dia, e desastres ocorrem, como os dois que citamos a seguir.1.7.1 FalhadosistemademsseisPatriotDuranteaGuerradoGolfo, em1991, oIraquelan couin umerosmsseisterra-terraScud(defabrica cao sovietica) contra Israel e Ar abia Saudita. A m de se protegerem contra esses ataques,astropasnorte-americanasinstalarambateriasdemsseisterra-arPatriot, osquaishaviamsido projetados no incio da decada de 70 para destrurem msseis cruzeiro e aeronaves sovieticas(voando a uma velocidade media de 2Mach), numa eventual guerra entre a OTAN e o Pacto deVarsovia.Umabateria demsseis Patriot consistedeumaunidadedecontrole computadorizada;deum radar de deteccao; e de ate 6 lan cadores qu adruplos de msseis. A unidade de controle dispoede um rel ogio que marca o tempo em decimos de segundo, armazenados em uma palavra inteira de24 bits; os calculos de determina cao das janelas de conrmacao e de engajamento (regi oes no ceudentro do qual o possvel alvo deve ser detectado pelo radar para que possa os msseis Patriotsejam lan cados) sao feitos em ponto xo, tambem com 24 bits.No dia 25 de fevereiro de 1991, uma bateria Patriot instalada em Dharan, na Ar abia Saudita,deixou de interceptar ummssil Scud quese aproximava. Como resultado, 28 soldados norte-americanos foram mortos devido ` a explosao do mssil Scud.Osresultadosdainvestiga cao, deacordocom[11], indicaramqueosmsseisPatriotnaoengajaram o Scud (apesar dos radares haverem detectado o mssil iraquiano) devido a um erronumerico de arredondamento.O tempo medido pelo rel ogio da unidade de controle e multiplicadopor 1/10 para representar o tempo em segundos, e armazenado em 23 bits, no formato de pontoxo; ocorre que 1/10 e um n umero que n ao tem representa cao nita em bin ario:(1/10)10 = (0, 0001 1001 1001 1001 1001 1001 1001 100 . . .)2Como a palavra usada para armazenar o rel ogio e de 24 bits, o erro e de aproximadamente(0, 0000 0000 0000 0000 0000 0001 100 . . .)2 (0, 0000 0009 5)10.E obvioque, ` amedidaqueotempodeopera cao dabateriademsseis aumenta, maiorseraoerro notempocalculado. Comoessetempo eusado parasecalcularasjanelasdedetec caoeengajamento de um alvo, isso ira causar um deslocamento da janela para baixo (i.e., a altitude naqual se espera que o alvo aparecer a na pr oxima varredura do radar da bateria ser a menordo quea que elese encontra). Com isso,o mssil continuar a trafegando em dire cao ao seu alvo, poremn ao ser a detectado e, portanto, os msseis Patriot n ao ser ao disparados.Ocorre que, devido ` as caractersticas do sistema Patriot, a m de se maximizar as chancesdederrubada domssil atacante,estedeveencontrar-se nomeio dajaneladeengajamento,aqual tem um comprimento de 274m. No dia 11 de fevereiro de 1991 (duas semanas antes da falhadosistema), vericou-se queapos 8hdeopera cao contnua,ajanelasofria umdeslocamentode55m; por extrapola cao, ap os 20h de uso contnuo, esse deslocamento seria de 137m, e a partir deent ao, n ao seria mais possvel detectar um mssil atacante.Quando ocorreu a falha, a bateria que deveria ter engajado o Scud estava operando conti-nuamente por 100h! Dessa forma, o erro acumulado era de0, 000000095 100 3600 10 = 0, 34s.Um mssil Scud viaja a uma velocidade terminal de 1.676m/s; em 0, 34s, ele percorre a dist anciade 569, 84m. O deslocamento da janela de deteccao, ap os 100h de opera cao, era de 687m. Logo,A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 27Introdu c ao ao Calculo Numerico Aritmetica no Computadorn ao havia comooScud serdetectado,j aqueajanelaencontrava-se aumaaltitudeinferior ` adele.1.7.2 ExplosaodofogueteAriane5No dia 4 de junho de 1996, o primeiro foguete Ariane 5, construdo pela Agencia Espacial Europeia,foi destrudopelosistemadecontroledefalhaapenas40sap osolan camentodasuabaseemKhourou, na Guiana Francesa. O Ariane 5 havia sido desenvolvido ap os 10 anos de trabalho, aum custo de 7 bilh oes de dolares. O custo do foguete, bem como da carga util transportada, erade 500 milh oes de dolares.Osresultadosdainvestigacao, ap osduassemanasdoincidente, indicaramqueoproblemaencontrava-senosoftwaredeguiageminercial. Umn umeroemponto-utuante, armazenadonuma palavra de 64 bits,e que representava a velocidade horizontal em rela cao `a plataforma delan camento, foi convertidoparaumn umerointeiro, noformatosinal-e-magnitude, de16bits.Como a velocidade era superior a 32.768 (o maior n umero representavel em 15 bits), ocorreu umafalha na conversao e o programa deixou de funcionar.1.8 ExercciosExerccio1.1Considere umsistemade ponto-utuante, F =(2, 24, 8), noqual os n umerosapresentamumamantissatal que1M 0, pode existir um n umero par (0, 2, 4, . . .) de razes reais em [a, b].3. Supondo quefe sua derivadaf

sejam contnuas em [a, b] e que o sinal def

seja constanteneste intervalo tem-se:(a) sef(a)f(b) < 0, entao existe uma unica raiz real em [a, b];(b) sef(a)f(b) > 0, entao nao existe raiz real em [a, b].Exemplo2.1Determinegracamenteosintervalosquecontemcadaumadasrazesreaisdasfun c oesf(x) = x3sen(x) ep(z) = z34z2+ 4z 1.Alem disso, e possvel que existam razes de multiplicidade maior do que um; nesse caso, aindaqueummetodoconsigadeterminarumadelas, adetermina caodamultiplicidadedeladeveserfeita posteriormente. Outro problema que afeta a determina cao numerica de uma raiz e quandoexistem duas razes tao pr oximas numericamente que a n ao se pode garantir para qual delas umdeterminado metodo ir a convergir.29Introdu c ao ao Calculo Numerico Razes de Fun coes Nao-LinearesNessecaptulo, apresentaremososmetodosdabissec cao, daposi caofalsa, deNewtonedasecante.2.2 MetododaBissec caoO metodo da bisseccao e um dos mais metodos mais simples para se obter uma raz de uma fun cao.Ele baseia-se na divisao sucessiva, ao meio, de um intervalo [a, b], no qual e garantido que h a umaraiz da equacaof(x) = 0, i.e., sign(f(a)) = sign(f(b)) - da o seu nome.Formalmente, sefe uma fun cao contnua no intervalo [a, b], esef(a)f(b)< 0,ent aoftemumzero em(a, b). Como h atrocadesinaldafun caofavaliada nos extremosdointervalo,emalgumpontoa0f(x + h) f(x)hA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 38Introdu c ao ao Calculo Numerico Razes de Fun coes Nao-LinearesFigura 2.11: O metodo de Newton-Raphson - interpreta cao geometrica.vemosqueocoecienteangulardaretatangenteeapr opriaderivadaf

(x0)(desdequef

(x)exista e seja contnua). Assim, podemos escrever (2.11) comox1 = x0f(x0)f

(x0)(2.13)e, generalizando para uma estimativaxk,xk+1 = xkf(xk)f

(xk), k = 0, 1, . . . (2.14)aqual eaequacao governante dometodo deNewton-Raphson. Ometodo deNewton-Raphsonpode ser expresso de forma algortmica como segue:Algoritmo2.4.1Newton-Raphsonproc newton raphson(input:x0,,,kmax; output:xk+1,k)fork = 0, 1, . . . , kmaxdoxk+1 xkf(xk)f

(xk)if | xk+1xk| < OR | f(xk+1) | 0,x0 = a; sef(b)f

(b) > 0,x0 = b. Caso contr ario,x0 =a+b2.Exemplo2.4Considere o polin omiop (x) = x35x2+ 8x 4. Calcule duas razes utilizando ometodo de Newton-Raphson.Solu cao:A f ormula iterativa neste caso exi+1 = xix3i 5x2i + 8xi43x2i 10xi + 8.Tomandox0 = 0, obtemos os valores mostrados na tabela 2.1.i xi0 01 0,52 0,83 0,954 0,9956521745 0,9999626796 0,9999999987 1,0Tabela 2.1: Valores aproximados da solu cao pelo metodo de Newton-Raphson.Para calcular a raiz seguinte, parte-se dex0 = 1, 75, com o qual obtemos os valores mostradosna tabela 2.2.Quest ao:Por que raz ao a convergencia para a segunda raiz e mais lenta?O metodo de Newton pode ser aplicado de maneira um pouco diferente quando a funcaof(x)e um polin omio, conforme veremos no captulo 3.2.5 Deriva caonumericaNo metodo de Newton, e necessario utilizar o valor de f

(xn) a m de se calcular a nova estimativaxn+1. No entanto, avaliarf

pode ser oneroso (a menos quefseja um polinomio) e, por isso, ` asvezes recorre-se a uma aproxima cao numerica def

.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 41Introdu c ao ao Calculo Numerico Razes de Fun coes Nao-LinearesFigura 2.13: Ausencia de razes reaisFigura 2.14: Segunda derivadaf

(x) = 0.Figura 2.15: Dist ancia inadequada entrex0ex.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 42Introdu c ao ao Calculo Numerico Razes de Fun coes Nao-Linearesi xi0 1,751 1,92 1,9529411763 1,9770662744 1,9886693185 1,9943672806 1,997191745...16 1,99999995617 1,999999956Tabela 2.2: Valores aproximados da solu cao pelo metodo de Newton-Raphson.Considere a deni cao da derivada def(x) em termos do limite,f

(x) f(x + h) f(x)h(2.17)Se f e linear (f(x) = ax+b), ent ao a Equa cao (2.17) e exata, i.e., para qualquer h = 0, ela nos d ao valor correto def

(x). Sef(x) n ao for linear, somente em casos muito especiais ela ser a exata;logo, h a um erro envolvido nessa aproxima cao, o qual pode ser mensurado usando o teorema deTaylor:f(x + h) = f(x) + hf

(x) +h22f

() (2.18)onde x 0, sex = 0, x V|| x|| = | | || x||, se IR, x V|| x + y || || x|| +|| y || sex, y V (desigualdadetriangular)Anorma deumvetor eoseucomprimento noespa covetorial V ;eumageneraliza cao dano cao de valor absoluto de um n umero real. Para o espa co vetorial IRn, a norma mais conhecidae a chamada norma Euclidiana, denida por|| x||2 =_n

i=1x2i_12(4.12)onde x = (x1, x2, . . . , xn)T. Particularmente, em IR2, temos || x||2 =_x21 + x22, que e a expressaopara a dist ancia de um ponto com coordenadas (x1, x2) em relacao `a origem do sistema de eixoscartesianos.Existem outras normas que sao bastante usadas em calculos numericos, como a norma-l|| x|| =nmaxi=1| xi| (4.13)e a norma-l1,|| x||1 =n

i=1| xi| (4.14)as quais sao bem mais simples e menos onerosas de se calcular do que a norma Euclidiana.4.4.2 NormasdematrizesUma vez especicada uma norma de um vetor, a norma matricial subordinada e denida como|| A|| = sup|| Au || : u IRn, || u || = 1 (4.15)para uma matrizAn n. Pode-se vericar que|| Ax|| || A|| || x||, x IRnPor exemplo, a norma matricial subordinada da norma vetorial || || e dada por|| A|| =nmaxi=1n

j=1| aij| (4.16)A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 74Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Lineares4.4.3 N umerodecondi caodeumamatrizNormasdevetores edematrizesnospermitemavaliar oquaosuscetvel aerros numericos ser auma computa cao empregando-se uma dada matrizA. Para tanto, suponha que se deseja resolvero sistemaAx = b, ondeA en n eA1existe.SeA1temseusvaloresperturbados(istoe, ligeiramentemodicados), gerandoumanovamatrizB, asolu caodosisteman ao emaisx=A1bmas x=Bb. Essaperturbacao podesermedida em termos do comprimento do vetorx x,|| x x|| = || x Bb || = || x BAx|| = || (I BA)x|| || I BA|| || x||ou|| x x|||| x|| || I BA||o que nos d a uma no cao do erro relativo entrex e x.Deforma an aloga, suponha queb foiperturbado,gerando umnovo vetor b. Sex e x sao assolu coes deAx = b eA x = b, podemos medir o erro absoluto entrex e x escrevendo|| x x|| = || A1b A1b || = || A1(b b) || || A1|| || b b ||e o erro relativo como|| x x|| || A1|| || b b || = || A1|| || Ax|||| b b |||| b || || A1|| || A|| || x|||| b b |||| b |||| x x|||| x|| || A1|| || A|||| b b |||| b ||o que nos diz que o erro relativo emx e limitado pelo n umero || A1|| || A||. Essa quantidade edenominada de n umero de condi cao deA, e e denotada por(A) = || A1|| || A|| (4.17)Vejamos um exemplo do uso de(A).Exemplo4.3Seja a matrizA e sua inversaA1,A =_1 1 + 1 1_A1=12_1 1 1 + 1_.Usando a norma-l, entao || A|| = 2 + e || A1|| =12(2 + ), de onde(A) =_2 + _2>42.Se 0, 01,entao(A) 40000. Issoquerdizerque,sebsofrerumapequenaperturba cao, aperturba c ao relativa na solu cao do sistemaAx = b ser a 40000 vezes maior!Umamatrizque tenhaumn umerode condi caomuitograndeeditamal-condicionada, epequenas varia coes nos valores de b induzir ao um grande erro relativo no vetor solu cao do sistema.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 75Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Lineares4.4.4 ErroscomputacionaisecondicionamentoQualquer solu cao de um sistema linear deve ser considerada uma solucao aproximada, em virtudede erros de arredondamento e outros. O metodo mais natural para determina cao da precisao deuma solu cao e vericar quao bem esta solucao satisfaz o sistema original, calculando o vetor resduo.Seasolu cao aproximada xfor umaboa aproxima cao, pode-seesperar quecada componenteder = b A x seja pequeno, pelo menos em um conceito relativo. Ha sistemas de equacoes, contudo,em que o resto nao proporciona uma boa medida da precis ao da solu cao. S ao sistemas nos quaispequenas alteracoes nos dados de entrada conduzem a mudan cas signicativas na solu cao. Estessao denominados sistemas inst aveisou mal-condicionados.Exemplo4.4A solu c ao exata do sistema_x1 + x2 = 21, 01x1 + x2 = 2, 01ex1 = x2 = 1. Supondo que, devido a erros, a solu cao calculada fossex1= 0x2= 2, 005o vetor resduo neste caso seriaRT= [0, 005; 0, 005].Entretanto, o erro em cada resposta,x1ex2, e de aproximadamenteuma unidade.Poroutrolado, oscoecientestambempodemcontererros. Supondoquealgumtipodeerrotenha mudado as equac oes acima para_x1 + x2 = 21, 0001x1 + x2 = 2, 007atemesmoumasolu c aobemdiferentedaanterior, comox1=100ex2= 98produziriaumresduo bem pequeno,RT= [0; 0, 003].Erros deste tipo, ao contr ario daqueles causados pela acumulacao de erros de arredondamento,n ao podem ser evitados por uma programa cao cuidadosa. Como, ent ao,determinar quando umproblema e mal-condicionado?Figura4.4: Ossistemasquedescrevemaintersec caodasretassao, daesquerdaparaadireita:bem-condicionado,mal-condicionadoe singular.Emgeral, tem-seasituacaomostradanagura4.4, paraocasodeduasretas. Quandoosistema e ordem maior, no entanto, deve-se recorrer a medidas algebricas para se estimar o mal-condicionamentodosistema. Umadessasmedidaseochamadodeterminantenormalizadodamatrizdoscoecientes. Para obte-lo,normaliza-se amatrizdecoecientesAdividindo-secadaA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 76Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Lineareslinha deA pela raiz quadrada da soma dos quadrados dos elementos de cada linha,norm| A| =a111a121 a1n1a212a222 a2n2............an1nan2n annn=| A|12 . . . n(4.18)onde i= _a2i1 + a2i2 + . . . + a2in. Diz-se, entao, queumamatriz Aemal-condicionadaseon umero norm| A| for pequeno, comparado com a unidade.Exemplo4.5Seja a matrizA =_1 11 1, 01_;verique se ela e mal-condicionada.Solu cao: Calculando-se1 =1 e2 =2, 0201, podemos obternorm| A| =1 11 1, 0112 0, 005ou seja,A e dita ser mal-condicionada.4.4.5 MetodositerativosDado um sistema nao-singular den equa coes linearesAx = b, um metodo iterativo para resolveresse sistema e denido por um conjunto de fun coes k(x0, x1, . . . , xk, A, b), ondex0 = 0(A, b) euma estimativa inicial para a solu caox = A1b ex1,x2,. . . sao as aproxima coes sucessivas paraa solu cao,x1= 1(x0, A, b)x2= 2(x0, x1, A, b)...xk= k(x0, x1, . . . , xk, A, b)Asfun coesknosdenemosmetodositerativos. Diz-sequeummetodo eestacionariose,para umm> 0, nn ao depende den para todon m,ou seja, = m= m+1=. . . Nessecaso, xn+1dependede,nomaximo, mvetores anteriores, xn, xn1, . . ., xnm+1. Por exemplo,param = 2, temosx0= 0(A, b) (4.19)x1= 1(x0, A, b) (4.20)xk= (xk2, xk1, A, b), k = 2, 3, . . . (4.21)O grau de um metodo estacionario e m (para m m) se, paran m 1, xn+1depende dexn,xn1,. . .,xn m+1mas nao parak < n m+. O grau de um metodo iterativo denido pelasequa coes (4.19)-(4.21) e 2.Ummetodo iterativo e ditolinearsetodas as fun coesisao fun coes lineares dex0, x1, . . .,xn1. Assim, um metodo iterativo estacionario linear de grau 1 pode ser expresso porxk+1 = Gxk + f (4.22)ondeG e uma matriz efum vetor, escolhidos adequadamente.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 77Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesPara um metodo como em (4.22), podemos nos referir a um sistema linear relacionado,(I G)x = f; (4.23)ondeIeamatrizidentidadedeordemn. Porexemplo, seG=I A, f b, ent ao(4.23)eequivalente aAx = b.A deni cao de um metodo iterativo pode tambem ser feita a partir de uma matriz separadora,Q. Podemos escrever o sistemaAx = b na forma equivalenteQx = (QA)x + b (4.24)isto e, x = (x, A, b), o que nos leva a escrever um processo iterativo, de aproximacoes sucessivas,comoQxk = (QA)xk1 + b, k = 0, 1, . . . (4.25)A matrizQ deve ser escolhida de tal forma que se possa calcular rapidamente osxke que aseq uenciax0,x1,. . . convirja rapidamente para a solu caox = A1b.A m de obter uma condi cao necess aria e suciente para que haja convergencia, reescrevemos(4.25) comoxk = (I Q1A)xk1 + Q1b (4.26)A solu caox satisfaz a equa caox = (I Q1A)xk1 + Q1b (4.27)i.e.,x e um ponto xo do mapax (I Q1A)xk1 + Q1b (4.28)Usando as equacoes (4.26) e (4.27), podemos obter uma express ao para o erroxk x comoxkx = (I Q1A)(xk1x) (4.29)e, aplicando normas, temos|| xkx|| || (I Q1A) || || (xk1x) || || (I Q1A) ||2|| (xk2x) ||... || (I Q1A) ||k|| (x0x) || (4.30)de ondelimk|| xkx|| = 0, se || (I Q1A) || < 1 (4.31)desde queA eQ sejam invertveis.4.4.6 RenamentoiterativoO primeiro metodo iterativo para a solu cao de um sistema de equacoes lineares Ax = b e o chamadorenamento iterativo. Para uma estimativa inicialx0, denimos o vetor erroe0 comoe0 = x x0(4.32)e o vetor resduor0comor0 = b Ax0. (4.33)O vetor erro nos diz o quantox0esta distante dex, e o vetor resduo nos diz o quantoAx0estadistante deb.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 78Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesMultiplicandor0porA1, temosA1r0 = A1b x0 = x x0 = e0(4.34)e, usando essa igualdade, podemos obter uma express ao parax:x = x0 + A1r0(4.35)Notequeaequa cao(4.35)envolveA1; mas, obviamente,n aopodemosutiliza-la, poisseacalcul assemos, a solucao do sistema seria imediata! Por outro lado, a equa cao (4.34) nos permiteescreverAe0 = r0(4.36)e, combinandoasequa coes(4.33), (4.36)e(4.35), podemosdescreverometododorenamentoiterativo como ___rk = b AxkresolveAek = rkxk+1 = xk + ek, k = 0, 1, . . . (4.37)Ometododorenamentoiterativo eutilizadoemconjuntocomummetododireto, comoaelimina caoGaussiana. TendofatoradoAnoprodutoLUeobtidoumasolu cao xparaAx=b(a qual pode n ao ser muito boa, devido a erros de arredondamento), fazemosx0 = x e renamosessa solucao, usando (4.37), ate quexkseja sucientemente bom. Note que a fatoracao LUpode,agora, ser utilizada para resolver Aek = (LU)ek = rk.Se consideramos que a solu cao obtida com a fatora cao LUdeA n ao foi exata, ent ao podemosdizer queU1L1= B A1. Usando a equa cao (4.35), escrevemosxk+1= xk + A1rk = xk + A1b A1Axk... B A1...= xk + B(b Axk) (4.38)De onde podemos mostrar que o metodo converge para uma solucao: subtraindox de ambos oslados da equacao (4.38), temosxk+1x = xkx + B(b Axk)... b = Ax...= xkx + B(Ax Axk) = (I BA)(xkx)e, tomando normas de ambos os lados da igualdade acima, vem, pela desigualdade triangular:|| xk+1x|| || I BA|| ||xkx||... || xkx|| = || I BA|| || xk1x||... || I BA||2||xk1x||... || I BA||k||x0x||o que nos diz que os erros convergem para 0 se || I BA|| < 1.Assimcomonos metodosdedeterminacaoderazesdefun coes, precisamosdenir algunscriteriosdeparadadoprocessoderenamento. Oprimeirodessescriterioseaestipula caodeumn umerom aximodeitera coes(kmax); osegundopodeserbaseadonanormadoresduork,devidamente escalonada por || b || (usando uma norma qualquer, previamente escolhida). Assim,as itera coes proceder ao enquanto|| rk|| < || b || (4.39)n ao for satisfeito; e um n umero real, escolhido de acordo com a exatidao requerida.Um algoritmo que expressa o metodo do renamento iterativo pode ser escrito como:A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 79Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesAlgoritmo4.4.1Renamento Iterativoproc renamento iterativo(input:A,L,U,b,x0,kmax,; output:xk)t || b ||fork = 0, 1, . . . , kmaxdork b Axkif || rk|| < t thenbreakendifresolve (LU)ek = rk, obtendoekxk+1 xk + ekendforendprocO exemplo abaixo [10] mostra o comportamento desse metodo.Exemplo4.6Seja o sistema__420 210 140 105210 140 105 84140 105 84 70105 84 70 60__x =__875539399319__cujasolu cao eovetorx = (1, 1, 1, 1)T. Utilizandoum computadorcomapenas6casasdecimaisde precis ao, obtemos como solucao inicial, atraves da eliminacao Gaussiana, com pivotamento, ovetorx = (0, 999988, 1, 000137, 0, 999670, 1, 000215)TAgora, dispondodosfatorestriangularesdafatora caoLU, podemosutilizaroalgoritmo4.4.1eobter:x = (0, 999994, 1, 000069, 0, 999831, 1, 000110)Tx = (0, 999996, 1, 000046, 0, 999891, 1, 000070)Tx = (0, 999993, 1, 000080, 0, 999812, 1, 000121)Tx = (1, 000000, 1, 000006, 0, 999984, 1, 000011)T4.4.7 MetodoiterativodeJacobiSuponha o sistema (4.1), comn = 3, sem perda de generalidade. Se os elementos da diagonal deA sao todos n ao-nulos, entao pode-se isolar cada vari avel x1,x2ex3atraves de___x1= c12x2+ c13x3+ d1x2= c21x1+ c23x3+ d2x3= c31x1+ c32x2+ d3ondecij=_ aijaii, i = j0, i = jdi= bi/aiiCom essa transformacao, o sistemaAx = b foi transformado em um sistema da forma(I C)x = dA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 80Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesondeC=D1(D A), d =D1b eD= diag(A)(isto e,amatriz formada pelos elementos dadiagonal deA). De forma equivalente, podemos escreverx = Cx + do que sugere uma correcao dex por aproxima coessucessivas,xk+1= Cxk + d = D1(D A)xk + D1b == (I D1A)xk + D1b, k = 0, 1, . . . (4.40)a qual dene o metodo iterativo de Jacobi.Amatrizseparadora,aqui,eD1; ometododeJacobiconverge seamatrizAfordiagonaldominante, i.e.,| aii| >n

j=1j=i| aij| (4.41)e, usando a norma-l,|| I D1A|| =max1inn

j=1j=iaijaiide onde pode-se vericar que a domin ancia diagonal e condi cao necess aria para a convergencia dometodo.Para obtermos a solu cao do sistemaAx = b via o metodo iterativo de Jacobi, podemos usar aforma equivalente a (4.40),xk+1 = xkD1Axk + D1b (4.42)Note que, do ponto de vista de eciencia do processo, deve-se efetuar as divis oes de cada linha deA e do elemento respectivo deb pelo elemento na diagonal deA antes de se iniciar as itera coes.Alemdisso, seocriterio deparada envolve oc alculodoresduork=b Axk, isso exigiriaumprodutomatriz-vetoramaisporitera cao, oquepodeserevitadoseusarmoscomocriteriodeparada rk+1 = D1rk+1,|| D1rk+1|| || D1|| || b || (4.43)pois rk+1 = D1rk+1 = D1b D1Axk+1e os termos no lado direito da equacao j a foram calculados, anteriormente, para se obterxk+1. Oalgoritmo a seguir utiliza essas ideias:A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 81Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesAlgoritmo4.4.2Metodo de Jacobiproc jacobi(input:A,b,x0,kmax,; output:xk+1)fori = 1, 2, . . . , n doqi a1iiendfort || q || || b ||fori = 1, 2, . . . , n doforj = 1, 2, . . . , n doaij aij qi% sobrescreveA comD1Aendforbi bi qi% sobrescreveb comD1bendforfork = 0, 1, . . . , kmaxdow Axkxk+1 = xkw + b rk+1 b wif || rk+1|| < t thenbreakendifendforendprocO exemplo abaixo ilustra o comportamento tpico do metodo de Jacobi:Exemplo4.7Resolva o sistema__4 1 1 01 4 0 11 0 4 10 1 1 4__x =__0110__cujasolu caoe x=(0, 1667, 0, 3333, 0, 3333, 0, 1667)T, usandoometododeJacobi comx0=(0, 0, 0, 0)Ta uma toler ancia = 102.Solu cao: Aplicando o metodo de Jacobi, obtemosx1= (0, 0, 25, 0, 25, 0)Tx2= (0, 125, 0, 25, 0, 25, 0, 125)Tx3= (0, 125, 0, 3125, 0, 3125, 0, 125)Tx4= (0, 1563, 0, 3125, 0, 3125, 0, 1563)Tx5= (0, 1563, 0, 3281, 0, 3281, 0, 1563)Tx6= (0, 1641, 0, 3281, 0, 3281, 0, 1641)Tx7= (0, 1641, 0, 3320, 0, 3320, 0, 1641)Tx8= (0, 1660, 0, 3320, 0, 3320, 0, 1660)Tou seja,com oito itera coes,obtemosuma aproximac aopara a solucao dentro da toleranciaespe-cicada.4.4.8 MetodoiterativodeGauss-SeidelAnalisandoometodo deJacobi,ve-se que, acadaiteracao,produzem-se todosos elementosdovetorxk+1, usandoapenasoselementosdovetorxk. Noentanto, nadaimpedeque, ` amedidaA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 82Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Linearesqueoselementosdexk+1saoproduzidos,elespossamserutilizadosparaproduzirospr oximoselementos do pr oprioxk+1. O metodo de Gauss-Seidel faz exatamente isso.De forma an aloga ao metodo de Jacobi, escrevemos, paran = 3,___xk+1,1= u12xk,2+ u13xk,3+ d1xk+1,2= l21xk+1,1+ u23xk,3+ d2xk+1,3= l31xk+1,1+ l32xk+1,2+ d3ou, em forma matricial,xk+1 = Lxk+1 + Uxk + d (4.44)ondeL= D1AL, U= D1AU, D=diag(A), d=D1beALeAUindicamaspor coesestritamente inferior e superior deA (isto e, sem a diagonal).NocasodometododeGauss-Seidel, podemosescreveracorrecao paraxk+1deformamaiscompacta; note que a expressao Lxk+1 + Uxkpode ser calculada comon

j=1j=i_aijaii_xj, i = 1, 2, . . . , nOcriterio deparada,noentanto,devesercalculado usandooresduork+1=b Axk+1, comomostra o algoritmo 4.4.3.Algoritmo4.4.3Metodo de Gauss-Seidelproc gauss seidel(input:A,b,x0,kmax,; output:xk+1)fori = 1, 2, . . . , n doqi a1iiendfort || q || || b ||fori = 1, 2, . . . , n doforj = 1, 2, . . . , n doaij aij qi% sobrescreveA comD1Aendforbi bi qi% sobrescreveb comD1bendforfork = 0, 1, . . . , kmaxdou xkfori = 1, 2, . . . , n doui b

nj=1j=iaijujendforxk+1 urk+1 b Axk+1if || rk+1|| < t thenbreakendifendforendprocO exemplo 4.8 ilustra o comportamento tpico do metodo de Gauss-Seidel:Exemplo4.8Resolva o sistema__4 1 1 01 4 0 11 0 4 10 1 1 4__x =__0110__A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 83Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Linearescujasolu caoe x=(0, 1667, 0, 3333, 0, 3333, 0, 1667)T, usandoometododeGauss-Seidel comx0 = (0, 0, 0, 0)Ta uma toler ancia = 102.Solu cao: Aplicando o metodo de Gauss-Seidel, obtemosx1= (0, 0, 25, 0, 25, 0, 125)Tx2= (0, 125, 0, 3125, 0, 3125, 0, 1563)Tx3= (0, 1563, 0, 3281, 0, 3281, 0, 1641)Tx4= (0, 1641, 0, 3320, 0, 3320, 0, 1660)Touseja, comquatroitera coes, obtemos umaaproximac aoparaasolucaodentrodatoleranciaespecicada.Da mesma forma que o metodo de Jacobi, uma condi cao necess aria e suciente para a conver-gencia do metodo de Gauss-Seidel e que a matrizA seja diagonal-dominante (ver equa cao 4.41).Existe um criterio de Sassenfeld que, se atendido, garante a convergencia do metodo. Para severicar se uma matriz de coecientes do sistema satisfaz a tal criterio, calcula-se os valoresS1,S2,. . .,Sn, denidos comoS1=1| a11 |(| a12| +| a13| + . . . +| a1n|S2=1| a22 |(| a21|S1 +| a23| + . . . +| a2n|...Sn=1| ann |(| an1|S1 + | an2|S2 + . . . +| ann1|Sn1(4.45)e, seSi< 1, 1 i nent ao o metodo de Gauss-Seidel ira convergir.Exemplo4.9Para a matriz do exemplo 4.8, verique se o criterio de Sassenfeld e atendido.Solu cao: Calculando os valores deSi, temos:S1=1| 4 |(| 1 | +| 1 |) = 0, 5 < 1S2=1| 4 |(| 1 |0, 5 +| 1 |) = 0, 375 < 1S3=1| 4 |(| 1 |0, 5 +| 1 |) = 0, 375 < 1S4=1| 4 |(| 1 |0, 375 +| 1 |0, 375) = 0, 1875 < 1e, comoSi< 1,1 i 4,ocriteriodeSassenfeldeatendidoe, porconseguinte,ometododeGauss-Seidel e convergente para um sistema com essa matriz de coecientes.OcriteriodeSassenfelde, noentanto, apenassuciente; umamatrizpoden aoatende-loe,mesmo assim, o metodo de Gauss-Seidel pode convergir, como mostra o exemplo abaixo.Exemplo4.10Resolva o sistema_1 11 3_x =_33_Solu cao: O criterio de Sassenfeld n ao e satisfeito pois, calculando os valores deSi, temosS1=1| 1 |(| 1 | = 1 < 1S2=1| 3 |(| 1 |1) = 0, 333 . . . < 1A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 84Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesNoentanto, o metodo de Gauss-Seidel converge para a solu cao x=(1, 5, 1, 5)Tem11itera c oes, a uma toler ancia de 104:x0= (0, 0000, 0, 0000)Tx1= (3, 0000, 2, 0000)Tx2= (1, 0000, 1, 3333)Tx3= (1, 6667, 1, 5556)Tx4= (1, 4444, 1, 4815)Tx5= (1, 5185, 1, 5062)Tx6= (1, 4938, 1, 4979)Tx7= (1, 5021, 1, 5007)Tx8= (1, 4993, 1, 4998)Tx9= (1, 5002, 1, 5001)Tx10= (1, 4999, 1, 5000)Tx11= (1, 5000, 1, 5000)TNote que a dominancia diagonal de uma matriz e relacionada com a ordem em que as equa coesse apresentam. Uma simples troca entre duas linhas pode ser desastrosa, como mostra o exemploa seguir.Exemplo4.11Seja o sistema apresentado no exemplo 4.10, com as linhas trocadas entre si, i.e._1 31 1_x =_ 33_Nessecaso, comoamatrizdosisteman aoediagonal-dominante, ometododeGauss-Seideldiverge, apresentando como primeiras estimativas os vetoresx0= (0, 0000, 0, 0000)Tx1= (3, 0000, 6, 0000)Tx2= (15, 0000, 12, 0000)Tx3= (39, 0000, 42, 0000)T...x8= (29523, 0000, 29526, 0000)Tx9= (88575, 000, 88572, 0000)T...x20= (5, 2302 109, 5, 2302 109)Tapesar do sistema ter a mesma soluc aox = (1, 5, 1, 5)T.4.4.9 Extrapola caodeummetodoiterativoUma das formas de garantir e/ou acelerar a convergencia de um metodo iterativo e utilizar umatecnica de extrapolacao,a qual consiste em se combinar a correcao da estimativaxk dada pelaequa cao governante do metodo iterativo com uma outra corre cao, semelhante. Em termos dasfun coes , isso pode ser expresso comoxk+1 = k(x0, x1, . . . , xk, A, b) + (1 )k(x0, x1, . . . , xk, A, b), IR (4.46)Note que, se = 1, temos o metodo iterativo original.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 85Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesPor exemplo, no caso do metodo de Jacobi, podemos usark(x0, x1, . . . , xk, A, b) = I (ou seja,a matriz identidade). Assim temos o metodo de relaxacao de Jacobi (JOR),xk+1 = (xkD1Axk + D1b) + (1 )xk, 0 < 1 (4.47)e,deforma an aloga, temoso metododasrelaxa coessucessivas(SOR),umavariante do metodode Gauss-Seidel,xk+1 = (Lxk+1 + Uxk + d) + (1 )xk, 0 < < 2 (4.48)Exemplo4.12Seja o sistema__2 1 0 11 2 1 00 1 2 11 0 1 2__x =__1331__cujasolu caoe(0, 1, 1, 0)T. Utilizando-seometodoJORpararesolve-lo, com=0, 65, aumatoler ancia = 102, obtemos:x0= (0, 0000, 0, 0000, 0, 0000, 0, 0000)Tx1= (0, 3250, 0, 9750, 0, 9750, 0, 3250)Tx2= (0, 0163, 0, 8937, 0, 8937, 0, 0163)Tx3= (0, 0349, 0, 9921, 0, 9921, 0, 0349)Tx4= (0, 0035, 0, 9884, 0, 9884, 0, 0035)Tx5= (0, 0038, 0, 9986, 0, 9986, 0, 0038)Tou seja, com cinco itera coes, obtemos uma aproximacao para a solucao dentro da tolerancia espe-cicada. Ometodo deJacobi,seutilizadopara resolvero mesmo metodo,n ao alcancaa solu caoapos 200 itera coes.Exemplo4.13Seja o sistema__4 1 1 01 4 0 11 0 4 10 1 1 4__x =__0110__cuja solu cao ex = (0, 1667, 0, 3333, 0, 3333, 0, 1667)T. Utilizando-se o metodo SOR para resolve-lo, com = 1, 1, a uma toler ancia = 102, obtemos:x0= (0, 0000, 0, 2750, 0, 2750, 0, 1513)Tx1= (0, 1513, 0, 3307, 0, 3307, 0, 1668)Tx2= (0, 1668, 0, 3336, 0, 3336, 0, 1668)Tou seja,comtres itera coes,obtemosuma aproximacaopara a solucaodentro da toleranciaespe-cicada (compare com o exemplo 4.8).4.5 MetododoGradienteO metodo do gradiente e indicado para resolver um SELA onde A e uma matriz simetrica, positivo-denida (SPD), i.e.xTAx > 0, x IRn(4.49)Umaoutracaractersticade matrizes SPDe que todos os seus autovalores saoestritamentepositivos.O metodo baseia-se na relacao existente entre a solucao de um SELA e a minimiza cao da formaquadr atica, quandoA for SPD. Assim, inicialmente veremos o que e a forma quadr atica e algunsexemplos da mesma.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 86Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes Lineares4.5.1 FormaQuadraticaA forma quadr atica e uma fun cao vetorial f: n dada porf(x) =12xTAx bTx + c (4.50)ondeA IRnn,x IRn,b IRnec IR.Por exemplo, considere o sistema_2 11 2_x =_11_(4.51)cuja solu cao ex= (1, 1)T. A matriz decoecientes e SPD e a gura 4.5 mostra o gr aco eascurvas denvel daforma quadr atica correspondente (comc = 0). Notequeo gr aco da fun caoe um parabol oide portanto, com apenas um ponto de mnimo e que aparentemente, o ponto(1, 1) (ou seja, a solu cao do sistema) e o ponto de mnimo da fun cao.J a as guras 4.6-4.8 mostram outras situacoes possveis, dependendo dos valores dos elementosdeA (todos os sistemas tiveram xada a sua solu cao em (1, 1) e os termos independentes foramcalculados adequadamente). Por exemplo, na gura 4.6, temos os gr acos para a matriz negativo-denida_ 2 11 2_i.e., xTAx < 0, x; note que a forma do gr aco e um parabol oide invertido, com apenas um pontode m aximo e a situacao oposta ` a de uma matriz SPD.Na gura 4.7, temos o caso em que a forma quadratica assume tanto valores negativos quantopositivos o gr aco da fun cao e a chamada sela. A matriz em questao e_1 44 5_Finalmente, ogr acoeascurvasdenvel paraaformaquadr aticaexibidos nagura4.8correspondem a uma matriz quase-singular,_1 22 3, 8_a qual apresenta innitas solucoes ao longo da reta na base do gr aco.Para vericarmos se isso e verdade, vamos calcularf

(x) = 0. Comofe uma fun cao vetorial,a sua derivada ou gradiente e dada porf

(x) =__x1f(x)x2f(x)...xnf(x)__(4.52)oqual representaumcampovetorial; paraumdadopontox, eleapontanadire caodemaiorvaria cao def(x). O gr aco def

(x) na gura 4.9 e tpico da situa cao em queA e SPD:Aplicando a equa cao (4.52) ` a (4.50), obtemosf

(x) =12ATx + 12Ax b (4.53)e, seA e simetrica, A = AT, de ondef

(x) = Ax b.Igualandof

(x)azero, obtemosAx=b, ouseja, osistemaquequeremosresolver. Portanto,podemos dizer queA.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 87Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.5: Gr aco def(x) e suas curvas de nvel paraA SPD.Figura 4.6: Gr aco def(x) e suas curvas de nvel paraA ND.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 88Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.7: Gr aco def(x) e suas curvas de nvel paraA indenida.Figura 4.8: Gr aco def(x) e suas curvas de nvel paraA singular.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 89Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.9: Gr aco def

(x) paraA SPD.minimizar a forma quadr atica f(x) =12xTAxbTx+c equivale a resolver o sistemaAx = b seA for simetrica.Se amatrizA e positivo-denida, alem desimetrica, ent ao a solu cao deAx =b e o mnimo( unico) de f(x); logo, para A SPD, a solu cao x = A1b e o ponto x que minimizaf(x). Isso podeser mostrado como segue.SuponhaAsimetrica, xumvetorquesatisfazAx=b, yumvetorsimilarax(emtermosgeometricos, y e um ponto pr oximo ax) ee = y x o vetor erro; ent ao,f(x + e) =12(x + e)TA(x + e) bT(x + e) + c=12xTAx + eTAx + 12eTAe bTx bTe + c... b = Ax; A = AT...=12xTAx bTx + c + eTb + 12eTAe bTe= f(x) + 12eTAe (4.54)ef(x + e) = f(x x + y) = f(y) = f(x) + 12(y x)TA(y x) (4.55)Agora, comoA e SPD, por hip otese, entao(y x)TA(y x) > 0, ye, portanto, f(y) > f(x). Isso mostra quex e o mnimo def(x), nesse caso.4.5.2 Descri caodometododoGradienteOgr acodaformaquadr atica, para ASPD, nos sugereumaestrategiaparalocalizarmos asolu cao do sistema: basta escorregar ao longo das paredes do paraboloide, pois isso nos levara,necessariamente, aopontodemnimo. Aquest aoquesecolocaagora e: qual dire caodevemostomar, a partir de umxk, para obtermos umxk+1que seja mais pr oximo da solu cao?A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 90Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.10: O vetor erroe0e o vetor resduor0.Lembramos queogradientef

(x)apontanadire cao demaioraumentodef(x), emsentidooposto ao fundo do parabol oide.E natural, portanto, que andemos ao longo da dire cao opostaao gradiente, isto e, f

(x) = b Ax. Ora, conforme visto anteriormente, rk = b Axk, de ondeestabelecemos as seguintes rela coes entre o vetor resduo e o gradiente def(x):rk= f

(x) (4.56)rk= b Axk... ek = xkx...= b AxAek... x = A1b...= Aek(4.57)A equa cao (4.56) nos diz que o resduo tem a mesma direcao do gradiente, porem sentido oposto;j aaequa cao(4.57)nosdizqueoresduo eovetorerro, transformadoporA(e, portanto, nomesmo espa co deb).Suponha, ent ao, que temos a seguinte situacao, conforme a gura 4.10. Como decidimos andarao longo do vetor resduo, a partir de x0, a nova estimativa x1 e um ponto sobre a reta r0, ou sejax1 = x0 + 0r0, 0 IR (4.58)O escalar0indica o deslocamento sobrer0. Para determinar o melhor0 ou seja, aquele parao qual || x1x|| e mnimo derivamosf(x1) em relacao a0e igualamo-la a zero:dd0f(x1) = f

(x1)Tdd0x1 = f

(x1)Tr0(4.59)Notequef

(x1)Tr0eoprodutoescalar entreos vetoresf

(x1)er0. Como oprodutoescalar edado poruTv = || u || || v || cos onde e o angulo formado entre os vetores u e v, ao igualarmos f

(x1)Tr0 a zero, estamos exigindoqueosvetoresf

(x1)er0sejamortogonaisentresi. Comof

(x1)= r1, issoimplicaquedoisresduos sucessivos sao ortogonais entre si; a gura 4.11 mostra tres situa coes tpicas para a solucaodo sistema (4.51), com as seq uencias de resduos (representados pelas retas) gerados pelo metododo Gradiente a partir de tres diferentes estimativas iniciais.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 91Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.11: Caminhostpicos no metodo do Gradiente.A partir da equa cao (4.59), pode-se o obter o valor de0:f

(x1)Tr0= 0... f

(x1) = r1...rT1 r0= 0(b Ax1)Tr0= 0(b Ax00Ar0)Tr0= 0(r00Ar0)Tr0= 0de onde0 =rT0 r0rT0 Ar0(4.60)O que signica minimizarf

(x1)?Os gracos mostrados na gura 4.12 mostram que, para0calculadoconformeaequa cao(4.60), anovaestimativax1corresponde aomnimodapar abolaobtida como se tivessemos cortado o parabol oide f(x) por um plano vertical ao plano xy quepassa pela retar0!Utilizandoasequacoes(4.58)e(4.60), alemdaexpressaoparaoresduo, devidamentege-neralizadasparaak-esimaiteracao, podemosescreverumalgoritmoquedescreveometododoGradiente. Antes, porem, note quer1 = r00Ar0conformeobtidonaderiva caodaequacao(4.60); essaexpress aonos permite economizar umproduto matriz-vetor da forma Axk, pois Ark j a tera sido calculado previamente para se obter k.O algoritmo pode ser, ent ao, escrito como segue.A.L. de Bortoli, C. Cardoso, M.P.G. Fachin, R.D. da Cunha 92Introdu c ao ao Calculo Numerico Resolucao de Sistemas de Equacoes LinearesFigura 4.12: e escolhido de tal forma quef(x1) e mnima.Algoritmo4.5.1Metodo do Gradienteproc gradiente(input:A,b,x0,kmax,; output:xk+1)t || b ||r0 b Ax0fork = 0, 1, . . . , kmaxdowk Arkk rTk rkrTk wkxk+1 xk + krkrk+1 rkkwkif || rk+1|| < t thenbreakendifendforendprocO exemplo seguinte ilustra o comportamento tpico do metodo do Gradiente:Exemplo4.14Resolva o s